Machine Learning for Infectious Disease Risk Prediction: A Survey

Some studies have used generalized linear models (GLMs) to predict disease risks in a single location or multiple locations. The general formulation of GLMs is given as follows [42]:

$$p(y|\bm{x},\bm{w},\sigma^{2})=\exp\left[\frac{y\bm{w}^{T}\bm{x}-A(\bm{w}^{T}\bm{x})}{\sigma^{2}}+\log h(y,\sigma^{2})\right],$$

(5)

where $y$ denotes the response variable to be predicted, $\bm{x}$ denotes the input feature vector used to predict $y$, $\bm{w}$ denotes the weighting parameter vector on $\bm{x}$, $\sigma^{2}$ represents the dispersion term, $A(\cdot)$ represents the log normalizer, $h(\cdot,\cdot)$ represents the base measure, and $\exp[\cdot]$ denotes a given probability distribution from the exponential family. Furthermore, there is a mean function $g^{-1}$ that maps $\bm{w}^{T}\bm{x}$ to the mean value of the response variable: $\mu=g^{-1}(\bm{w}^{T}\bm{x})$ [42]. In different instances, different probability distributions are used to model disease risks. For example, in [22], Zhang et al. used the Poisson distribution to model case numbers as integer values. In [23, 24], Pei et al. used the Gaussian distribution to model the disease case numbers as continuous values; they also used the Bernoulli distribution to model the status of getting infected. Based on the distribution selected for modeling the data, the mean function varies accordingly, such as $\mu=\exp(\bm{w}^{T}\bm{x})$ for Poisson distribution, $\mu=\bm{w}^{T}\bm{x}$ for Gaussian distribution, and $\mu=\frac{1}{1+e^{-\bm{w}^{T}\bm{x}}}$ for Bernoulli distribution. Although GLMs have a similar regression framework to generate predictions, their specific structures and the correpsonding inference methods in different works are different, as they are elaborately designed based on various assumptions to reflect specific disease transmission processes. For example, Zhang et al. incorporated the effects of intra-regional, inter-regional, and external factors on disease risk into a unified Poisson regression-based framework to model epidemic diffusion between multiple locations [22]. In their approach, information from prior knowledge about disease transmission is taken into consideration to specify what heterogeneous risk-related factors are involved in three aspects of disease transmission: intra-transmission, inter-transmission, and external transmission. Specifically, in the intra-transmission part of their framework, climate data (e.g., temperature, rainfall), geographical data (e.g., elevation), and demographic data (population) are combined to predict self-infections within a region; in the inter-transmission part of their framework, a diffusion matrix with a prior structure constrained by the transportation network within a region is used to describe disease transmission between locations; and in the external-transmission part of their framework, a quadratic function that has unimodal patterns is employed to model the effect of seasonal imported cases of the disease. Similarly, Pei et al. have utilized a multivariate regression model denoted the group sparse Bayesian learning model (GSBL), which is based on a transmission network, to predict disease dynamics [23, 24]. In contrast to [22], they focused on the sentinel selection problem, which arises due to limited disease surveillance resources. Sentinel selection, which is closely related to the concept of active surveillance in the domain of public health, is the selection of representative locations from all targeted locations at which to conduct disease surveillance. Pei et al. have formulated this problem as a learning process of a row-sparse disease transmission network, in which sentinels are indicated by the non-zero rows [23, 24]. Thus, their model uses the disease data of these sentinels and an inferred network to recover or predict the global dynamics of all target locations. In addition, the model does not use other prior knowledge about disease transmission and only uses historical case number data to infer a transmission network, so it can be more easily applied than the above Poisson regression model to various diseases and other domains.

In contrast to GLMs, which assign a pre-defined prior distribution to the parameters of a model (i.e., a distribution that is irrelevant to the observational data), empirical Bayesian models usually estimate prior distributions from historical observations. A typical example is the semiparametric empirical Bayes framework for epidemic modeling proposed by Brooks et al. [25]. This framework first estimates the prior, i.e., the shape of the influenza-like illness (ILI) curve, the noise, the peak height, the peak week, and the pacing, using a set of uniform distributions over the historical observations. It then generates the underlying ILI curve of the current ILI season by linearly adjusting the piecewise quadratic curves of historical seasons using the current year’s CDC baseline weekly ILI level.

To predict disease risks, Gaussian process (GP) models assume that the random variable $f(\bm{x}_{i})$ in continuous domains (e.g., time or space) follow a Gaussian distribution with the mean $\mu=m(\bm{x}_{i})$ and the variance $\sigma_{i}$, and that the joint distribution of a finite set of these variables $\bm{f}=[f(\bm{x}_{1}),\cdots,f(\bm{x}_{M})]$ will follow the multivariate Gaussian distribution with the mean $\bm{\mu}=[m(\bm{x}_{1}),\cdots,m(\bm{x}_{M})]$ and the covariance $\bm{\Sigma}_{ij}=\mathcal{K}(\bm{x}_{i},\bm{x}_{j})(i,j=1,\cdots,M$), where the $\mathcal{K}$ is the kernel function and $M$ is the number of observations [43, 42]. In GP models, given the training set $\mathbf{X}=\{{\bm{x}_{1},\cdots,\bm{x}_{M}}\}$, we have $\bm{f}_{X}\sim\mathcal{N}(\bm{\mu}_{X},\mathbf{K}_{X,X})$, where $\mathbf{K}_{X,X}$ is the $M\times M$ covariance matrix of the data set $\mathbf{X}$. Given a test set $\mathbf{X}_{*}$, the joint distribution $p(\bm{f}_{X},\bm{f}_{X_{*}}|\mathbf{X},\mathbf{X}_{*})$ is represented as follows:

$$\left(\begin{array}[]{l}\bm{f}_{X}\\ \bm{f}_{X_{*}}\end{array}\right)\sim\mathcal{N}\left(\left(\begin{array}[]{l}\bm{\mu}_{X}\\ \bm{\mu}_{X_{*}}\end{array}\right),\left(\begin{array}[]{cc}\mathbf{K}_{X,X}&\mathbf{K}_{X,X_{*}}\\ \mathbf{K}_{X,X_{*}}^{\top}&\mathbf{K}_{X_{*},X_{*}}\end{array}\right)\right).$$

(6)

