Inference Description
Methods for fitting model to data
flepiMoP can be used to conduct forward simulations of a model with user-defined parameter values, or, it can be used to iteratively run a model with some unknown parameters, compare the model output to ground truth data, and find parameter values that optimize the fit of the model to data (i.e., conduct model "inference"). We have developed a custom model inference method that is based on standard Markov Chain Monte Carlo (MCMC)-based approaches to Bayesian inference for dynamic models, but is adapted to deal with some of the particular challenges of large-scale epidemic models, including i) long times and high computational resources required to simulate single model runs, ii) multiple subpopulations with location-specific parameters but inter-location transmission, iii) a high-dimensional parameter space, iv) the need to produce real-time epidemic projections, and v) the availability of parallel computing resources.
Notation
Background
This section can be skipped by those familiar with Markov Chain Monte Carlo approaches to Bayesian inference.
Bayesian inference
Our model fitting framework is based on the principles of Bayesian inference. Instead of estimating a single "best-fit" value of the unknown model parameters, our goal is to evaluate the consistency of every possible parameter value with the observed data, or in other words, to construct a distribution that describes the probability that a parameter has a certain value given the observations. This output is referred to as the posterior probability. This framework assumes that the model structure accurately describes the underlying generative process which created the data, but that the underlying parameters are unknown and that there can be some error in the observation of the data.
The likelihood function can be defined for a model/data combination based on an understanding of both a) the distribution of model outcomes for a given set of input parameters (if output is stochastic), and b) the nature of the measurement error in observing the data (if relevant).
For complex models with many parameters like those used to simulate epidemic spread, it is generally impossible to construct the full posterior distribution either analytically or numerically. Instead, we rely on a class of methods called "Markov Chain Monte Carlo" (MCMC) that allows us to draw a random sample of parameters from the posterior distribution. Ideally, the statistics of the parameters drawn from this sample should be an unbiased estimate of those from the complete posterior.
Markov Chain Monte Carlo methods
In many Bayesian inference problems that arise in scientific model fitting, it is impossible to directly evaluate the full posterior distribution, since there are many parameters to be inferred (high dimensionality) and it is computationally costly to evaluate the model at any individual parameter set. Instead, it is common to employ Markov Chain Monte Carlo (MCMC) methods, which provide a way to iteratively construct a sequence of values that when taken together represent a sample from a desired probability distribution. In the limit of infinitely long sequences ("chains") of values, these methods are mathematically proven to converge to an unbiased sample from the distribution. There are many different MCMC algorithms, but each of them relies on some type of rule for generating a new "sampled" parameter set from an existing one. Our parameter inference method is based on the popular Metropolis-Hastings algorithm. Briefly, at every step of this iterative algorithm, a new set of parameters is jointly proposed, the model is evaluated at that proposed set, the value of the posterior (e.g., likelihood and prior) is evaluated at the proposed set, and if the posterior is improved compared to the previous step, the proposed parameters are "accepted" and become the next entry in the sequence, whereas if the value of the posterior is decreased, the proposed parameters are only accepted with some probability and otherwise rejected (in which case the next entry in the sequences becomes a repeat of the previous parameter set).
The full algorithm for Metropolis-Hastings Markov Chain Monte Carlo is:
Evaluate the likelihood and prior at the proposed parameter set
Inference algorithm
Likelihood
In our algorithm, model fitting involves comparing timeseries of variables produced by the model (either transmission model state variables or observable outcomes constructed from those variables) to timeseries of observed "ground truth" data with the same time points. For timeseries data that arises from a deterministic, dynamic model, then the overall likelihood can be calculated as the product of the likelihood of the model output at each timepoint (since we assume the data at each timepoint was measured independently). If there are multiple observed datastreams corresponding to multiple model outputs (e.g., cases and deaths),
And the overall likelihood taking into account all subpopulations, is the product of the individual likelihoods
Fitting algorithm
The method we use for estimating model parameters is based on the Metropolis-Hastings algorithm, which is a class of Markov Chain Monte Carlo (MCMC) methods for obtaining samples from a posterior probability distribution. We developed a custom version of this algorithm to deal with some of the particular mathematical properties and computational challenges of fitting large disease transmission models.
There are to major unique features of our adapted algorithm:
Parallelization – Generally MCMC methods starting from a single initial parameter set and generating an extremely long sequence of parameter samples such that the Markov process is acceptably close to a stationary state where it represents an unbiased sample from the posterior. Instead, we simulate multiple shorter chains in parallel, starting from different initial conditions, and pool the results. Due to the computational time required to simulate the epidemic model, and the timescale on which forecasts of epidemic trajectories are often needed (~weeks), it is not possible to sequentially simulate the model millions of times. However, modern supercomputers allow massively parallel computation. The hope of this algorithm is that the parallel chains sample different subspaces of the posterior distribution, and together represent a reasonable sample from the full posterior. To maximize the chance of at least local stationarity of these subsamples, we pool only the final values of each of the parallel chains.
Multi-level – Our pipeline, and the fitting algorithm in particular, were designed to be able to simulate disease dynamics in a collection of linked subpopulations. This population structure creates challenges for model fitting. We want the model to be able to recreate the dynamics in each subpopulation, not just the overall summed dynamics. Each subpopulation has unique parameters, but due to the coupling between them (), the model outcomes in one subpopulation also depend on the parameter values in other subpopulations. For some subpopulations (), this coupling may effectively be weak and have little impact on dynamics, but for others (), spillover from another closely connected subpopulation may be the primary driver of the local dynamics. Thus, the model cannot be separately fit to each subpopulation, but must consider the combined likelihood. However, such an algorithm may be very slow to find parameters that optimize fits in all locations simultaneously, and may be predominantly drawn to fitting to the largest/most connected subpopulations. The avoid these issues, we simultaneously generate two communicating parameter chains: a "chimeric" chain that allows the parameters for each subpopulation to evolve quasi-independently based on local fit quality, and a "global" chain that evolves only based on the overall fit quality (for all subpopulations combined).
Note that while the traditional Metropolis-Hastings algorithm for MCMC will provably converge to a stationary distribution where the sequence of parameters represents a sample from the posterior distribution, no such claim has been mathematically proven for our method.
Generate initial state
Make "global" decision about proposed parameters
Update the recorded subpopulation-specific likelihood values (chimeric and global) with the likelihoods calculated using the proposed parameters
Else
: REJECT the proposed parameters for the global chain and make subpopulation-specific decisions for the chimeric chainMake "chimeric" decision:
Else
: REJECT the proposed parameters for the chimeric parameter chain for this locationEnd if
End making chimeric decisions
End if
End making global decision
There are a few important notes/limitations about our method currently:
All parameters to be fit must be location-specific. There is currently no way to fit a parameter that has the identical value across all locations
The pipeline currently does not allow for fitting of the basic parameters of the compartmental epidemic model. Instead, these must be fixed, and the value of location-specific "interventions" acting to increase/reduce these parameters can be fit. All parameters related to the observational/outcomes model can be fit, as well as "interventions" acting to increase or reduce them.
At no point is the parameter fitting optimizing the fit of the summed total population data to total population model predictions. The "overall" likelihood function used to make "global" parameter acceptance decisions is the product of the individual subpopulations likelihoods (which are based on comparing location-specific data to location-specific model output), which is not equivalent to likelihood for the total population. For example, if overestimates of the model in some subpopulations were exactly balanced by underestimates in others, the total population estimate could be very accurate and the total population likelihood high, but the overall likelihood we use here would still be low.
It is currently not possible to infer parameters of the measurement process encoded in the likelihood function. For example, if the likelihood is chosen to be a normal distribution, which implies an assumption that the observed data is generated from the underlying truth according to a normal distribution with mean zero, then the standard deviation must be specified, and cannot be inferred along with the other model parameters.
There is an option to use a slightly different version of our algorithm, in which globally accepted parameter values are not pushed back into the chimeric likelihood, but the chimeric likelihood is instead allowed to continue to evolve independently. In this variation, the chimeric acceptance decision is always made, not only if a global rejection happens.
Hierarchical parameters
The baseline likelihood function used in the fitting algorithm described above allows for parameter values to differ arbitrarily between different subpopulations. However, it may be desired to instead impose constraints on the best-fit parameters, such that subpopulations that are similar in some way, or belong to some pre-defined group, have parameters that are close to one another. Formally, this is typically done with group-level or hierarchical models that fit meta-parameters from which individual subpopulation parameters are assumed to draw. Here, we instead impose this group-level structure by adding an additional term to the likelihood that describes the probability that the set of parameters proposed for a group of subpopulations comes from a normal distribution. This term of the likelihood will be larger when the variance of this parameter set is smaller. Formally
Last updated