Limits of epidemic prediction using SIR models

12/13/2021
by   Omar Melikechi, et al.
Duke University
0

The Susceptible-Infectious-Recovered (SIR) equations and their extensions comprise a commonly utilized set of models for understanding and predicting the course of an epidemic. In practice, it is of substantial interest to estimate the model parameters based on noisy observations early in the outbreak, well before the epidemic reaches its peak. This allows prediction of the subsequent course of the epidemic and design of appropriate interventions. However, accurately inferring SIR model parameters in such scenarios is problematic. This article provides novel, theoretical insight on this issue of practical identifiability of the SIR model. Our theory provides new understanding of the inferential limits of routinely used epidemic models and provides a valuable addition to current simulate-and-check methods. We illustrate some practical implications through application to a real-world epidemic data set.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 16

11/09/2021

Probabilistic predictions of SIS epidemics on networks based on population-level observations

We predict the future course of ongoing susceptible-infected-susceptible...
06/11/2020

The Limits to Learning an SIR Process: Granular Forecasting for Covid-19

A multitude of forecasting efforts have arisen to support management of ...
09/23/2021

On the parametrization of epidemiologic models – lessons from modelling COVID-19 epidemic

A plethora of prediction models of SARS-CoV-2 pandemic were proposed in ...
05/13/2021

Calibration of a SEIR-SEI epidemic model to describe the Zika virus outbreak in Brazil

Multiple instances of Zika virus epidemic have been reported around the ...
08/05/2020

Modelling collective decision-making during epidemics

The outcome of an epidemic outbreak can be critically shaped by the coll...
02/01/2020

Online Bayesian phylodynamic inference in BEAST with application to epidemic reconstruction

Reconstructing pathogen dynamics from genetic data as they become availa...
05/12/2020

Remarks on a data-driven model for predicting the course of COVID-19 epidemic

Norden E. Huang, Fangli Qiao and Ka Kit Tung presented a data-driven mod...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The Susceptible-Infectious-Recovered (SIR) model, first introduced in the early twentieth century, is a mathematical model describing the spread of a novel pathogen through a population [15, 23, 22, 24]

. This model is governed by the ordinary differential equations

(1)

According to this model, the population is divided into three groups or “compartments,” each of which represents a proportion of the population. The compartment consists of the proportion of individuals who have never been infected with the pathogen and are therefore susceptible to infection. The compartment consists of the proportion who are currently infected and contagious. The third compartment consists of the proportion of individuals who have either recovered from the pathogen and are immune or have died, and are therefore removed from the population of susceptible individuals. Since the SIR model assumes all recovered individuals are permanently immune to the pathogen, the value of can be obtained from and via the identity .

During the last century the SIR equations have been modified and extended to model a diverse range of epidemics including Ebola, cholera, H1N1, tuberculosis, HIV/AIDS, influenza, malaria, Dengue fever, Zika, and most recently SARS-CoV-2 [5, 10, 12, 16, 17, 19, 21, 30]. In many of these examples, additional terms are added to account for pathogen specific characteristics of transmission. Additional compartments may also be added to model different subpopulations. One such example is the SEIR model, which includes a subpopulation of exposed (E) but non-infectious cases [25]. Collectively, the SIR model and its extensions and variations provide epidemiologists with a vast array of interpretable and highly expressive models to understand and predict the behavior of outbreaks. However, incorporating too many features can have subtle but important drawbacks including limited or unreliable inference of model parameters early in an epidemic. The main contribution of this article is insight on the inferential limits of epidemics (as captured by estimated parameters) which can be obtained from noisy, real-time observations of an outbreak.

Our work is motivated by the application of SIR and related compartmental models to real-time analysis of epidemics of human disease such as the 2014 Ebola outbreak in West Africa and the ongoing SARS-CoV-2 (“coronavirus”) pandemic. In such outbreaks, the initial aim of the public health response is to extinguish the epidemic while the number of infected individuals is still small, or at least to significantly slow the rate of infection to allow time for the pathogen to be better understood and effective therapeutics or vaccines to be developed. The stay-at-home orders instituted by many countries due to SARS-CoV-2 are one recent example which has had profound global economic impacts. As such, mathematical models employed in the real-time analysis of epidemics must provide accurate inferences about properties of the epidemic – encapsulated by the parameters of the model – early in the epidemic, when only a small fraction of the population has been infected. Hereafter, we refer to the estimation of unknown model parameters from observations as the inverse problem.