The covariance of these variables is calculated by choosing the appropriate kernel function $\mathcal{K}(\cdot,\cdot)$ and is used to describe the characteristics of processes. Due to the inherent ability of a covariance matrix to model the similarity between data points, conventional GP models are generally used as interpolation models. However, some recent studies have extended their use by applying them to epidemic prediction tasks. For instance, Senanayake et al. proposed a model based on GP regression that predicts influenza cases by capturing the spatiotemporal dependency of data [26]. They constructed a non-linear kernel with both spatial and temporal components, and spatiotemporal covariance components, to address the challenges associated with the complicated characteristics of disease dynamics, such as temporal characteristics (i.e., periodicity, non-stationarity, and short- and long-term dependency) and spatial characteristics (i.e., the distance between locations and morphology of a region). Zimmer and Yaesoubi proposed a GP-based framework to forecast seasonal epidemics [27]. In contrast to Senanayake et al. [26], they did not design kernels to represent the dependencies between spatial locations but rather focused on exploring the temporal dependency between within-seasonal and between-seasonal time series.

These methods, which are popular in the field of recommender systems [44, 45], are also used to predict disease risks. For instance, Chakraborty et al. proposed matrix factorization with nearest-neighbor regression (MFN), which incorporates MF regression and nearest-neighbor-based regression, for ILI count prediction [28]. In their MFN model, they integrated disease-related features, historical disease dynamics, and the disease dynamics to be predicted across time into a prediction matrix. Then, they factorized the prediction matrix as a factor-feature matrix and a factor-prediction matrix, such that the prediction matrix could be reconstructed by multiplying the factor-feature matrix and factor-prediction matrix. Subsequently, they incorporated nearest neighbor regression to correct the reconstructed prediction matrix with the $K$ nearest samples.

Some studies have improved the robustness of predictions by generating them from ensemble models rather than from a single model. For instance, [28] used an ensemble model derived from the fusion of outputs of models trained on data from various sources. Specifically, they trained multiple MFN models on data whose effects on disease dynamics had been previously studied and combined their results to give the final prediction. Similarly, [27] trained various GP regression models on different features of disease-related data and aggregated these models’ results to generate the final predictions.

RNNs are widely used to model the temporal dependency of time series data, such as voice or text data. RNN modules are formed based on the assumption that the current output is not only related to the current input but also depends on the previous inputs. Thus, as infectious disease dynamics are a type of time series data, they can also be modeled by RNNs. For example, the interactively and integratively connected deep recurrent neural network model (I²DRNN) [29] uses stacked RNN modules to capture spatiotemporal dependencies from heterogeneous and multiple-scale risk-related data. The model structure contains three components: (1) an input module, which is used to integrate heterogeneous (i.e., fine-, coarse-, and same-scale) data; (2) a hidden module, which is designed as a hierarchical structure to extract dependencies from different locations and the heterogeneous factors of different scales; (3) and an output module, which is used to generate final predictions based on extracted hidden features. RNN architectures based on a gating mechanism, such as a long short-term memory (LSTM) network [46], have also been used in recent studies for disease risk prediction, due to their ability to preserve the long-term information of data sequences. Venna et al. proposed an LSTM-based deep learning model that consists of multiple LSTM cells sequentially stacked over time [30]. In the sequential structure, every LSTM cell takes two inputs (except the first cell, which only takes the dynamic data as input): the dynamic data at a single time point and the output of the previous cell, to generate a prediction for the next time step. That study also examined the effects of climate variables by applying the symbolic time-series approach and the effects of regions with geographical proximity by applying weighted summation to adjust the output of an LSTM to generate the final prediction. Volkova et al. adopted two LSTM layers to learn temporal dependencies from the ILI case number and social media data, respectively, and merged their outputs via a fully connected layer to generate predicted ILI proportions [31].

In contrast to RNN models, which capture the temporal dependency of sequential data, GNN models can deal with data with graphical structures [47]. That is, a GNN model is based on a graph with a given structure and encodes structural information by passing messages between nodes of the graph. Due to this ability of GNN models to capture characteristics in graph structures, they are naturally used to capture and represent spatial patterns of disease dynamics, which can be regarded as driven by a disease transmission network. For example, the spatio-temporal graph neural network (STGNN) [32] utilizes daily mobility data from Google to construct the structure of time-varying disease transmission networks. Based on their constructed network, Kapoor et al. designed two types of edges, i.e. edges between nodes within the network at the same time, and edges between nodes within the network at the current time and nodes within the network at the previous time, to characterize varying spatiotemporal dependencies driven by cross-regional human mobility and the effect of historical risk trends, respectively [32]. Moreover, GNNs are not limited to modeling intuitive spatial relationships by delineating the network structure between locations; they can also be used to model the dependency between extracted features. For instance, the spectral temporal graph neural network (StemGNN) [33] uses a graph convolutional network (GCN) structure to model temporal dependency to predict newly confirmed COVID-19 cases. Specifically, instead of modeling the time series in the time domain, it utilizes the graph Fourier transform (GFT) to model inter-series correlations within the spectral domain and the discrete Fourier transform (DFT) to model intra-series temporal correlations within the frequency domain, and then feeds the representation of correlations into a GNN.

