Tracking the transmission dynamics of COVID-19 with a time-varying coefficient state-space model.

The spread of COVID-19 has been greatly impacted by regulatory policies and behavior patterns that vary across counties, states, and countries. Population-level dynamics of COVID-19 can generally be described using a set of ordinary differential equations, but these deterministic equations are insufficient for modeling the observed case rates, which can vary due to local testing and case reporting policies and nonhomogeneous behavior among individuals. To assess the impact of population mobility on the spread of COVID-19, we have developed a novel Bayesian time-varying coefficient state-space model for infectious disease transmission. The foundation of this model is a time-varying coefficient compartment model to recapitulate the dynamics among susceptible, exposed, undetected infectious, detected infectious, undetected removed, hospitalized, detected recovered, and detected deceased individuals. The infectiousness and detection parameters are modeled to vary by time, and the infectiousness component in the model incorporates information on multiple sources of population mobility. Along with this compartment model, a multiplicative process model is introduced to allow for deviation from the deterministic dynamics. We apply this model to observed COVID-19 cases and deaths in several U.S. states and Colorado counties. We find that population mobility measures are highly correlated with transmission rates and can explain complicated temporal variation in infectiousness in these regions. Additionally, the inferred connections between mobility and epidemiological parameters, varying across locations, have revealed the heterogeneous effects of different policies on the dynamics of COVID-19.

INTRODUCTION

Since its spread at the beginning of 2020, the COVID‐19 pandemic has made clear the impact that human mobility and government policies can have on the spread of novel respiratory diseases. Throughout the world, different approaches have been taken to mitigate the spread of COVID‐19, including travel bans, mask mandates, quarantine policies, and capacity restrictions. The diversity of actions taken presents a unique opportunity and an urgent demand to assess the impact of these approaches quantitatively and systematically.

Statistical approaches offer opportunities to explore the relationship between different policies and actions and disease dynamics while allowing for variation in the outcome from both measured and unmeasured (or even unmeasurable) factors. However, a purely marginal comparison of factors—as might be obtained using a regression‐based approach—can miss important facets of the dynamics of disease spread. That is why the foundation for most infectious disease models is a compartmental model for the disease states, which captures these dynamics. The traditional example of this is a compartmental “SIR” model that employs ordinary differential equations to model transition from “susceptible” to “infected” to “recovered” states. , , , However, the traditional SIR model lacks the flexibility to account for the sophisticated and rapidly changing transmission dynamics of COVID‐19. To capture the key characteristics of the COVID‐19 pandemic, a number of alterations to the traditional SIR model are needed, including changes to the structure of the compartmental model, incorporation of time‐varying factors that may influence the dynamics of the disease spread over the course of the pandemic, and introduction of stochasticity that results from limitations in model assumptions as well as noise in the observed data.

First, there are characteristics of COVID‐19 spread that require an adaptation in the structure of the compartmental model. A key characteristic of COVID‐19 is that many of those infected are at their most infectious before they are diagnosed, with viral load typically highest at the onset of symptoms. Thus, there can be a misalignment if a compartmental model assumes people can only be infectious following diagnosis through a positive test. Further, there are many who have milder or asymptomatic cases and are never officially diagnosed as “infected,” but who may still spread the disease. The framework of a dynamic model should address these facets of how the data (from testing) lines up with the principles of infectious spread.

Second, since the emergence of SARS CoV‐2, the dynamics of COVID‐19 detection and spread have evolved, and the model should allow for an associated evolution in some of its parameters. Public health guidance has changed throughout 2020 and 2021, with changing guidance and regulations on whether to wear masks, travel, quarantine, gather in groups, attend in‐person meetings and school, etc., which can change the rate of transmission from those who are infected to those who are still susceptible. Testing has improved, become more easily accessible, and in some cases evolved to include regular testing even without symptoms, all of which can influence the probability that someone infected with the disease is detected. Treatment has improved, including through the adaptation of corticosteriods and remdesivir among patients with severe disease and care techniques like prone positioning, use of high‐flow oxygen therapy, and intubation timing. , These strategies may help in reducing mortality rate of patients with severe disease since the start of the pandemic; indeed, there is evidence of a decrease in risk of mortality among those hospitalized for COVID‐19 over 2020, both in the United States and the United Kingdom. New strains have evolved and gained prevalence during the pandemic, with different infectivity and severity. , Because of these changing dynamics, parameters of the compartmental model—including rate of detection, probability of infection given a contact with someone infected, and rates of mortality versus recovery among the infected—have changed over the pandemic, suggesting the need for time‐varying coefficients for modeling.

Finally, there is stochasticity within the observed data compared to the dynamics of the model. This stochasticity results both from imperfect compliance with the assumptions of the model and from practical constraints in collecting and reporting the data. For example, the standard compartment model assumes that a community's population is well‐mixed, with an equal chance of contact among any pair of members of that population. For COVID‐19, as with most infectious diseases, this assumption is overly simplistic and requires allowance for stochasticity, as local spikes in cases were sometimes the result of outbreaks within an institution or organization, including significant outbreaks in prisons, nursing homes, and meat‐processing plants, , , suggesting the pandemic dynamics were in part driven by more frequent contacts within population subsets, rather than across a well‐mixed community population. Further, the model assumes homogeneity in individual transmissability, while in fact some spread is driven by superspreaders. , , , Also, measurement error is introduced through data collection. Public health officials have made an enormous effort to collect and publish counts of cases and deaths during the pandemic, but understandably there were occasional patterns in the data related to data collection and reporting rather than dynamics of the virus' spread. For example, Colorado included death counts for deaths that occurred earlier but had not yet been reported on April 24, 2020. Further, on weekends and holidays, reporting rates can be lower than usual, with reporting higher following the break to incorporate the backlog.

Here, we develop a time‐varying coefficient state‐space model that uses a structure appropriate for COVID‐19 and allows for stochasticity and measurement error in data, as well as evolution of some model parameters over time, with the aim of investigating how a specific factor was associated with virus spread. This model has an advantage over regression‐based models of factors that may affect COVID case counts over time, since under a compartmental modeling framework we are modeling the process of disease spread, rather than correlating two time series. Further, by incorporating elements that address time‐variability in model parameters and stochasticity inherent in the relationship between the available COVID‐19 data and the compartmental model, the model addresses limitations in a classic SIR‐style model for disease spread. Specifically, we introduce a multiplicative process model and a negative binomial data model to account for fluctuations in the rates of disease detection and transmission and variability in data reporting. Focusing on retrospective estimation and inference, our model offers a framework to explore and draw inference on the effect of different human mobility behavior and related policies on the pandemic parameters in the model. It therefore meets a critical need to understand how policy choices might affect the dynamics of spread.

