A flexible framework for local-level estimation of the effective reproductive number in geographic regions with sparse data

Abstract

Background

Our research focuses on local-level estimation of the effective reproductive number, which describes the transmissibility of an infectious disease and represents the average number of individuals one infectious person infects at a given time. The ability to accurately estimate the infectious disease reproductive number in geographically granular regions is critical for disaster planning and resource allocation. However, not all regions have sufficient infectious disease outcome data; this lack of data presents a significant challenge for accurate estimation.

Methods

To overcome this challenge, we propose a two-step approach that incorporates existing $\:{R}_{t}$ estimation procedures (EpiEstim, EpiFilter, EpiNow2) using data from geographic regions with sufficient data (step 1), into a covariate-adjusted Bayesian Integrated Nested Laplace Approximation (INLA) spatial model to predict $\:{R}_{t}$ in regions with sparse or missing data (step 2). Our flexible framework effectively allows us to implement any existing estimation procedure for $\:{R}_{t}$ in regions with coarse or entirely missing data. We perform external validation and a simulation study to evaluate the proposed method and assess its predictive performance.

Results

We applied our method to estimate $\:{R}_{t}\:$using data from South Carolina (SC) counties and ZIP codes during the first COVID-19 wave (‘Wave 1’, June 16, 2020 – August 31, 2020) and the second wave (‘Wave 2’, December 16, 2020 – March 02, 2021). Among the three methods used in the first step, EpiNow2 yielded the highest accuracy of $\:{R}_{t}$ prediction in the regions with entirely missing data. Median county-level percentage agreement (PA) was 90.9% (Interquartile Range, IQR: 89.9–92.0%) and 92.5% (IQR: 91.6–93.4%) for Wave 1 and 2, respectively. Median zip code-level PA was 95.2% (IQR: 94.4–95.7%) and 96.5% (IQR: 95.8–97.1%) for Wave 1 and 2, respectively. Using EpiEstim, EpiFilter, and an ensemble-based approach yielded median PA ranging from 81.9 to 90.0%, 87.2-92.1%, and 88.4-90.9%, respectively, across both waves and geographic granularities.

Conclusion

These findings demonstrate that the proposed methodology is a useful tool for small-area estimation of $\:{R}_{t}$, as our flexible framework yields high prediction accuracy for regions with coarse or missing data.

Introduction

Infectious disease outbreaks pose significant threats, impacting health, operations, and the economy. Effective outbreak management relies on the ability to estimate the intensity of local transmission and strategize appropriately, which is vital for controlling disease transmission, enhancing health outcomes, and reducing economic impacts. However, current approaches often target broad geographic areas, which may not yield actionable insights for local-level or institutional decision-making [1]. This lack of granularity can be attributed to variations in community characteristics such as age distribution, prevalence of comorbidities, and social behaviors [2, 3], all of which influence disease transmission and susceptibility to severe outcomes. Moreover, policies derived from coarse data may lead to suboptimal local disease control [4, 5]. Therefore, enhancing local-level transmission estimation and forecasting is essential for tailoring policies and optimizing interventions.

The effective reproductive number is a key epidemiological metric used to describe the transmissibility of infectious diseases for a specific time and setting; it is critical to understand the intensity of transmission. $\:{R}_{t}$ represents the average number of people infected by one infectious individual at a given time t. An $\:{R}_{t}$ value greater than $\:1$ signifies that the disease is spreading, while an $\:{R}_{t}$ less than $\:1$ indicates that the transmission is declining [6, 7]. $\:{R}_{t}$ changes with the progression of an epidemic and the implementation of control measures [8,9,10]. It also varies spatially due to the differences in population density, mobility, social interactions, and public health policies [11, 12]. Unfortunately, local-level estimation of $\:{R}_{t}$ is challenging due to data sparsity or unavailability of infectious disease testing or encounter records. For example, insufficient testing and changes in case definitions or reporting can substantially impact estimates of R_t [13, 14]. Due to substantial lags in producing estimates of disease epidemiological metrics from traditional surveillance systems [15], electronic health records (EHR) obtained directly from health systems or hospitals are useful for real-time estimation of such metrics, including $\:{R}_{t}$. However, healthcare systems do not serve all populations, which inhibits the ability to forecast disease spread in communities that are not represented in the healthcare system. Additionally, estimates of $\:{R}_{t}$ may be unstable for communities that are underrepresented in the healthcare system and hence offer few data points for accurate estimation.