More recent studies have applied combinations of multiple neural network structures to model complex spatiotemporal patterns of disease transmission. These models with mixed neural network structures make full use of the aforementioned common neural network structures, i.e., RNNs, GNNs, and CNNs, to design new composite architectures serving various purposes. Usually, these architectures contain two separate modules—i.e., a spatial module and a temporal module—that are connected to form an integrated model that is subsequently optimized in an end-to-end manner to capture and model spatial and temporal dependencies simultaneously. For example, Wu et al. proposed a model which incorporates CNN, RNN, and residual structures, named CNNRNN-Res, to capture spatiotemporal dependencies in historical disease dynamics [34]. The CNN module uses an adjacency CNN filter to represent the adjacent graph of different regions, which is employed to integrate the information from neighbors. The RNN module uses a gated recurrent unit (GRU) to capture the temporal correlation in data. To solve the problem of overfitting, a sparse residual link is used to skip connections with some previous layers. The Cross-location attention-based GNN (ColaGNN) [35] is designed for long-term ILI prediction. It uses location-aware attention to infer the spatial influence between different regions from learned hidden features (denoting temporal dependencies), which are extracted from an RNN module. In addition, it employs a dilated convolution module to learn attributes for each node from historical disease trends and thereby capture multiple-scale local temporal dependency. Based on the above-mentioned network structure and node attributes, a graph message-passing mechanism is used to integrate the spatiotemporal information, which is then used to generate ILI predictions. The Hierarchical spatial-temporal framework (HierST) [36] includes a temporal module that combines two time-series architectures—the long- and short-term time series network (LSTNet) and the neural basis expansion analysis for time series (N-BEATS)—to model temporal dependency; and a spatial module that contains the gated edgeGNN, which adaptively adjusts the connections of edges, and the nodeGNN, which learns the representation of node features. The novelty of this approach is also reflected by the introduction of prior knowledge of common sense to constrain the model inference. Specifically, given that the predictions for different administrative levels (i.e., country, state, and county) should be close to each other, [36] designed a consistency optimization objective that includes items representing the difference between predictions at different spatial scales in addition to the difference between ground truth and predictions. The epidemic forecasting model based on functional neural process (EPIFNP) proposed by Kamarthi et al. [37] also includes temporal and spatial modules, which are implemented by a probabilistic neural sequence encoder and a stochastic correlation graph, respectively. Instead of generating point estimates of forecast value, the EPIFNP model generates the probability distribution of prediction via a probabilistic generative process model to evaluate the uncertainty of prediction. The population-level disease prediction model (named PopNet) proposed by Gao et al. assumes that an undirected disease transmission network drives disease dynamics [38]. PopNet learns the connection structure of this network by using population and geographical distance to calculate the similarity of each pair of locations. Then, based on the learned network structure, PopNet uses two graph attention networks (GATs) to obtain node embedding from real-time disease data and updated disease data (i.e., revised data released after the initial release), respectively. Subsequently, PopNet fuses these two kinds of node embeddings by spatial latency-aware attention (S-LAtt) and temporal latency-aware attention (T-LAtt), sequentially. S-LAtt uses a feature similarity-based attention mechanism and considers the marginal effects of time latency on final predictions to learn the edge weights between pairs of nodes together to update node embedding. T-LAtt uses GRU networks to learn temporal dependency. Finally, PopNet concatenates the learned node embeddings to generate final predictions.

The encoder–decoder framework is a general framework that consists of different deep modules to manage sequential data and was initially used in machine translation [48]. As trends in disease dynamics exhibit a sequential dependency that is similar to that of sentences of text, an encoder–decoder framework with an RNN structure can also be used to predict epidemic trends within a given time period [39, 40]. Adhikari et al. proposed EpiDeep to predict weighted ILI (wILI) using an encoder–decoder framework, together with deep clustering components [39]. EpiDeep uses an LSTM-based encoder to encode an input influenza sequence as latent variables that contain temporal information, and a deep clustering component (an improved deep embedded clustering (IDEC) module [49]) to learn the embedding of the existing observed epidemic trend in the current season whose trend is to be predicted, and then clusters this embedding with the most similar epidemic trends in historical seasons. EpiDeep also uses this approach to learn and cluster the embedding of full-length historical trends. Next, it learns a mapping function to map the embedding of the incomplete sequence to the space of the full-length sequence. Finally, EpiDeep uses a decoder to predict the future sequence of the epidemic trend in the current season by taking the mapped clustering embedding and the encoded trend (both are in the current season) as inputs. Cui et al. also used an encoder–decoder framework to predict the dynamics of COVID-19 pandemic [40]. The encoder component employs several CNN modules with different kernel sizes to extract temporal features within multiple time ranges from the data of case numbers and regional visitor counts. They designed a graph-based module that characterizes spatial patterns by modeling human mobility and infection processes in each range in which graph structure is learned from data by the attention mechanism, and then fuses features learned from each range into one feature via a multi-headed self-attention mechanism. The decoder component employs a temporal embedding module to embed the case numbers and death numbers, and the obtained embedding is fed into the multi-head attention layers with the output of the encoder. Finally, the features obtained from the encoder and decoder are passed through a multilayer perceptron (MLP) to generate a final prediction.

However, knowledge on the properties and transmission of emerging diseases may be ambiguous, and there may be a paucity of observational data on epidemic trends. Additionally, the data for regions with inadequate systems for surveillance and reporting of infectious diseases may be incomplete and noisy. Therefore, it is of interest to determine how to use abundant and high-quality data on similar diseases to the disease of interest or for regions with similar characteristics to the region of interest to facilitate the prediction of a given disease epidemic. One approach used is transfer learning (TL) architecture, which consists of one source task and one target task, and is designed to transfer knowledge learned from the previous model to enable the learning of the target task [50]. Similarly, the COVID augmented ILI deep network (CALI-NET) is a heterogeneous transfer learning (HTL) framework for COVID-ILI forecasting [41] that applies the EpiDeep [39] model as the source model to learn representations of temporal dependency from historical wILI data. [41] designed the COVID-augmented exogenous model (CAEM) to encode representations of spatiotemporal features of exogenous data signals of COVID-19 by Laplacian regularization of a geographical adjacent matrix and a GRU module, and used these encoded representations in the target model. They also designed a knowledge distillation (KD) loss, which consists of the hint loss between the mapped representations from source and target models, and the imitation loss between source predictions and ground truth to ensure effective knowledge transfer.

The above-described studies demonstrate that previously developed deep learning-based models use spatial and temporal modules, which are designed to capture disease patterns over space and time, respectively. The relationship between spatial information is typically modeled by a GNN module. However, using a GNN module is challenging, due to its unknown network structure. Some approaches pre-build connections in a network by using proxy data, such as human mobility data [32], population and geographical distance data [38], and regional demarcation data [41], and then feed the network into a GNN module or constrain model learning by applying graph Laplacian regularization. If there is a lack of related data, a graph structure can be constructed by using the attention mechanism to infer the weights of edges in a network structure during model optimization [35, 36, 38, 40]. Temporal patterns are modeled by mining the temporal dependency of time series data using various RNN structures, such as a classic RNN [35], a GRU [34, 41, 37, 38], LSTM [39], or an LSTNet/N-BEATS [36]. However, in addition to using an RNN and GNN to characterize temporal and spatial patterns, respectively, these neural network modules can be flexibly designed to serve these purposes. For instance, in addition to RNN modules being used to extract temporal dependencies, CNN modules (such as temporal convolution modules [40] or dilated CNNs [35, 38]) and GNN modules (such as spectral temporal GNNs [33]) can also be used for this purpose; and in addition to GNN modules being used to capture spatial dependencies, a modified CNN module can also be used for this purpose [34].