As mentioned, a factor of particular interest is the influence of local mobility on COVID‐19 spread within a community. Recent studies have suggested strong ties between COVID‐19 infection rates and human mobility. Initially in China, there was a strong relationship between the number of COVID‐19 cases and human mobility. This finding is consistent with the theory of infectious disease spread in highly coupled metapopulations. , This relationship weakened after control measures were put into place to restrict the movement of people in and out of Wuhan province. This previous study used real‐time mobility data as well as travel history data to explore the relationship to spread of the disease. They concluded that the drastic control measures implemented in China substantially mitigated the spread of COVID‐19. Decreased mobility was also shown to have a protective effect against COVID‐19 transmission in the United States, a result that agreed with other findings that mobility inflow into a county early in the pandemic was associated with increases in early case counts. While the official response in the United States to COVID‐19 has been heterogeneous in terms of lock‐downs, this work showed that social distancing helps to reduce the spread of the disease, and should remain part of personal and institutional responses to the pandemic until a vaccine is widely adopted and should continue to be considered as a key protective public health policy in future pandemics of respiratory diseases.

While the relationship between mobility and rates of COVID‐19 infection has previously been explored, , , , , analysis of this relationship has largely relied on regression‐based comparisons of time series data, rather than through incorporating observed mobility within the dynamics of an epidemiological model. For example, Kraemer et al used generalized linear models with daily new/cumulative cases as the outcome and mobility as one of the covariates. On the other hand, compartmental models with time‐varying coefficients have been employed to study the transmission dynamics of COVID‐19, , , , , , , , , , but a large portion of these works were not specifically aimed at investigating the effect of mobility on the spread of COVID‐19. For example, Wu et al estimated the basic reproductive number of COVID‐19 under a compartmental modeling framework but did not attempt to associate it with mobility. The compartmental models in several existing works were overly simplified, and thus might miss key features of the dynamics of COVID‐19. For example, Wang et al and Chen and Qiu did not explicitly consider a compartment for infectious but undetected individuals; Zhou and Ji did not include a compartment for exposed yet noninfectious individuals. In addition, some existing works focused on obtaining point estimates of the epidemiological parameters but did not offer uncertainty quantification. For example, Giordano et al fitted their compartmental model by minimizing the sum of the squares of the errors but did not provide confidence intervals for the parameters. Lastly, many existing works (such as Li et al, Hao et al, and Yan et al ) assumed a Poisson model for the case counts, which can be vulnerable to overdispersion and large variability in the observed data. Aiming to overcome these limitations, we develop a time‐varying coefficient state‐space epidemiological model that allows us to explore the relationship between mobility and COVID‐19 spread. The structure of this model incorporates mobility as a factor that influences the dynamics of the epidemic, while also accounting for variation over time in some model parameters and stochasticity within the observed case and fatality data compared to the model dynamics. This approach allows us to explore how mobility influenced the dynamics of COVID‐19 spread, while the model development provides a structure that can be extended to explore how other factors (such as environmental factors and mitigation policies) influence the dynamics of COVID‐19 pandemic and spread of other diseases in the future. We consider a partition of the population into eight compartments to realistically reflect key features of COVID‐19 spread. By incorporating a multiplicative process model and a negative binomial data model, our model allows for stochastic deviations from the trajectories given by a system of differential equations and better captures the large variability in data reporting. Through Bayesian inference, we provide coherent uncertainty quantification for the epidemiological parameters. Our model is applied to several U.S. states and Colorado counties.

This article is organized as follows. In Section 2, we present a deterministic compartmental model for COVID‐19 with time‐varying coefficients, characterized by a system of differential equations. In Section 3, by extending the compartmental model outlined in Section 2, we develop a state‐space model for COVID‐19, which better accounts for stochasticity in disease spread. We apply our proposed method to analyze county‐level COVID‐19 data in the U.S. state of Colorado and state‐level data in the United States in Section 4. We conclude this article with some remarks in Section 5. Extra results are deferred to the online supporting materials.

A TIME‐VARYING COEFFICIENT COMPARTMENTAL MODEL FOR COVID‐19 DYNAMICS

Given its flexibility in structure and natural connection with dynamical systems, the compartmental modeling approach has been extensively employed in epidemiology , , , , and pharmacokinetics. , This approach has been widely adopted to understand the dynamics of COVID‐19. , , , , , , , , , We first outline a basic compartmental model for recapitulating the transmission dynamics of COVID‐19, which is characterized by a system of differential equations. It will serve as the cornerstone for our statistical model in Section 3. We partition the population in a region (which can be a county, state, or country) into the following eight compartments, each representing a specific stage of COVID‐19:

Susceptible individuals (): those who have not been infected by the disease and are at risk of becoming infected;

Exposed individuals (): those who have been exposed to the disease but are not yet capable of infecting others;

Undetected infectious individuals (): those who have the disease, are able to infect others, but have not yet been detected as having the disease;

Detected infectious individuals (): those who have the disease (either with or without symptoms), are able to infect others, and have been detected as having the disease;

Individuals removed from without being detected (): those who had the disease but are then removed from the possibility of being infected again or spreading the disease for any reason (eg, significant reduction in viral load, recovery, or death), without ever being diagnosed with the disease;

Individuals who are hospitalized : those who have been identified as cases and are hospitalized.

Individuals removed after having been detected (): those who have been through the state (and possibly state as well) but are then removed from the possibility of being infected again or spreading the disease for any reason except death. They are assumed to eventually recover from the disease;

Individuals who died (): those who have been through the state and are eventually deceased.

