key: cord-0220017-zos2d0qf authors: Bayer, Damon; Fintzi, Jonathan; Goldstein, Isaac; Lumbard, Keith; Ricotta, Emily; Warner, Sarah; Busch, Lindsay M.; Strich, Jeffrey R.; Chertow, Daniel S.; Parker, Daniel M.; Boden-Albala, Bernadette; Dratch, Alissa; Chhuon, Richard; Quick, Nichole; Zahn, Matthew; Minin, Vladimir N. title: Semi-parametric modeling of SARS-CoV-2 transmission in Orange County, California using tests, cases, deaths, and seroprevalence data date: 2020-09-06 journal: nan DOI: nan sha: 6eed788673a3e045d562f14cce0c54a524492bbc doc_id: 220017 cord_uid: zos2d0qf Mechanistic modeling of SARS-CoV-2 transmission dynamics and frequently estimating model parameters using streaming surveillance data are important components of the pandemic response toolbox. However, transmission model parameter estimation can be imprecise, and sometimes even impossible, because surveillance data are noisy and not informative about all aspects of the mechanistic model. To partially overcome this obstacle, we propose a Bayesian modeling framework that integrates multiple surveillance data streams. Our model uses both SARS-CoV-2 diagnostics test and mortality time series to estimate our model parameters, while also explicitly integrating seroprevalence data from cross-sectional studies. Importantly, our data generating model for incidence data takes into account changes in the total number of tests performed. We model transmission rate, infection-to-fatality ratio, and a parameter controlling a functional relationship between the true case incidence and the fraction of positive tests as time-varying quantities and estimate changes of these parameters nonparameterically. We apply our Bayesian data integration method to COVID-19 surveillance data collected in Orange County, California between March, 2020 and March, 2021 and find that 33-62% of the Orange County residents experienced SARS-CoV-2 infection by the end of February, 2021. Despite this high number of infections, our results show that the abrupt end of the winter surge in January, 2021, was due to both behavioral changes and a high level of accumulated natural immunity. SARS-CoV-2 is a novel human coronavirus that is associated with high morbidity and mortality (Cummings et al., 2020; Wu and McGoogan, 2020; Song et al., 2020) . Like other human coronaviruses, SARS-CoV-2 is transmitted person to person through close contact via respiratory droplets and airborne particles, with the latter being a major concern in crowded indoor settings and around activities that generate aerosols (WHO, 2021) . To suppress viral transmission, many geographic/administrative regions have implemented mitigation measures (business and educational institutions closures, physical distancing, requiring facial masks, etc.). To estimate the effect of these policies and their relaxations, we develop a method that can estimate time varying parameters of the SARS-CoV-2 transmission dynamics from time series of positive and negative diagnostic tests detecting presence of SARS-CoV-2, deaths due to COVID-19, and seroprevalence data that quantifies abundance of antibodies developed in response to a SARS-CoV-2 infection. 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) . At later pandemic stages, transmission dynamics modeling helped evaluate intervention effectiveness (Knock et al., 2021) and to quantify transmission advantages of genetic variants (Davies et al., 2021) . 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 (Anderson et al., 2020; Miller et al., 2020; Morozova et al., 2021; Irons and Raftery, 2021) . However, neither national nor subnational/regional modelers fully integrate all surveillance data available to them, because inclusion of each additional data source leads to an increase in model complexity, which complicates statistical inference and reduces computational efficiency of this inference. Incorporating case incidence data into inference proved particularly problematic, because a data generating model for cases needs to take into account preferential sampling of symptomatic individuals and dependence of case counts on 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. 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 first innovation is the model for cases, where we use a flexible betabinomial 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. Our second innovation is nonparametric estimation of time-varying parameters that control both the transmission model and surveillance model. By allowing our model to adapt to temporal changes in transmission and surveillance, we are able to identify how changes in policy affected the near-term progression of the outbreak. To benchmark the ability of our model to capture temporal trends in transmission dynamics we use simulated data to compare our estimation of changes in the effective reproductive number to analogous estimates produced by simpler semi-parametric methods that do not attempt to model multiple compartments and use only cases as their data source (Cori, 2019; Bhatt et al., 2020) . This comparison shows that failing to account for fluctuations in the number of tests performed can result in misleading inferences about effective reproduction numbers, a critical quantity in infectious disease epidemiology. Our data integration approach allows our model to provide a much more detailed picture of the spread of an infectious disease beyond the effective reproduction number, including estimating the time varying infection fatality ratio, the total number of infected individuals, and changes in testing strategy. Data integration further allows us to produce reasonable short term forecasts of deaths and testing positivity. We demonstrate these enhanced capabilities by fitting our compartmental model to surveillance data collected in Orange County, California -the sixth most populous county in the United States of America (U.S.A.), with an estimated 3.2 million inhabitants as of 2019 (United States Census Bureau, 2020). We analyze Orange County surveillance data collected between March 3, 2020 and February 28, 2021, prior to the start of widespread vaccine availability. We find both basic and effective reproductive numbers varied widely during the first year of the pandemic, which is expected in light of implementations and subsequent relaxations of mitigation measures. The infection-to-fatality ratio was also changing through time, increasing notably during the fall-winter surge. On a positive note, changes in the parameter controlling our surveillance model for cases reflects improvements in tests availability, resulting in lower preference for testing symptomatic individuals. Our use of an Orange County seroprevalence survey allows us to estimate that 33-62% of Orange County residents had experienced SARS-CoV-2 infection by March of 2021. Finally, we demonstrate that our model produces reasonable short term (up to 4 weeks ahead) probabilistic forecasts of deaths and test positivity. We start with time series of daily numbers of SARS-CoV-2 diagnostic 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 weekly intervals. Figure 1 shows such a collection of aggregated time series for Orange County, CA, corresponding to the observation period spanning days between March 30, 2020 and February 28, 2021. 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 OCHCA. 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, because most of these repeat positive tests do not represent new infections and stem from follow up testing during patient treatment. 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 the weeks are indexed by l and L = 48 is the total number of weeks. 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. Additionally, we use seroprevalence data from Bruckner et al. (2021) , which consists of 343 seropositive cases (Z l * ) among 2979 (U l * ) tests conducted between July 10 and August 16, 2020. For simplicity, we lump all seroprevalence test dates to a single time point -August 16, 2020 -corresponding to week l * = 20. To formulate the surveillance model for cases Y, deaths M, and seropositive cases Z l * , 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 5 compartments: S = susceptible individuals, E = infected, but not yet infectious individuals, I = infectious individuals, 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(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 ) = D(t 0 ) = 0 and normalize X(t 0 ) so that X(t 0 ) T 1 = 1, where 1 = (1, 1, 1, 1, 1) T . 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 (t), N IR (t), N ID (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 β(t) is the transmission rate, which varies over time, 1/γ is the mean latent period duration, 1/ν is the mean infectious period duration, and η(t) is the infection-to-fatality ratio (IFR), which varies over time. Note that the time varying transmission rate allows our model to capture effects of interventions and changes in human behavior. Time-varying IFR should capture changes in age profiles of infected individuals over time and other factors (e.g., stress on healthcare providers during surges). 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 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 , and I 0 as free model parameters. Hence, our model describes the forward time evolution of the population still at risk. 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) . 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 seroprevalence data and two time series: numbers of new cases and deaths reported during some pre-specified time periods (e.g., weeks), as well as cumulative cases. 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 l observed in time interval (t l−1 , t l ], where l = 1, . . . , L. Since our ODEs track the latent cumulative fraction of deaths N ID (t l ), we can compute ∆N ID (t l ) = N ID (t l ) − N ID (t l−1 ) -the latent fraction of the population that died in the interval (t l−1 , t l ]. We model the observed death count M l as a realization from the following negative binomial distribution: where N is the population size, µ D l and σ 2 l are the mean and variance of the negative binomial distribution, ρ ∈ [0, 1] is the mean overall death detection probability, δ l ∈ [0, 1] is the death reporting probability based on the time of data collection, and φ > 0 is an overdispersion parameter. Informally, our mortality model says that, on average, the observed number of deaths, M l , is a fraction of the true death count estimated by the model, N × ∆N I p D (t l ), with some noise due to underreporting, delayed reporting, and sampling variability. Next, we develop a model for the number of positive tests (cases), Y l , observed in time interval (t l−1 , t l ]. We start with a simple binomial model with per-test positivity probability ψ l : where T l is the number of COVID-19 diagnostic tests administered during time interval (t l−1 , t l ]. We could proceed to modeling ψ l 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 l−1 , t l ] follows the following Beta distribution: where κ is an overdispersion parameter and µ C l is the mean test positivity probability. We assume that mean test positivity odds is proportional to the unobserved odds of transitioning from exposed to infectious, where α l > 0. This functional form ensures that on average probability of detecting a SARS-CoV-2 infection grows with the population incidence. Parameter α l can be thought of as an effect of testing guidelines/practices. A model with α l ≈ 0 (i.e., µ C l ≈ ∆N I e I p (t l )) says that in interval (t l−1 , t l ] testing is done approximately by sampling individuals uniformly at random, so that the positivity probability over time interval l is equal to the fraction of the population that transitions from the latent to infectious state. As we increase α l above 0, 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 l ) = T l × µ C l . 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 l and ∆N EI (t l ), 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 where the variance under an analogous pure binomial model would be µ C l (1 − µ C l )/T l . 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 SARS-CoV-2 incidence and implicitly allows for heterogeneity in the mean test positivity probability. Finally, we model the number of observed seropositive cases Z l * among U l * tests with a binomial distribution: This simple model assumes the seroprevalence data comes from a high quality study based on random sampling, which does not exhibit the problems observed in the testing data. We are now ready to describe our inferential Bayesian procedure. First, we slightly re-parameterize our model by replacing β(t) with a basic reproductive number R 0 (t) = β(t)/ν. We parameterize each of our time-varying parameters, R 0 (t), η(t), α(t) as piecewise constant functions, where each vector defining the constants a priori follows a Gaussian Markov random field (GMRF). More precisely, we define the auxiliary vectors: which follow the Gaussian Makrov random field priors: and define the piecewise constant functions: Initial proportion of R 0 + D 0 + E 0 who are recovered Logit-normal(-7.9, 9.85e-05) 3.71e-04 (3.71e-04, 3.72e-04) D 0 Initial proportion of R 0 + D 0 who are deceased Logit-normal (-7.9, 9. In addition, we parameterize initial compartment fractions as S 0 , . This construction allows us to specify independent prior distributions for S 0 , I 0 , R 0 , and D 0 , preserving the sum-to-one constraint on the original initial compartmental fractions. Next, we collect all our model parameters into a vector θ = (S 0 ,Ĩ 0 ,R 0 ,D 0 ,R 0 , γ, ν,η, ρ, φ,α, κ, σ R 0 , σ η , σ α ). Our probabilistic construction described above implies that the likelihood function -probability of observing incidence, mortality, and seroprevalence data -can be written in the following way: where Pr(M l | θ), Pr(Y l | θ), and Pr(Z l * | θ) are Negative-binomial, Beta-binomial, and Binomial probability mass functions defined by equations (3), (6), and (7) respectively. We encode available information about our model parameters in a prior distribution with density π(θ). We assume that all univariate non-GMRF distributed 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 slice sampling implemented in the stemr R package (Fintzi and Bayer, 2021) . We used ten Markov chains run in parallel to draw a total of 35,570 posterior samples. 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/damonbayer/stemr_oc. To validate our model, we simulated one dataset with parameters given in Supplementary Table A1 . The resulting dataset is presented in Figure A1 . The number of tests at each time point is the exact same as in the Orange County data set, and the parameters were deliberately chosen to produce data similar to the Orange County data. Priors used for this model fit are the same as those used in the Orange County model (see Table 1 ). The posterior distribution of this model fit is presented in Supplementary Figures A2 and A3 . The posterior distribution appropriately captures each of the true parameters, with the exception of σ α , which is largely unchanged from the prior distribution, resulting in an underestimate. However, this does not negatively affect the model's ability to capture the true α l . Our simulation results do not reveal statistical or implementation problems. Fitting the simulated data took approximately five days to generate 35,570 posterior samples among 10 chains, totalling around 1200 CPU hours. This number of chains and posterior samples resulted in satisfactory effective sample sizes for posterior inference. MCMC details are presented in Appendix D. We decided against repeating the above simulation multiple times in order to examine frequentist properties of our Bayesian procedure. Although we usually find such studies informative and useful, computational cost and carbon footprint of such a study would be substantial. We applied two state-of-the-art methods for effective reproduction number estimation, EpiEstim and epidemia, to the same simulated dataset. These semi-mechanistic methods do not attempt to estimate the unobserved number of susceptible and recovered individuals and do not take changes in testing volume into account. See Supplementary Section B for a brief description of EpiEstim and epidemia statistical models. EpiEstim and epidemia estimated effective reproduction number trajectories, shown in Figure B4 , exaggerate fluctuations of the true effective reproduction number, overestimating it during surges and underestimating this parameter during transmission lulls. These results underscore the importance of accounting for fluctuations in testing volume even when using semi-mechanistic statistical modeling aimed only at the time-varying effective reproduction number estimation. Next, we apply our Bayesian inferential procedure to COVID-19 surveillance data collected in Orange County, California between March 3, 2020 and February 28, 2021. By the end of the modeling period, approximately 14% of Orange County residents were at least partially vaccinated. Because our model does not incorporate vaccination directly, we felt it would be inappropriate to continue modeling beyond February 2021. Throughout the modeling period, a variety of non-pharmaceutical interventions were enacted and sometimes lifted at the state, county, and city level. Notably, in-person school closures, indoor dining bans, and mask mandates were in effect for most or all of this period. Fitting the simulated data took approximately five days to generate 35,570 posterior samples from 10 chains, totaling around 1400 CPU hours. Although computationally expensive, these run times show that our model can be fit frequently enough to be useful for a real-time policy response. This number of chains and posterior samples resulted in a satisfactory effective sample size for posterior inference. Details about the posterior samples are presented in Appendix D. Our main interest is in understanding differences in transmission dynamics and surveillance efforts throughout this period. Since our model is highly parametric, we used existing knowledge of SARS-CoV-2 transmission dynamics to formulate informative priors for all model parameters that we list in Table 1 . We briefly highlight some of our assumptions. Our priors for the initial compartment sizes reflect our belief that the number of infections was small, but potentially underreported by a factor of 10, at the beginning of the pandemic. Since our observation period starts close to the the stay-at-home order taking effect in California, we assume that March 2020 basic reproduction number should be around 1.0 to reflect reduced contacts during this period. Lengths of latent and infectious periods a priori assumed to be 0.5 and 1 week respectively, with substantial variance. Based on Orange County, CA seroprevalence study, we assume initial infection-to-fatality ratio to be around 0.3% (Bruckner et al., 2021) . We compare reported deaths in Orange County, CA with estimates of U.S. county-specific excess mortality to set the prior for death reporting probability to be around 0.9 (Stokes et al., 2021) . Prior and posterior distributional summaries of all model parameters are available in Supplementary Figures C.1-C.8. These figures also contain results of our sensitivity analysis, where we examined the effects of our prior assumptions on our inference. Our main conclusion is our results are not too sensitive to reasonable prior perturbations. The upper-left plot of Figure 3 presents the posterior distribution of the basic reproductive number (R 0 ) for Orange County. Throughout the late spring and summer, the basic reproductive number is estimated to be slightly above 1.0, with some probability of being below 1.0 in the early fall. Beginning in October, the basic reproductive number begins to rise and nears 2.0 at the peak of the winter wave. This rise in the fall may be associated with the school re-openings that occurred around this time. Despite the high basic reproductive number throughout the modeling period, the upper-right plot of Figure 3 shows that the effective reproductive fell below 1.0 for much of the summer and again in late January, following the winter surge, allowing us to separate the effects of reducing the average community contact rate and accumulated infection-induced immunity. We take a short detour here and compare our inference of effective reproductive number with EpiEstim and epidemia results, presented in Figure B5 . As in our simulation study, EpiEstim and epidemia broadly agree with our estimated R t trend, but with a couple of important differences. First, EpiEstim and epidemia appear to be confused by the increase in testing availability during Spring and early Summer of 2020, producing unrealistically high effective reproductive number estimates. Second, both EpiEstim and epidemia results suggest that R t started increasing in mid Summer, while our model points to October as the start of the increase in transmission. We speculate that increase in testing volume is a culprit again here, because testing super-sites have opened in Orange County in July 2020, significantly improving access to SARS-CoV-2 diagnostic testing. Overall, EpiEstim and epidemia produces R t estimates that are higher than our model counterparts, especially at the peaks of disease transmission. We proceed with describing inference results for the other two time-varying parameters in our model: infectionto-fatality ratio and parameter α that governs the relationship between testing positivity and the true proportion of newly infected individuals in the population. The posterior distribution of the infection-to-fatality ratio is presented in the middle-left plot of Figure 3 . The IFR is estimated to be consistent over time, hovering around 0.3%, but our estimates are less certain near the end of the modeling period. This rise could have been caused by a combination of the overwhelmed healthcare system and the increasing prevalence of the Alpha variant at this time, which infections have been tied to more severe outcomes (Grint et al., 2021) . The middle-right plot and bottom plots of Figure 3 present three perspectives on testing policy and case detection: the posterior α, the weekly latent:case ratio, and the cumulative latent:case ratio. Generally, α, drifts lower over time, indicating that testing policy became less preferential toward selecting infected individuals as testing became more accessible. This trend is reversed slightly during the summer and winter waves, which is reflected in the decreasing weekly latent:case ratio during these times. The cumulative latent:case ratio also drifts lower over time, eventually arriving and a final cumulative latent:case ratio of 4:1 -8:1. We plot posterior medians and Bayesian credible intervals of the latent cumulative death counts (N ID (t)) between March 2020 and February 2021, using three credibility levels in left plot of Figure 4 . Reported death counts are shown as black circles in the same plot. The plot reflects and overall death reporting rate of 86% -95%. The center plot of Figure 4 shows the posterior distributions of cumulative number of infections (N S E (t)) occurred in Orange County, with the cumulative observed cases displayed as black circles. As in Figure 3 , this shows a cumulative latent:case ratio of 4:1 -8:1, with 1/3 -2/3 of all Orange County residents having been infected by the end of February 2021. From this plot, we also note that our posterior estimate of seroprevalence in mid-August of 2020 (11.0%-13.4%) closely matches the 11.5% estimate from (Bruckner et al., 2021 ) -this is not surprising since we explicitly used this seroprevalence data in our inference. The right plot of Figure 3 shows prevalence of SARS-CoV-2 infected individuals at a particular time, (E(t) + I(t)). At the peak of the winter wave, we estimate that 1.8% -6.1% of Orange County residents had an active infection at the same time. 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 left column of Figure 5 . Conceptually, one 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 the other figures, we use three credibility levels in Figure 5 , and observed values are displayed as black circles. 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. Next, we retrospectively produce one and four week ahead probabilistic mortality forecasts, both updated weekly as we march forward in time. Our one week forecasts are precise and well calibrated. The four week ahead forecasts predict the data well in times of relative stability, but exhibit poor performance when model time-varying parameters are changing rapidly, such as during the winter wave. 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. For both posterior predictive and forecast distributions, we use model parameters, drawn from the posterior distribution, to generate cases from the betabinomial distribution shown in equation (6) and then divide them by the total number of tests. In practice, during forecasting we would not know the future total numbers of weekly tests, but to get reasonably calibrated forecasts one only needs to get this number approximately right to predict positivity. The posterior predictive and forecasts for positivity in the right column of Figure 5 exhibit similar behavior to the deaths posterior predictive and forecasts: excellent fit in the posterior predictive, good forecasting performance in the one week ahead forecasts and lackluster performance in the four week ahead forecasts when some time-varying parameters are rapidly changing. However, we note that if we were performing this forecasting in real time, we would draw incorrect conclusions about the direction of transmission dynamics in some weeks, but would be back on track after at most a couple of weeks after these forecasting blunders. In practice, this would be acceptable because policy makers would not want to act on just a couple of alarming forecasts and would want to see forecasting consistency for some time before enacting a new policy aimed at reducing community transmission. In collaboration with the Orange County Health Care Agency, we have been producing weekly situation reports based on the model described in this paper. These situation reports continued until May 2021 when we decided to stop producing them due to lack of detailed vaccination data that was needed to extend our model to account for accumulation of vaccine induced immunity. The most recent reports are available at https://www.stat.uci.edu/ oc_covid_model/. The companion github repository contains the source code for producing the reports and all our archived reports dating back to March 2020: https://github.com/damonbayer/uci_covid_modeling2. We developed a Bayesian SARS-CoV-2 transmission model that integrates information from incidence, mortality, and seroprevalence data. Our approach combines an ODE-based SEIR compartmental model of SARS-CoV-2 transmission dynamics and a carefully constructed surveillance model for cases, deaths, and seroprevalence. 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. Another distinguishing feature of our approach is nonparametric modeling of changes in key transmission and surveillance model parameters. Since we are integrating multiple sources of information, we can afford to be fairly ambitious and to include three such parameters into our model. Changes in one of these parameters accounts for changing the strength of preferentially testing SARS-CoV-2 infected individuals, which helps us avoid an important source of potential bias when inferring transmission model parameters. We reconstruct latent dynamics of the two pre-delta COVID-19 waves in Orange County, CA and estimate that 33-62% of Orange County residents experienced SARS-CoV-2 infection by March 2021. Retrospective analysis shows that our model produces accurate and well calibrated one week ahead and reasonable four week ahead forecasts, but the latter lack accuracy during periods when SARS-CoV-2 transmission dynamics and mitigation policies change 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 other characteristics (Kim et al., 2021; Bhargava et al., 2020; Petrilli et al., 2020) . Similarly, all transmission model parameters may depend on the vaccination status of an individual. Fortunately, compartmental models can be extended to account for these complexities. For example, we can stratify each model compartment by age, vaccination status, and geographical location, as is commonly done in epidemiological modeling (Li and Brauer, 2008; Van den Driessche, 2008) . We have addressed changes in control/mitigation measures and in human behavior by nonparametrically modeling variability of some of the SARS-CoV-2 transmission model parameters across time. Anderson et al. (2020) and Miller et al. (2020) use parametric approaches to model effects of mitigation measures on R 0 . It would be interesting to try a semi-parametric approach that combines together parametric and non-parametric components, which would allow us to include indicators of human behavior (e.g., mobility data as in (Miller et al., 2020) ) into our inference and forecasting. 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, 2020) 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 simulationbased 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; . 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 et al. (2021) . We hope that our methodology and other works in this spirit, along with better quality of surveillance data, will provide us with better predictive analytics tools when the next pandemic strikes. Supplementary Materials for "Semi-parametric modeling of SARS-CoV-2 transmission in Orange County, California using tests, cases, deaths, and seroprevalence data" We performed a small simulation study on one data set in order to validate our models. We purposely chose parameter values that resulted in data similar to the Orange Country data used in the main text. Exact values for these parameters are presented in Table A1 . The resulting data is plotted in Figure A1 . Figures A2 and A3 summarize the posterior distribution of model parameters and their functions (e.g., ODE trajectories). Commentary on these results is presented in Section 3 of the main text. As a comparison, we apply the widely used EpiEstim (version 2.2-3) method (Cori, 2019) and the newly developed epidemia (Bhatt et al., 2020 ) (version 0.7.0) method to both simulated and real data used in this study. In contrast to the compartmental model used in this paper, these methods have related the mean of current incidence to a weighted sum of previous incidence and the effective reproduction number R t . Let I t be the incidence at time t, R t be the effective reproduction number at time t, and g(t) be the pdf of the generation time distribution (the time between an individual becoming infected and infecting another individual; under the compartmental model framework this is usually taken to be equivalent to the sum of the latent period and the infectious period). Then the mean relationship used is: The model for EpiEstim is then EpiEstim then allows the user to specify a time window during which R t is assumed to be constant, for this analysis, we follow Gostic et al. (2020) and use single day windows. Generation time distributions are assumed to be gamma distributions with shape parameter 2 and means equal to the true mean generation time (in the case of the simulated data) or 11.5 days for the real data analysis. Posterior samples are drawn from the conjugate posterior derived from the above model. R package epidemia builds on the EpiEstim framework by incorporating an observation model, modeling the effective reproduction number as a random walk, and modeling unobserved incidence as a random variable (log-normal for this analysis). τ ∼ exp(λ) Hyperprior for unobserved incidence I ν ∼ exp(τ) Prior on unobserved incidence ν days before observation Prior on variance parameter for incidence I t |I ν , . . . , I t−1 ∼ Log-Normal(R t s