As noted in the recent review by Hamelin et al, many disease models proposed in the literature follow a similar structure: (1) a model is proposed, (2) a subset of model parameters are inferred from the literature, and (3) the remaining parameters are fit from data using least squares or maximum likelihood estimation [14]. In order for these parameter estimates to be reliable the parameters must be statistically identifiable, ruling out settings in which multiple parameter values are equally consistent with observed data. Such issues were first considered in the context of compartmental models by Bellman and Aström in 1970 [3]. Specific details relevant to the SIR model may be found in [14]. Here we provide a brief overview of the well-posedness of the inverse problem.

A model is structurally identifiable when there is a single value of the parameters consistent with noise-free data observations. A comprehensive review of analytic methods for assessing structural identifiability is given in [9]; alternatively, software packages such as DAISY can be used [4]. There are many examples in the literature [7, 8, 11, 12, 20, 27, 28, 29]. In particular, structural identifiability of the SIR parameters, and , is well understood with strong theoretical support. See [14] for specific cases based on different observations of the compartments.

In practice, data observed during an epidemic tend to be very noisy, so we are far from the idealized noise-free case. Practical identifiability is the ability to discern different parameter values based on noisy observations. Despite considerable recent attention [2, 1, 9, 26], far less is known about practical identifiability. Present theoretical methods rely on sensitivity analysis and the computation of the Fisher information matrix, which is analytically intractable in the SIR model and its extensions. Instead, it is common to see Monte Carlo methods employed, wherein the model is simulated for a set value of the parameters, noise is added to the simulated observations, and a fitting procedure is conducted on the noisy data [9, 14, 17, 28]. The fidelity of parameter estimates relative to the known values is summarized using the average relative estimation error, which is then plotted as a function of the noise intensity.

Interestingly, the lack of practical identifiability manifests in a remarkably similar manner across multiple, different model formulations even in cases where the parameters are known to be structurally identifiable. As the magnitude of the noise is increased, Monte Carlo parameter estimates concentrate along a curve stretched throughout parameter space indicating a functional relationship between model parameters [6, 12, 20, 27, 28]. Importantly, there are often great disparities in parameter values along this curve and hence huge uncertainty in the parameters. See Figure 1 for a representative example in the specific case considered herein.

The goal of this article is to provide theoretical tools for understanding practical identifiability in the context of the SIR model. We propose a formulation based on realistic observations early in an outbreak. Then, using linearizations similar to those of [25], we construct analytically tractable approximations to the SIR dynamics from which theoretical guarantees of the performance of the inverse problem are developed. We begin by introducing the model under consideration, reemphasizing ideas discussed previously to provide overt examples of the challenges of practical identifiability.

2 Statistical model

The data available to infer the parameters of an SIR model are usually noisy, biased measurements of the rate of change in the size of the susceptible compartment, discretized to unit time intervals . For simplicity, we take the time unit to be one day. Here, represents the total population size in the jurisdiction under study. The quantity is the number of newly infected individuals between day and day . Data on daily confirmed cases, hospitalizations, or deaths are all examples of observable data that depend on the underlying value of . Specifically, all are discrete convolutions of of the form , where

is the probability that an infected person goes on to be diagnosed, hospitalized, or die, and

is the conditional probability that a person tests positive, is hospitalized, or dies days after becoming infected given that the corresponding outcome will eventually occur. It is likely that the parameters and, to a lesser extent, change over the course of an epidemic. However, changing values of these parameters can only make inference more difficult, and since our main focus is on studying limitations of inference, as a starting point we assume that and are fixed and known.

While the inverse problem with known initial conditions but unknown parameters is well-posed when even a partial trajectory of is observed, in reality we observe corrupted with noise, and we always have to work with finitely many discrete-time observations. In epidemic modeling, unlike some other inverse problems, we do not even have control of the sampling rate and are generally stuck with at best daily monitoring data. To simplify exposition, we focus on a simple but flexible noise model in which the observed data

are realizations of a random variable satisfying

for some known . In this case, and for . While our results apply to many noise models, to fix ideas we begin with Gaussian noise

(2)

Initially, suppose that the variances

are known. A simple procedure for solving the inverse problem from data is maximum likelihood. The gradient of the log-likelihood can be obtained by numerically solving an extended ODE system [13] which allows for easy fitting via gradient-based optimization methods. It can be shown that, even when the trajectory is observed only at discrete time intervals and the peak of infections has not yet occurred, the maximum likelihood estimator (MLE) exists and is unique, and so the model is structurally identifiable [14].

