Limits of epidemic prediction using SIR models

Abstract

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.

Keywords: Epidemic prediction; Hypothesis testing; Identifiability; Nonlinear dynamics; Parameter inference; SIR model.

Fig. 1

The trajectory of the SIR model used in the simulation (left). Plots of $\widehat{\beta }$ vs $\widehat{\gamma }$ from 1000 realizations from the sampling distribution of their MLE (right)

Fig. 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 $\Vert {\phi }_t(x,{\theta }_{ϵ}(\omega ))-{\phi }_t(x,\theta )\Vert$ between a perturbed trajectory and the true trajectory at time t. The gray, green, and black curves correspond to 90 perturbed trajectories, one for each of 90 equally spaced angles $\omega$ in $[0,2\pi )$. The black curves correspond to the angles $\pi /4$ and $5\pi /4$. The green curves correspond to the remaining angles in the intervals $[\pi /4-\pi /12,\pi /4+\pi /12)$ and $[5\pi /4-\pi /12,5\pi /4+\pi /12)$, i.e. in intervals of width $\pi /6$ centered at $\pi /4$ and $5\pi /4$, respectively. The gray curves correspond to those angles in $[0,2\pi )$ outside these two intervals. Note the distances corresponding to angles close to $\pi /4$ and $5\pi /4$ (the green and black curves) are smaller than those distances corresponding to angles farther away from $\pi /4$ and $5\pi /4$ (the gray curves), which supports the claim that the inverse problem is least practically identifiable for parameter perturbations approximately along a line of slope 1. The approximate lower bound of Approximation 1 is in red. The peak time of the trajectory corresponding to $\theta$ is indicated by the vertical blue line, and 80% of it by the vertical orange line. The first through fourth columns have population sizes ${10}^4,{10}^5,{10}^6$, and ${10}^7$, respectively, with only one initial infection in each case. The perturbation sizes for the first through fourth rows are $ϵ=.03,.03,.06$, and .1, respectively. The SIR paramaters for the first through fourth rows are $(\beta ,\gamma )=(.21,.14),(.21,.07),(.42,.07)$, and (1.68, .14), which give respective $R_0$ values of 1.5, 3, 6, and 12. Note the approximate lower bound holds roughly up to 80% of the peak time in all cases despite the wide range of parameters. Finally, we remark that the two seemingly “distinct" classes of gray curves in each plot correspond to different subsets of the 90 distinct angles. This as well as the multimodality of certain curves (which becomes more apparent when our graphs are extended further beyond the peak time) are consequences of the nonlinearity of the SIR model and are not directly relevant to our analysis

Fig. 3

Logarithm of distance between perturbed and true trajectories for different parameter values and population sizes. Everything is the same as in Fig. 2 except now we plot $log\Vert {\phi }_t(x,{\theta }_{ϵ}(\omega ))-{\phi }_t(x,\theta )\Vert$ instead of $\Vert {\phi }_t(x,{\theta }_{ϵ}(\omega ))-{\phi }_t(x,\theta )\Vert$. This gives a better view of the approximate lower bound early in the epidemic. Note the vertical axis is now a log scale

Fig. 4

Type II error as a function of perturbation size and noise level. The left panel shows the empirical and theoretical type II errors for the angles $\omega =0,\pi /4$, and $\pi$ as a function of perturbation size $ϵ$ with fixed noise level $\sigma =0.3$. The right panel shows the empirical and theoretical type II errors for the same angles as a function of noise level $\sigma$ with fixed perturbation size $ϵ=.03$. In each case the SIR parameters are those from Sect. 2, namely $(\beta ,\gamma )=(.21,.07)$, $N={10}^7$, and initial condition $i_0=1/N$. The time horizon T is 60 days into the epidemic, which in this case is 60 days prior to the peak time. The significance level is $\alpha =.05$. Here, theoretical refers to the first (red) and second (black) approximations of type II error ${\mathcal{E}}_2(\omega )$ in Approximation 2, i.e. Eqs. (8) and (9), respectively. Empirical refers to the type II error obtained by performing 1000 simulations of the noisy SIR model (2) followed by a likelihood ratio test of the hypothesis in (7) for each set of parameters. More specifically, the red and black curves lying over the blue line are the type II error approximations (8) and (9) when $\omega =\pi /4$, those lying over the purple line are when $\omega =\pi$, and those lying over the green line are when $\omega =0$, with the blue, green, and purple curves corresponding to the empirically computed type II error rates when $\omega =\pi /4,0$, and $\pi$, respectively. In each case both theoretical results closely align with the empirical ones, with the first approximation being slightly better than the second as expected. Also as predicted, the empirical type II errors all approach $1-\alpha =.95$ both as perturbation size goes to 0 and as the noise level gets large, and this approach is most rapid when $\omega =\pi /4$. In each case the noise model is Case 2, ${\sigma }_t=N\sigma i_t$. For the simulated blue, green, and purple curves, we used a numerical integrator to obtain the $i_t$ values, while for the red and black curves we used the pre-peak approximation $i_t\approx e^{\delta t}{i}_0$