Data assimilation techniques, which are widely applied in atmospheric and oceanic sciences and in numerical weather forecasting [90], aim to utilize observations to optimize mechanism-based models. Thus, they have also been applied in disease dynamic prediction [11, 51, 52, 53]. In a set of data assimilation models used for epidemic prediction, the Kalman filter (KF) and its variants [91, 92] and particle filter (PF) [93] methods have been used to estimate model statuses. For instance, Shaman and Karspeck applied data-assimilation techniques to the problem of influenza forecasting and generated retrospective ensemble forecasts of influenza seasons from 2003 to 2008 in New York City, USA [11]. They proposed the SIRS–EAKF framework, which uses the ensemble adjustment Kalman filter (EAKF) [92]) and a PF [93] to assimilate the observations of infections (i.e., estimates of influenza infections from Google Flu Trends) into the susceptible–infectious–recovered–susceptible (SIRS) model [94]), which is a humidity-forced compartmental model. The SIRS–EAKF framework can estimate the posterior of probabilistic distributions of system states (i.e., susceptible populations $S_{t}$ and infected populations $I_{t}$) and epidemiological parameters (e.g., the mean infectious period $D$, the average duration of immunity $L$, and the maximum and minimum of daily basic reproductive number $R_{0max}$ and $R_{0max}$) in the used SIRS model. In [94], Shama et al. represented model states and epidemiological parameters by a set of variables $Z_{t}=(S_{t},I_{t},R_{0max},R_{0min},L,D)$. Then the posterior of $Z_{t}$ can be represented as follows:

$$p(Z_{t}|y_{t},y_{t-1},\cdots)\propto p(y_{t}|Z_{t})p(Z_{t}|y_{t-1},\cdots),$$

(9)

where the first term on the right-hand side is the likelihood of observational disease risk given states and parameters, while the second term is the prior distribution of the states and parameters. For KF, these two terms are assumed to be Gaussian distributions; in contrast, for PF, these two terms are not under these assumptions. Subsequently, Shaman et al. used similar data-assimilation techniques to generate weekly influenza forecasts for the influenza season in 2012 and 2013 in 108 cities in the USA [51]. Yang et al. tested four types of compartmental models and two types of filter models in their model-data assimilation framework and analyzed the epidemiological characteristics of influenza dynamics from the 2003–2004 season to the 2012–2013 season in 115 cities in the USA [52]. Pei et al. developed a model-data assimilation framework based on a metapopulation compartmental model to accurately predict the spatial spread of influenza [53]. In this metapopulation compartmental model, which is based on a humidity-driven SIRS model [94] that they had used in their previous studies [11, 51, 52], they divided a population into different groups in terms of geographical locations (i.e., different states), and incorporated two types of human mobility (i.e., fixed commuting flows and irregular movement of visitors).

Balcan et al. proposed the global epidemic and mobility (GLEaM) computational model for simulating infectious disease transmission [95, 96]. The GLEaM model is a global model based on a stochastic compartmental model at the meta-population level (i.e., a SLIR model considering both symptomatic and asymptomatic infections) and incorporates multiscale human mobility (short-range commuting and long-range airline flows) to effectively capture disease transmission patterns. Tizzoni et al. used the GLEaM model to model the disease transmission of 2009 H1N1 influenza and utilized the Monte Carlo maximum likelihood method to estimate some parameters [54]. Similarly, Zhang et al. proposed an epidemic computational framework based on the GLEaM model [55]. The computation is performed via three steps: (1) microblogging data from Twitter and surveillance data are used to estimate initial infections; (2) epidemiological parameters are searched in four-dimensional space by running Monte Carlo simulations with selected sampling points, and the GLEaM model is used to generate simulations; (3) a set of best-fit models is selected by using a multi-model information approach, which minimizes the loss of information (which is calculated by the Akaike information criterion (AIC)).

Aside from data-assimilation methods and simulation-based methods, some machine learning approaches are proposed to estimate the model states and epidemiological parameters. In these works, the loss function is generally formulated as the difference between states simulated using epidemiological models and the ground truth of these states. The general formulation of such loss function can be represented as follows:

$$\mathcal{L}_{e}(\bm{\theta}_{e})=\ell(\hat{\bm{y}}_{e},\bm{y}_{e})=\ell(f_{e}(\bm{\theta}_{e},\bm{s}_{0}),\bm{y}_{e}),$$

(10)

where $\bm{\theta}_{e}$ (in the parameter space $\bm{\Theta}_{e}$) denotes epidemiological parameters (e.g., contact rate and recovery rate) in a given epidemiological model, $\bm{y}_{e}$ denotes the ground truth of the target variable (usually are model states with records, e.g., infected case number and death number), $\hat{\bm{y}}_{e}$ denotes predictions on the target variable, $\bm{s}_{0}$ denotes the initial value of model states, and $f_{e}$ denote the given epidemiological model, which is generally described by an ODE. With such an ODE representation, $\hat{\bm{y}}_{e}$ can be calculated using model parameters and initial values. Given the above loss function, the optimal model parameters can be inferred by minimizing the loss:

$$\hat{\bm{\theta}}_{e}=\arg\min_{\bm{\theta}_{e}\in\bm{\Theta}_{e}}\mathcal{L}_{e}(\bm{\theta}_{e}).$$

(11)

For example, Zou et al. formulated a loss function with a logarithmic-type mean square error (MSE) [56]. Based on this loss function, parameters can be optimized by the general gradient-based optimizer. Moreover, Zou et al. also developed a novel compartmental model, named the SuEIR model—an improved SEIR model that considers a scenario of untested or unreported cases of COVID-19—and trained it with their machine learning approach. Wang et al. formulated a similar loss based on MSE to estimate the parameters of epidemiological models. Based on this loss function, they formalized the learning procedure as the AutoODE algorithm, which infers the model parameters of mechanism-based models by an automatic differentiation method. In addition, based on a case study on the forecasting of COVID-19 dynamics, they proposed the spatiotemporal SuEIR model, which is an extension of the SuEIR model [56] that better models spatiotemporal patterns of COVID-19 spread. The power-law degree and data priori jointly regularized non-negative network inference ($D^{2}PRI$) approach of Wang et al. [58] is based on a SIR model at the meta-population level. As this model regards the infectious interactions between individuals at different locations as a transmission process in a disease propagation network, Wang et al. [58] formulated the parameter inference of edge weights in the network and disease transmission rate in the SIR model as an integrated network inference problem. Moreover, according to prior knowledge of network structure, i.e., the power-law distribution of node degree and the features extracted from mobility-related data, they designed corresponding regularization items to constrain the parameter inference.