The eight compartments, together with the possible transitions among compartments, are illustrated in Figure 1 and are designed to capture important features of the dynamics of COVID‐19. It has been well recognized that asymptomatic individuals, presymptomatic individuals, and those who are symptomatic but have not yet been diagnosed with the disease (collectively denoted by ), contribute significantly to the spread of COVID‐19. , , Among this group of undetected cases, transmission may differ between those who are symptomatic vs. asymptomatic—for example, symptomatic cases have been found to have a higher chance of passing COVID‐19 to close contacts, although those who are symptomatic may be less likely to mix with the general population, if they stay home while sick. However, since the observed data (to be explained in Section 3.2 and Section C in Appendix S1) do not provide any information for this partition, we do not separately model symptomatic and asymptomatic cases. As time progresses, undetected infectious individuals become detected or get removed without ever being detected. Detected infectious individuals become hospitalized or get removed without ever being hospitalized. Individuals who are hospitalized may recover or die, and we assume individuals who die must go through the hospitalization compartment. To reflect this, the proposed model specifically forces individuals in to pass through first. Nevertheless, detected individuals can recover without having gone to the hospital, and so individuals in can, but do not have to, pass through first.

FIGURE 1

The conceptual dynamics for our eight‐compartment time‐varying coefficient dynamical model. Square nodes represent observed states and elliptical nodes represent unobserved states. Compartment transitions are labeled with their rates

We use the notation of a compartment followed by a time index to denote the size of the compartment at a specific time point. For example, denotes the number of susceptible individuals at time . The transmission of COVID‐19 can be characterized by the flow of individuals through these compartments over time, given by the following system of differential equations with time‐varying coefficients: Here, denotes the population of active individuals. A simplification in model (1) to (8) is that we ignore births and deaths, an approximation that is appropriate for a fast‐spreading pandemic like COVID‐19. For simplicity, we also do not consider immigration and emigration. Given the initial condition (ie, the initial sizes of the compartments) and parameter values, the trajectory of the epidemic process can be deterministically obtained by solving the system of differential equations. In Sections 3.1 and 3.2, we will discuss how to allow for additional flexibility and stochasticity in the time‐varying coefficient compartmental model (1) to (8), and how to link the observed data on COVID‐19 cases and deaths to the epidemic process. The initial conditions of the model are not trivial and will be discussed in Section 3.3.

The time‐varying coefficients in (1) to (8) bring extra flexibility to the compartmental model to capture the complex and rapidly varying dynamics of COVID‐19. Specifically, models the disease transmission rate between the undetected infectious and susceptible individuals at time . The rate can be understood as the number of effective contacts (ie, contacts that are sufficient for disease transmission) made by an average member of the undetected infectious individuals per unit time. The probability of each of these contacts being with a susceptible individual is . Therefore, undetected infectious individuals lead to a rate of new infections . For the detected infectious individuals, the rate of effective contacts is reduced by a factor of , since those diagnosed with COVID‐19 are likely to have reduced contact with others due to potential quarantine or self‐isolation, leading to a rate of new infections . Lastly, we assume there is no meaningful amount of transmission among hospitalized people () due to precautions undertaken in the hospital setting. To incorporate additional covariate information such as mobility and changes in policy, a detailed model of will be provided in Section 3.4.

After being infected, a susceptible individual first goes through an “exposed” state, meaning that the individual has been exposed to the disease but is not immediately able to infect others. The exposed individuals enter the undetected infectious state at a rate of , and naturally represents the latent period, which is the time interval between when an individual is exposed to the disease and when the individual becomes capable of infecting other susceptible individuals. Then, the undetected infectious individuals are either diagnosed with the disease at a detection rate of or removed from the infectious state at a rate of without ever being diagnosed. Similar to , the detection rate is expected to vary along with the development of COVID‐19 pandemic. Estimation of plays an indispensable role in understanding the epidemiological mechanics of COVID‐19. In our model, we further assume that an infectious individual always needs to go through the “undetected” state before the individual is detected as having the disease. It is possible that an exposed individual can be directly diagnosed with the disease through contact tracing, but we treat this as a special case that the individual spends zero time at the undetected infectious state. Infectious individuals that have been detected become hospitalized at rate or recover (without ever being hospitalized) at rate . Hospitalized individuals are removed either through recovery at a rate of or death at a rate of . Naturally, deceased individuals do not contribute to the dynamics of disease transmission. We also assume that recovery from COVID‐19 confers immunity to reinfection, and thus the recovered individuals can neither spread the disease nor be infected again. This is a simplifying assumption motivated by the rate of reinfection being so low during the time period studied that it does not meaningfully alter the transmission dynamics. In Section 3.5 we discuss the modeling and interpretation of each of these parameters in greater detail.

For epidemiological models, the reproductive ratio is a fundamental quantity to track the pandemic and it serves as a threshold that predicts the spread of an infection. , For models that are more sophisticated than the simple SIR model, the reproductive ratio is usually computed based on the equilibrium reproduction number, which characterizes the persistence of the pandemic. For the proposed model (1) to (8), the equilibrium reproduction number is which is derived using the endemic equilibrium argument in Section A in Appendix S1. If after time , the number of infectious individuals will decrease after time and lead to the disease‐free equilibrium. Therefore, an indicates containment of the disease at time . On the contrary, if after time , the pandemic will persist with a nontrivial equilibrium and the larger the equilibrium reproduction number, the larger the population of infectious individuals at the equilibrium will be. Due to its important role in characterizing disease spread, inference on is one of our major interests. In addition, differentiating against , we observe that is monotonically decreasing in if and only if . By restricting and modeling the relationship between and (accomplished via prior specification, see Section 3.5), it is guaranteed that . As , the difficulty of detection increases and monotonically converges to . That is, the pandemic dynamics are dominated by the undetected infectious individuals. On the other hand, as increases, the deviation between and shrinks. That is, the pandemic dynamics will be governed more by the detected infectious individuals.

STATE‐SPACE MODEL FOR COVID‐19 WITH TIME‐VARYING COEFFICIENTS

While the compartmental model in Section 2 describes the underlying deterministic and smooth disease dynamics, there is stochasticity in the actual disease spread. This is due to factors such as nonuniform mixing of the population and heterogeneity in day‐to‐day activity patterns. Additionally, data on COVID‐19 cases and deaths are always observed with substantial statistical noise, such as reporting error and reporting delay, which cannot be easily captured by the differential equations. To mitigate these issues, we cast the compartmental model in Section 2 in a state‐space modeling framework. Following an approach similar to Dukic et al and Osthus et al, we construct a process model which allows the epidemic process to stochastically deviate from the solution given by the compartmental model (1) to (8). Moreover, we build a data model for the observed data, which further takes into account the measurement error of COVID‐19 cases and death counts.

Process model