There are several techniques for estimating the effective reproductive number; e.g., Wallinga and Tenius method [16], Cori et al. method [17], EpiFilter [6], and EpiNow2 [18]. EpiEstim is a tool integrated into widely used software platforms for estimation of the reproductive number from incidence time series, accounting for uncertainties in the serial interval distribution [17]. EpiFilter utilizes all available incidence information and does not rely on window size assumptions, making it more statistically robust, particularly during periods of low incidence [6]. The modeling approach of EpiNow2, an R software package [18], is based on the methods developed by Cori et al. and Thompson et al. [19], and it incorporates the delays from infection to reporting to estimate the effective reproductive number utilizing a range of open-source tools [20]. Although there are several methods for estimating the effective reproductive number, there is a notable gap in accurately estimating $\:{R}_{t}$ for small areas with coarse or missing data. To address this gap, we develop a flexible framework that effectively allows us to implement any estimation procedure for $\:{R}_{t}$ (e.g., EpiEstim, EpiFilter, EpiNow2, etc.) in regions with coarse or entirely missing data.

The motivation for this study, funded by the National Library of Medicine of the National Institutes of Health (NIH) and the Center for Forecasting and Outbreak analytics of the Centers for Disease Control and Prevention (CDC), is to inform allocation of mobile health clinics for infectious disease interventions to medically underserved communities at high risk of infectious disease outbreaks. We address the challenge of local level $\:{R}_{t}$ estimation in areas with coarse or entirely missing data by proposing a two-step estimation procedure which incorporates existing R_t estimation techniques along with an ensemble estimation using data from geographic regions with sufficient data (step 1) into a spatial modeling framework to predict $\:{R}_{t}$ in regions with sparse or missing data (step 2).

The remainder of the paper is organized as follows: the Materials and Methods section outlines the framework of the study, detailing the statistical techniques and data utilized. In the Results section, we present our findings and interpret the results. The Discussion section explores the implications of the results, acknowledges the limitations of our study, and suggests directions for future research. The paper concludes with the Conclusion section, summarizing the main points.

Materials and methods

To estimate the region-specific effective reproductive number $\:\left({R}_{t}\right)$ for an infectious disease in small areas with coarse or entirely missing data, we develop a two-step procedure. In the first step of our analyses, we estimate $\:{R}_{t,i}$, using the existing techniques EpiEstim, EpiFilter, and EpiNow2 for each geographic region $\:i$ [6, 17, 18]. We denote the resulting estimates by $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, and $\:{\widehat{R}}_{t,Now,i}$, respectively. We also perform an ensemble-based estimation ($\:{\widehat{R}}_{t,Ensemble,i})$ by considering equal weights for each estimation technique. In the second step, we employ a Bayesian Integrated Laplace Approximation (INLA) model incorporating spatial and sociodemographic information. For each geographic region $\:i$ with available outcome data (e.g., Covid-19 cases), we denote these estimates by $\:{\widehat{R}}_{t,Estim,i}^{S}$, $\:{\widehat{R}}_{t,Filter,i}^{S}$, $\:{\widehat{R}}_{t,Now,i}^{S}$, $\:{\widehat{R}}_{t,Ensemble,i}^{S}$, which are referred to as spatially and covariates smoothed estimates. When geographic region $\:i$ is not included in the spatial (and covariate-adjusted) INLA model fitting in step 2 (e.g., geographic region $\:i$ has entirely missing data), we denote the prediction in this region by $\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Filter,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Now,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Ensemble,i}^{{S}{\left(i\right)}}$, where $\:{S}{\left(i\right)}$ indicates information from region $\:i$ is not used in the model fitting in step 2. More details of initial estimation, smoothing and predictions of $\:{R}_{t}$ in regions with entirely missing data, along with their evaluations, are discussed in the respective sections.

Data description

The number of cumulative confirmed COVID-19 case data for US counties is publicly available on the New York Times (NYT) GitHub repository [21]. The daily reported cases are calculated from the county-level cumulative data. We obtained the ZIP code level daily COVID-19 case data from Prisma Health’s (SC’s largest not-for-profit health care provider) COVID-19 registry. We utilize data from SC counties and ZIP codes during the first COVID-19 wave (Wave 1, between June 16, 2020 – August 31, 2020) and the second wave (Wave 2, between December 16, 2020 – March 02, 2021). Since many ZIP codes in the Prisma Health system data lack sufficient case numbers, we restrict our analysis to those ZIP codes with at least 50 cases in each period to ensure the reliability of our results (Wave 1: N = 30 ZIP codes; Wave 2: N = 45 ZIP codes). Population demographics are characterized by age, sex, race/ethnicity, employment, and insurance coverage. These data are sourced from the United States Census Bureau website [22]. The social vulnerability index (SVI) data is obtained from the Agency for Toxic Substances and Disease Registry website [23].

Initial estimation of R _t

