We are interested in studying and forecasting transmission dynamics of severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) -a novel human coronavirus that is associated with high morbidity and mortality. Community spread of the virus was first reported in December 2019 in Wuhan, China, and by September 2020 infected more than 25 million individuals globally, resulting in more than 800,000 deaths (WHO, 2020a) . It is closely genetically related to SARS-CoV, which in 2003 caused the first known severe outbreak associated with a human coronavirus (Chen et al., 2020) . Infection with SARS-CoV-2 causes coronavirus disease , which ranges from asymptomatic or mild clinical illness in the majority of people, to Acute Respiratory Distress Syndrome (ARDS) and extrapulmonary complications in 5-30% of patients (Cummings and et. al, 2020; Wu and McGoogan, 2020; Song and et. al, 2020) . SARS-CoV-2, like other human coronaviruses, is thought to be transmitted primarily by respiratory droplets. Transmission may also occur through fomites and aerosols, with the latter being a major concern in crowded indoor settings with activities that generate aerosols (WHO, 2020b) . To suppress viral transmission, many geographic/administrative regions have implemented mitigation measures (business closures, physical distancing, requiring facial masks, etc.). To estimate effect of these policies and their relaxations, we develop a method that can estimate parameters of the SARS-CoV-2 transmission dynamics from relatively short time series of reported cases -positive diagnostic tests detecting presence of SARS-CoV-2 -and deaths due to Modeling and forecasting dynamics of SARS-CoV-2 transmission allows government officials to make informed decisions when managing the spread of the virus. In the early stages of the pandemic, such modeling played an important role in alerting the public about potential dangers of unmitigated virus spread (Prem and et. al, 2020; Ferguson and et. al, 2020; Davies et al., 2020) . Differences in mitigation strategies, surveillance efforts, and population characteristics across countries and even across different regions within one country prompted development of regional modeling of SARS-CoV-2 transmission. For example, Anderson and et. al (2020) , Miller and et. al (2020) , and Morozova et al. (2020) provide great examples of rigorous data-driven regional modeling. However, neither national nor subnational/regional modelers fully integrate all surveillance data available to them. In particular, incorporating case incidence data into inference proved particularly problematic, because the number of cases depends not only on the unobserved number of infections, but also on the number of diagnostic tests performed, which varies significantly temporally and spatially. However, even with delayed reporting, positive diagnostic tests (cases) are the earliest indicators of changing disease dynamics, so taking advantage of this source of information is important for producing timely forecasts and for policy decision making.
In this paper, we show how to fit a mechanistic model of SARS-CoV-2 spread to incidence and mortality time series, while accounting for time-varying number of diagnostic tests performed. The mechanistic model is a fairly standard ordinary differential equation (ODE) model that describes changes in the proportions of the population residing in model compartments (e.g., susceptible and infectious compartments). Death counts are modeled with a negative binomial distribution that allows for overdispersion often observed in surveillance data -a standard practice in infectious disease epidemiology. Our main innovation is the model for cases, where we use a flexible beta-binomial distribution, whose mean is a product of the total number of tests performed and a non-linear function of unobserved infections modeled by the ODE model. The "beta" part of the beta-binomial distribution ensures that our estimates are not unduly influenced by large fluctuations of COVID-19 diagnostic test positivity fractions.
We first validate our model using simulated data, demonstrating importance of adjusting for the number of tests performed for transmission parameter inference. Next, we fit our compartmental model to surveillance data collected in Orange County, California. Orange County is the sixth most populous county in the United States of America (U.S.A.), with an estimated 3.18 million inhabitants as of 2019 (United States Census Bureau, 2020) . The third confirmed COVID-19 case in U.S.A. was identified in Orange County on January 26, 2020 (NBC Bay Area, 2020), though community transmission following this case in a traveler from Wuhan, China was limited. In mid-March cases in OC began to surge, triggering major public health interventions, including a mandatory shelter in place order (Orange County Health Care Agency, 2020). A state-wide shelter in place order was put in place the following day (Gavin Newson, Governor of California, 2020a) . Cases then flattened and remained at relatively low levels until the public health and executive orders were relaxed beginning May 4 (Gavin Newson, Governor of California, 2020b) . In June cases locally, statewide, and nation-wide surged. In late June the Governor ordered indoor dining and bars to close again in an attempt to curb the growth of cases (NBC Los Angeles, 2020) .
To account for changes in SARS-CoV-2 transmission dynamics due to implementation and subsequent relaxation of mitigation measures, we analyze three 36 day time periods separately: March 31 -May 05, May 06 -June 10, and June 11 -July 16. The first period overlaps with the state-wide stay-at-home order. Stay-at-home order termination and fast re-opening of most businesses happened in the second and third time periods. Our analysis shows that the basic reproduction number has been steadily increasing after the stay-athome order ended. On a positive note, changes in parameters controlling our surveillance model for observed numbers of positive COVID-19 diagnostic tests seem to reflect changing testing guidelines that moved from testing mostly symptomatic individuals to casting a wider testing "net" over the Orange County population. In addition to estimating transmission and observation model parameters, we report the inferred latent incidence and prevalence and retrospectively evaluate ability of our model to produce short-term forecasts.
We start with time series of daily tests (positive and negative), case counts (positive tests), and deaths observed over some time period of interest. Daily resolution is too noisy due to reporting delays, weekend effect, etc., so we aggregate/bin the three types of counts in 3 day intervals. We use 3 day aggregation, because it allows us to analyze shorter time periods than more standard weekly aggregation, while still smoothing over decline of surveillance reporting during weekends and other data artifacts. Figure 1 shows such a collection of aggregated time series for Orange County, CA, corresponding to the observation period spanning days between March 31, 2020 and July 16, 2020. The data was compiled from anonymized individual test results provided by the Orange County Health Care Agency (OCHCA). We define cases as either confirmed or presumed COVID-19 diagnoses that have been officially reported to the state public health authorities. We used specimen collection dates and dates of deaths to tabulate test, case, and death counts. If an individual was tested more than once, we used only the first positive test of the individual for construction of test and case time series. We denote the vector of binned tests by T = (T 1 , . . . , T L ), the vector of case counts by Y = (Y 1 , . . . , Y L ), and the vector of deaths by M = (M 1 , . . . , M L ), where counts T l , Y l , and M l are aggregated over a 3 day interval indexed by l and L is the total number of these intervals. We do not model changes in the numbers of tests performed. Rather, we condition on test counts in the specification of the sampling model for the vector of case counts, which describes the probability of the observed case count given the observed number of tests and unobserved/latent incidence of cases over each time interval. To formulate this data generating model and its companion for the vector of death counts M, we first need a model for latent trajectories of incidence and prevalence of SARS-CoV-2 infections.
To model latent incidence and prevalence trajectories, we divide all individuals in a population of interest (e.g., population of Orange County, CA) into 6 compartments: S = susceptible individuals, E = infected, but not yet infectious individuals, I e = infectious individuals at early stages of infection, I p = infectious individuals who progressed in their infection stage, R = recovered individuals, D = individuals who died due to COVID-19. Possible progressions of an individual through the above compartments are depicted in Figure 2 . We model the time-evolution of the proportions of individuals occupying the above compartments with a set of deterministic ordinary differential equations (ODEs). For simplicity, we assume a closed and homogeneously mixing population, although it is possible to relax these assumptions, and we also assume that recovery confers immunity to subsequent infection over the duration of the modeling period. Let X(t) = S (t), E(t), I e (t), I p (t), R(t), D(t) T denote the population proportions in each compartment at time t, and let X(t 0 ) = x 0 denote the population proportions at time t 0 , the start of the modeling period. By convention, we model the population at risk, i.e., those individuals who may still move throughout the model compartments. Hence, we take R(t 0 ) = R 0 = D(t 0 ) = R 0 = 0 and normalize X(t) so that X(t) T 1 = 1. Since we want to fit this model to incidence data, it is convenient to also keep track of cumulative fractions of the population that experience transitions between compartments from t 0 to t: N(t) = N S E (t), N EI e (t), N I e I p (t), N I e R (t), N I p R (t), N I p D (t) T . To describe mathematically how vectors X(t) and N(t) change through time, we first define rates of transitions between compartments, with possible transitions corresponding to the arrows in Figure 2 :
where β is the infection rate, δ is the fraction by which progressed infectious individuals become less infectious than individuals at the early stage of infection, 1/γ is the mean latent period duration, 1/ν e is the mean early infectious period duration, 1/ν p is the mean progressed infectious period duration, and η is the probability of transitioning from I p to death, rather than to recovery. We do not estimate δ and fix it to 0.8 in all our analyses.
Equipped with the population-level transition rates, we are ready to define ODEs for our model:
subject to initial conditions X(t 0 ) = x 0 and N(t 0 ) = 0, where x 0 = (S 0 , E 0 , I e,0 , I p,0 , R 0 , D 0 ) are initial compartment proportions. We set R 0 = 0 and D 0 = 0, because these proportions do not play a role in future dynamics of the epidemic, leaving S 0 , E 0 , I e,0 , and I p,0 as free model parameters. Hence, our model describes the forward time evolution of the population still at risk. Setting R 0 = 0 and D 0 = 0 works well at the start of the epidemic, but complicates parameter interpretation as the sizes of R and D compartments become non-negligible. We will make the initial concentrations R 0 and D 0 free parameters in the future revisions of the model. The above equations are redundant, and typically only the prevalence ODEs in the left column are used in mathematical modeling. However, the cumulative incidence/transition representation of the model, shown by the ODEs in the right column, is useful for statistical modeling of infectious disease dynamics (Bretó and Ionides, 2011; Ho et al., 2018) . In practice, we solve the subset of the above ODEs that are needed to track X(t) and the parts of N(t) that "connect" our transmission model to data. We proceed to make this connection in the next subsection. 
Recall that we would like to fit our transmission model to two time series: numbers of new cases and deaths reported during some pre-specified time periods (e.g., 3 day periods or weeks). First, we assume that conditional on the fractions of latent transitions between compartments N(t), case and death counts are independent of each other, because they are just noisy realizations of information encoded by N(t). Similarly, conditional on N(t), case counts are independent across time intervals and death counts are independent across time intervals. This leaves us with formulating models for cases and deaths in each individual observation interval.
Consider the number of deaths M observed in time interval (t −1 , t ], where = 1, . . . , L. Since our ODEs track the latent cumulative fraction of deaths N I p D (t ), we can compute ∆N I p D (t ) = N I p D (t ) − N I p D (t −1 ) -the latent fraction of the population that died in the interval (t −1 , t ]. We model the observed death count M as a realization from the following negative binomial distribution:
where N is the population size, µ D and σ 2 are the mean and variance of the negative binomial distribution, ρ ∈ [0, 1] is the mean death detection probability, and φ > 0 is an overdispersion parameter. Informally, our mortality model says that, on average, the observed number of deaths, M , is a fraction of the true death count estimated by the model, N × ∆N I p D (t ), with some noise due to underreporting and sampling variability. Next, we develop a model for the number of positive tests (cases), Y , observed in time interval (t −1 , t ]. We start with a simple binomial model with per-test positivity probability ψ :
where T is the number of COVID-19 diagnostic tests administered during time interval (t −1 , t ]. We could proceed to modeling ψ as a deterministic function of latent incidence, but observed variation in estimated positivity probabilities cannot be explained by changes in incidence alone. Variable testing guidelines, test shortages, and other factors must be at play. We use another layer of randomness to account for these unobserved factors affecting positivity probabilities and assume that the positivity probability in interval (t −1 , t ] follows the following Beta distribution:
where κ is an overdispersion parameter and µ C is the mean test positivity probability. We assume that mean test positivity log-odds is a linear function of the unobserved log-odds-transformed fraction of the population that transitioned from early to progressed infectiousness, ∆N I e I p (t ) = N I e I p (t ) − N I e I p (t −1 ) in interval (t −1 , t ]:
where α 0 > 0 and 0 < α 1 < 1. This somewhat arbitrary functional form nonetheless ensures that on average probability of detecting a COVID-19 case grows with the population incidence, as shown in Figure 3 . Parameters α 0 and α 1 can be thought of as reflecting testing guidelines. A model with α 0 = 0 and α 1 = 1.0 (i.e., µ C = ∆N I e I p (t )) says that in interval (t −1 , t ] testing is done approximately by sampling individuals uniformly at random, so that the positivity probability over time interval is equal to the fraction of the population that transitions from early to progressed infection. As we increase α 0 above 0 and decrease α 1 below 1, the model mimics preferential testing of individuals who are more likely to have severe infection (e.g., testing only individuals with certain symptoms). We can streamline our surveillance model for case counts by integrating over positivity probabilities and arriving at the following Beta-binomial distribution:
Properties of the Beta-binomial distribution imply that E(Y ) = T ×µ C . This means that our model predicts that on average cases grow linearly with the number of tests administered. Keeping in mind our assumed relationship between µ C and ∆N I e I p (t ), the average number of cases also grows with accumulation of new infections. Furthermore, the variance of the fraction of tests that are positive under the beta-binomial distribution is Figure 3 : Relationship between COVID-19 diagnostic test positivity probability and latent incidence in the population defined by equation (5).
where the variance under an analogous pure binomial model would be µ C (1−µ C )/T . Hence, the overdispersion parameter, κ, can be interpreted in terms of the excess variance of the beta-binomial model relative to a pure binomial distribution. In summary, our Beta-binomial distribution for observed case counts ensures that we do not confuse increase in testing for increase in COVID-19 incidence and implicitly allows for heterogeneity in the mean test positivity probability.
We are now ready to describe our inferential Bayesian procedure. First, we slightly re-parameterize our model by replacing β with a basic reproductive number R 0 = β 1/ν e + δ/ν p , obtained as the dominant eigenvalue of the next-generation matrix for our ODEs linearized about the disease-free equilibrium state (Diekmann et al., 2010) . In addition, we parameterize initial compartment fractions as S 0 , E 0 = (1 − I 0 )(1 − S 0 ), I e,0 = I e,0 I 0 (1−S 0 ), and I p,0 = (1− I e,0 ) I 0 (1−S 0 ), where I 0 = I e,0 + I p,0 / (1 − S 0 ) is a fraction of infected individuals who are infectious and I e,0 = I e,0 / I e,0 + I p,0 is a fraction of early stage infectious individuals out of all infectious individuals in both early and progressed stages. Next, we collect all our model parameters into a vector θ = (S 0 ,Ĩ 0 ,Ĩ e,0 , R 0 , γ, ν e , ν p , η, ρ, φ, α 0 , α 1 , κ). Our probabilistic construction described above implies that the likelihood function -probability of observing incidence and mortality data -can be written in the following way:
where Pr(M l | θ) and Pr(Y l | θ) are Negative-binomial and Beta-binomial probability mass functions defined by equations (3) and (6) respectively. We encode available information about our model parameters in a prior distribution with density π(θ). We assume that all parameters are a priori independent and list our prior assumptions in Table 1 . Since our model is highly parametric, we rely on informative prior distributions that we parameterize using existing scientific studies. We base all our inferences and predictions on the posterior distribution of all model parameters:
We approximate this posterior using Hamiltonian Monte Carlo implemented in the Stan probabilistic programming language for fitting Bayesian models, version 2.21.2 (Carpenter et al., 2017 ; Stan Development Team, These priors are used only for the initial model fit at the beginning of the epidemic. In models fit to later stages of the epidemic, method of moments estimators are used to create a prior based on the posterior of the relevant compartments on the relevant day in the previous model. 2020). We used four Markov chains in parallel to draw a total of 8,000 posterior samples with the first 4,000 discarded as warmup. Convergence and mixing was assessed using potential scale reduction factors, effective posterior sample sizes, and traceplots of model parameters. Model code and data are available at the following GitHub repository: https://github.com/vnminin/uci_covid_modeling/.
We used two simulation scenarios to validate our modeling approach and to examine whether accounting for time-varying testing makes a difference (see Supplementary Section A). We first simulated 1,000 datasets from our data generating model and obtained posterior distributions of our model parameters based on these data sets. In the second simulation study, we examined the effect of ignoring variability in the number of tests conducted over time. The results of both simulation studies are summarized in Supplementary Figures A2-A3  and Supplementary Table A4 . In the first simulation study, low mean absolute errors and coverage of credible intervals close to the nominal 95% demonstrate reasonable performance of our inferential procedure. The second simulation study shows that ignoring fluctuations in diagnostic tests indeed biases parameter estimation. In particular, 95% credible intervals for basic reproductive number R 0 contains the true value of this parameter only in 8% of simulations ( Table A4 ). This result highlights the importance of appropriately adjusting for the number of tests performed when analyzing case counts. We show summaries of the priors and posteriors for both models fit to a single simulated data set, depicted in Figure A1 , in Figures A4 and A5. These summaries suggest some of our parameters may be weakly identifiable, if at all. For instance, the prior and posterior summaries look nearly identical for mean durations of stay in E and I e compartments. Therefore, close attention needs to be paid to specifying meaningful informative prior distributions for these parameters.
Next, we apply our Bayesian inferential procedure to COVID-19 surveillance data collected in Orange County, California. To account for time-varying disease dynamics and surveillance effort, we analyze separately three 36 day long time periods: March 31-May 05, May 06-Jun 10, and June 11-July 16. The first period overlaps with the state-level stay-at-home order -the most stringent virus transmission containment effort implemented in California. The second period corresponds to a series of restriction relaxations implemented by the California four-stage reopening roadmap. During the third chosen period a mix of relaxations and restrictions occurred, with some businesses opening, universal masking orders, re-closures of in-door dining, and a re-instatement of the stay-at-home order. Our main interest is in quantifying differences in transmission dynamics and surveillance efforts among the three chosen periods. Three sets of prior and posterior distributional summaries of all model parameters are available in Appendix Figures B1-B3. Below we walk the reader through the most important parameter estimation results in detail.
The upper-left plot of Figure 4 shows prior and posterior distributions of the basic reproduction number (R 0 ) for Orange County in the three time intervals. The basic reproduction number is estimated to be near or under 1.0 in April and May, around 1.0 in May and June, and above 1.0 in June and July. Note that since the effective reproductive number can be computed as R e (t) = S (t)R 0 and the proportion of susceptibles in Orange County was still close to 1.0 before June 2020, we expect R e (t) to be only slightly smaller than our calculated R 0 in the first two periods. Indeed, the middle plot in the top row of Figure 4 demonstrates that the posterior distributions of S (t) at the end of these periods are close to 1.0, while the third period ended with 2-8% drop in susceptibles. The posterior distributions of R e (t) are shown in the upper-right plot of Figure 4 . The leftmost plot in the middle row shows percent change of R e (t) between the first and second and second and third periods. Posterior medians of this percent change say that transmission increased by 18% in the second period and by 24% in the third period. It is of interest if R e > 1.0, because that would indicate that the number of new infections should be increasing in the future. We compute a Bayes factor in favor of the hypothesis that R e > 1 -BF R e >1 = Pr(M, Y | R e > 1)/Pr(M, Y | R e ≤ 1) -the probability of observing the data in hand under the assumption that R e > 1 divided by the probability of observing the data in hand under the assumption that R e ≤ 1. The Bayes factors for May 5, June 10, and July 16 are 0.22, 0.63, and 30, respectively. This indicates that in the latest time period the cases and deaths data are much more likely under the assumptions that R e > 1 than under the alternative of R e ≤ 1, providing strong evidence in favor of R e > 1 ( Kass and Raftery, 1995) . For the earlier time periods, the Bayes factors provide compelling evidence neither in favor of R e > 1 nor in favor of R e ≤ 1.
Prior and posterior distributions of the infection-to-fatality ratio (IFR), represented by the parameter η in our model, are depicted in the middle plot of the second row in Figure 4 . The most recent 95% Bayesian credible interval is (0.002-0.013), meaning that somewhere between 0.2% and 1.3% of all SARS-CoV-2 infections result in death. The estimated IFR for the most recent period is lower than the previous two periods, but all 3 posterior distributions largely agree with previously reported IFR plausible ranges (Verity and et. al, 2020; Meyerowitz-Katz and Merone, 2020) . The rightmost plot in the second row and the leftmost plot in the third row of Figure 4 depict the prior and posterior distributions of parameters related to testing guidelines. Recall that a setting with α 0 = 0 and α = 1.0 corresponds to testing individuals in the population uniformly at random. The posterior distributions of these parameters indicate that, over time, α 0 is decreasing, while α 1 is increasing, suggesting that testing is becoming closer to sampling uniformly at random. Still, even in the latest period, the value of α 0 is far from zero, which means that infected individuals are much more likely to get tested than susceptibles. Prior and posterior distributions of the death reporting probability -the model parameter we call ρ -are shown in the middle plot of the last row in Figure 4 . Our results show that 52% to 97% of all deaths are reported in each observation period, similarly to our prior assumption. We assess the under-reporting factor of cases by taking the ratios of the observed to estimated cumulative incidence counts at the end of each time period. The resulting posterior distributions are displayed in the lower-right plot of Figure 4 , showing that reported cases underestimate the actual number of infections by at most 10 fold.
We now turn to the estimated latent dynamics of the SARS-CoV-2 transmission in Orange County. First, we plot posterior medians and Bayesian credible intervals of the latent cumulative death counts (N I pD (t)) between March and July 2020, using three credibility levels: 50%, 80%, and 95% (top row of Figure 5 ). We also extend the time horizon by 4 weeks into the future in all plots of this figure to provide our predictions about the SARS-CoV-2 spread. Reported death counts are shown as black bars in the same plots. When available, the reported counts for days in the prediction region are shown as light-gray bars. These bars show our predictions to be quite accurate in all the time periods, but with a tendency to underestimate in the initial time period, likely due to the ending of the stay at home order. Note that we expect that predicted latent deaths to be greater than reported deaths, due to underreporting.
The middle row of Figure 5 shows the posterior distributions of cumulative number of infections (N S E (t)) occurred in Orange County between March and July 2020. We show observed cumulative incidences as black , infection-to-fatality ratio (IFR), parameters of the beta-binomial observational model for cases, α 0 , and α 1 , and death reporting probability ρ, and case underreporting factor. For each parameter, we report three posterior distributions corresponding to the three observation periods. The shaded areas are mirrored posterior density plots. Black circles correspond to prior/posterior medians. Black vertical lines indicate 50%, 80%, and 95% Bayesian credible intervals. bars in the same plots. As before, the light-gray bars indicate observed incidence in the predicted period. The bottom row of Figure 5 shows progression of SARS-CoV-2 prevalence (total number of infectious individuals at one particular time, I e (t) + I p (t)) over the time period under study. These plots show that SARS-CoV-2 prevalence increased approximately 10 fold between April and July. The three prevalence plots disagree about the level of prevalence for overlapping intervals, by a factor of roughly 2. This disagreement is mostly likely a result of a model misspecification. For example, our unrealistic assumption of a constant infection rate, or equivalently R 0 , in each observation period may be at fault here.
Finally, we turn to model diagnostics and model-based forecasting of observable quantities. Let us first look at our model-based predictions of reported deaths in the top row of Figure 6 . Conceptually, you can think about these predictions as counts that are generated by first sampling all parameters of our model from the posterior distribution, plugging them into the mortality negative binomial distribution, and sampling from this distribution at each time point of interest. As in Figure 5 , we use three credibility levels in Figure 6 . The fact that during the observation time interval the reported death counts fall within 95% credible interval of their model-based predictive distributions indicates that our model fits the Orange County mortality data reasonably well. The credible intervals outside of the observation interval summarize our four week ahead probabilistic forecasts of mortality. When available, the observed deaths from the forecasting period are presented as lightgray bars. In the first time period, our forecast underestimated the actual number of observed deaths; however, our accuracy was much greater in the second and third periods.
Since forecasting cases is impossible without knowing how many tests will be conducted in the future, we focus on positivity fraction (cases divided by the total number of tests) instead. To generate predictive distribution of the positivity fraction we use model parameters, drawn from the posterior distribution, to generate realizations from the beta distribution shown in equation (4). Again, the observed positivity fractions fall within 95% credible intervals of their posterior predictive distribution. Our forecasts also perform well in the first time period, but are quite inaccurate in the second time period, where a major rise in testing positivity began just after the model training period. Our model also did not forecast decline in test positivity in the second half of July and beginning of August. Although our model did not forecast well test positivity in the second and third periods due to our constant R 0 assumption, moving the 36 day long observation period into the future would allow our model to pick up on increases and decreases in transmission.
We developed a Bayesian model that integrates information from incidence and mortality data. Our approach combines a realistic, but not overly complex, ODE-based compartmental model of SARS-CoV-2 transmission dynamics and a carefully constructed surveillance model for cases and deaths. Importantly, our method accounts for variability in the number of SARS-CoV-2 diagnostics tests across time, thus ensuring that we do not confuse increases in testing with increases in incidence. In supplementary simulations, we have also shown that failing to account for changes in testing can lead to poor estimation of the basic reproduction number R 0 . In our application of the model to Orange County, CA surveillance data shows that the SARS-CoV-2 transmission steadily increased after the California state-wide stay-at-home order ended. The model allowed us to estimate IFR, magnitude of incidence and death underreporting, and changes in prevalence in Orange County, CA. Retrospective analysis showed that our model produced reasonable forecasts, but they were off during periods when SARS-CoV-2 transmission dynamics and mitigation policies changed rapidly.
Our primary focus in this work was on developing a framework for integrating multiple data streams into a transmission model. However, there are a number of extensions we could pursue to improve the realism of the assumed transmission dynamics and strengthen the model's forecasting skill. Our model assumes that the population of interest is well mixed and that all individuals in the population infect others and get infected at the same per capita rate. In fact, the actual SARS-CoV-2 transmission process is much more complex because individuals come into contact with each other based on their geographical and social network proximity. Furthermore, it is well established that COVID-19 disease progression process depends on the individual's age and possibly other characteristics (Kim and et. al, 2020; Bhargava and et. al, 2020; Petrilli and et. al, 2020) . Fortunately, compartmental models can be extended to account for these complexities. For example, we can stratify each model compartment by age and geographical location, as is commonly done in epidemiological modeling (Li and Brauer, 2008; Van den Driessche, 2008) . Fitting models with complex transmission heterogeneities is challenging because of a large number of model parameters we need to estimate. Hence, fitting increasingly granular models requires richer surveillance data, e.g., cases, negative tests, and deaths stratified by age and geographical location. For this reason, it is often preferable to work with line lists of positive and negative tests, which often include the geographic and demographic information needed to fit such models.
Changes in control/mitigation measures and in human behavior make SARS-CoV-2 transmission dynamics time-varying. In our Orange County, CA surveillance data analysis, we have addressed this complication by breaking the observation time interval into smaller sub-intervals and by allowing all our model parameters to vary across these sub-intervals. However, a more optimal solution would keep most of the parameters constant across time, while allowing a small number of parameters to be time-varying. To address changes in human behavior and resulting changes in contact rates, it is important to allow for temporal variability in R 0 . Anderson and et. al (2020) and Miller and et. al (2020) use parametric approaches to model effects of mitigation measures on R 0 . It would be interesting to try Bayesian nonparametric approaches within our framework (Xu et al., 2016; Tang et al., 2019) , or a semi-parametric approach that combines together parametric and non-parametric components. Other aspects of the model where time-varying dynamics could reasonably be incorporated are the IFR and the mean test positivity and death detection rate in the surveillance model. Additional extensions to time-varying dynamics should be implemented thoughtfully, with particular attention paid to controlling model complexity and to model evaluation.
Our general modeling framework could be made useful for anticipating and preparing for a surge in COVID-19 hospitalizations by adding compartments for hospitalization and intensive care admission. Daily counts of the numbers of people who are hospitalized or receiving intensive care are available in Orange County, CA and in many other geographic locations. Adding additional compartments corresponding to hospitalization and intensive care would allow us to incorporate these counts into our surveillance model.
In this paper, we have sidestepped the thorny issue of reporting delays by restricting our analyses to time periods in which the data have stabilized. Hence, our analyses should be robust to reporting delays so long as we have either "run out the clock" on the extent of the delays or there are reporting delays do not differ between positive and negative COVID-19 diagnostic tests. A useful set of extensions that would make our model more useful for real-time surveillance involve estimating the reporting delay distribution (Höhle and an der Heiden, 2014; Stoner and Economou, 2019) and using this distribution in our surveillance model.
Finally, we would like to point out that our deterministic representation of the latent epidemic process could be substituted for a fully stochastic model where the latent epidemic is represented as a Markov jump process, albeit with some loss of computational efficiency. In our large population setting, this could be achieved via simulation-based methods (Bretó et al., 2009; Andrieu et al., 2010; Dukic et al., 2012) , data augmentation (Pooley et al., 2015; Nguyen-Van-Yen et al., 2020) , or a variety of approximations of the latent stochastic epidemic process (Lekone and Finkenstädt, 2006; Cauchemez and Ferguson, 2008; Fintzi et al., 2020) . Scaling our model to state or national level could be done by analyzing multiple counties independently or by building a Bayesian hierarchical model that would allow borrowing information among counties. An even more ambitious undertaking would be allowing importation/exportation events across county lines, as was done by Pei and Shaman (2020) . We hope that our methodology and other works in this spirit, along with better quality of surveillance data, will contribute to statistically rigorous COVID-19 epidemic modeling.
Miller, A. and et. al (2020), "Mobility trends provide a leading indicator of changes in SARS-CoV-2 transmission," medRxiv. Supplementary Materials for "Using multiple data streams to estimate and forecast SARS-CoV-2 transmission dynamics, with application to the virus spread in Orange County, California"
To validate the model we generated one thousand simulated sets of observed cases and tests based on actual tests performed in Orange County, CA during April 14th to May 19th. This period was chosen because there was a large increase in the number of tests during this time period. An example of one of the one thousand data sets we simulated is displayed in Figure A1 . We then fit two separate versions of the model to each of the simulated data sets. The first version of the model is the model described in the main text, to which we refer to as the "Test-Aware" model. The second model used a negative binomial distribution to model cases, much as a negative binomial distribution is used to model deaths, and did not account for number of tests. This second model will be referred to as the "No Tests" model. The parameters used to simulate the data sets, as well as the priors used in fitting the models are described in the tables below (Tables A1, A2 , and A3). Summaries of priors and posteriors of all model parameters from both models fit to the sample data set in Figure A1 are displayed in Figures A4 and A5 .
For each parameter and model we calculated three metrics to assess model performance. To examine bias we used relative absolute difference ( Figure A2 ), calculated as the absolute difference between the posterior median and the true parameter value, divided by the true parameter value. We also calculated relative credible interval width ( Figure A3 ). Relative credible interval width was calculated as the difference between the 97.5 and 2.5 quantiles of the posterior distribution divided by the true parameter value. Finally, we summarized coverage of model parameters, defined as the total number of simulations for which the 95% credible interval contained the true parameter value divided by the total number of simulations times one hundred, in Table A -4. Overall the Test-Aware model has lower median relative absolute error ( Figure A2 ), with median credible interval widths occasionally wider than the No Tests model ( Figure A3 ). Across the board, the Test-Aware model has coverage that is as good as or better than the No Tests model (Table A4) . Of particular interest is the coverage for the basic reproductive number R 0 . The Test-Aware model covered the true value of R 0 in 98% of simulations, while the No Tests model coverage of R 0 was only 8%, supporting our decision to incorporate testing data into our surveillance model. These priors are used only for the initial model fit at the beginning of the epidemic. In models fit to later stages of the epidemic, method of moments estimators are used to create a prior based on the posterior of the relevant compartments on the relevant day in the previous model. 
In the main analysis we report some summaries of parameter posteriors, but did not have room in the main results section to summarize every model parameter. We provide visual summaries of posteriors for all thirteen model parameters for each time period in Figures B1, B2 and B3.
We conducted four sensitivity analyses in order to see how our results changed depending on the specified priors. In each additional analysis we changed only one prior from the main analysis. In two analyses we changed the prior for R 0 to be lognormal(log(2.5), .5) or lognormal(log(.53), .78). In two analyses we changed the prior on the fraction of people who were initially infected so that the prior mean was either twice as large or half as large as the mean used in the main analysis. In our model construction fraction of people initially infected is equal to 1 minus fraction of people who are initially susceptible. Overall, we found that changing the prior on R 0 changed the magnitude of our inference on R 0 somewhat, but generally didn't affect key results. , infection-to-fatality ratio (IFR), and parameters of the beta-binomial observational model for cases, α 0 , and α 1 . For each parameter, we report three posterior distributions corresponding to the three observation periods. The shaded areas are mirrored posterior density plots. Black circles correspond to prior/posterior medians. Black vertical lines indicate 50%, 80%, and 95% Bayesian credible intervals. Sensitivity analysis for Figure 4 where R 0 ∼ lognormal(log(2.5), .5).
Changing the prior on the fraction of people who were initially infected changed our results on prevalence noticeably, and highlights the importance of this prior in our model. The results of the sensitivity analyses are shown in the figures below. , infection-to-fatality ratio (IFR), and parameters of the beta-binomial observational model for cases, α 0 , and α 1 . For each parameter, we report three posterior distributions corresponding to the three observation periods. The shaded areas are mirrored posterior density plots. Black circles correspond to prior/posterior medians. Black vertical lines indicate 50%, 80%, and 95% Bayesian credible intervals. Sensitivity analysis for Figure 4 where mean initial infected is halved. Figure B1 where R 0 ∼ lognormal(log(2.5), .5). Figure B1 where mean initial infected is doubled. Figure B1 where mean initial infected is halved. distributions for all model parameters for the model fitted to data from 5/6/2020-6/10/2020. Black dots represent medians, and black lines are 95% credible intervals. Sensitivity analysis for Figure B2 where R 0 ∼ lognormal(log(2.5), .5). distributions for all model parameters for the model fitted to data from 5/6/2020-6/10/2020. Black dots represent medians, and black lines are 95% credible intervals. Sensitivity analysis for Figure B2 where mean initial infected is doubled. distributions for all model parameters for the model fitted to data from 5/6/2020-6/10/2020. Black dots represent medians, and black lines are 95% credible intervals. Sensitivity analysis for Figure B2 where mean initial infected is halved. Figure B3 where R 0 ∼ lognormal(log(2.5), .5). Figure B3 where where mean initial infected is doubled. Figure B3 where where mean initial infected is doubled.