Since COVID‐19 data are reported daily, it is natural to consider a discretized version of the compartmental model in Section 2 with a time step of 1 day. Denote the true but unobservable populations of the eight compartments on day . Define as the solution to differential equations (1) to (8), starting from the state of the prior day, .

Denote by the unobservable number of individuals diagnosed with COVID‐19 on day . Similarly, let and denote the unobservable number of individuals who are hospitalized or died from COVID‐19 on day , respectively. In terms of states and , represents the number of individuals who moved from to , represents the number of individuals who moved from to , and represents the number of individuals who moved from to . As discussed before, it is unrealistic to assume that the true yet unobservable numbers of the diagnosed and reported deceased individuals perfectly agree with the solution to the compartmental model, especially for a complicated pandemics such as COVID‐19. Similarly, the number of reported hospital admissions will likely not perfectly agree with the true number of admitted individuals with COVID‐19. To introduce additional flexibility, motivated by Davis and Wu, we model , , and via multiplicative processes where , , and refer to the numbers of new detections, hospitalizations, and deaths on day , according to the solution to our compartmental model. For simplicity, we model the multiplicative processes , , and as random variables independent in with so that , , , , , and . The unit mean assumption, which mimics the zero mean assumption for the additive models, leads to , , and . Compared to the traditional additive models, the multiplicative process model (10) is more natural to introduce deviations between the nonnegative solution to differential equation models and unobservable process with nonnegative observable space. In summary, the underlying epidemic process is centered at the solution to differential equations (1) to (8) but is also allowed to be stochastically different from the deterministic process specified by the compartmental model.

The solution to the differential equations (1) to (8) is not available in closed form. A numerical approximation method, such as the Runge‐Kutta solver, must be employed. Here we use a simple and computationally efficient first‐order difference calculation, which does not compromise the model performance. From (5), a first‐order difference update to the states means that . Thus, , and . We assume that in (10) only affects the transition from state to state , only affects the transition from to , and only affects the transition from to . Thus, we have In essence, the error process therefore offers additional flexibility in the detection rate from the compartmental model. Accordingly, the empirical calculation of will be adjusted by replacing in (9) with and with . As both and individuals are unable to spread the disease, does not influence the equilibrium reproduction number.

Lastly, based on our assumptions, the error terms do not have an impact on the states , , , and , and these states are consistent with the corresponding solutions to the compartmental model , , , and , calculated using the first‐order approximation.

Data model

Let be the number of newly reported positive cases on day , where . Similarly, denote the number of newly reported hospitalizations and the number of newly reported deaths on day . In the event that only cumulative numbers of cases, hospitalizations, and deaths are available, the values of , , and can be computed using a first‐order difference. We assume negative binomial models for the daily reported case, hospitalization, and death counts, where and , respectively. Parameter allows for overdispersion, relative to a Poisson distribution, and grants more flexibility to the model for fitting the COVID‐19 data. We place a Gamma prior on , , while and are treated similarly.

Model for the initial condition

The initial condition of the infection dynamics refers to the initial population sizes of the compartments. The calendar time corresponding to can be chosen for each modeling context. We assume that at , the total number of removed individuals is 0, that is, . The number of detected infectious individuals at time 0, , is observed and assumed to be non‐zero. In the analyses of Section 4, we choose to be a time early in the pandemic when there are multiple detected cases.

The numbers of susceptible, exposed, and undetected infectious individuals at time 0 are not observed and need to be estimated. We assume the initial population size of the exposed and undetected infectious individuals is times the initial size of the detected infectious individuals, We place a gamma distribution prior on , , with a prior mean of 5. This choice is based on the findings by Li et al that 86% of all infections were undocumented at the beginning of the epidemic in China.

Next, we place a uniform distribution prior on the proportion of exposed individuals among . That is, we assume Lastly, we have , where is the (known) population size.

Model for the time‐varying disease transmission rate with covariates

We now turn to the modeling of the epidemiological parameters. We start with the model construction for the time‐varying disease transmission rate between the undetected infectious and susceptible individuals, , which is an important parameter that characterizes the speed of disease transmission. We model the transmission rate as time‐varying to account for changes in COVID infection rates due to human behavior, , enactment of government policies, , evolution of new viral strains, and other time‐sensitive factors. Since , it is easier to work with the log transformation of . A simple and flexible way of modeling is through temporal splines, where is the regression intercept term, is a vector of basis functions evaluated at time (excluding the intercept), and is the vector of regression coefficients. We primarily consider a piecewise‐linear spline for with knots selected by quantiles, but other forms (piecewise constant, natural cubic splines, etc.) are possible in this framework. The intercept term and regression coefficients are further modeled through normal priors, and with . Here, denotes a normal distribution with mean 0.5 and standard deviation 0.1 restricted to .

The framework in (12) makes it straightforward to incorporate additional covariates for modeling . Let denote covariates. Motivated by the hybrid of single index models and additive models, we further generalize (12) as In Section 4, we will include the information on human mobility via the terms .

The transmission rate between the detected infectious and susceptible individuals is assumed to be reduced by a factor of due to potential quarantine and self‐isolation. We place a uniform distribution prior on , .

Model for the other epidemiological parameters

In this section, we will detail the model and prior specification for the remaining epidemiological parameters. Most parameters in model (1) to (8) have practical implications corresponding to clinical characteristics of COVID‐19. Therefore, we elicit informative priors for these parameters by summarizing the findings in the literature. , , , , , , , , We acknowledge that there is not yet a consensus on all clinical characteristics of COVID‐19. For example, Li et al estimated a mean incubation period of 5.2 days, while the same quantity was estimated to be 7.75 and 9.1 days in Jiang et al and Deng et al, respectively. As a result, our priors are chosen to be consistent with the majority of the literature. Our general methodology works for any choices of priors. We also note that the timing of COVID‐19 events is highly variable at the individual level, but the parameters in (1) to (8) are defined at the population level, which are less variable. The informative priors improve the interpretability of the results by assigning larger prior mass to the parameters around clinically meaningful values. In addition, since our model is an eight‐compartment model with only three compartments observed, some model parameters are nonidentifiable. Such nonidentifiablility is mitigated through the use of informative priors. ,

In the following paragraphs, we detail the interpretation and modeling choices for the remaining parameters, which are summarized in Tables 1 and 2.

TABLE 1

Model parameters, interpretations, and prior distributions