We utilize the techniques EpiEstim, EpiFilter, and EpiNow2 to obtain the initial estimates of the effective reproductive number. EpiEstim and EpiFilter assume the incidence $\:{I}_{t}$ at time $\:t$ is Poisson distributed with mean $\:{R}_{t}\sum\:_{s=1}^{t-1}{I}_{t-s}{w}_{s}$, where $\:{R}_{t}$ is the effective reproductive number at time $\:t$, $\:s$ is the time since infection of the case, $\:{w}_{s}$ is the probability that a primary case takes time between $\:s-1$ and $\:s$ days to generate the secondary infection, and $\:{w}_{s}$ is obtained from the serial interval distribution [6, 17]. EpiEstim uses an expectation-maximization algorithm to reconstruct daily incidence from aggregated data for estimating $\:{R}_{t}$ [24]. It computes estimates of the effective reproductive number, $\:{\widehat{R}}_{t,Estim,i}$, using a sliding time window; the default seven days window setting is used. EpiFilter employs a recursive Bayesian smoothing technique to estimate the effective reproductive number, $\:{\widehat{R}}_{t,Filter,i}$. It unifies and extends two popular techniques EpiEstim and the Wallinga-Teunis method [6]. EpiNow2 employs a Bayesian latent variable model that incorporates a probabilistic programming language, Stan [20, 25]. Detailed description of the EpiEstim, EpiFilter, and EpiNow2 methods are available in the literature [6, 17, 18, 20]. Each of the three methods estimates $\:{R}_{t}$ using different parameters or data assumptions. All three methods depend on the generation time distribution. Since infection times are often unobserved, the generation time is usually approximated using the serial interval distribution [17]. In all three methods, the parameters for both the serial interval distribution and incubation period are utilized.

From the literature, we obtain a mean of 5 days (95% CI: 4.94–5.06) and a standard deviation (SD) of 2.4 for the incubation period [26], a mean of 4.7 days and an SD of 2.9 for the serial interval distribution [27]. The ‘estimate_R’ function of the EpiEstim, R software package version 2.2.4 is used to obtain $\:{\widehat{R}}_{t,Estim,i}$. We calculate the estimate, $\:{\widehat{R}}_{t,Filter,i}$, by using the R software function ‘epiSmoother’, available at https://github.com/kpzoo/EpiFilter. The ‘epinow’ function from the R software package EpiNow2, version 1.4.0, is utilized to get $\:{\widehat{R}}_{t,Now,i}$. We develop an ensemble method to obtain $\:{\widehat{R}}_{t,Ensemble,i}$ by averaging $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, and $\:{\widehat{R}}_{t,Now,i}$. The ensemble method aims to reduce the bias and variance in epidemiological forecasting by averaging the various model outputs [28].

Proposed covariate-adjusted spatial smoothing via INLA

We propose a sociodemographic covariate adjusted, two-step Integrated Nested Laplace Approximation (INLA) model to refine the initial estimation of the effective reproductive number at the local-level and to provide estimates for the regions with coarse or missing data. The use of INLA for spatial data prediction is supported by extensive research, demonstrating its effectiveness in handling spatial dependency and missing data challenges [29, 30]. The INLA method is designed to derive the posterior marginal distribution for latent Gaussian models and serves as a computationally efficient alternative to traditional Markov Chain Monte Carlo (MCMC) methods for Bayesian inference [29, 31]. INLA utilizes a series of nested Laplace approximations to calculate the marginal posterior and hyperparameters of latent Gaussian models. The INLA framework is employed to integrate the spatial dependencies of the effective reproductive number estimates into our model [32]. To fit the INLA model, we use the R software package called “R-INLA” [29, 30, 33]. We incorporate the effect of the spatial location (counties or ZIP codes) in the INLA function through Besag’s model [34, 35]. The Besag model employs spatial random effects for each area to capture the correlation among neighboring areas, adjusting the estimated $\:{R}_{t}$ based on the estimates of the neighboring areas. Neighborhood structures are defined using the ‘poly2nb()’ function, identifying areas as neighbors if their boundary polygons share at least one vertex. We then convert the neighborhood structure into a binary adjacency matrix using the ‘nb2mat()’ function. Each element of the adjacency matrix represents whether two locations are neighbors (1 if they share a boundary, 0 if they do not).

The model equation in INLA is:

$$\:{\widehat{R}}_{t,m,i}={\beta\:}_{0}+\beta\:{X}_{i}+{u}_{i}+{\epsilon\:}_{i}$$