Fig. 5

Consequences of type II error. The top panels show the total number of infections 10 days past the time of peak infection as a percentage of the total population. The bottom panels show the duration of the epidemic, which is defined to be the first day past the peak when less than 10 individuals are infectious. The left panels correspond to the parameter ${\theta }_0=(0.21,0.14)$ and the right panels to ${\theta }_0=(.21,.07)$. The red lines give the total percent infected or duration of the epidemic for the true parameter ${\theta }_0$ in each of their respective plots, while the blue and green curves give these values for ${\theta }_{ϵ}(\pi /4)$ and ${\theta }_{ϵ}(5\pi /4)$ over a range of $ϵ$ values, respectively. In all cases $N={10}^7$

Fig. 6

Practical identifiability of $\delta$. The left and center panels use the first type II error approximation in Approximation 2 to graph type II error as a function of perturbation size, $ϵ$, and noise level, $\sigma$, respectively. Each of the rainbow colored curves in both panels correspond to one of 150 different values of $\omega$ spread uniformly across $[0,2\pi )$. The color chart in the right panel indicates the colors corresponding to different angles $\omega$: The light red curves correspond to the $\omega$ closest to $\pi /4$ and $5\pi /4$, the yellow are those a bit farther away, the blue still farther, and the purple are those farthest from $\pi /4$ and $5\pi /4$, i.e. closest to $3\pi /4$ and $7\pi /4$. Finally, the dark red curve in each of the two panels corresponds to $\omega =\pi /4$ and $5\pi /4$, which have the same type II error. Note the rapid fall off in type II error as angles get farther from $\pi /4$ and $5\pi /4$, especially as a function of $ϵ$. This agrees with the empirical observation in Fig. 1 that MLE favors parameters lying along a line of slope 1. In particular, the inverse problem for $\delta$ is practically identifiable

Fig. 7

Going from hypothesis test (12) to (7)

Fig. 8

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 for each reporting rate

Fig. 9

Type II error rate as a function of $ϵ$ and $\sigma$ at significance level $\alpha =0.1$, reporting rate $p=0.15$, and $T=14$ days of new infection observations

Fig. 10

Error analysis corresponding to Figs. 2, 3. As in Figs. 2, 3, the horizontal axis in each graph is the number of days since the start of the epidemic, the vertical blue line is the peak time, and the orange line is 80% of the peak time. The vertical axis is the logarithm of the relative error, $log|E_t^{ϵ}|-log\Vert {\phi }_t^{ϵ}-{\phi }_t\Vert$, at time t. There are 50 green curves in each plot which are the relative errors corresponding to 50 angles in the intervals $[\pi /4-\pi /12,\pi /4+\pi /12)$ and $[5\pi /4-\pi /12,5\pi /4+\pi /12)$. The red line in each plot is the average linear approximation of the 50 green curves. Also as in Figs. 2, 3, the first through fourth columns have population sizes ${10}^4,{10}^5,{10}^6$, and ${10}^7$, respectively with one initial infection in each case, and the first through fourth rows have parameters $(\beta ,\gamma ,ϵ)=(.21,.14,.03),(.21,.07,.03),(.42,.07,.06)$, and (1.68, .14, .1) which give $R_0$ values of 1.5, 3, 6, and 12. The main point is that for every combination of parameters the relative error is exponentially small until about 80% of the peak time at which point it becomes $\mathcal{O}(1)$ and subsequently blows up exponentially