Checking Stochastic SIR Models Against Epidemic Data

What Should a Stochastic SIR Model Match?

A stochastic SIR model does not predict one fixed epidemic curve. It specifies a probability law for possible epidemics, so a comparison with data must distinguish the model's typical behavior from the variation among its possible outcomes. A useful check examines several features together: incidence through time, the timing and height of the peak, cumulative infections, and variation in these quantities across outbreaks. A discrepancy in overall scale may be consistent with under-reporting; a discrepancy in timing or curve shape generally cannot be removed by a single constant reporting fraction.

Consider a finite population of size \(N\), with state \((S_t,I_t,R_t)\) at time \(t\ge0\), on a probability space \((\Omega,\mathcal{F},\mathbb{P})\). The state space is \(\{(s,i,r):s,i,r\text{ are nonnegative integers and }s+i+r=N\}\). In the continuous-time stochastic SIR model, infections occur at rate \(\beta S_tI_t/N\), changing the state by \((-1,+1,0)\); recoveries occur at rate \(\gamma I_t\), changing it by \((0,-1,+1)\). Here \(\beta,\gamma\ge0\). These rates define the model's law given its parameters and initial state; they do not say that every realized epidemic follows its expected curve.

For bins \((t_{k-1},t_k]\), let \(C_k\) be the number of new infections in the model during bin \(k\). The incidence vector is \((C_1,\ldots,C_K)\), and the final number of new infections is \(C_{\mathrm{tot}}=\sum_{k=1}^K C_k\), when the bins cover the epidemic being assessed. If the initial number infected is \(I_0\), the final number ever infected is \(I_0+C_{\mathrm{tot}}\). For observed data, write \(Y_k\) for reported cases assigned to bin \(k\). Define peak timing as the earliest bin index attaining the largest incidence; this convention handles ties.

Adding an Observation Model

Latent infections are often not observed one for one. A simple under-reporting model assumes each new infection is reported independently with a constant probability \(r\), where \(0<r\le1\), and that reports are assigned to the correct time bin. Conditional on the latent counts, \(Y_k\) has a binomial distribution with \(C_k\) trials and success probability \(r\). This assumption excludes reporting delays, changing ascertainment, and false-positive reports. If those mechanisms matter, a constant reporting fraction is not an adequate observation model.

Proof. Conditional on \(C_k\), a binomial count has conditional mean \(rC_k\) and conditional variance \(r(1-r)C_k\). Taking expectations of the conditional mean gives \(\mathbb{E}[Y_k]=r\mathbb{E}[C_k]\). The conditional-variance identity gives \(\operatorname{Var}(Y_k)=\mathbb{E}[\operatorname{Var}(Y_k\mid C_k)]+\operatorname{Var}(\mathbb{E}[Y_k\mid C_k])\), which equals \(r(1-r)\mathbb{E}[C_k]+r^2\operatorname{Var}(C_k)\). For any two bins \(a,b\), if \(\mathbb{E}[C_a]>\mathbb{E}[C_b]\), multiplying by \(r>0\) preserves the strict inequality, so \(\mathbb{E}[Y_a]>\mathbb{E}[Y_b]\). Thus the unique maximizing bin is unchanged. This is a statement about expected curves, not a guarantee that the peak of one noisy reported epidemic occurs in the same bin as the peak of one latent realization. \(\square\)

The variance formula separates two sources of variation. The term \(r^2\operatorname{Var}(C_k)\) comes from variation in epidemic incidence; the term \(r(1-r)\mathbb{E}[C_k]\) is additional reporting variation. Consequently, comparing the variance of observed counts directly with the variance of latent model counts is not a like-for-like check. The reporting mechanism must be included.

Fit Scale, Then Inspect Shape

When there are several bins, ratios \(Y_k/\mathbb{E}[C_k]\) can be unstable or undefined if a model expectation is zero. A simple alternative is to fit one common scale by least squares. Let \(x_k=\mathbb{E}[C_k]\) be the model-implied latent mean incidence and \(y_k\) the observed count. The estimate below is a descriptive scale fit to these summaries; it is not, by itself, a likelihood estimate or a full uncertainty analysis.

