Early-Phase Local-Area Model for Pandemics Using Limited Data: A SARS-CoV-2 Application

Building on Equation (1.1), for counties $i=1,\cdots,n$, and time points $t=1,\cdots,N$, we posit moment constraints on the number of new infections on day $t$ in location $i$ given previously reported existing infections prior to day $t$, as:

where $\mathcal{F}_{i,t}$ is the filtration of past incident cases information, $\phi_{i}$ is a dispersion parameter, $g(\cdot)$ is a known variance function, and $\mu_{it}=R_{it}\Lambda_{it}$ when $R_{it}\in\mathcal{F}_{i,t-1}$. When $g(\cdot)$ is unknown, smoothing techniques, e.g. local polynomial fitting, can be used to estimate the variance function (Chiou and Müller, 1998). More discussions on unknown variance functions can be found in Section A of the Supplementary Materials.

As previously noted, $\Lambda_{it}\triangleq\sum_{s=1}^{t}I_{i,t-s}\omega_{s}$ represents the infection potential, determined by the current number of infectious individuals and the infectiousness function $\omega_{s}$. We assume that the $\omega_{s}$’s represent the probability distribution of the serial interval or generation time, and set $\omega_{s}=0$ when $s=0$ or $s>\eta$, with $\eta$ denoting the duration from infection to either recovery or mandatory quarantine. We use $\omega=\{\omega_{s},1\leq s\leq\eta\}$ to denote the vector of $\omega_{s}$. In practice, $\omega$ is a vector of positive values, often obtained from epidemiological studies of infectious pathogens (He et al., 2020). In cases where multiple studies on serial intervals exist, a meta-analysis can be conducted to obtain a pooled estimate.

To estimate $R_{it}$ and investigate the impact of local area factors on disease transmission, we further assume a time series model as follows:

where $D_{1,it}=\big{\{}X_{ij}\big{\}}_{0\leq j\leq t}$, $D_{2,it}=\big{\{}I_{ij}\big{\}}_{0\leq j\leq t-1}$, $D_{3,it}=\big{\{}R_{ij}\big{\}}_{0\leq j\leq t-1}$, and $\mathcal{F}_{i,t-1}=\sigma(D_{1,it}\cup D_{2,it})$. The term $\sum_{m=1}^{q}\theta_{m}f_{m}(D_{1,it},D_{2,it},D_{3,it})$ describes the time-series dependency of the disease transmission rate. Specifically, $f_{m}(\cdot)$ is a function of past data and $\theta_{m}\geq 0$ for $1\leq m<q$. $f_{q}\equiv 1$ is set as the constant function and $\theta_{q}$ is the intercept term. The $h(\cdot)$ is a known link function. Moreover, $W_{it}\in\mathbb{R}^{p_{1}}$ and $X_{it}\in\mathbb{R}^{p-p_{1}}$ are vectors of local-area factors for location $i$ on day $t$ associated with disease transmission. The factors in $W_{it}^{T}$ have fixed location-specific effects, while the factors in $X_{it}^{T}$ have common effects across locations. To simplify notion, we consolidate the presentation of all local-area factors as $Z_{it}=(X_{it}^{T},\overbrace{0,\cdots,0}^{p(i-1)},W_{it}^{T},\overbrace{0,\cdots,0}^{p(n-i)})^{T}$, present the vector of all regression coefficients of these factor as $\beta=(\check{\beta}_{0}^{T},\check{\beta}_{1}^{T},\cdots,\check{\beta}_{n}^{T})^{T}\in\mathbb{R}^{p}$, and denote the vector of all other regression coefficients as $\theta=(\theta_{1},\cdots,\theta_{q})^{T}\in\mathbb{R}^{q}$.

In practice, a choice for Equation (2.2) can be a non-stationary log-transmission model with an autoregressive process as follows:

In the above models, to counter potential biases arising from data errors and incorrectly specified model assumptions, we incorporate $\epsilon_{it}$ as a measurement error term for $R_{it}$. This term can also be viewed as a random effect term of $R_{it}$, capturing unexplained deviations in $R_{it}$ that could result from inaccuracies in reported infection data or model mis-specifications. At this point, our proposed model breaks down the serial dependency of incidence cases into two components: Equation (2.1) encapsulates the serial dependency inherent in the generative nature of disease transmission given $R_{it}$, and Equation (2.2) portrays the serial dependency arising from the $R_{it}$ interdependence between adjacent time points, elucidated by covariates and time series structures.