where $\:{\beta\:}_{0}$ is the intercept of the model, $\:\beta\:$ is the coefficients vector for sociodemographic community-level covariates $\:{X}_{i}$, $\:{\widehat{R}}_{t,m,i}$ corresponds to the estimates $\:{\widehat{R}}_{t,Estim,i}$,$\:\:{\widehat{R}}_{t,Filter,i}$, $\:{\widehat{R}}_{t,Now,i}$, or $\:{\widehat{R}}_{t,Ensemmble,i}$ at area $\:i$.The term $\:{u}_{i}$ specifies the spatial random effect for location $\:i$, which depends on the neighboring locations defined by the adjacency matrix, and the random Gaussian error term is $\:{\epsilon\:}_{i}$.

The sociodemographic covariates used for the analyses at the county and ZIP code levels are the percentage of the population for different age groups (“Age 0–19”, “Age 20–44”, “Age 45–64”, “Age 65 and over”), the percentage of the population for different races and ethnicity (“White”, “Hispanic”, “Black or African American” and “others”), the percentage of the population “Employed” and “Unemployed”, the percentage of “Health Insurance Coverage”, the percentage of the population of “Male” and “Female”, the Social Vulnerability Index (SVI), and household median income. Household median income is normalized from 0 to 1, and SVI values range from 0 to 1, where higher values indicate greater social vulnerability. We consider the covariates “Age 0–19”, “White”, “Unemployed”, “Uninsured”, and “Female” population groups as the reference categories.

Assessment of prediction of $\:{\widehat{R}}_{t,m,i}$ for regions with entirely missing data

We employ the two-step spatial (covariate-adjusted) INLA model using the initial estimate, $\:{\widehat{R}}_{t,m,i}$, for each method m (i.e., EpiEstim, EpiFilter, EpiNow2, and Ensemble) at location $\:i$, as the response variable. We validate the model using two approaches. In the first approach, we randomly select 90% of the regions to hold out for the training set and 10% for the test set. We then apply the proposed two-step spatial (covariate-adjusted) INLA model to the training set and use this model to predict the effective reproductive number for the geographic regions in the test set, which we denote by $\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}$. We then compare the agreement between $\:{\widehat{R}}_{t,m,i}$, and $\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}$ for the geographic regions that are left out from the model fitting in the second step. The prediction accuracy measurement metrics are calculated using the formulas in Eqs. 1–3. This process is repeated 200 times.

In the second validation approach, we hold out geographic regions with sufficient data for the test set – specifically, Charleston, Greenville, Horry, and Richland counties and ZIP codes 29605, 29642, 29680, and 29681. Since there is no gold standard for comparing $\:{R}_{t}$ estimates, we select the above regions to predict $\:{R}_{t}$ for test cases with entirely missing data, assuming the initial estimates from regions with sufficient data will serve as a baseline for comparison. The prediction in this approach is performed by leaving out one area at a time.

The prediction accuracy measurement metrics are calculated using the following formulas:

Mean Absolute Error,

$$\:MAE=\frac{1}{T}\sum\:_{t=1}^{T}|\:{\widehat{R}}_{t,m,i}-{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}|$$

(1)

Root Mean Squared Error

$$\:RMSE=\sqrt{\frac{1}{T}\sum\:_{t=1}^{T}{\left(\:\:{\widehat{R}}_{t,m,i}-{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}\right)}^{2}}$$

(2)

Percentage Agreement

$$\:PA=\frac{1}{T}\sum\:_{t=1}^{T}\frac{\text{m}\text{i}\text{n}(\:{\widehat{R}}_{t,m,i}\:,\:\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}})\:\:}{\text{m}\text{a}\text{x}(\:{\widehat{R}}_{t,m,i}\:,\:\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}})}$$

(3)

where, $\:{\widehat{R}}_{t,m,i}$ is the initial estimates, $\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}$ is the prediction of effective reproductive number for each method $\:m$, at time $\:t$ and location $\:i$, and $\:T$ is the total number of time points (days).

Simulation study

We conduct a simulation study to validate our two-step spatial (covariate-adjusted) INLA framework under controlled conditions. To generate realistic effective reproductive number $\:({R}_{t}$) values, we capture the time-varying nature of disease transmission using a sinusoidal function. The parameters for the sinusoidal function are drawn from a multivariate normal (MVN) distribution, with a covariance structure derived from the adjacency matrix of the counties in SC. This structure ensures that spatial correlations are incorporated, allowing neighboring counties to have similar transmission patterns. The parameters are then scaled to fit epidemiologically meaningful ranges, representing the baseline reproductive number, amplitude, frequency of disease waves, and linear adjustments. Using these parameters, $\:{R}_{t}$​ values are generated and used to simulate daily COVID-19 case counts at the county level. The number of new infections at time $\:t$ is based on total infectiousness, which depends on past cases weighted by a discretized gamma-distributed serial interval.