In contrast to models that infer values or probabilistic distributions of model parameters from observations, other models estimate the variation of epidemiological parameters and formulate them as functions of covariates. For instance, Arik et al. proposed the use of time-varying functions to model parameters [59]. That is, they used an improved compartmental model that is based on the SEIR model: instead of using the static epidemiological parameters in the traditional compartmental model, they used learnable functions to estimate parameter values from various covariates, which enable parameter values to vary over time. Specifically, they used the generalized additive model to encode the effects of covariates on epidemiological parameters. Baek et al. predicted the disease dynamics of multiple regions by using a stochastic SIR model [60]. This stochastic model employs a mixed-effects model that incorporates a random-effects term within each region and a fixed-effects term between different regions to encode the effects of static and time-varying covariates on the disease transmission rate.

Some studies have utilized the epidemiological concept to guide the construction of model structures. For instance, rather than directly describing and predicting malaria dynamics by using ODEs, Shi et al. constructed a nonlinear stochastic model based on the formula for vectorial capacity (VCAP) and entomological inoculation rate (EIR), which are epidemiological concepts derived from corresponding ODEs that describe malaria transmission [61]. VCAP is defined as the daily rate of future inoculations from mosquitoes to humans caused by a currently infected human case [97], whereas EIR is defined as the number of infectious bites received from mosquitoes per day by a human [97]. Based on the epidemiological meaning of local transmission, the local infections at time $t$ can be formulated as a function involving the VCAP/EIR and the infections at time $t-1$. Then, by considering the effects of cross-regional transmission, Shi et al. used periodic function modeling to depict the periodic transmission patterns [61]. Thus, their nonlinear stochastic model consists of the items of local infections and imported infections. Liu et al. developed a multivariate regression model based on a next-generation matrix of a meta-population vector–human compartmental model to predict the malaria risk in multiple locations[62]. The next-generation matrix [98, 99], consisting of the non-linear relationships between epidemiological parameters that can be derived from a compartmental model, represents the change in model variables from one time step to the next time step, thereby enabling the prediction of disease dynamics.

Some studies have added epidemiological constraints and regularizations, which are derived from compartmental models, to standard objective functions of supervised machine learning models to aid model parameter optimization. Hua et al. proposed the social media based simulation (SMS) model for influenza dynamics prediction [63]. This model incorporates two learning spaces: the social media space, which is designed to identify individuals’ health statuses from social media posts; and the epidemiological simulation space, in which a transmission network is built to simulate disease propagation between individuals. These two spaces are linked by minimizing the loss in terms of the inconsistency between the health status at the population level, which is obtained from the social media space and the simulation space. Kargas et al. applied epidemiological constraints in tensor factorization approaches to predict disease dynamics [64] by devising spatio-temporal tensor factorization with epidemiological regularization (STELAR). A tensor is an intuitive and natural structure used to represent and preserve the complex structure of high-dimensional data, especially spatiotemporal data with multiple risk-related factors. Tensor factorization is usually employed for dimensionality reduction and data decomposition rather than to predict disease transmission dynamics. STELAR enables the prediction of long-term epidemic trends by the addition of the latent epidemiological regularization of the SIR model into a standard tensor factorization method, i.e. canonical polyadic decomposition (CPD). Osthus et al. proposed a dynamic Bayesian (DB) influenza forecasting approach that models discrepancies between mechanistic model-generated simulations and observations [65]. This approach assumes that the uncertainty of prediction cannot be fully explained by observational noise and therefore models a wILI as the sum of three items: the logit of the infections that are described by the SIR model, a common discrepancy item for all influenza seasons, and a specific discrepancy item for each influenza season.

Some deep learning models, such as GNNs [66, 67] and RNNs [68, 69], also incorporate mechanistic models to constrain the learning of model structures and parameters, such that they effectively fit the realistic situation of disease transmission. In addition to the prediction loss, some of the aforementioned methods further introduce the epidemiological loss to ensure that the prediction of deep learning models is consistent with the dynamics described by the epidemiological model. The general formulation of the loss function is shown as follows:

$$\mathcal{L}=\mathcal{L}_{d}(\bm{\theta}_{d})+\mathcal{L}_{e}(\bm{\theta}_{e})=\ell(\hat{y}_{d},y_{d})+\ell(\hat{\bm{y}}_{e},\bm{y}_{e}),$$

(12)

where $\mathcal{L}_{d}$ denotes the prediction loss to ensure the prediction accuracy of the deep learning model, and $\mathcal{L}_{e}$ denotes the epidemiological-constrained loss. The approaches to introduce $\mathcal{L}_{e}$ vary in different works. For example, spatio-temporal attention network (STAN) is a GAT model with epidemiological constraints that is designed for long-term pandemic prediction [66]. In the network, nodes denote different locations and have both static features (latitude, longitude, population size, and population density) and dynamic features (case number, and the information on hospitalizations across all timestamps), whereas the edges denote disease transmission. The weights of edges in the STAN are calculated in terms of geographical proximity and population size. After constructing the network structure of the STAN, the graph attention mechanism is used to update node attributes and feed them to a GRU to extract temporal features. Epidemiological constraints are incorporated into STAN learning and prediction, as it generates two types of output by multitasking prediction settings: (1) epidemiological parameter predictions (i.e., the transmission rate and recovery rate); and (2) disease dynamics predictions (i.e., the increases in the infected and recovered case numbers). Gao et al. [66] also designed the loss function for model optimization based on the above-mentioned two kinds of outputs. In their loss function, the first item is the prediction loss which captures short-term trends by calculating errors between the dynamics predicted by the deep modules and real case numbers, the second item is the epidemiological loss which captures long-term trends by calculating the errors between the disease dynamics simulated with the SIR model and real case numbers. The causal-based graph neural network (CausalGNN) model proposed by Wang et al. is another framework that constrains the dynamic attention-based GNN module with an epidemiological model (i.e., the SIRD model) [67]. Similar to [66], the CausalGNN model generates two types of outputs: data-driven predictions of case numbers; and the epidemiological parameters in the SIRD model, which are used to simulate predictions of case numbers. The loss function in this model consists of two $l_{1}$-norm items, which are errors between the ground truth of the case number and the two types of predictions. Unlike the model developed by Gao et al. [66], the CausalGNN model also feeds the simulations obtained from the SIRD model together with the input features into the data-driven model to generate the model outputs. Wang et al. [68, 69] proposed an epidemic prediction framework, named deep learning based epidemic forecasting with synthetic information (DEFSI), to conduct short-term and high-resolution ILI incidence prediction. The novelty of DEFSI is that it generates fine-scale ILI incidence data from an agent-based simulator (EpiFast) of an SEIR model, whose transmission parameters are estimated from the surveillance data. The obtained synthetic data are used to train a two-branch LSTM model to capture the within-season and between-season temporal dependencies of the incidence trends, and the model outputs are merged to generate final predictions.