With the measurement error term in Equation (2.2), the estimation of $R_{it}$ and parameters in the proposed model becomes difficult. In a simple scenario when $\epsilon_{it}$ is assumed to be 0, Equation (2.2) reduces to

In this situation, an estimator of the parameters $\gamma=(\beta^{T},\theta^{T})^{T}$, denoted as $\hat{\gamma}$, can be obtained by solving $U_{N}({\gamma})=0$, where $U_{N}(\gamma)$ is the quasi-score function

However, when $\epsilon_{it}$ is a non-degenerate random variable, calculating the expectation in $\mu_{it}=\mathbbm{E}[R_{it}|\mathcal{F}_{i,t-1}]\Lambda_{it}$ and solving $U_{N}(\gamma)=0$ are generally difficult. Similar models were studied by Davis, Dunsmuir and Streett (2003) to model count data using the log link function, where the measurement error was assumed to be a stationary process and calculated using an autoregressive moving average recursion with a distributed lag structure. Here, given that $\{R_{it}\}_{1\leq t\leq N}^{1\leq i\leq n}$ are not directly observed, we describe an iterative algorithm for estimating parameters in (2.2) without requiring any stationary structure and distribution assumption on the measurement error term. Briefly, the proposed iterative algorithm consists of two steps. The first step uses the second layer of the model, the time series model, to construct a semi-parametric, locally efficient estimator of $\beta$ based on a location-shift regression model. In this step, $\beta$ and $R_{it}$ can be written as functions of $\theta$ using an initial estimate of $\theta$ or an estimate from the last iteration. In the second step, the parameters were updated using the first layer of the model, the quasi-score function. Specifically, for $t=q,\cdots,N$, $i=1,\cdots,n$, we first rearrange Equation (2.2) to

We use the block bootstrap method introduced by Hall (1985) and Künsch (1989) to quantify the uncertainty of the proposed online estimator. We will illustrate this calculation for $\hat{\beta}$. The confidence band of $\{\hat{R}_{it}\}_{1\leq t\leq N}^{1\leq i\leq n}$ can be obtained similarly. Given the time series data of daily case incidence and local area factors, denoted as$\{(I_{it},Z_{it})\}_{1\leq t\leq N,1\leq i\leq n}$, where $Z_{it}$ is a $p$-dimensional vector of covariates at time $t$ in location $i$, we assume that $\beta$ can be estimated using a block of the time series data as long as the length of the block exceeds $\ell$. Specifically, we require that the expression $\left(\sum_{i}\sum_{t\in\mathcal{S}_{i}}(Z_{it}-\bar{Z}_{\mathcal{S}})(Z_{it}-\bar{Z}_{\mathcal{S}})^{T}\right)^{-1}$ be well defined for any arbitrary set $\mathcal{S}_{i}\subset\{1,\cdots,N\}$ with $|\mathcal{S}_{i}|\geq\ell$, $i=1,\cdots,n$ and $\bar{Z}_{\mathcal{S}}$ is defined as an analogy to $\bar{Z}$.

For $B$ bootstrapping samples, we assume the sampled block $\mathcal{S}_{b}$ is the same for all counties in one bootstrapping observation for simplicity, $b=1,\cdots,B$. That is, we denote a block of data up to time $t$ with length $\ell$ as

for $t=\ell,\cdots,N$. We adopt the suggested $\ell=O(N^{1/3})$ from Bühlmann and Künsch (1999) and independently resample $B$ blocks of data with replacement, for $B=O(\lfloor(N-p-\tau_{0})/\ell\rfloor)$ (Bühlmann, 2002) to obtain $B$ bootstrap samples $\xi_{t_{1}},\xi_{t_{2}},\ldots,\xi_{t_{B}}$.

For the effect of $j$-th covariate, $\beta_{j}$, where $j=1,\cdots,p$, let’s denote the estimator based on the $b$-th sample as $\hat{\beta}^{*}_{t_{b},j}$ and the estimator using the whole data as $\hat{\beta_{j}}$. One may use $\{\hat{\beta}^{*}_{t_{b},j}-\hat{\beta}_{j}:b=1,\cdots,B\}$ to approximate the empirical distribution of $(\hat{\beta}_{j}-\beta_{j})$ or use $\{\hat{\beta}^{*}_{t_{b},j}:b=1,\cdots,B\}$ as an approximation of the empirical distribution of $\hat{\beta}_{j}$. Thus, the level $1-\alpha$ bootstrap confidence interval of $\beta_{j}$ based on the $B$ samples, denoted as $(L_{j,\alpha/2,B}^{*},U_{j,\alpha/2,B}^{*})$, can be constructed by