To evaluate the performance of different estimation methods, we generate 50 independent sets of daily COVID-19 case data. The analysis follows the same workflow as for real data, including initial $\:{R}_{t}$​ estimation, two-step spatial (covariate-adjusted) INLA smoothing, and prediction of $\:{R}_{t}$​ for counties with entirely missing data. For initial estimation, we apply EpiEstim, EpiFilter, and an ensemble-based approach, excluding EpiNow2 due to its high computational cost. To assess the predictive performance of our framework, we systematically remove selected counties from the INLA model fitting process, making their data unavailable. The model is then trained using data from the remaining counties, and $\:{R}_{t}$​ is predicted for the excluded counties. Unlike real data, where the true $\:{R}_{t}$ is unknown, the simulation study generates true $\:{R}_{t}$ values, allowing to benchmark and compare the performance of different estimation methods. Further details on the simulation methodology are provided in the Supporting Information 1 (pages 14–17).

Results

Initial estimation and spatial (covariate-adjusted) INLA smoothing of $\:{R}_{t}$

For the $\:{R}_{t}$ estimation in this section, we utilized the county and ZIP code level COVID-19 daily case data for the first wave (June 16, 2020, through August 31, 2020) and second wave (December 16, 2020, through March 02, 2021) of the disease in SC. We measure the spatial (covariate-adjusted) INLA model at each time point using DIC (Deviance Information Criterion), Watanabe-Akaike Information Criterion (WAIC), and Log Marginal Likelihood. The minimum and maximum values of these diagnostics across time points are reported in the supporting information (Table S1).

Figure 1 displays the estimates of effective reproductive number for Charleston, Greenville, Horry, and Richland counties during the second COVID-19 wave in SC. In general, the proposed two-step spatial INLA estimates, $\:{\widehat{R}}_{t,Estim,i}^{S}$, $\:{\widehat{R}}_{t,Filter,i}^{S}$, $\:{\widehat{R}}_{t,Now,i}^{S}$, $\:{\widehat{R}}_{t,Ensemble,i}^{S}$, demonstrated reduced peaks compared to the initial estimates, $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, $\:{\widehat{R}}_{t,Now,i}$, $\:{\widehat{R}}_{t,Ensemble,i}$.

Fig. 1

Comparison of two-step spatial (covariate-adjusted) INLA smoothing with initial $\:{R}_{t}$ estimates at the county level in Wave 2. $\:{\widehat{R}}_{t,m,i}^{S}$ (blue line) represents the two-step spatial (covariate-adjusted) INLA smoothed estimates, $\:{\widehat{R}}_{t,Estim,i}^{S}$, $\:{\widehat{R}}_{t,Filter,i}^{S}$, $\:{\widehat{R}}_{t,Now,i}^{S}$, and $\:{\widehat{R}}_{t,Ensemble,i}^{S}$, while $\:\:{\widehat{R}}_{t,m,i}$ (red line) represents the initial estimates, $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, $\:{\widehat{R}}_{t,Now,i}$, and $\:{\widehat{R}}_{t,Ensemble,i}$, for select Charleston, Greenville, Horry, and Richland counties during COVID-19 Wave 2 (between December 16, 2020 – March 02, 2021) in SC. The state level initial estimate of $\:{R}_{t}$ is presented with dashed gray line (SC State). The plots in the rightmost panel present the average daily cases for the respective counties over the same period

S1-S3 Figs in the supplementary materials (Supporting Information 1) provide a detailed comparison of the effective reproductive number estimates, $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, $\:{\widehat{R}}_{t,Now,i}$, $\:{\widehat{R}}_{t,Ensemble,i}$ and $\:{\widehat{R}}_{t,Estim,i}^{S}$, $\:{\widehat{R}}_{t,Filter,i}^{S}$, $\:{\widehat{R}}_{t,Now,i}^{S}$, $\:{\widehat{R}}_{t,Ensemble,i}^{S}$, across various counties and ZIP codes. Results are presented for Charleston, Greenville, Horry, and Richland counties and ZIP codes 29605, 29642, 29680, and 29681, as these results serve as benchmarks when evaluating the prediction accuracy. The bar graphs displaying daily COVID-19 cases reveal spikes and declines that correspond to the fluctuations in the $\:{R}_{t}$ estimates. For instance, peaks in case counts align with the increases in both sets of estimates, $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, $\:{\widehat{R}}_{t,Now,i}$, $\:{\widehat{R}}_{t,Ensemble,i}$ and $\:{\widehat{R}}_{t,Estim,i}^{S}$, $\:{\widehat{R}}_{t,Filter,i}^{S}$, $\:{\widehat{R}}_{t,Now,i}^{S}$, $\:{\widehat{R}}_{t,Ensemble,i}^{S}$. The figures show that when daily cases increase over a period, corresponding $\:{R}_{t}$ estimates exceed 1.0, indicating infection growth. For example, in Fig. 1, daily cases for Charleston, Greenville, and Horry counties increase between time points 10 and 25, with $\:{R}_{t}$ estimates staying above 1.0. After time points 25, as cases decline almost everywhere, $\:{R}_{t}$ estimates fall to 1.0 or lower, particularly with our two-step spatial INLA method (blue lines).

