0
$\begingroup$

1) Some background
I'm currently learning about compartmental ODE modeling for epidemiology (inspired by the current pandemic!) and I've been exploring parameter estimation of SIR-like models using optimization tools in Julia. One of the many challenges of modeling COVID-19, I have learned, is that the available time series data sets are probably quite noisy. Infection prevalence and incidence data in particular is probably very noisy, and subject to both under- and over-counting. In my experimentation with parameter estimation in Julia, I've found that small to moderate changes in the data point values can sometimes lead to significant changes in the parameter estimates. Consequently, I've become interested in modeling the error structure of the observed data so that I can get a better sense of the uncertainty in the parameter values. This leads me to

2) My Question: How can I model reporting error/noise in COVID-19 prevalence data?

By "prevalence", I mean the following definition from Wikipedia:

"Prevalence is a measurement of all individuals affected by the disease at a particular time."

This differs from "incidence" data which, to quote Wikipedia again, is "a measurement of the number of new individuals who contract a disease during a particular period of time".

3) More detailed background
As a simple example, consider the basic SIR (Susceptible-Infected-Recovered) model:

$$\frac{dS}{dt} = -\beta \frac{SI}{N}, \qquad \frac{dI}{dt} = \beta \frac{SI}{N} - \gamma I, \qquad \frac{dR}{dt} = \gamma I $$

where $N = S(t) + I(t) + R(t) = \text{const}$. Let's say that $t$ is in days. The prevalence on day $t$ would be $I(t)$; that is, $I(t)$ is the number of individuals that are actively infected on day $t$. The daily incidence for day $t$, which I'll denote by $\Delta C_t$, would be given by $\Delta C_t = C(t) - C(t-1)$, where $C(t)$ denotes the number of cumulative infections that have occurred by day $t$ (starting from some day $t_0$). (Note that $C(t) = I(t) + R(t)$.) So in other words, $\Delta C_t$ is the number of people that became infected in the one day period from $t-1$ to $t$.

When incidence data is available, there are some reasonable ways to model the error structure. For example, letting $\Delta C_1^{\text{obs}},\ldots,\Delta C_n^{\text{obs}}$ be the observed incidence data and $\Delta C_1^{\text{true}},\ldots,\Delta C_n^{\text{true}}$ be the "true" (but unobservable) incidence data, one reasonable (IMO) error model would be

$$\hspace{2cm} \frac{\Delta C_t^{\text{true}} - \Delta C_t^{\text{obs}} }{\Delta C_t^{\text{obs}}} = 1 + \epsilon_t, \quad \text{where } \epsilon_t \overset{\text{iid}}{\sim} N(0,\sigma^2) \qquad (1)$$

for a chosen value of $\sigma$. (Perhaps a truncated normal should actually be used--truncating $\epsilon_t$ to $[0,\infty)$ would ensure that $\Delta C_t^{\text{true}} \geq 0$, which we obviously want.) I find that the above model makes intuitive sense: It says that the relative error in the number of new cases reported on day $t$ is normally distributed with mean $0$ and variance $\sigma^2$. I've tested out the above model by simulating many sets $\{\Delta C_t^{\text{true}} \}_{t=0}^{n}$, fitting the SIR model to the simulated data sets, and then examining the distributions of the parameter estimates for $\beta$ and $\gamma$. The results I'm getting seem reasonable.

Now I'd like to repeat the procedure I just described, but using prevalence data. In the case of COVID-19, incidence data seems to be the most common type of data reported, but I have a data set I'm interested in that only contains prevalence data. (And I would just like to note: The relevance and importance of my question goes beyond my particular data set, and beyond just COVID. For example, the CDC collects and reports influenza prevalence data as part of its Epidemic Prediction Initiative.) When it comes to modeling error in prevalence data, things are trickier because the daily change in $I(t)$, call it $\Delta I_t := I(t) - I(t-1)$, is given by
\begin{align*} \Delta I_t &= \Delta C_t - \Delta R_t, \end{align*}

where $\Delta R_t = R(t) - R(t-1)$. In other words: $$\Delta I_t = (\# \text{ of new infections on day } t) - (\# \text{ of new recoveries on day } t).$$

Thus, $\Delta I_t$ depends on both the incidence of infection and the incidence of recovery. So my thinking is that an error model for $\Delta I_t$ should account for over- and under-reporting of both $\Delta C_t$ and $\Delta R_t$. (So it therefore probably does not make sense to simply substitute $\Delta C_t$ with $\Delta I_t$ into (1)...or could such a model be justified?)

Herein lies my dilemma: $\Delta I_t$ depends on $\Delta R_t$, and often there is not available or reliable data for $\Delta R_t$. By not "reliable" I mean that in the case of many COVID data sets, people who are 2 weeks post-infection are automatically classified as "recovered" (unless they've died, but for this toy model I'm ignoring deaths). Thus, I don't know if I would be able to simulate "noisy" $\Delta R_t$ data.

So to summarize... Is there a reasonable way to model the error in $\Delta I_t$ when prevalence data is the only data we have? If so, what are some error models that I could try? (I would also be interested to hear feedback and/or thoughts on error models for incidence data as well.)

$\endgroup$
1
+150
$\begingroup$

You’ve taken a little bit of an idiosyncratic approach to this problem, but it points to the two main ways people deal with this in the epi modeling community.

Probably the most common approach is to frame the problem as a partially observed Markov process—a model with an assumed Markov structure to the state space and an observation model that represents how (some) state variables are sampled. In this approach, implemented generically in the pomp package in R and described here for example http://www2.uaem.mx/r-mirror/web/packages/pomp/vignettes/intro_to_pomp.pdf, one doesn’t explicitly model the error structure between time point as long as one can simulate it under your model. This is necessary because most of the time for interesting models, the likelihood is intractable. Perhaps this framework will be enlightening to you.

The other approach less often used in practice is joint Bayesian estimation of the model parameters and the missing data using MCMC. Here’s an idealized tutorial https://si.biostat.washington.edu/sites/default/files/modules/2017_SISMID_8_Lab7.pdf, and a much deeper discussion https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6275108/. In this approach, you’ll see familiarity with your question about the missing “delta-R” data. Here, one simulates the missing data under the model while also estimating parameters, so that you can work with the tractable full data likelihood.

A third approach, much more similar to what you’ve laid out in notation and ideas—using a Gaussian approximate to describe the complete data likelihood above— is also described here https://iazpvnewgrp01.blob.core.windows.net/source/archived/Sustained_reductions_in_transmission_have_led_to_declining_COVID_19_prevalence_in_King_County_WA.pdf. It’s new enough we don’t have the formal publication yet.

In comparing your post to the approaches here, the challenge you’re facing is that you’re trying to write out something approximately equivalent to the likelihood of the process and observation model with the hidden states “integrated out.” This is intractable. Digging into either of the approaches above and the literature around them will get you where you need to go! Enjoy!

$\endgroup$
1
  • 1
    $\begingroup$ Thanks for your answer and for the useful links! $\endgroup$
    – Leonidas
    Jul 21 at 17:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.