Problems become apparent when one seeks to study uncertainty in the estimated parameters. Figure 1 gives a stark indication of the challenges. We simulate data from an SIR model with parameters and initial conditions for . These parameters were selected to roughly approximate the dynamics of the coronavirus epidemic in New York City prior to the lockdown of March 16, 2020. The trajectories for are shown in the left panel. By , about 1 percent of the population has been infected, and the peak size of the infected compartment occurs around . The right panel of Figure 1 is obtained by repeatedly simulating data from (2) using the trajectory in the left panel, with and . For each replicate simulation, the model is fit by maximum likelihood. The resulting estimates of are shown in Figure 1, which plots against . These are samples from the sampling distribution of the maximum likelihood estimator for these parameters. The estimates exhibit very tight concentration along a line of slope . The variation in observed for these values is large, ranging from to . This high degree of uncertainty occurs despite the fact that we have observed data up through the time when over half the population has been infected.

Figure 1: The trajectory of the SIR model used in the simulation (left). Plots of vs from 1000 realizations from the sampling distribution of their MLE (right).

The linear shape of the plot in Figure 1 suggests a practical identifiability problem in this model. That is, while the MLE exists and is unique, the curvature of the log-likelihood in the neighborhood of the MLE is very small in the direction where lie along a line. We are not the first to notice this phenomenon. Previous works include [9, 14, 17, 28], which experience qualitatively similar issues despite notable differences in the formulation of the likelihood in those settings.

While various empirical studies exist, our main contribution is a theoretical analysis of this phenomenon and the resulting limitations for solving the inverse problem from noisy observations. We take a two-step approach to the analysis. First, we characterize sensitivity of trajectories to perturbations of the parameters , and show that perturbations of in the directions and (equivalently, along the line of slope through ), closely approximate the smallest variation in the trajectory among all perturbations for which . We then give a sharp, computable lower bound on for times prior to the peak infection time. Taken together, these results provide a rigorous explanation for the phenomenon in Figure 1.

In the second part of the analysis, we relate the problem of uncertainty quantification to hypothesis tests of the form

for

. We use the result of the first part of our analysis to approximate the type II error of the test, which in turn allows for both theoretical and empirical analysis of the limits of epidemic prediction using SIR models.

3 Results

3.1 Perturbation bound for SIR trajectories

Informally, the phenomenon in Figure 1 is a manifestation of the fact that very different values of can lead to SIR model trajectories that are very close. To formalize this, let be the -trajectory of the SIR model starting from with parameters . To aid the reader, all relevent notation is summarized in Table 1.

Notation Description
Shorthand for initial conditions with
Shorthand for the parameters of the SIR model
The reproductive number
An important combination of the model parameters appearing in later analysis
Perturbation of such that
Perturbation of in the direction such that
Solution of the SIR equation with initial condition and parameter
Observed data on days 1 through
Likelihood of given observed data
Table 1: Summary of notation used throughout this article.

For , let denote the circle of radius about . That is,

where . We set and assume throughout that ; if not, then the reproductive number is nonpositive and the epidemic does not grow even at time zero. Similarly, we assume . This ensures values of the perturbed parameters are also strictly greater than one. The main result of this subsection is the following proposition.

Proposition 1.

Let denote the Euclidean norm on and let be the time of peak infection. Then for all ,

(3)

Furthermore the infimum is approximately achieved when or .

The derivation of (3) is in the Appendix. Proposition 1 says that for any perturbation of , the distance between the perturbed trajectory and the true trajectory at time is approximately bounded below by the left side of (3). The “” in (3) indicates the bound is subject to error. Specifically, our derivation of Proposition 1 involves two approximations: First, we approximate the SIR model by a differential equation (11) whose solution is given by (12) (see Section A in the Appendix). Second, we use a first-order Taylor expansion to approximate perturbations of resulting from perturbations in parameter space. Despite these approximations, numerical simulations indicate (3) is sharp for a significant amount of time prior to and holds for a wide range of parameter values and population sizes (see Figure 2).

In addition to giving a sharp lower bound, Proposition 1 successfully predicts the directions in parameter space, namely and , corresponding to the most uncertainty about parameters even when data are observed up to the time of peak infection, as in Figure 1. Moreover, we find that the lower bound in (3) approximately holds for the trajectory alone. That is, if and are the trajectories corresponding to and , respectively,

(4)

and the infimum is again achieved when and . This observation will be used for the hypothesis testing in Section 3.2 since our statistical model depends crucially on , which in turn depends only on rather than on and together. The approximation error implicit in the symbol in (3) and (4) is the one quantity we do not have rigorous control over.