The proposed method is a general approach for various counting processes and is robust to mis-specification of distribution assumptions. When the measurement error can be ignored, i.e., $\epsilon_{it}=0$, the Poisson process used in (Cori et al., 2013) becomes a special case of Equation (2.1), with $\nu_{it}=\mu_{it}$. Since $U_{t}(\gamma)\in\mathcal{F}_{t}\triangleq\sigma(\cup_{i=1}^{n}\mathcal{F}_{i,t})$ and $\partial\mu_{it}/\partial\gamma,\mu_{it},\nu_{it}\in\mathcal{F}_{i,t-1}$, then ${(U_{s}(\gamma),\mathcal{F}_{s}),s\geq q}$ forms a mean zero martingale sequence, indicating that $U_{N}(\gamma)$ is an unbiased estimating function. For asymptotic properties of the estimators, we consider the least favorable case where we only have one county’s data and discard subscript $i$ in all notations. Denote the true parameter value as $\gamma_{0}$, and according to Theorem 3 of Kaufmann (1987), we have:

When the measurement error cannot be ignored, i.e., $\epsilon_{it}\neq 0$, the proposed iterative algorithm consists of two major steps: one that utilizes a semi-parametric locally efficient estimator to express $\beta$ and $R_{it}$ as functions of $\theta$, and another that updates the estimates by solving the score equation. However, the latter step is subject to estimation bias due to the non-zero measurement error term in equation (2.2), which can cause bias in the score equation (Cameron and Trivedi, 2013). Correcting this bias in the general case represented by (2.2) is challenging, but it can be tackled in the special case described by (2.3) by requiring $\mathbbm{E}\big{(}e^{\epsilon_{t}}|\mathcal{F}_{t-1}\big{)}=c_{0}$ for some constant $c_{0}$ (Zeger, 1988). Therefore, we can correct the bias by modifying the quasi-score equation (2.9) into an unbiased estimating equation as follows:

In practice, $c_{0}$ can be treated as a nuisance parameter and estimated by solving the quasi-score equation. Moreover, under Equation (2.3), solving the modified quasi-score Equation (2.9) is equivalent to obtaining

where $\hat{\theta}^{(k)}$ is the MLE based on conditional profile likelihood of Poisson counts. Therefore, the proposed online estimation procedure guarantees the following properties. The proof of the theorems is presented in Section C of the Supplementary Materials.

Theorem 3.3 shows that the bound of the step-wise difference of the online estimator between iterations decreases as the number of observation time increases. The consistency of the estimator remains an open question as discussed in Davis, Dunsmuir and Streett (2003). The difficulty arises from the measurement error, $\epsilon_{it}$, and the series dependency of disease incidence, $\Lambda_{it}$, in the model. Both elements are essential in modeling the dynamics of disease transmission. Similar but simplified models, which do not model the measurement errors and series dependency of disease incidence were studied for different contexts, e.g. studies of economics or disease without series dependency of disease incidence (Davis, Dunsmuir and Wang, 1999, 2000; Davis, Dunsmuir and Streett, 2003; Fokianos, Rahbek and Tjøstheim, 2009; Neumann et al., 2011; Doukhan, Fokianos and Li, 2012; Doukhan, Fokianos and Tjøstheim, 2012). Inspired by the derivation and discussion in Davis, Dunsmuir and Streett (2003), we conjecture that conditions requiring stationary and uniformly ergodic of $\log(R_{t}\Lambda_{t})$ are needed to establish the consistency of estimators of model (2.3), although these conditions may be difficult to verify in practice.

We generated data on the daily numbers of SARS-CoV-2 infections, along with associated local-area factors influencing the disease transmission rate. We generated data in a single location with time-varying factors and suppressed the location subscript $i$ in this section. We fit the proposed model, as described by Equation (2.3), assuming that the $R_{t}$ values exhibited an AR(1) dependency over time and were linked with an intercept term and two time-varying factors, i.e. $q=2$ and $p=2$. Each scenario was replicated 1,000 times to calculate evaluation statistics. The infectiousness function was assumed to be known and was modeled using the probability density function of a gamma distribution, in a manner mirroring the results of a previous epidemiological study of SARS-CoV-2 cases in Wuhan, China (Li et al., 2020). The estimation algorithm was applied with a pre-estimation period of 5 days, i.e., $\tau_{0}=5$.