Prediction of $\:{R}_{t}$ in regions with entirely missing data

In this section, we evaluated the accuracy of our two-step spatial INLA model in estimating $\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Filter,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Now,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Ensemble,i}^{{S}{\left(i\right)}}$ for regions with entirely missing data. We utilized daily COVID-19 case data at both county and ZIP code levels for the first wave (June 16, 2020, through August 31, 2020) and second wave (December 16, 2020, through March 02, 2021) of the pandemic in SC. Table 1 presents the prediction accuracy achieved by the two-step spatial INLA model through comparison with the one-step estimates $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, $\:{\widehat{R}}_{t,Now,i}$, and $\:{\widehat{R}}_{t,Ensemble,i}$. That is, we used the two-step spatial INLA model based on a subset of geographic regions (training set) to predict estimates of $\:{R}_{t}$ for the remaining geographic region (i.e., validation set). Since there is no gold standard for comparison, we compared the predicted effective reproductive number, $\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}$, to the initial estimate $\:{\widehat{R}}_{t,m,i}$.      

At both county and ZIP code levels during the second wave of COVID-19, predicted $\:{R}_{t}$ values, $\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}$, based on the two-step spatial INLA estimators achieved high level of accuracy (Table 1). For counties, the median percentage agreement was 90.0% (IQR: 89.0–90.9%) for $\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$, 87.2% (IQR: 86.3–88.1%) for $\:{\widehat{R}}_{t,Filter,i}^{{S}{\left(i\right)}}$, 92.5% (IQR: 91.6–93.4%) for $\:{\widehat{R}}_{t,Now,i}^{{S}{\left(i\right)}}$, and 90.5% (IQR: 89.6–91.2%) for $\:{\widehat{R}}_{t,Ensemble,i}^{{S}{\left(i\right)}}$. The prediction, $\:{\widehat{R}}_{t,\:Now,i}^{{S}{\left(i\right)}}$, achieved the highest accuracy, with a median PA of 92.5% at the county level. Similarly, $\:{\widehat{R}}_{t,Now,i}^{{S}{\left(i\right)}}$ exhibited the highest median PA of 96.5% for SC ZIP codes. In the first wave, the accuracy was also high, particularly for $\:{\widehat{R}}_{t,Now,i}^{{S}{\left(i\right)}}$, which obtained a median PA of 90.9% at the county level and 95.2% at the ZIP code level (see Table S2). However, $\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$ showed the lowest accuracy, with a median PA of 81.9% at the ZIP code level, and a median PA of 85.7% at the county level. Overall, the two-step spatial (covariate-adjusted) INLA procedure showed high performance accuracy in computing $\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}$ for areas with missing data, as further evaluated through root mean squared error (RMSE) and mean absolute percentage error (MAPE).

Table 1 Comparison of two-step Spatial (covariate-adjusted) INLA prediction, $\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}$, with initial estimates, $\:{\widehat{R}}_{t,m,i}$, in Wave 2

In the second validation approach, we selected Charleston, Greenville, Horry, and Richland counties (county level analysis), and ZIP codes 29605, 29642, 29680, and 29681 (zip code level analysis) as the test set. Figure 2 provides a detailed comparison of the predictions $\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Filter,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Now,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Ensemble,i}^{{S}{\left(i\right)}}$ to initial estimates $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, $\:{\widehat{R}}_{t,Now,i}$, $\:{\widehat{R}}_{t,Ensemble,i}$ across counties–Charleston, Greenville, Horry, and Richland and ZIP codes – 29605, 29642, 29680, and 29681 (also see S4-S8 Figs). The results comparing prediction accuracy for this validation setting are presented in Table S3 for the first wave and Table S4 for the second wave. During the second wave, $\:{\widehat{R}}_{t,Now,i}^{{S}{\left(i\right)}}$ demonstrated the highest accuracy, with a percentage agreement (PA) range of 91.1 − 94.0% across the four counties and 96.1 – 98.7% across the four ZIP codes. $\:{\widehat{R}}_{t,Filter,i}^{{S}{\left(i\right)}}$ achieved the lowest PA range for counties (86.7 – 90.8%) and $\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$ achieved the lowest PA range for ZIP codes (86.9–91.5%). Similar results were obtained during the first wave.

Fig. 2