Figure 2: Distance between perturbed and true trajectories for different parameter values and population sizes. In each graph the horizontal axis is the number of days since the start of the epidemic and the vertical axis is the distance between a perturbed trajectory and the true trajectory at time . In each case . The light colored curves correspond to perturbed trajectories, one for each of equally spaced points on the circle . The curves corresponding to the angles and are in black. The lower bound of Proposition 1 is in red, and the peak time is indicated by the vertical blue line. The first column has population size million and the second column population thousand with only one initial infection in each case. The first row has parameters for an of , and the second row has parameters for an of .

3.2 Hypothesis testing for the inverse problem

In this subsection we revisit the inverse problem in light of the perturbation bounds (3) and (4). To give context and motivate the main result of this subsection, namely Proposition 2 and its subsequent discussion, we first give a brief overview of simple hypothesis testing and the Neyman-Pearson Lemma.

Suppose we observe data taking values in a space

and that these data are drawn from an unknown probability distribution belonging to a parametrized family of probability distributions

. Given two parameters and , a natural question is whether the observed data came from or . This is a simple hypothesis test, denoted by

(5)

where and are the null and alternative hypotheses, respectively. Simple here refers to the fact that both and correspond to single values which completely determine the distributions and . The aim is to decide whether to reject in favor of , which is done by choosing a subset of called the rejection region. This choice of completely determines the test: If , then reject in favor of ; if , then do not reject . Type I error occurs when is true but is rejected, and type II error occurs when is false but not rejected. This is quantified111Type I and II error rates are commonly denoted by and , but since is already used as an SIR parameter we adopt the unconventional notation and . as

Ideally one would find a rejection region that simultaneously minimizes type I and type II error rates, but this is generally impossible. Instead, a common statistical paradigm is to fix a significance level and minimize subject to the constraint . For such an , a region is called a most powerful level- rejection region if and for all satisfying . That is, minimizes type II error over all rejection regions with type I error equal to . The Neyman-Pearson Lemma gives the most powerful rejection region in the case of a simple hypothesis test.

Lemma 1.

(Neyman-Pearson) Let denote the likelihood function for data and a parameter , and fix . Then there exists an such that

(6)

is a most powerful level- rejection region for the hypothesis test (5).

is called the likelihood ratio and the decision to reject or not reject based on the rejection region is called the likelihood ratio test. Since the Neyman-Pearson Lemma guarantees the likelihood ratio test is most powerful in our setting, we henceforth consider only the rejection region and set and .

Reject Type I error Success
Do not reject Success Type II error
Table 2: Simple hypothesis test

Returning now to the SIR model, our hypothesis test of interest is

(7)

where, as before, is a perturbation of of size in the direction . The observed data are for any time before the time of peak infection, with each as in (2). For the rest of this paper will denote the type II error rate of (7) for the likelihood ratio test with angle . As such, the likelihood ratio test (6) minimizes thereby providing the most powerful technique for detecting differences of order in the SIR model parameters. We also set where, recall, , and let

denote the standard normal cumulative distribution function. With this notation we now present the main result of this subsection.

Proposition 2.

For any , , and significance level ,

where and . Moreover,

Figure 3: Type II error as a function of perturbation size and noise level. The red and black curves are the first and second approximations of in Proposition 2, respectively. The blue curve is the empirical type II error rate when , which was obtained by performing simulations of the noisy SIR model (2). Similarly, the green and purple curves correspond to the empirical type II error rates for and . In each case the SIR parameters are as in the example from Section 2, namely and with a population size of and initial condition . The time horizon is days into the epidemic, which in this case is days prior to the time of peak infection. The significance level is . As predicted, the type II errors all approach both as perturbation size goes to and as the noise level gets large, and this approach is most rapid when . In each case the noise model is Case 2, . For the simulated blue, green, and purple curves, we used a numerical integrator to obtain the values, while for the red and black curves we used the pre-peak approximation . In the left panel ; in the right panel .

The derivation of Proposition 2 is in the Appendix. The first and second approximations of correspond to the red and black curves in Figure 3, respectively. Comparing these to the empirical type II error rates (the blue, green, and purple curves) we see these approximations are sound. In particular, the last part of Proposition 2 indicates the angles and give rise to the largest type II error rate for the hypothesis test (7) with perturbation size and significance level . To quantify the magnitude of type II error in these cases, we substitute into the second approximation to get

(8)

Note that as the noise level goes to , the sum under the square root goes to infinity and the entire expression goes to . That is, if there is no noise then the type II error rate of the likelihood ratio test will vanish. If there is any noise at all however, the sum is finite and (8) becomes arbitrarily close to as either , the probability of detecting an infected individual, or , the perturbation size, go to . For example, if we set the type I error rate to then as either or go to , the probability of making a type II error will approach . Similarly, type II error will go to as goes to infinity. This limit is unrealistic though since is the variance of observed data and as such should be less than the population size. This leads us to consider two cases for noise.

  1. Noise proportional to population size, i.e. for some .

  2. Noise proportional to new infections, i.e. for some .