Parameter	Interpretation	Prior
E(0)	Initial number of exposed individuals	E(0)E(0)+I^u(0)∼Unif(0,1)
τ	Reduction in transmission rate for detected infectious individuals	Unif(0,1)
ϕ_C	Overdispersion parameter for daily cases	Ga(50,100)
ϕ_H	Overdispersion parameter for daily hospitalizations	Ga(50,100)
ϕ_D	Overdispersion parameter for daily deaths	Ga(50,100)
ζ_β	Intercept for disease transmission rate (log scale)	N(0,1)
ξ_β, ξ_βθ	Regression coefficients for disease transmission rate (log scale)	N(0,σβ2I)
σ_β	Standard deviation of regression coefficients for disease transmission rate	N_+(0.5,0.1^2)
a_ψ	Stabilized detection fraction	Be(55.5,18.5)
b_ψ	a_ψ−b_ψ is the initial detection fraction	b_ψ/a_ψ∼Unif(0,1)
c_ψ	Speed of detection fraction increase	Ga(5,100)
ζ_ν	Intercept for hospitalization fraction (logit scale)	N(−2,1)
ξ_ν	Regression coefficients for hospitalization fraction (logit scale)	N(0,σν2I)
σ_ν	Standard deviation of regression coefficients for hospitalization fraction	N_+(0.5,0.1^2)
ζ_ω	Intercept for death fraction (logit scale)	N(−1,1)
ξ_ω	Regression coefficients for death fraction (logit scale)	N(0,σω2I)
σ_ω	Standard deviation of regression coefficients for death fraction	N_+(0.5,0.1^2)

TABLE 2

Model parameters with informative priors elicited based on clinical characteristics of COVID‐19

Parameter	Interpretation	Prior
κ	Multiplicative factor of initial exposed and undetected individuals relative to detected	Ga(25,5)
α	Inverse of latent period	Be(31.5,58.5)
η_0	Inverse of infectious period before detection (for detected individuals)	Be(32.4,58.4)
λ_0	Inverse of infectious period between detection and hospitalization (for detected individuals)	Be(82.0, 163.9)
γ_I0	Inverse of infectious period after detection (for detected individuals)	Be(19.3, 115.8)
γ_H0	Inverse of infectious period after hospitalization (for detected individuals)	γH0−1+λ0−1=γI0−1
ρ_0	Inverse of whole infectious period (for undetected individuals)	ρ0−1=η0−1+γI0−1
δ_0	Inverse of time to death after the end of infectious period (for detected individuals)	Be(29.1,494.6)

Latent period

The parameter denotes the rate of exposed individuals becoming infectious, and represents the latent/preinfectious period, that is, the time period between exposure to the disease and being able to infect others. Note that the latent period in our paper is different from the incubation period, where the latter refers to the time period between exposure and symptom onset (see Figure 2). For COVID‐19, it is well recognized that patients usually become infectious before the onset of symptoms. We place a beta distribution prior on , , with a prior mean of and standard deviation of 0.05.

FIGURE 2

Natural history of COVID‐19, definitions of terms, and corresponding parameters

Detection rate

The parameter represents the rate of undetected infectious individuals being detected as having the disease at time . We model by where represents the proportion of undetected infectious individuals on day that are eventually diagnosed with the disease, and is a parameter such that represents the time interval between when an individual is capable of infecting others and when the individual is detected as having the disease. In practice, detection may occur after the onset of symptoms (eg, a person feels sick a day or two and then goes to get a test) or earlier than the start of symptoms (eg, a group like a college athletic organization getting mandatory regular tests that identify presymptomatic or asymptomatic cases). Without additional information, we assume that the average time to detection for undetected infectious individuals is similar to the time to symptoms onset. In other words, on average, undetected infectious individuals (who are eventually diagnosed) are diagnosed with the disease when their symptoms start to appear. Therefore, is roughly equal to the incubation period. We place a beta distribution prior on , , with a prior mean of and standard deviation of 0.05. Next, we model with a three‐parameter curve, where , and . As a result, starts from at time , monotonically increases with increasing , and converges to as . The underlying assumption is that at the beginning of the pandemic, due to the limited testing capacity, the proportion of detection starts from a low level . As time progresses, testing capacity increases monotonically, and eventually, anyone who wants a test can get it, making the proportion of detection stabilize at . The parameter characterizes the speed of testing capacity increase. We assume with a prior mean of 0.75 and standard deviation of 0.05. Further, we model the ratio between and and place a uniform distribution prior on it, . Lastly, we assume .

Infectious period and hospitalization rate.

The parameter denotes the rate of detected infectious individuals becoming hospitalized at time . We model by , where controls the proportion of detected infectious individuals on day who will eventually be hospitalized, and represents the time interval between when an individual is detected as having the disease and when the individual is hospitalized. We place a beta distribution prior with mean 1/3 and standard deviation 0.03 on : . Next, we consider the logit transformation of , . We model the transformed with a linear B‐spline, Here, is the intercept term, is a vector of basis functions evaluated at time (excluding the intercept), and is the vector of regression coefficients. The intercept term and regression coefficients are given normal distribution priors, and with .

The parameter denotes the rate of detected infectious individuals becoming noninfectious without being hospitalized. We have , where is a parameter such that represents the time interval between disease detection and when an individual is no longer able to infect others (Figure 2). We place a beta distribution prior with mean 1/7 and standard deviation 0.03 on : .

The parameter denotes the rate of undetected infectious individuals becoming noninfectious without ever being detected as having the disease. We have Here, is the proportion of undetected infectious individuals on day that are never diagnosed with the disease, and is the time interval during which a never‐detected infectious individual is capable of infecting others, that is, the whole infectious period.

By the definitions of and , we have In other words, the whole infectious period equals to the time interval between when an individual is capable of infecting others and when the individual is diagnosed with the disease, plus the time interval between diagnosis and when the individual is no longer able to infect others.

Death rates

The parameter is the rate of death for hospitalized individuals at time (recall that we assume all individuals who die will pass through the hospitalized compartment). We assume that hospitalized individuals who will not die become no longer infectious with rate at time . We have where represents the fraction of deaths among hospitalized individuals at time , and accordingly, represents the corresponding fraction of recoveries. By constructing and as a function of the time‐invariant and and the time‐varying in (17), we can easily estimate a time‐varying recovery rate and time‐varying death rate without introducing additional identifiability concerns.