Comparison of two-step spatial (covariate-adjusted) INLA prediction with the initial $\:{R}_{t}$ estimates in Wave 2. $\:{\widehat{R}}_{t,m,i}^{{S}{\left(i\right)}}$ (blue line) represents the two-step spatial (covariate-adjusted) INLA prediction of effective reproductive number, $\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Filter,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Now,i}^{{S}{\left(i\right)}}$, $\:{\widehat{R}}_{t,Ensemble,i}^{{S}{\left(i\right)}}$, while $\:{\widehat{R}}_{t,m,i}$ (red line) represents the initial estimates, $\:{\widehat{R}}_{t,Estim,i}$, $\:{\widehat{R}}_{t,Filter,i}$, $\:{\widehat{R}}_{t,Now,i}$, $\:{\widehat{R}}_{t,Ensemble,i}$, for select Charleston, Greenville, Horry, and Richland counties during COVID-19 Wave 2 (between December 16, 2020 – March 02, 2021) in SC. Here, $\:{S}^{\left(i\right)}\:in$$\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$ indicates that geographic region $\:i$ was not included in the INLA model fitting. For example, if $\:i$=Greenville (GVL), $\:{\widehat{R}}_{t,Estim,GVL}^{{S}{\left(GVL\right)}}$ means Greenville County was not used in the spatial (covariate-adjusted) INLA model fitting in step 2. The plots in the rightmost panel present the average daily cases for the respective counties over the same period

Simulation results

We evaluated the performance of our two-step spatial (covariate-adjusted) INLA smoothing and prediction framework for estimating $\:{R}_{t}$ at the county level. The analysis was conducted using simulated COVID-19 case counts data for counties in SC and focused on assessing the accuracy of initial estimates, improvements achieved through spatial INLA smoothing, and predictive ability of the framework for counties with completely missing data.

Figure S9 presents a comparison between the initial $\:{R}_{t}$​ estimates obtained using EpiEstim, EpiFilter, and the ensemble approach, the smoothed $\:{R}_{t}$​ estimates generated through the two-step spatial (covariate-adjusted) INLA smoothing, and the true generated $\:{R}_{t}$​ values used in the simulation. The results demonstrate that the initial estimates provide a reasonable estimation of $\:{R}_{t}$​ across counties. However, after applying spatial INLA smoothing, the estimates align more closely with the true generated $\:{R}_{t}$ values. The percentage agreement (PA) values in Table S6 also indicate that the smoothed $\:{R}_{t}$ estimates align well with the true values. In this scenario, all available data were used, meaning that the initial estimates were available for all counties before applying the spatial INLA model.

Figure S10 evaluates the predictive capability of our approach for estimating $\:{R}_{t}$​ in counties with completely missing data. To simulate this scenario, we excluded each county one by one from the model fitting process. For example, Charleston County was left out to replicate a situation where no data was available for that county. The model then predicted $\:{R}_{t}$​ for Charleston by incorporating spatial information from neighboring counties along with sociodemographic covariate data specific to this county. This process was repeated for Greenville, Horry, and Richland counties to evaluate the model’s robustness. The results show that the $\:{R}_{t}$​ values predicted using our two-step spatial (covariate-adjusted) INLA model closely align with the true generated $\:{R}_{t}$​ values. The percentage agreement values in Table S6 further validate the robustness of the proposed framework, showing that it maintains high prediction accuracy even for the counties with missing data.

For the initial estimates, the median PA was 92.55% (IQR: 92.63–93.74%) for $\:{\widehat{R}}_{t,Estim,i}$, 91.52% (IQR: 91.12–91.83%) for $\:{\widehat{R}}_{t,Filter,i}$, and 92.09% (IQR: 91.72–92.37%) for $\:{\widehat{R}}_{t,Ensemble,i}$. Our framework showed higher performance, the PA was 93.10% (IQR: 92.43–93.75%) for $\:{\widehat{R}}_{t,Filter,i}^{S}$ and 93.10% (IQR: 92.47–93.91%) for $\:{\widehat{R}}_{t,Ensemble,i}^{S}$. For region with entirely missing data that makes impossible to estimate $\:{R}_{t}$ with the existing methods, the PA was 90.76% (IQR: 89.67–91.38%) for $\:{\widehat{R}}_{t,Estim,i}^{{S}{\left(i\right)}}$, 91.30% (IQR: 90.79–91.57%) for $\:{\widehat{R}}_{t,Filter,i}^{{S}{\left(i\right)}}$. These findings highlight the effectiveness of our two-step spatial (covariate-adjusted) INLA framework in estimating $\:{R}_{t}$ at the local level with coarse or missing data.

Discussion