In both cases is constant and independent of . Case 2 involves which is not expressible in closed-form. However, we can use Proposition 1 and its derivation, specifically the approximate solution (12), to circumvent this issue by replacing with . As discussed in Section 3.1, this approximation is appropriate early in the epidemic. In Case 1, Equation (8) becomes

In addition to the aforementioned limits, we see in this case that the expression, and hence the type II error, approaches as the population goes to infinity (so that goes to ). In Case 2, Equation (8) becomes

(9)

The above expression does not depend on population size, , nor on the SIR parameters and , while the asymptotic results for , , and still apply. Since Case 1 has noise proportional only to , it implicitly assumes relative noise is larger earlier in the outbreak which may not be realistic. Case 2 avoids this since relative noise will be small whenever the reported number of infected individuals is small, e.g. early in an epidemic. For this reason and its invariance under different model parameters and population sizes, we consider only Case 2 moving forward.

3.3 Empirical analysis: Implications of Proposition 2

Proposition 2 says that given an SIR parameter , the probability of failing to reject the hypothesis when the alternative is true can be very high, especially when the angle of perturbation is or . In this subsection we take a closer look at what this means for epidemic prediction.

Consider for concreteness the familiar setting , , and . Figure 4 shows the total number of infections days past the time of peak infection (left panel) as well as the duration of the epidemic (right panel), which is taken to be the first day past the peak with less than infected individuals. The red line gives these values for while the blue and green curves give these values for and over a range of values.

Figure 4: Consequences of type II error.

Setting or and letting noise be as in Case 2, the first approximation in Proposition 2 can be rearranged to obtain

From this we can compute the consequences of type II error. For example, suppose we are days into an epidemic and wish to test the hypothesis versus as above. Moreover, suppose (perfect diagnostics),

(infection standard deviation of

of new cases), and we set and . Then the above equation gives . Thus, reading off the right panels in Figure 4, we see that under these fairly generous conditions a type II error – which has a chance of occurring – will result in underestimating the total number of infections of an epidemic by over of the total population and the duration of an epidemic by over of the predicted duration. For a type II error in this setting will result in overestimating the total infections by nearly of the total population and the duration by approximately of the predicted one. A similar computation shows when with all other parameters the same, and again from the left panels in Figure 4

we observe significantly different predicted outcomes depending on whether or not the null hypothesis

is rejected.

Figure 5: Practical identifiability of . Each curve is the first approximation of type II error in Proposition 2 for different values of spread evenly across . The dark red curve in each panel corresponds to and . The light red curves are angles closest to and , the yellow are those a bit farther away, the blue still farther, and the purple (which overlap the blue) are those farthest from and , i.e. closest to and . Note the rapid fall off in type II error as angles get farther from and , especially as a function of . This agrees with the empirical observation in Figure 1 that MLE favors parameters lying along a line of slope . In particular, the inverse problem for is practically identifiable.

While our empirical and theoretical results indicate the inverse problem of finding is prone to error, they also show inference of is robust and reliable (see for instance Figures 1 and 5). In particular, since is completely determined by222We choose to focus on because, unlike which depends on the average number of people an infectious person will come in contact with, depends only on the pathogen, not on human social behavior, and therefore tends to be more stable and better approximated in practice. and , knowledge of reduces the inverse problem to finding , the reciprocal of the average number of days an individual is infectious. In this case our new hypothesis test becomes

(10)

for some real number . Furthermore, knowing implies lies on the line of slope with vertical intercept . So by a simple geometric argument (see Figure 6), the above hypothesis test is equivalent to the hypothesis test (7) with and angle if or if . So by (9) the type II error of (10) is

Thus rather than consider the original hypothesis test, one can first infer , then consider the hypothesis test (10) with type II error rate as above.

Figure 6: Going from hypothesis test (10) to (7).

3.4 Empirical analysis: NYC Covid cases, March 2020

In this section we discuss the extension of the theoretical results on parametric non-identifiability to a real world dataset. Consider the Spring 2020 COVID-19 outbreak in New York City. The New York City Health Department keeps a repository of all public COVID-19 data online [18]. Using their daily case data as a proxy for new infections, we directly apply equation (2) to the noisy data. We assume the proportion of positive cases which are detected by testing, , is 15%, and the population size To focus on estimation early in the pandemic, we focus on reported daily cases from February 29, 2020 through March 14, 2020, which approximately represent the first two weeks of the pandemic in New York City.