We place beta distribution priors on . Specifically, with a prior mean of and standard deviation of 0.01. We set such that . The fraction of deaths , and therefore, we consider the logit transformation of , . We model the transformed with a linear B‐spline, Here, is the intercept term, is a vector of basis functions evaluated at time (excluding the intercept), and is the vector of regression coefficients. The intercept term and regression coefficients are given normal distribution priors, and with .

Implementation

We use the Stan modeling framework to sample from the posterior distribution of all parameters. Stan uses a Hamiltonian Monte Carlo procedure to efficiently sample from the posterior distribution. A major advantage to using the Stan framework is that it can easily accommodate nonconjugate and truncated prior distributions, which allows us for more flexibility in choosing a prior distribution that is scientifically relevant for each parameter. We interface with the model via the rstan package in R.

ANALYSIS OF THE COVID‐19 DATA

Data collection

We apply this model to county‐level data in the U.S. state of Colorado (Section 4.2) and state‐level data in the United States (Section 4.3) in the regions highlighted in Figure 3. Daily counts of cases and deaths were obtained from The New York Times.

FIGURE 3

Study locations (in blue) for analysis of county‐level data in Colorado (A) and state‐level analysis in the United States (B)

We collected mobility data from SafeGraph via the covidcast R package, which uses anonymized location data from mobile phones to generate different views of mobility over time. We use three mobility metrics from this data source—the fraction of mobile devices that did not leave the immediate area of their home in each day (“completely home”), the fraction of mobile devices that spent more than 6 hours at a location other than their home during the daytime (“full‐time work”), and the number of daily visits made by those with SafeGraph's apps to restaurants in a county. Each of these metrics are available daily, and we applied a kernel smoother to obtained smoothed values of each metric. To better capture changes in movement patterns, the smoothing for completely home and full‐time work variables used only weekday data and the smoothing for restaurants used only data from Friday through Sunday. An example of these mobility data for Colorado is shown in Figure 4, along with the dates of some key statewide policies related to COVID‐19 control. The mobility data for all regions is shown in Figures S1 to S3 and S19 to S21.

FIGURE 4

Example of mobility data used in analysis. Measurements from SafeGraph were used to characterize mobility patterns, including the fraction of mobile devices that were completely at home each day, the fraction that spent 6 hours or more away from home (“full‐time work”), and the number of daily visits to restaurants. This figure shows these data for the state of Colorado. The gray curves represent the daily values and the red curve represents the smoothed values used for modeling. Dates of several key statewide policies, as well as weekly new cases of COVID‐19 in the state, are also shown

The number of total (adult and pediatric) COVID‐related hospitalizations was obtained from U.S. Health and Human Services data via the covidcast R package. For states, these data are available one a daily basis beginning on July 14, 2021. For counties, these data are available on a weekly basis beginning July 31, 2021. The weekly hospitalization were disaggregated into seven integer‐valued daily values that sum to the weekly count.

Models for counties in Colorado

Model setup

We first apply this model to data from counties in Colorado. The first case of COVID‐19 in Colorado was reported on March 5, 2020. On March 10, 2020, the Colorado Governor issued a Declaration of Disaster Emergency and on March 27, 2020, a statewide “Stay at Home” order was given (Figure 4). Although some of the early outbreaks in Colorado were located in resort areas that have a small resident population, COVID‐19 quickly spread throughout the state. On April 27, 2020, the state transitioned to a “Safer at Home” order that was later amended to county‐specific regulations related to measures of local pandemic risk.

For modeling, we selected the ten largest counties (by population) in Colorado (Figure 3 and Table S1). This excludes some of the counties with ski resort communities where the earliest outbreaks occurred, but represents the vast majority of the state's population and communities varying from urban to rural. For each county, we select as the first day with six or more detected cases in that county and use each day as the model time‐step. The number of initial detected cases () was taken as the total number of reported cases in the five days before the modeling start date. We model the dynamics through January 31, 2021, at which point the widespread introduction of vaccines fundamentally changed the underlying dynamics of disease spread.

All models are fit with the priors specified in Tables 1, and we set . In practice, we found that the prior distributions in Table 2 were not always informative enough for regions with a large number of cases and may lead to clinically implausible parameter estimates. Therefore, we scaled the priors in Table 2 based on the cumulative number of detected cases. Specifically, we scale the prior variances of , , , , , and based on the ratio of the total number of cases of the county and that of Pueblo county. If the ratio is , the prior standard deviation is scaled to of the default prior standard deviation in Table 2.

For , we used linear B‐splines with 10 degrees of freedom (df) as the basis functions . In addition, we included the three smoothed mobility measures in the regression component as in (13). To improve interpretability of the mobility measure coefficients, we projected the splines for onto the orthogonal complement of the space spanned by the mobility time series. This follows approaches common in the spatial literature , for avoiding confounding from collinear trends. For and , we used linear B‐splines with 5 degrees of freedom. In model fitting for all counties, we used four chains each with 1000 warmup iterations and 1500 postwarmup iterations for a total of 6000 postwarmup iterations. We set the maximum tree depth to 14 and adapt_delta to 0.95. Convergence was assessed using the metric provided by rstan, specifically achieving values less than or equal to 1.01, and through visual inspection of trace plots for key parameters.

As a sensitivity analysis, we additionally fit models with two or four times larger prior standard deviation for the clinical parameters (, , , , and ). We also fit models with 5 and 20 degrees of freedom for the spline, instead of 10 degree of freedom.

County model results

The time‐varying for each Colorado county is displayed in Figure 5. There is a similar tri‐modal trend across all counties: peaks in in late March (at the start of modeling), July, and November. These trends match the rise and fall of new cases, shown in Figure 6, with the largest peak in cases happening in November and December in all counties. The dynamics of some counties are slightly different, such as additional peaks in in early Fall 2020.

FIGURE 5

Posterior mean (solid line) and 95% credible interval (dashed line) of for Colorado county models. Horizontal line represents

FIGURE 6

New cases by date in Colorado counties. Points show observed counts, line shows average model predictions (ie, mean of the posterior predictive distribution)

The impact of the stay‐at‐home orders is clearly evident in the decline in throughout April and the resulting downturn of new cases in late April and early May. The equilibrium reproductive number remained near 1 until a small rise in early July. In early September, the began to increase in the state overall, with a particularly large rise in Boulder county, where there were several high‐profile outbreaks among college‐aged individuals (see corresponding panel in Figure 6). As the number of new daily cases decreased in Boulder county, the estimated declined again.