Data were simulated under various scenarios, including settings where the fitted models were correctly specified or misspecified, as shown in Table 1. For the scenarios with correctly specified models (scenarios 1-7), data were generated using the parameters $(\theta_{1},\theta_{2},\beta_{1},\beta_{2})=(.7,.5,-.02,-.125)$. Here, the simulated values $\{R_{t},t\geq 1\}$ were designed to mimic the trend of the $R_{it}$ reported in Pan et al. (2020). We varied factors such as the days of observation ($T$), initial incident cases ($I_{0}$), and utilized different distributions of error terms ${\epsilon_{t}\sim_{i.i.d}\epsilon,t=1,\ldots,T}$ for data generation. In the settings where the model was mis-specified (scenarios 8-9), data were simulated with the assumption that the $R_{t}$ values followed an AR(2) dependency (i.e., $q=3$) or were associated with an intercept term and three covariates (i.e., $p=3$). The intercept term is $\theta_{q}$, while the two time-varying covariates ($Z_{t,1},Z_{t,2}$) were simulated independently to mimic the real data of temperature in Philadelphia and data of social distancing from daily cellular telephone movement, provided by Unacast (Unacast, 2021), which measured the percent change in visits to nonessential businesses, e.g., restaurants and hair salons during March 1st and June 30th, 2020. Specifically, daily temperatures were generated from a shifted normal distribution, i.e., $Z_{t1}\overset{d}{=}5+\Big{(}t-\frac{T}{2}\Big{)}\big{/}8+\mathcal{N}(0,9)$, for $t=1,\cdots,T$. The social distancing data were generated from a uniform distribution and applied a logit transformation afterwards, i.e., $Z_{t2}\overset{i.i.d}{\sim}2+logit(Uniform(0,1)),$ for $t=1,\cdots,T$. Additionally, in the misspecified scenarios with three covariates, the third covariate was generated from $Z_{t3}\overset{d}{=}\Big{(}t-\frac{T}{2}\Big{)}\big{/}9+\mathcal{N}(0,5)$, for $t=1,\cdots,T$. The variations of error terms were assumed to be small as the variation of $\log(R)$ was small with $R$ ranging from 0.5 to 3 and has a comparable size of the variation of covariates times their effects.

Table 2 presents the bias and coefficient of variation of the estimators for each scenario. In scenarios where the models were correctly specified, both the bias and the coefficient of variation decreased with an increase in either the days of observation (scenarios 1-3) or the initial case count (scenarios 4-5). When the variance of the error term increased to a magnitude smaller than the product of the time-varying local-area factors and their corresponding effect sizes (comparing scenarios 6 and 2), the proposed estimator still demonstrated robust performance, with less than a 5% increase in relative bias and a 0.2 increase in the coefficient of variation. Even when the error term exhibited a heavier tail (comparing scenarios 7 and 2), the relative bias increased by less than 4%, and the coefficient of variation increased by less than 0.01.

Two types of mis-specification were considered. In scenario 8, an important time-varying covariate was omitted from the model. Due to the correlation between the generated local-area factors over time, this scenario differs from merely increasing the error term’s variance. In both mis-specified cases, the estimator $\hat{\theta}$ is closer to the corresponding partial correlation coefficients rather than to the true values. Although the mis-specification led to an increase in estimation bias, the algorithm still provided a relatively accurate estimator for the correctly specified parts, showing a small bias and coefficient of variation (comparing scenarios 2, 8, and 9). Most importantly, as illustrated in Figure 1 from a random replication of the simulation, the estimated $R_{t}$ values still captured the tendency of the true $\{R_{t}\}_{t\geq 0}$, especially when the time-series sample size increased.

Two types of mis-specification were considered. In scenario 8, an important time-varying covariate was omitted from the model. Owing to the correlation between the generated local-area factors over time, this scenario represents more than just an increase in the error term’s variance. In both cases of mis-specification, the estimator $\hat{\theta}$ was closer to the corresponding partial correlation coefficients than to the true values. While the mis-specification led to increased estimation bias, the algorithm still provided a relatively accurate estimator for the correctly specified parts, exhibiting a small bias and coefficient of variation (as seen when comparing scenarios 2, 8, and 9). Most importantly, as illustrated in Figure 1 from a random replication of the simulation, the estimated $R_{t}$ values still accurately captured the trend of the true $\{R_{t}\}_{t\geq 0}$, particularly as the time-series sample size increased.