To connect with the earlier analysis, we are following Case 2 as discussed in section 3.2, in which the noise is proportional to new infections, i.e. for some which is also inferred via maximum likelihood. We thus model daily infections by

from which we obtain the log-likelihood function

Maximizing the above likelihood gives estimates

and a corresponding estimate of . The SIR curve generated by the maximum likelihood estimates of and is shown in Figure 7 with corresponding 95% confidence regions based on the maximum likelihood estimate of

Figure 7: New York City public testing results for COVID-19 from the first known case on February 29, 2020 to March 15, 2020. We have used maximum likelihood estimation to generate an SIR trajectory through the noisy data.

Returning to the testing framework, Figure 8 provides type II error estimates based on the approximation of Equation (9) with significance level , reporting rate , and days of new infection counts. At the MLE , we see the type II error is greater than 80% for all values of such that or with corresponding .

Figure 8: Type II error rate as a function of and at significance level , reporting rate , and days of new infection observations.

Thus, while there may be a large disparity between the true SIR parameters and our maximum likelihood estimates – hence large differences in the estimate of – the hypothesis testing framework has very low power to detect such differences. This result is based on the most difficult to detect perturbations in . However, it provides pessimistic but important lower bounds on the extent to which one can rely on parameter estimates from noisy, early pandemic data.

One final note on the preceding example. The MLEs of and in the previous analysis are quite sensitive to the reporting rate. For reference, Table 3 provides corresponding MLE estimates for , , and as a function of .

0.02 15.59 14.99 1.29 1.04
0.05 4.82 4.22 1.37 1.14
0.1 2.20 1.58 1.36 1.39
0.15 1.44 0.82 1.35 1.7
0.2 1.07 0.46 1.35 2.35
0.25 0.86 0.24 1.36 3.58
Table 3: Maximum likelihood estimates of , , and are shown for different choices of reporting rate . The corresponding estimates of using the MLEs is also provided.

However, the type II error plot in Figure 8 is largely unchanged for the range of in the preceding table. Since is fairly robust to different choices of , our conclusion about the limited power of testing holds true for the range of considered. Therefore, one has limited statistical power to detect large differences in SIR model parameters in the worst case scenario, regardless of the choice of reporting rate.

There is an important distinction to make. Having low power to detect a difference is not equivalent to being unable to tell that there is a difference. Certain parameter values are essentially impossible given natural assumptions about the dynamics of a pandemic. For example, is the average time an infected individual can spread the disease before they are either recovered or removed from the population by quarantine. Extremely large values of and hence small values of are likely unrealistic. Thus, the inclusion of side or prior information on and/or akin to the analysis in Section 3.3 can greatly improve one’s ability to disambiguate different SIR parameters.

4 Discussion

The preceding analysis was based on a simple implementation of the SIR model. Practitioners studying future outbreaks may consider a multitude of modifications to our model construction which result is very different likelihood functions. Thus, we have decided to conclude this article with a short discussion of how one may adapt our techniques to these different settings to better understand issues of practical identifiability with noisy or sparse observations.

To construct an analytically tractable approximation to the type II error, we relied heavily on a linear system, namely (11), that approximates the SIR equations. Such approximations are suitable locally in time and are therefore appropriate when one is focused on the early stages of an outbreak. Importantly, a similar approach can be used to construct analytic approximations to the dynamics of any compartment. Such expressions will depend on unknown SIR parameters, fixed parameters such as population size, and other parameters such as reporting rate. One can use these approximations in different model formulations to construct analytic expressions for type II error and therefore better understand potential limitations of a novel model for pandemic dynamics.

Declaration of interests

The authors declare no competing interests.

Acknowledgments

This work was partially funded by NIH R01-ES028804 from the National Institute of Environmental Health Sciences, and the Duke University DOMath summer program.