The predicted and observed number of daily hospitalizations and deaths are provided in Figures S4 and S5, respectively. On most days, no deaths were reported, but the predicted number of new deaths clearly shows the April and November peaks in mortality. The predicted hospitalizations followed the observed counts in most counties, with some notable deviations in Douglas county.

Posterior means for other parameters are shown in Figure S6 to S16. The proportion of detected infected individuals for each day is shown in Figure S17. In all counties, this proportion starts small but increases through April and May and remains elevated throughout the rest of the pandemic, reflecting increased testing capacity.

Results from sensitivity analyses are shown in Section B.6 of Appendix S1. Increasing the prior variances for clinical parameters had minimal impact on model fit. In the Denver county model, changing the degrees of freedom did not have a meaningful impact on the predicted cases (Figure S40) or (Figure S42). But for Boulder county, the model with increasing flexibility yielded large spikes in (Figure S43) while the model with less flexibility in the splines greatly oversmoothed the trend in predicted cases (Figure S41). These results suggest that the choice of 10 df for the splines in provides a good amount of flexibility while reducing overfitting.

Impact of mobility on within‐county transmission

County‐wide mobility played an important role in modeling the transmission rate () in most counties. The effects of mobility on transmission are summarized in four ways. Figure 7 shows the correlations between the mobility data time series and , that is, . Figure S38 shows the posterior means (and 95% credible intervals) for the coefficients for each mobility term. Figure S39 shows the posterior correlation between the mobility coefficients () and , that is, . Figure S18 shows the combined time‐varying effect of mobility on infectiousness (ie, ).

FIGURE 7

Correlation of mobility data time series with the time‐varying transmission parameter . Points and horizontal error bars represent posterior means and 95% credible intervals, respectively

In all counties, there is a clear positive correlation between the proportion of people working full time and disease transmission (Figure 7). Furthermore, the coefficient for the proportion of people working full time is positive in almost all counties (Figure S38), indicating a positive association between this measure of mobility and COVID spread, and in all except Boulder, Larimer, Pueblo, and Weld counties the credible interval excludes zero. This clearly demonstrates how changing patterns in work trends (data shown in Figure S2) are associated with transmission of COVID‐19. This reflects not only the sharp reduction in the proportion of people working full time after the stay‐at‐home orders at the start of the pandemic in March 2020, but also the gradual increase in mobility due to more people working away from home after the April 27 “Safer‐at‐Home” order that relaxed some restrictions through the summer and fall (Figures 4 and S2). The counties with the largest correlations and coefficients are all in the metropolitan areas surrounding the economic centers of Denver, Boulder, and Colorado Springs. In addition, for most counties, we observe strong positive posterior correlations between the coefficient for full time work and .

There is a positive, but somewhat weaker, correlation between the number of people visiting restaurants and transmission. This is consistent with the trend seen with full time work: more people being active outside the home is correlated with greater transmission. However, the coefficient for the number of people visiting restaurants is heterogenous across counties (Figure S38). This must be viewed in light of concurrent adjustment for the proportion of individuals working from home, which is positively correlated with restaurant visits.

When concurrently adjusting for full‐time work and restaurant visits, we did not see meaningful relationships between the proportion of individuals staying at home and transmission (credible intervals in Figure 7 include zero), except Jefferson County, which has a negative correlation. It is important to note that this does not mean that stay‐at‐home orders in these counties were ineffective, but rather the mobility measures of people outside the home (at work and restaurants) were more predictive of transmission than the numbers of people staying home. Similar to restaurant visits, the coefficient for the proportion of people staying at home is heterogenous across counties, due to its correlation with the other mobility measures (Figure S38).

Models for U.S. states

Data and model setup

We also fit the model to nine U.S. states that are representative for their different trends in the number of cases throughout the pandemic. The selected states were (Figure 3): Colorado, South Dakota, and Wisconsin, which had small early waves and a large late wave; California, Florida, and Texas, which had large numbers of cases mid‐summer and a late wave; Iowa, which had moderate numbers followed by a late wave; Georgia, which had relatively small summer and fall peaks; and New York, which had a very large early peak in the spring and a late wave in November.

For these models, we chose the initial time point to be the first day after the state‐level case count exceeded 100 cases (Table S2) and modeled dynamics through January 31, 2021. Prior specification and model structure was the same as in the county‐level models described in Section 4.2.1. For the state‐level models, we included the same three measures of mobility, but calculated at the intra‐state level (Figures S19‐S21), in the term. For the clinical parameters, the prior standard deviations were scaled against the number of cases in South Dakota.

State model results

The value of for each state is displayed in Figure 8 and the observed and predicted numbers of new cases for each state are plotted in Figure S22. There is a similar tri‐modal trend across all states: peaks in in April, July, and November to December. New York began with a very large , that precipitously declined as the epidemic in New York City was brought under control in April 2020. After the initial spread of COVID‐19 spurred a round of lockdown measures, most states had values around 1 before later local peaks in lead to additional waves of COVID‐19. California, Florida, and Texas all had increases in in late June that preceded large numbers of cases in July (Figure S22). Meanwhile, South Dakota had an above 1 from July through October, which led to continuous rate of growth in new cases during that time (Figure S22). Consistent with the county‐level models, the full‐state Colorado model had its largest peak in November.

FIGURE 8

Posterior mean (solid line) and 95% credible interval (dashed line) of in state models. Horizontal line represents

Similar to the county results, the sensitivity analyses showed minimal impact of increasing prior variance for clinical parameters and considerable oversmoothing when the degrees of freedom in was reduced (Figures S44‐S47).

Impact of within‐state mobility

For some states, including the overall Colorado model, there is a strong relationship between one or more of the mobility factors and the transmission parameter (Figure 7). The models for California, Colorado, and New York have correlation for the full‐time work mobility that are positive, indicating that the transmission rates are positively associated with the proportion of individuals being away from home more than six hours per day. The large correlation and coefficient for New York reflects the link between the drop in mobility in March and April 2020 (see Figure S37) and the drop in cases after an initial large peak during that time. Colorado's positive correlation and coefficient value also reflects a link between a drop in mobility and a drop in cases in April 2020. Although the credible interval for the coefficient for full‐time work in California includes zero (Figure S38), the correlation between that measure of mobility and transmission is strongly positive (Figure 7). Similar to the Colorado county models, the coefficient for full‐time mobility in these three states was highly positively correlated with (Figure S39). The impact of people working away from home full time was limited in the other state models, as evidenced by the credible intervals including zero for correlations with transmission (Figure 7). It is possible that the mobility data provide a better reflection of potential for contacts in some states compared to others. In a state with heavy tourist traffic, for example, the mobility data for working full‐time or at home would not capture tourist patterns, although the restaurant visits would to some degree. This may partially explain why a weaker association was observed in states like Florida and South Dakota.