One of the most significant features of the dynamics of diseases such as influenza and malaria, which are triggered or influenced by the climate, is their obvious correlation with trends in climate conditions. For example, Jeffrey and Melvin reanalyzed data obtained in laboratory experiments on guinea pigs [100] and demonstrated that absolute humidity (AH) greatly affected the influenza virus transmission (IVT) and influenza virus survival (IVS) in temperate regions [101]. On this basis, Shaman et al. further examined these findings at the human population level and proposed an AH-forced SIRS model, which quantitatively defines the relationship between AH and basic reproduction number, to simulate the seasonal patterns of influenza-related deaths [94]. In their subsequent studies, Shaman et al. used an AH-driven SIRS model that has a similar structure to their SIRS model [94] to incorporate AH data into their data-assimilation framework and generate influenza predictions for cities in the USA [11, 51]. In regions with non-temperate climate patterns, such as subtropical and tropical regions, the relationship between AH and influenza can be quite different. In contrast to the above-mentioned studies, which have used compartmental models to incorporate risk-related factors, Tamerius et al. collected data from numerous sites across a global range and explored the relationship between several climate variables (i.e., temperature, solar radiation, specific humidity, and precipitation) and the influenza peaks by rank order analysis. They also established univariate and multivariate logistic regression models to assess this relationship [102]. Venna et al. used a symbolic time-series approach to model the nonlinear relationship between climate variables (i.e., precipitation, temperature, and sun exposure) and dynamics of influenza [30].

Vector-borne diseases, such as malaria, are greatly affected by climate conditions, as these influence the biological features of vectors and pathogens, such as the survival rate of mosquitoes and the incubation period of Plasmodium. The quantitative relationship between temperature or rainfall and VCAP was defined by Ceccato et al. [103]. Shi et al. utilized temperature and rainfall data to estimate some epidemiological parameters in their model, which incorporates VCAP, to assist with the prediction of malaria dynamics which show obvious seasonal fluctuations [61]. Zhang et al. used temperature and rainfall and other disease-related data to form the feature vector of locations to learn their impact on intra-regional transmission risk with the Poisson regression model [22].

Human behavior plays an important role in infectious disease transmission [104], and human mobility is a crucial factor affecting the range and distribution of disease risk [105]. Several studies have examined the effects of different human mobility patterns on disease transmission. For instance, Balcan et al. studied two kinds of mobility patterns—short-range commuting flows and long-range flights—and explored their effects on the spatiotemporal dynamics of infectious disease by integrating a mobility network into the GLEaM model [95]. Pei et al. used commuter data from the US census survey²²2https://www.census.gov/topics/employment/commuting.html. Cannot be accessed on March 25, 2023 to support the prediction of influenza spread [53]. In their model, recurrent commuters and random visitors are considered to compose the connections between locations in the meta-population compartmental model. Kapoor et al. utilized human mobility data from two Google mobility datasets (the Google COVID-19 Aggregated Mobility Research Dataset and Google Community Mobility Reports³³3https://www.google.com/covid19/mobility/. Accessed March 25, 2023) to construct a human mobility network at the daily scale [32] for studying the COVID-19 pandemic. In their network, nodes represent locations, and spatial and temporal edges denote human movement between neighborhoods and connections between historical days and the current day, respectively. Cui et al. used visitor counts in the census block group provided by SafeGraph⁴⁴4https://docs.safegraph.com/docs/places#section-patterns. Accessed March 25, 2023 to model human mobility [40]. Unlike the method of Kapoor et al. [32], which directly constructs the edges of a human mobility network from movement data, the method of Cui et al. [40] treats visitor counts and other disease-related factors as features of each location and uses the attention-based mechanism to learn the structure of a disease transmission network.

In addition to physical movement, the online behavior of individuals can partly reflect the patterns of disease risk. For example, some researchers have assumed that there is a correlation between changes in ILI levels and online search activities. For example, Ginsberg et al. developed Google Flu Trends, which uses the proportion of ILI-related search queries to overall search queries from Google as an explanatory variable for predicting ILI physician visits (the outcome) by fitting a simple linear model [106]. They focused on designing an automated method to identify the ILI-related search queries from a huge number of search records, which contribute most to the model when fitting with ILI physician-visit data. However, maintenance of the website for generating Google Flu Trends stopped in August 2015⁵⁵5https://ai.googleblog.com/2015/08/the-next-chapter-for-flu-trends.html. Accessed March 25, 2023. Nevertheless, the Google Extended Trends (GET) application programming interface remains accessible and provides the statistics of online search trends at various temporal and geographical granularities that researchers can use to train their models [107]. Yang et al. utilized GET data to devise the autoregression with Google search data (ARGO) method for influenza epidemic estimation [107]. They obtained the search terms that are strongly correlated with the ILI from Google Correlate⁶⁶6www.google.com/trends/correlate/. Cannot be accessed on March 25, 2023 (which stopped providing data after March 28, 2015) and search trends from Google Trends⁷⁷7https://trends.google.com/trends/. Accessed March 25, 2023, and then used these data as input for an ARX model. Kandula et al. also used Web-based search activity data from GET to nowcast ILI dynamics at subregional geographic scales (i.e., state-level) [108]. They formalized strongly correlated query terms for a specific region in a period as explanatory variables and used them to train an ARIMA model to predict the response variable (i.e., ILI observations). Then, they treated the forecast of the ARIMA model as an additional explanatory variable and used it with the original explanatory variables to train a random forest model for making final predictions of the response variable. Similarly, search query data from Baidu⁸⁸8https://index.baidu.com/. Accessed March 25, 2023, one of the biggest search engines in China, were used by Yuan et al. to monitor influenza case counts in China [109]. They also used a multiple regression model to predict case counts based on the previous case count and a composite index of searches.