References

  • [1] E. Balsa-Canto, A. A. Alonso, and J. R. Banga (2008) Computational procedures for optimal experimental design in biological systems. IET Syst. Biol. 2 (4), pp. 163–172. External Links: Document, ISSN 17518849, Link Cited by: §1.
  • [2] E. Balsa-Canto, A. Alonso, and J. Banga (2009) An iterative identification procedure for dynamic modeling of biochemical networks. BMC Systems Biology 4, pp. 11 – 11. Cited by: §1.
  • [3] R. Bellman and K.J. Åström (1970) On structural identifiability. Mathematical Biosciences 7 (3), pp. 329–339. External Links: ISSN 0025-5564, Document, Link Cited by: §1.
  • [4] G. Bellu, M. P. Saccomani, S. Audoly, and L. D’Angiò (2007-10) DAISY: A new software tool to test global identifiability of biological and physiological systems. Comput. Methods Programs Biomed. 88 (1), pp. 52–61. External Links: Document, ISSN 01692607 Cited by: §1.
  • [5] F. Brauer, C. Castillo-Chavez, and Z. Feng (2019) Mathematical Models in Epidemiology. Texts in Applied Mathematics, Vol. 69, Springer New York, New York, NY. External Links: Document, ISBN 978-1-4939-9826-5, Link Cited by: §1.
  • [6] A. P. Browning, D. J. Warne, K. Burrage, R. E. Baker, and M. J. Simpson (2020) Listen to the noise: identifiability analysis for stochastic differential equation models in systems biology. bioRxiv (October). External Links: Document, ISSN 26928205 Cited by: §1.
  • [7] N. J.B. Brunel (2008) Parameter estimation of ODE’s via nonparametric estimators. Electron. J. Stat. 2 (March 2007), pp. 1242–1267. External Links: Document, ISSN 19357524 Cited by: §1.
  • [8] J. D. Chapman and N. D. Evans (2009-10) The structural identifiability of susceptible-infective-recovered type epidemic models with incomplete immunity and birth targeted vaccination. Biomed. Signal Process. Control 4 (4), pp. 278–284. External Links: Document, ISSN 17468094 Cited by: §1.
  • [9] O. Chis, J. R. Banga, and E. Balsa-Canto (2011) Structural Identifiability of Systems Biology Models: A Critical Comparison of Methods. PLoS One 6 (11), pp. 27755. External Links: Document, Link Cited by: §1, §1, §2.
  • [10] B. J. Coburn, B. G. Wagner, and S. Blower (2009-06) Modeling influenza epidemics and pandemics: Insights into the future of swine flu (H1N1). BMC Med. 7, pp. 30. External Links: Document, ISSN 17417015, Link Cited by: §1.
  • [11] A. C. Daly, D. Gavaghan, J. Cooper, and S. Tavener (2018) Inference-based assessment of parameter identifiability in nonlinear biological models. J. R. Soc. Interface 15 (144). External Links: Document, ISBN 0000000294997, ISSN 17425662 Cited by: §1.
  • [12] M. C. Eisenberg, S. L. Robertson, and J. H. Tien (2013-05) Identifiability and estimation of multiple transmission pathways in cholera and waterborne disease. J. Theor. Biol. 324, pp. 84–102. External Links: Document, ISSN 00225193 Cited by: §1, §1, §1.
  • [13] T. H. Gronwall (1919-07) Note on the Derivatives with Respect to a Parameter of the Solutions of a System of Differential Equations. Ann. Math. 20 (4), pp. 292. External Links: Document, ISSN 0003486X Cited by: §2.
  • [14] F. Hamelin, A. Iggidr, A. Rapaport, G. Sallet, G. Sallet Observability, F. Hamelin, A. Iggidr, A. Rapaport, and G. Sallet (2020) Observability, Identifiability and Epidemiology A survey. Technical report External Links: Link Cited by: §1, §1, §1, §2, §2.
  • [15] W. Kermack and A. G. Mckendrick (1927-08) A contribution to the mathematical theory of epidemics. Proc. R. Soc. London. Ser. A, Contain. Pap. a Math. Phys. Character 115 (772), pp. 700–721. External Links: Document, ISSN 0950-1207, Link Cited by: §1.
  • [16] A. Khaleque and P. Sen (2017-02) An empirical analysis of the Ebola outbreak in West Africa. Sci. Rep. 7. External Links: Document, ISSN 20452322, Link Cited by: §1.
  • [17] C. Lee, Y. Li, and J. Kim (2020-10) The susceptible-unidentified infected-confirmed (SUC) epidemic model for estimating unidentified infected population for COVID-19. Chaos, Solitons and Fractals 139, pp. 110090. External Links: Document, ISSN 09600779, Link Cited by: §1, §1, §2.
  • [18] NYC Health Department (2021) NYC Coronavirus Disease 2019 (COVID-19) Data. GitHub. External Links: Link Cited by: §3.4.
  • [19] S. Pasquali, A. Pievatolo, A. Bodini, and F. Ruggeri (2021) A stochastic SIR model for the analysis of the COVID-19 Italian epidemic. External Links: 2102.07566, Link Cited by: §1.
  • [20] C. Piazzola, L. Tamellini, and R. Tempone (2020) A note on tools for prediction under uncertainty and identifiability of SIR-like dynamical systems for epidemiology. Math. Biosci.. External Links: Document, 2008.01400, ISSN 18793134 Cited by: §1, §1.
  • [21] A. Rachah and D. F.M. Torres (2015) Mathematical modelling, simulation, and optimal control of the 2014 ebola outbreak in West Africa. Discret. Dyn. Nat. Soc. 2015. External Links: Document, 1503.07396, ISSN 1607887X Cited by: §1.
  • [22] L. S. R. Ross and H. P. Hudson (1917-05) An application of the theory of probabilities to the study of a priori pathometry.—Part III. Proc. R. Soc. London. Ser. A, Contain. Pap. a Math. Phys. Character 93 (650), pp. 225–240. External Links: Document, ISSN 0950-1207, Link Cited by: §1.
  • [23] L. S. R. Ross (1916-02) An application of the theory of probabilities to the study of a priori pathometry.—Part I. Proc. R. Soc. London. Ser. A, Contain. Pap. a Math. Phys. Character 92 (638), pp. 204–230. External Links: Document, ISSN 0950-1207, Link Cited by: §1.
  • [24] R. Ross and H. P. Hudson (1917-05) An application of the theory of probabilities to the study of a priori pathometry.—Part II. Proc. R. Soc. London. Ser. A, Contain. Pap. a Math. Phys. Character 93 (650), pp. 212–225. External Links: Document, ISSN 0950-1207, Link Cited by: §1.
  • [25] T. Sauer, T. Berry, D. Ebeigbe, M. Norton, A. Whalen, and S. Schiff (2020-06) Identifiability of infection model parameters early in an epidemic. medRxiv, pp. 2020.06.15.20132217. External Links: Document, Link Cited by: §1, §1.
  • [26] S. Srinath and R. Gunawan (2010-09) Parameter identifiability of power-law biochemical system models. J. Biotechnol. 149 (3), pp. 132–140. External Links: Document, ISSN 01681656, Link Cited by: §1.
  • [27] N. Tuncer, H. Gulbudak, V. L. Cannataro, and M. Martcheva (2016-09)

    Structural and Practical Identifiability Issues of Immuno-Epidemiological Vector–Host Models with Application to Rift Valley Fever

    .
    Bull. Math. Biol. 78 (9), pp. 1796–1827. External Links: Document, ISSN 15229602 Cited by: §1, §1.
  • [28] N. Tuncer and T. T. Le (2018-05) Structural and practical identifiability analysis of outbreak models. Math. Biosci. 299 (February), pp. 1–18. External Links: Document, ISSN 18793134, Link Cited by: §1, §1, §1, §2.
  • [29] A. F. Villaverde (2018) Observability and structural identifiability of nonlinear biological systems. arXiv 2019. External Links: 1812.04525, ISSN 23318422 Cited by: §1.
  • [30] X. Yang, S. Wang, Y. Xing, L. Li, R. Y. Da Xu, K. J. Friston, and Y. Guo (2020) Revealing the Transmission Dynamics of COVID-19: A Bayesian Framework for $R_t$ Estimation. External Links: 2101.01532, Link Cited by: §1.