The impact of restaurant mobility was heterogeneous across states. In California and Colorado, the coefficient was not meaningfully different from zero (Figure S38) but restaurant mobility was positively correlated with transmission (Figure 7). Restaurant mobility was also positively correlated with transmission in New York, where the coefficient was negative. The models for Florida, Georgia, and Texas also yielded negative posterior means for the coefficients of restaurant mobility. However, this must be viewed in conjunction with the negative posterior means for the coefficient of the proportion of individuals completely home (for which a negative value means less disease spread when more people are staying home). The net result was a minimal effect of mobility on disease spread in those particular states (Figure S37).

Overall, the correlations between the time‐varying transmission rate and the proportion of people completely home were stronger among states than in the Colorado county models. This was particularly true in California and New York, where the increase in the proportion of individuals working from home was particularly large (Figure S19).

In Iowa and Texas, the correlation between transmission and mobility measures was opposite from what was seen in the majority of models. In these two states, more people working from home was positively correlated with transmission and more people visiting restaurants was negatively correlated with transmission (Figure 7). These surprising correlations may be driven in part by the reduction in mobility around the time of the December transmission waves (compare Figure 8 and Figure S20). Or it may be that differences in policies (such as mask‐wearing) mean that mobility alone is not representative of transmission risk in these states.

CONCLUSIONS AND DISCUSSIONS

We have presented a flexible model for a time‐varying infectious disease that includes: temporal variation in key model parameters through the use of splines and a regression term, a multiplicative state‐space process model for allowing for day‐to‐day heterogeneity in disease spread, and overdispersion in the observed number of cases. Together, these elements allow our model to capture important features of data observed in a pandemic.

The inclusion of the mobility terms in the infectiousness parameter was of particular interest to us. We observed that the mobility measures were correlated with transmission in many counties and regions that we modeled (Figure 7). However, it did not play a role in the modeled transmission in all regions. The heterogeneous impact of mobility data may reflect the importance of other factors, such as compliance with mask‐wearing mandates. While derived from individual‐level information on movement, the mobility data we incorporated into the models was averaged at the county and state levels. This relatively coarse spatial scale means that we cannot differentially model the connection between disease transmission and mobility in at‐risk or highly active subpopulations. Individual‐level information on movement could improve the strength of the evidence between mobility and COVID transmission. Furthermore, a potential downside to using cell phone movement as a proxy for human mobility is that there are disparities in cell phone ownership and usage among groups of people. In particular, older Americans are less likely to own a cell phone, and thus may be poorly represented by this form of mobility data. Similarly, cell phone use may also be affected by socioeconomic status and not capture institutionalized populations well. Coston et al have shown the cell phone mobility data to have low coverage of these populations, which leads to bias in the mobility metrics. Despite these limitations of the mobility data, we saw a clear effect of reduced mobility in the sharp reduction in cases early in the pandemic (Figure 5). Furthermore, higher rates of working outside the home were associated with increased transmission in many counties and states.

The multiplicative process model introduces extra flexibility to the compartmental model, which better reflects fluctuations in the rates of detection and transmission. In addition to nonrandom mixing of the population, a key violation of the traditional SIR model assumptions is that each individual is equally infectious. There is now considerable evidence , , , that superspreaders play an important role in the spread of COVID‐19, and the flexibility afforded by the process model can incorporate this heterogeneity in spread by allowing for a multiplicative shift to the solution to differential equations of the compartmental model at each time point. This has the additional advantage of allowing the process to explicitly account for heterogeneity in case reporting, such as the tendency of many states to report fewer cases on weekends.

Our approach also presents advantages over approaches that only compare time series of cases and other factors to identify marginal relationships. , , By incorporating information on mobility into a model for disease spread, our model can more accurately represent a potentially causal role of factors in disease spread. This further allows for comparison of mobility information with the latent rate of disease transmission, rather than just the number of reported cases.

One general challenge in the fitting and interpretation of the compartment model is the limited data for the number of individuals at each stage of the dynamics. For our model, although the framework presented in Section 2 contains eight compartments, we only observe data from three of those compartments (, , and ). Several parameters influence the number of new cases , which can lead to some identifiability issues. For example, a moderate number of identified cases could occur due to a low detection (small ) among a large number of undetected infectious individuals (large , arising from large ), or it could be due to high detection (large ) among a small number of infectious individuals (smaller , arising from smaller ). To mitigate this, we introduced a nondecreasing structure for . Nonetheless, it remains nontrivial to accurately estimate the true number of undetected cases at any given time and the results of Figure S17 should be interpreted cautiously. Inclusion of information on testing rates could be used to partially address this by providing a basis for , but quality and quantity of information on testing rates and strategies vary widely between and within jurisdictions, which deserves more modeling efforts in the future.

Several extensions are possible to the models we presented here. The additive model framework for the time‐varying parameters could be expanded to include environmental factors such as temperature, humidity, and air pollution. There has been suggestive evidence , , , , , , , that these factors may influence spread and severity of symptoms. As noted above, more detailed mobility information could be incorporated to represent small‐scale movement patterns. Although we used the same prior distributions and model settings across all counties and states, a natural extension of this approach would be to fit a joint model that includes all regions (counties or states) together. This would allow for information on the non‐spatial parameters (eg, ) to be estimated using shared information. If available, area‐specific migration information could be included to account for movement between counties or states, such as college students moving at the start of an academic year. However, the major universities in Colorado were primarily conducting remote classes during the time period of modeling, and so this likely would not have a major impact on the results for the counties modeled here. Direct movement of susceptible and infectious cases between regions could then be included as well. The primary drawback to such an approach is the computational complexity that arises from having separate spline parameters in each region.

Overall, the proposed model provides a rigorous, process‐driven framework for modeling the impact of time‐varying factors on infectious disease spread. Our analyses showed an important role of mobility in the spread of COVID‐19 in several Colorado counties and U.S. states.

Appendix S1. Supporting information

Click here for additional data file.