Other important online behaviors of individuals that reflect the health status and contact patterns of humans are social media posts and virtual interactions. Although these data do not directly quantitatively reveal the extent and patterns of disease transmission, they can be mined to obtain some useful information. For example, Achrekar et al. [110] calculated Pearson correlation coefficients and fitted a regression model, finding a linear correlation between the number of Twitter users who tweeted influenza-related posts and the percentage of ILI physician visits. Based on this evidence, they proposed the social network enabled flu trends (SNEFT) framework, which is based on an ARX model, to predict ILI case numbers from ILI physician visits and collected Twitter data. Similarly, Zhang et al. counted ILI-related tweets to estimate the number of influenza infections in a given week and in given locations, but without filtering retweets and successive posts [55].

However, in [110, 55], only the number of users who post influenza-related information or the number of tweets relevant to the influenza keywords were used to estimate the number of influenza infections; other features of textual information were not utilized. In contrast, Volkova et al. extracted detailed linguistic features and communication patterns as latent embeddings from social media posts of a Twitter dataset and fed these data, together with ILI data, into a joint neural network model based on LSTM modules to predict ILI dynamics [31]. Hua et al. used a Bayesian model to identify individual health states (i.e., healthy, exposed, and infectious) from Twitter posts, aggregated these states at the population level, and then employed these states to inform the epidemiological parameters of a simulation model [63].

In addition to data that describe disease severity and transmission in a whole population (i.e., incidence rates or case numbers), there are abundant individual-patient data that can provide insights into disease severity and medical resource utilization. Gao et al. extracted daily data on hospitalizations, intensive care unit stays, and frequency of diagnosis codes from IQVIA’s medical claim data⁹⁹9https://www.iqvia.com/solutions/real-world-evidence/real-world-data-and-insights. Accessed March 25, 2023 [66] and combined these data with the number of COVID-19 cases (active, confirmed, and death cases) provided by Johns Hopkins University¹⁰¹⁰10https://github.com/CSSEGISandData/COVID-19. Accessed March 25, 2023 to serve as the dynamic features of nodes in a graph neural network to facilitate the characterization of spatiotemporal patterns. In [38], Gao et al. extracted the statistics of disease-related features from the same IQVIA claim dataset to facilitate the training and prediction of a downstream deep learning model.

Naturally, current approaches assume that disease dynamics are related to historical dynamics, and this assumption is obviously consistent with the course of disease transmission. Some models assume that the dynamics at the current time step develop from the dynamics at the immediately previous time step, such as the ODEs of the SIR model and other similar compartmental models [71, 73], or that they are related to the dynamics of multiple time steps ahead, due to the existence of disease incubation, such as autoregressive-type models [19, 20], and many deep learning models based on RNN modules [31, 30, 34, 36, 37, 41, 38], CNN modules [35, 40], or GNN modules [33]. Specifically, unlike traditional RNN modules, which can naturally consider multiple time steps, some CNN-based deep learning models capture the temporal patterns at different time scales by conducting dilated convolution across time [35] or encode temporal patterns in various time windows by using kernels of different sizes [40]. There is also a GNN-based model [33] that uses the GFT to capture the inter-series temporal dependency in the spectral domain, and uses the DFT to learn the intra-series in the frequency domain.

Similarly, the spatial pattern is also a necessary characteristic in the description of disease propagation. Some studies have focused on the modeling of disease transmission in a specific location, such as compartmental models at the population level [71, 73] and machine learning models for single locations [31, 30, 27], whereas many other studies have predicted disease dynamics at multiples locations simultaneously. Some of these studies have assumed that spatial patterns can be reflected by the statistical correlation of data between different locations; that is, locations with similar feature conditions have similar disease severities. For instance, the GP regression model [26] uses designed kernels to calculate the covariance matrix, which reflects spatial similarity. Other studies have modeled disease transmission between locations by assuming that a disease transmission network exists, such as compartmental models at the meta-population level [86] and the individual level [89], a multivariate regression model [22, 23], and deep learning models based on a CNN module [34] or a GNN module [32, 35, 36, 38]. The nodes of disease transmission networks typically consist of locations or individuals of interest, whereas edges consist of the interactions or relationships between locations or individuals. Although network-based methods aim to achieve the same goal—i.e., generate epidemic predictions for multiple locations by considering their relationships—they have different focuses on adopted network hypotheses, based on which they can be distinguished. However, although such methods assume that a disease transmission network exists, the interactions in a transmission network cannot be observed directly in the real world, which means there is no ground truth for a disease transmission network. Therefore, some ODE-based models construct a transmission network based on physical contact patterns, and some machine learning models infer network structure from data during model optimization. In doing so, various hypotheses are made regarding network properties to approximate the real situation. In the following, we introduce three types of network structure hypotheses that have often been used in previous studies.

The first common type of setting of network properties is whether edges in a network are undirected or directed, where an undirected edge and a directed edge denote a symmetrical and an asymmetrical correlation or dependency between different locations, respectively. In studies that have assumed that a network is undirected (e.g., [38]), some symmetrical similarity measurements of features have been considered as the quantitative relations between different locations. This implies that the disease dynamic of locations with similar conditions can be used to improve the prediction in a target location. In other studies, a network has been assumed to be directed [22, 23, 34, 32, 35, 36] which implies that there are asymmetrical interactions between different locations. Some of these studies have incorporated prior knowledge and data on the drivers of disease transmission, such as human mobility, which causes cross-regional disease transmission, to assist with network inference[22, 32]. Other studies have directly inferred a transmission network from disease dynamic data [23, 34, 24, 35, 36].

The second common type of setting of network properties is whether the network structure is static or dynamic. Some studies have used a fixed network to predict disease risk [34, 23, 22], which means that after a model is trained on historical data, an inferred transmission network remains constant during the whole prediction. However, in the real world, the structure of a disease transmission network changes over time, due to changes in risk-related factors (e.g., climate, socio-economic indicators, and human behaviors) and interactions between these factors. Therefore, many studies have tried to characterize dynamic disease transmission networks [32, 62].

