A tight fit of the SIR dynamic epidemic model to daily cases of COVID-19 reported during the 2021-2022 Omicron surge in New York City: A novel approach

Abstract

We describe a novel approach for recovering the underlying parameters of the SIR dynamic epidemic model from observed data on case incidence. We formulate a discrete-time approximation of the original continuous-time model and search for the parameter vector that minimizes the standard least squares criterion function. We show that the gradient vector and matrix of second-order derivatives of the criterion function with respect to the parameters adhere to their own systems of difference equations and thus can be exactly calculated iteratively. Applying our new approach, we estimated a four-parameter SIR model from daily reported cases of COVID-19 during the SARS-CoV-2 Omicron/BA.1 surge of December 2021-March 2022 in New York City. The estimated SIR model showed a tight fit to the observed data, but less so when we excluded residual cases attributable to the Delta variant during the initial upswing of the wave in December. Our analyses of both the real-world COVID-19 data and simulated case incidence data revealed an important problem of weak parameter identification. While our methods permitted for the separate estimation of the infection transmission parameter and the infection persistence parameter, only a linear combination of these two key parameters could be estimated with precision. The SIR model appears to be an adequate reduced-form description of the Omicron surge, but it is not necessarily the correct structural model. Prior information above and beyond case incidence data may be required to sharply identify the parameters and thus distinguish between alternative epidemic models.

Keywords: COVID-19; Newton–Raphson algorithm; Omicron variant; SARS-CoV-2; Susceptible-Infected-Removed model; heterogeneous mixing; inversion problem; nonlinear least squares; parameter identification; reduced form models; sloppy models; structural models; underreporting.

Figure 1.

Daily reported cases of COVID-19, 12/1/2021–3/15/2022, in New York City. The connected gray datapoints show the raw case counts $(c_t)$ . The red datapoints, covering 12/4/2021–3/12/2022, show the centered 7-day moving averages $(y_t)$ . For color figure, see online version.

Figure 2.

Estimates of the proportions of the Delta, Omicron BA.1, and Omicron BA.2 Variants in New York City During 11/24/2021–3/16/2022. Weekly averages (solid datapoints) were derived from a compilation maintained by the New York City Health Department. The Stata pchipolate interpolation routine was then employed to estimate the intervening days (connecting curves). The Omicron BA.2 category included strains identified as BA.2, BA.2.12.1, and BA.2.75. Not shown are the proportions of all other strains, which represented no more than 0.7 percent of total samples in any one week. For color figure, see online version.

Figure 3.

Daily reported and predicted cases of COVID-19, 12/1/2021–3/15/2022, in New York City. Estimates Based Upon All Variants (Blue) and Omicron BA.1 Variant Only (Orange). Blue datapoints represent all reported cases $y_t$ , while orange datapoints represent Omicron BA.1 cases ${\hat{y}}_t=h_t{y}_t$ . The respective blue and orange curves connect the predicted values of the output variable $X_t$ . The displayed estimates are based upon the 4-dimensional Newton–Raphson algorithm described in Propositions 2 and 3. The estimates based upon the alternative EM procedure, described in Propositions 4 and 5, were identical. For color figure, see online version.

Figure 4.

Panel A: Least squares criterion $\mathit{V}$ (left axis) and first partial derivative $\mathrm{\partial }\mathit{V}/\mathrm{\partial }\mathit{\beta }$ (right axis) as functions of the parameter $\mathit{\beta }$ . Panel B: Least squares criterion $\mathit{V}$ (left axis) and first partial derivative $\mathrm{\partial }\mathit{V}/\mathrm{\partial }\mathit{\alpha }$ (right axis) as functions of the parameter $\mathit{\alpha }$ . For color figure, see online version.

Figure 5.

Least squares criterion $\mathit{V}$ projected onto the $(\mathit{\beta },\mathit{\alpha })$ plane. The center point identifies the minimum where $(\beta *,\alpha *)$  =  (0.233, 0.941). The parameters $i_0$ and N were held constant at their optimum values given in Table 1. For color figure, see online version.

Figure 6.

Contours of $\mathit{V}$ running through the optimum point ${\mathbf{\Theta }}^{\mathbf{\ast }}$ projected onto the $(\mathit{\beta },\mathit{\alpha })$ plane. All plotted contours run through the optimum $(\beta *,\alpha *)=$ (0.233, 0.941), with the parameters $i_0$ and N held constant at their optimum values. With $\beta$ also held constant at its optimum value (green contour), the second-order derivative was ${\mathrm{\partial }}^{2}V({\mathit{\theta }}^{*})/\mathrm{\partial }{\alpha }^2$  = 4.912 $\times$ 10⁶. With $\alpha$ also held constant at its optimum value (magenta contour), the second-order derivative was ${\mathrm{\partial }}^{2}V({\mathit{\theta }}^{*})/\mathrm{\partial }{\beta }^2$  = 6.075 $\times$ 10⁶. Along the ray defined by $w*$  = 1.073 (orange contour), the second-order directional derivative was $(d^{2}V/d{u}^2){|}_{u=0}$  = 0.418 $\times$ 10⁶. For color figure, see online version.

Figure 7.

Estimated parameters based upon a simulated SIR model with Gaussian Error. The SIR model of (2), (3) and (4) was simulated with assumed known parameters $\mathbf{\Theta }=(\beta ,\alpha ,i_o,N)=(0.6,0.7,0.04,{10}^5)$ over $T=100$ time periods. The observed case incidence was then computed as $y_t=X_t+{\epsilon }_t$ , where the independent Gaussian errors ${\epsilon }_t$ were randomly drawn with zero mean and standard deviation $\sigma =100$ . Panel A: Plots of parameter estimates ${\beta }^{(j)}$ versus ${\alpha }^{(j)}$ with all parameters unconstrained (gold solid points) and conditional upon a known population N (dark blue circles). Panel B: Box-and-whisker plots of the corresponding estimates ${\mathcal{R}}_0^{(j)}$ of the basic reproduction number. The boxes show the 25th, 50th, and 75th percentiles while the capped lines show the bounds of interquartile range. For color figure, see online version.

Figure B1.

Path of $({\mathit{\beta }}^{(\mathit{k})},{\mathit{\alpha }}^{(\mathit{k})},{\mathit{V}}^{(\mathit{k})})$ through successive iterations of the 4-dimensional Newton–Raphson algorithm: New York City Omicron wave. The algorithm began at point A, where $({\beta }^{(0)},{\alpha }^{(0)},V^{(0)})=(0.4,0.8,3907)$ . The step size for successive iterations, defined in equation (16), was $q=\mathrm{0.5.}$ At point B, where $({\beta }^{(10)},{\alpha }^{(10)},V^{(10)})=(0.446,0.730,101)$ , the path of the algorithm had entered the ravine described in Figure 6. The path continued along the ravine until the stopping point at $({\beta }^{(40)},{\alpha }^{(40)},V^{(40)})=(0.233,0.941,76.7)$ . The stopping criterion was $|V^{(k)}-V^{(k-1)}|<{10}^{-4}$ . The execution time was 0.436 s.