Proof. Put \(A=\sum_kx_k^2>0\) and \(B=\sum_kx_ky_k\). Expanding the objective gives \(\sum_k y_k^2-2rB+r^2A\). Completing the square, this is \(\sum_k y_k^2-B^2/A+A(r-B/A)^2\). Since \(A>0\), minimizing over \(r\in[0,1]\) is equivalent to choosing the point in that interval closest to \(B/A\). If \(B/A<0\), the closest point is \(0\); if \(0\le B/A\le1\), it is \(B/A\); and if \(B/A>1\), it is \(1\). These cases give the stated formula. The assumption \(A>0\) is necessary: if every \(x_k=0\), the objective does not depend on \(r\), so no unique scale is determined. \(\square\)

After fitting the scale, inspect the residuals \(y_k-\widehat r x_k\), the observed and predicted peak bins, and cumulative totals. A small overall error can still hide a systematic early-versus-late pattern. In addition, a scale discrepancy is not uniquely attributable to under-reporting: a wrong initial number infected, a wrong population denominator, or an incorrect incidence definition can also alter the scale.

Final Size and Outbreak Variability

Incidence by bin describes the course of an outbreak; final size describes its cumulative reach. Compare the observed cumulative count with the model's distribution of final size, not only with its mean. In a stochastic epidemic, some realizations may end quickly with few infections while others grow much larger. A single observed final size can therefore be unusual without being impossible under the model. Conversely, a model whose average final size is close to the observation may still assign too much probability to very small or very large outbreaks.

Under independent constant-probability reporting, the expected reported total is \(r\mathbb{E}[C_{\mathrm{tot}}]\), since expectation is additive. Its variance also includes reporting noise. If reports for all cases are conditionally independent and \(C_{\mathrm{tot}}\) is the latent total, the same variance calculation gives \(r(1-r)\mathbb{E}[C_{\mathrm{tot}}]+r^2\operatorname{Var}(C_{\mathrm{tot}})\). This permits a comparison of observed final-size variation with a model that includes ascertainment. If only one epidemic is observed, its value alone cannot estimate the across-outbreak variance reliably; the model's variability is a distributional prediction, not an estimate obtained from a single realization.

Interpreting a Mismatch

A constant reporting fraction changes the expected level of incidence and final size, but it does not change the expected peak bin under the theorem's conditions. A peak that is systematically too early or too late, or an incidence curve with a different rise-and-fall pattern, therefore calls for more than a constant scale correction. Transmission-related parameters can affect epidemic growth and timing, while recovery dynamics and susceptible depletion also shape the trajectory. However, a shape mismatch is not proof that transmission alone is wrong: changing ascertainment over time, reporting delays, errors in the observation window, or a misspecified initial state can create similar patterns.

Check model fit in a sequence that keeps these interpretations separate:

• Confirm that observed cases and model infections use the same time bins and case definition.
• Compare incidence curves, including the peak bin and whether the observed shape can be a constant multiple of the model mean.
• Compare cumulative incidence and final-size distributions, allowing for the reporting mechanism.
• Compare variability using model-generated distributions with reporting variation included, rather than treating the mean curve as a deterministic prediction.
• Interpret discrepancies as evidence against a particular combination of transmission, initial conditions, and observation assumptions—not as identification of one cause from one summary.

Common mistake. Suppose model-implied expected incidence is \((2,8,10,4)\), but an analyst multiplies it by \(1/2\) and claims the resulting curve matches data \((1,9,7,2)\). The calculation gives \((1,4,5,2)\), not \((1,9,7,2)\); in particular, the observed peak is in bin 2 while the scaled expected peak remains in bin 3. A reporting fraction can lower the expected level, but it cannot move the peak of the expected curve under constant-probability, no-delay reporting. The fitted-scale calculation and residual inspection are needed to reveal this distinction.