In addition to basic network properties, more advanced network characteristics such as the power law of degree distribution are ubiquitous in real-world networks, such as the World Wide Web, worldwide airline networks, and interurban commuting networks [111, 112, 113]. Some studies have used this power law to constrain network structure learning for disease transmission. Specifically, Wang et al. [58] adopted a metapopulation SIR model to model disease propagation between different locations and thereby construct an infection network. That is, by assuming that the node degrees of the disease transmission network followed the power-law distribution, they formalized the prior probability of the whole network and incorporated it into their Bayesian framework of network inference. However, in recent studies that have inferred disease transmission networks, this aspect has not been fully explored.

The first category of model is based on probabilistic modeling, which takes uncertainty into account by incorporating the probabilistic distribution of parameters or functions. Usually, these models first assign a prior distribution for targeted variables, and then use Bayes’ theorem to calculate its posterior distribution. This category can be subdivided into two classes. The first class assumes that model parameters follow a probabilistic distribution. For instance, the empirical Bayes framework proposed by Brooks et al. to predict ILI trends uses historical data to estimate the prior distribution of model parameters and produces the posterior distribution of epidemic curves [25]. Other Bayesian inference methods, including the Kalman filtering method and its variants [92], and PFs [93], have also been used in disease dynamic prediction [11, 54, 51, 52, 55, 53]. Rather than assuming that the model parameters follow a probability distribution, stochastic process-based models (e.g., the GP model) define the probability distribution over functions [43]. Some representative studies have used the GP model to predict disease dynamics and provide the uncertainty of results. For instance, Senanayake et al. used the Gaussian process regression model with designed spatiotemporal kernels to model the spatiotemporal characteristics of influenza transmission and used variational inference to optimize the model to adapt to a large dataset [26]; Zimmer and Yaesoubi trained an independent GP model for each location to capture the correlation between historical data in the previous season and predict the dynamics $t$ weeks ahead based on available data for the current year [27]. Furthermore, some studies have used deep learning models to perform stochastic processes; this type of model is called neural processes (NPs) [118]. Its extensions, such as the functional neural process (FNP) [119] and the recurrent neural process (RNP) [120, 121], have also been developed to capture complex dependency. For instance, the EPIFNP is a neural process model that incorporates RNP and FNP modules for disease dynamics prediction [37]. It uses a probabilistic deep-sequence model to encode the sequence data of influenza outbreaks as latent variables, and uses a stochastic correlation graph to learn the correlations between these variables. Finally, it uses an MLP module to parameterize the predictive distribution of model outcomes.

Another way is to use ensemble approaches, which collect a set of predicted values from multiple trained models and calculate their probabilistic distributions to enhance the robustness of predictions. For instance, the COVID-19 Forecast Hub ¹²¹²12https://covid19forecasthub.org. Accessed March 25, 2023 ensembles the forecasts from different models from various institutions to generate US COVID-19 death data [122]; and [123] evaluates the uncertainty from data by using bootstrapping, which is a data sampling technique, to re-sample an entire training dataset to multiple subsets and then uses the subsets to train a set of models. Based on the trained models, multiple prediction results are produced and can be used to calculate the confidence interval.

The most common methods used to evaluate the accuracy of point estimation include the root-mean-square error (RMSE), as shown in Eq.(13); the mean absolute error (MAE), as shown in Eq.(14); the mean absolute percentage error (MAPE), as shown in Eq.(15); and the root mean squared percent error (RMSPE), as shown in Eq.(16). These indicators all calculate the deviation of predicted values $Y^{*}$ from ground truth $Y$, where $Y_{t}^{*}$ and $Y_{t}$ denote the predicted value and the ground truth at time step $t$, respectively. The Pearson correlation coefficient (CORR), as shown in Eq.(17), is used to evaluate the correlation between predicted trend $Y^{*}$ and real trend $Y$; where $\bar{Y^{*}}$ and $\bar{Y}$ denote the mean value of the predicted trend and the real trend, respectively, during the time slot from $1$ to $T$.

$$RMSE=\sqrt{\frac{1}{T}\sum_{t=1}^{T}\left(Y_{t}-Y^{*}_{t}\right)^{2}}$$

(13)

$$MAE=\frac{1}{T}\sum_{t=1}^{T}|Y_{t}-Y^{*}_{t}|$$

(14)

$$MAPE=\left(\max_{t}\frac{|Y_{t}-Y^{*}_{t}|}{Y_{t}}\right)*100$$

(15)

$$RMSPE=\sqrt{\frac{1}{T}\sum^{T}_{t=1}\left(\frac{Y_{t}-Y^{*}_{t}}{Y_{t}}\right)^{2}}$$

(16)

$$CORR=\frac{\left.\sum_{t=1}^{T}\left(Y_{t}-\bar{Y}\right)\left(Y^{*}_{t}-\bar{Y^{*}}\right)\right)}{\sqrt{\sum_{t=1}^{T}\left(Y_{t}-\bar{Y}\right)^{2}}\sqrt{\sum_{t=1}^{T}\left(Y^{*}_{t}-\bar{Y^{*}}\right)^{2}}}$$

(17)

Some indicators are widely used for interval estimation [127]. For instance, prediction interval coverage, denoted as $k_{M}(c)$, calculates the percentage of observed values falling into the $c$ (i.e., 50% or 95%) confidence interval of predicted distributions of $M$ [122]; calibration score is the integral of $\|k_{M}(c)-c\|$ over $c$ from $0,\cdots,1$, as shown in Eq.(18) [37], and a calibration plot shows the relationship between $c$ and $k_{M}(c)$ [37]. Another indicator is the logarithmic score, also called the log-score, which is used in the CDC’s influenza prediction Challenge in the US. It is calculated as follows: given the predicted distribution of the outcome, first calculate the sum of probability of bins within a given interval around the true value, and then determine the natural logarithm of the calculated sum, which represents the final score [128, 27].

	$\displaystyle CS(M)$	$\displaystyle=\int_{0}^{1}\|k_{M}(c)-c\|$		(18)
		$\displaystyle\approx 0.01\sum_{c\in\{0,0.01,\cdots,1\}}\|k_{M}(c)-c\|$