Appendix A Derivation of Proposition 1

The main observation leading to (3) is that remains close to early in the epidemic. Motivated by this, we replace with in the SIR model to obtain

(11)

The corresponding solution starting from is

(12)

Fix and set , , , and . The expression of interest, namely , can be written

where is the error that results from approximating by . As mentioned in the main text, we do not have rigorous control over but numerical simulations indicate it is negligible compared to early in the epidemic. In particular, the approximation

(13)

is a good one prior to the time of peak infection. Thus we turn attention to the difference . Taylor expanding about gives, to first order,

where denotes the derivative in and . The time-varying matrix

has two eigenvalues, a nonpositive one,

, and a nonnegative one, . The nonpositive eigenvalue

has unit eigenvectors

for all . This suggests the perturbations that yield the smallest separation between and are those corresponding to the directions or, equivalently, the angles and . Now

(14)

where and . Plugging in gives

(15)

and similarly for , only negative. Thus, in combination with (13),

for all , which is precisely (3).

Appendix B Derivation of Proposition 2

Fix and set , , and . The likelihood for observed data is

and similarly for . So the log-likelihood ratio between and is

If the satisfy (2) with parameter , i.e. , then

So, letting be the value in (6) for the likelihood ratio test,

where

(16)

Now implies and hence

So, letting denote a standard normal random variable,

By an entirely similar computation (note the symmetry between and ),

and so . Therefore, since ,