We also constructed the bootstrap confidence interval according to (2.12), and the coverage probability of the bootstrap confidence interval for $\theta_{2}$ is shown in Table 3 as an illustration. Table 3 confirms that the proposed bootstrap confidence interval provides the expected coverage probability when choosing $\ell=45=\tau_{0}+O(T^{1/3})$ for $T=300$, following the recommendations in Bühlmann and Künsch (1999), and setting the number of replications to $B=200=O(T^{2/3})$, as per Bühlmann (2002). When a smaller block length $\ell=30$ is chosen for $T=300$, the accuracy of the bootstrap estimator is reduced, leading to a conservative confidence interval. In contrast, when a larger block length $\ell=60$ is chosen, the bootstrap estimator’s accuracy may improve, but the coverage probability becomes less than ideal due to the high correlation between bootstrap subsamples. Furthermore, the optimal block length is also related to $T$, as $\ell=30$ and $\ell=45$ would only result in conservative confidence intervals for $T=450$. Practically, we suggest choosing the smaller value of $\ell$ that satisfies $\ell=O(N^{1/3})$ as recommended in Bühlmann and Künsch (1999).

The estimated covariate effects on $R_{it}$ with bootstrap intervals for the three time periods are shown in Table 4. Among the four local-area covariates we examined, social distancing appeared to be the most influential factor related to disease transmission, particularly during the first wave of the pandemic. It exhibited a significant positive association with disease transmission. During the first wave, a 50% reduction in the frequency of visiting non-essential businesses was estimated to reduce $R_{it}$ by an average of 11%. This finding aligns with the effects of social distancing reported in other studies from the first wave of the pandemic, but our study encompasses more counties across the US (Rubin et al., 2020; Courtemanche et al., 2020). Over time, the impact of social distancing diminished. This change could be attributed to multiple factors, including possible changes in the virus and variations in population behavior, which might have prevented our social distancing measure from capturing all aspects of population mobility. Furthermore, we observed a gradual decrease in social distancing over time, from a median of 31% reduction in the first wave to 11% reduction in the third wave. It’s also possible that the effect of social distancing is nonlinear, and its impact was minimal when only a small amount of social distancing was implemented.

We also observed a significant heterogeneity in disease transmission across counties with varying population densities during the first wave. The reproduction number was higher in more populated areas compared to rural areas. Simultaneously, we found that wet-bulb temperature was negatively associated with $R_{it}$, suggesting that transmission was generally higher in colder weather compared to warmer days. Although the regression coefficient of temperature was small, its impact on policy-making could be substantial due to the wide range of temperatures a county might experience within a year. During the first study period, the largest increment in wet-bulb temperatures observed within a county was 21^∘C, which could reduce COVID-19 transmission by 15%, with other covariates held constant. However, temperature could also affect social distancing values, as outdoor activities tend to increase during spring and summer, potentially offsetting the reduction in transmission due to increased temperatures. Both the heterogeneity associated with varying population densities and the impact of temperature decreased over time, indicating that the entire US began to experience a more uniform pandemic situation with less seasonality as we moved into the middle or later stages of the pandemic. Additionally, the effect of vaccination coverage on disease transmission during the first half of 2021 was not found to be significant. This lack of significance could be due to factors such as relatively lower vaccine coverage during early-stage vaccination rollout, behavioral changes post-vaccination, and challenges in achieving herd immunity.

We compared the projected daily case counts with the actual observed numbers for the three study periods, and the accuracy of the projection was assessed using the weekly Percentage Absolute Error (PAE). Mean and median PAE are shown in Table 4 by week and study period. The mean PAE ranges from 20.2% to 27.9% for the first two weeks, and 22.5% to 40.5% for the second two weeks, respectively. We also compared our projection results with the CDC COVID-19 Cases Ensemble Forecast, which combines forecasts submitted by a large and variable number of contributing teams using different modeling techniques and data sources (Ray et al., 2020; Cramer et al., 2022). During the three study periods, our model showed similar performance to the CDC model in the 1-week projection and was superior in the 3-4 week projections. The CDC model’s average PAE across all states and weeks was 22%, 32%, 44%, and 57% for the 1 to 4-week forecast windows, respectively (Du et al., 2022). It’s worth noting that the CDC ensemble model is designed for state-level projection while our model focuses on local-level, specifically county-level, projections. Projecting at the county level is more challenging due to the smaller case incidence numbers and higher noise-to-signal ratio in county-level data.