Estimation of time-varying effective reproductive number $\:{R}_{t}$ is necessary for understanding the transmission dynamics of infectious disease and implementing the necessary public health interventions. We developed a two-step procedure for estimating the effective reproductive number for infectious diseases in granular geographic regions. Our approach incorporates existing R_t estimation procedures (EpiEstim, EpiFilter, and EpiNow2), along with an ensemble-based estimate using data from geographic regions with sufficient data (step 1), into a spatial modeling framework to predict $\:{R}_{t}$ in regions with sparse or missing data (step 2). The choice of $\:{R}_{t}$ estimation methods is not exhaustive. Our flexible framework allows us to implement any existing estimation procedure for $\:{R}_{t}$ in regions with coarse or entirely missing data. Our proposed methodology improves the estimation of the effective reproductive number $\:{R}_{t}$ in small areas, especially in areas where the disease data are coarse or missing. By integrating the Bayesian spatial INLA model, this method efficiently utilizes the available data and borrows information from neighboring areas and similar regions through sociodemographic covariates, enhancing the robustness of $\:{R}_{t}$ estimates. Decisions regarding the implementation, removal, or adjustment of infection control measures are guided by a combination of epidemiological, social, and economic factors [36, 37]. Effective reproductive number estimates assists in this process by providing critical information about the potential future course of an outbreak. Through small area estimation of the effective reproductive number, we can better identify and prioritize high risk areas, which allows for optimal resource allocation [2]. For example, this current NIH and CDC-funded study is intended to better identify and prioritize medically underserved communities at high-risk of infectious disease outbreaks for priority delivery of mobile health clinics for infectious disease screening, treatment, and vaccination.

Our study has several limitations. One primary limitation is that our approach focuseson the retroactive estimation of $\:{R}_{t}$. This is particularly useful for understanding past outbreaks and forforecasting outbreaks when the effective reproductive number in each region is similar between outbreaks (e.g., first two waves of Covid-19 in SC; Table S5). While this can be a reasonable assumption for endemic viruses, this is far from guaranteed for emerging viruses where disease transmission dynamics and/or human behavior changes can substantially change with the introduction of a new variant. Therefore, a natural extension is to extend our approach for forecasting future values of $\:{R}_{t}$, particularly in granular geographic regions with sparse data.

We evaluated model performance using three commonly used techniques (EpiEstim, EpiFilter, EpiNow2) along with the ensemble-based estimation. Similar to most $\:{R}_{t}$ estimation methods, EpiEstim and EpiFilter rely on the reported incidence data [6, 17]. Inaccuracies or delays in reporting can lead to biased estimates. EpiNow2 also makes certain assumptions about the underlying transmission dynamics, such as the generation time distribution and reporting delays. These assumptions can introduce uncertainty into the estimates, particularly during the periods of low incidence [18]. Additionally, our study did not include the other $\:{R}_{t}$ estimation procedures, such as Epidemia [38], EpiInvert [39], EpiRegress [40], estimateR [41], GrowthPredict [42], and Extended Kalman Filter [43].

Due to the limitations of using existing $\:{R}_{t}$ estimators as a benchmark to evaluate the performance of our proposed estimator, there is no gold standard for comparison. To address this, we performed a simulation study to compare the $\:{R}_{t}$ estimates with ground truth $\:{R}_{t}$ values. While the simulation study demonstrated our proposed estimator was accurate and closely resembled ground truth estimates, it was conducted in a single simulation setting, with fixed parameter ranges and temporal dynamics, which may not fully capture the diversity of real-world scenarios.

Future research may extend our framework by incorporating the forecasting of $\:{R}_{t}$ values and subsequent case projections, providing critical insights into future disease dynamics. This would help policymakers with decision-making, support outbreak detection, and strengthen preparedness for future outbreaks. In our current approach, we utilized three existing $\:{R}_{t}$estimation methods and combined them using an ensemble approach with equal weights. Future effort also could explore incorporating additional estimation methods and developing an ensemble estimation technique with optimal weights to improve the accuracy and robustness across various scenarios. Furthermore, simulations under different settings could be conducted to capture a wider range of epidemic scenarios.

Conclusion

This study presented a two-step spatial INLA approach to estimate the effective reproductive number in small areas with coarse or missing data. Our approach incorporated existing estimators into a Bayesian spatial model that adjusted for spatial variation by borrowing information from neighboring regions and incorporating sociodemographic covariates. This approach provided predictions of the effective reproductive number for regions with missing data. Importantly, our proposed framework is flexible, allowing for implementation of any existing estimation procedure for $\:{R}_{t}$ in geographically granular regions. While this study focused on Covid-19, our methodology is translatable for estimation of the effective reproductive number for other respiratory infectious diseases.