Confidence in the dynamic spread of epidemics under biased sampling conditions
CConfidence in the dynamic spread of epidemicsunder biased sampling conditions.
James D. Brunner, PhD and Nicholas Chia, PhD Mayo Clinic Center for Individualized Medicine, Department of Surgery
Corresponding author:James D. Brunner, PhD Email address: [email protected]
ABSTRACT
The interpretation of sampling data plays a crucial role in policy response to the spread of a disease during an epidemic, suchas the COVID-19 epidemic of 2020. However, this is a non-trivial endeavor due to the complexity of real world conditions andlimits to the availability of diagnostic tests, which necessitate a bias in testing favoring symptomatic individuals. A thoroughunderstanding of sampling confidence and bias is necessary in order make accurate conclusions. In this manuscript, weprovide a stochastic model of sampling for assessing confidence in disease metrics such as trend detection, peak detection,and disease spread estimation. Our model simulates testing for a disease in an epidemic with known dynamics, allowing us touse Monte-Carlo sampling to assess metric confidence. This model can provide realistic simulated data which can be used inthe design and calibration of data analysis and prediction methods. As an example, we use this method to show that trends inthe disease may be identified using under biased samples each day, and an estimate of disease spread can be madewith additional − unbiased samples each day. We also demonstrate that the model can be used to assess moreadvanced metrics by finding the precision and recall of a strategy for finding peaks in the dynamics. INTRODUCTION
Policy decisions in the face of an epidemic rely on perceptions of the population dynamics of an infectious disease, for examplewhether cases are growing or shrinking (Centers for Disease Control and Prevention, 2020; Occupational Safety and HealthAdministration, 2020; Lee et al., 2020). If everyone were tested every day, then this would be a matter of looking at the trend inthe total number of positive tests per day. However, this scenario is unrealistic. In actuality, the population sampling neededto track an epidemic in a community will depend on the nature of the question we would like to answer. Sometimes, thesequestions are in conflict with each other. For instance, the primary goal of healthcare providers is to identify infected patientsin the hospital or clinical setting so that appropriate treatment and protective measures may be prescribed, while at the otherextreme the epidemiologist concerned with infection prevention within the population may be interested in determining thenumber of infected individuals in the population in order to focus efforts at limiting disease spread. Clearly, the former isvery targeted towards patients showing up at a clinic with certain symptoms while the latter requires broad testing of bothsymptomatic and asymptomatic populations.The reality is that testing needs to serve multiple purposes with a finite number of tests. In particular, the COVID-19pandemic and response from world leaders has shed light on the need for a better understanding of community infection dataand how to use it, both for decision making and media reporting. Some of the particular challenges are the significant proportionof asymptomatic carriers of the disease (Bai et al., 2020; Mizumoto et al., 2020; Poll´an et al., 2020) and the changing availabilityof the testing. Both have resulted in testing that is strongly biased towards infected individuals and not representative of theproportion of cases in the population. In this manuscript, we focus on two intermediate use cases — individuals and businessesthat need to estimate risk, i.e., the probability that an infected person will be present in a given situation, and public policymakers that need to understand changing trends in the spread of the disease. We show that this can be done with a combinationof biased and unbiased sampling that requires many fewer tests to be performed every day, but importantly must include thenumber of negative tests in addition to the number of positive cases that is more widely reported. Notably, this is in line withthe World Health Organization’s global surveillance guidelines, which include reporting of total tests so positive percentage canbe determined (World Health Organization, 2020).The purpose of this manuscript is to introduce a method for assessing confidence in conclusions made from biased samplingof the spread of an epidemic, and therefore providing a tool in calibration of data analysis and prediction methods. We beginwith a calculation of the number of tests needed to identify a significant portion of infected individuals in a given day. We then a r X i v : . [ q - b i o . P E ] J u l escribe a stochastic dynamical model that simulates testing over the course of an epidemic with known dynamics. We thenshow how we can use a given model of epidemic dynamics to investigate the amount of testing needed to estimate diseasespread and trends in disease. We further use our approach to simulate testing with variable bias and error and investigate theroles of bias and error in testing. We find that amount of testing needed to identify most infected people in a population of300 million (approximately the population of the United States), is extremely high. On the other hand, we show that trendsin the spread of the disease can be accurately identified by sign (i.e., positive or negative) with less than 10000 biased testsper day. We show that approximately 1000 additional unbiased tests can be used to estimate the bias in testing. This can beused to estimate the extent of disease spread in the community on a daily basis. Our approach can also be used to assess thereliability of many data analysis techniques. We demonstrate this by assessing a strategy for finding peaks in the dynamics ofthe outbreak by using a smoothed numerical derivative of the data. Finally, we show the importance of understanding biasunder the conditions of limited testing by examining COVID-19 data from within the United States of America.As a special note, we emphasize that while many important efforts are being made to model the spread of COVID-19 anddetermine how testing can be used to reduce that spread (Reich et al., 2020; Piguillem and Shi, 2020; Alvarez et al., 2020;Chowdhury et al., 2020), our approach does not attempt to predict disease spread. Instead, we are testing confidence in dataanalysis sampled from known dynamics. In other words, rather than trying to predict the future, our work focuses on estimatingconfidence in current trends under non-ideal conditions. METHODS
Sampling all infected
Sampling the population in order to identify all or some large proportion of the infected individuals will require a large amountof testing. If we assume that testing is done in a single day, the proportion of infected in the population is roughly constantin that day, and no person is tested more than once, the number of positive cases will follow a hypergeometric distributionbased on the number of cases in the population, the bias in the testing, and the number of tests performed. We can compute thecumulative distribution function, and so compute the number of tests needed so that the probability that some number of casesare found.We take testing to be a process of sampling with replacement in a population of size T . To do this, we must compute thenumber of infected people in the population, which is given by rT . In principle, we may use a hypergeometric distributionwith population T = rT + ( − r ) T and infected rT . However, this approach assumes no bias in the testing, meaning that theprobability of each individual (symptomatic or not) being tested is uniform. In the real world, symptomatic individuals aremuch more likely to be tested for a disease. We therefore introduce non-dimensional a bias parameter B and take the apparentnumber of infected to be rBT , so that the apparent total population is rBT + ( − r ) T . For B >
1, this biases the testing towardsinfected individuals, representing the fact the symptomatic individuals are more likely to be tested, and also more likely tobe infected. The rest of this paper considers evaluating the infected population under biased testing conditions as a means tounderstand trends in real data.In our simulations of biased testing, we use a bias on the order of B =
10. However, this is done only to illustrate themethod and should not be taken as an estimate of the bias from real data. We discuss below how bias may be estimated from acombination of biased and unbiased samples.
A stochastic model for disease sampling.
We developed a stochastic model to simulate the sampling of a population that is undergoing an epidemic with known dynamics.That is, given an underlying set of dynamics tracking asymptomatic infected individuals, symptomatic infected individuals, andnon-infected individuals, we simulate testing members of this population for the disease. Let I ( t ) , I ( t ) , and H ( t ) represent thenumber of asymptomatic infected individuals, symptomatic infected individuals, and non-infected individuals in the population,respectively, for t ∈ [ , T ] . Our model does assume that the dynamics of the disease spread are some known functions of threecompartments. In practice, such dynamics are often simulated by simple systems of ordinary differential equations. However,this need not be the case, and dynamics can even be given directly as functions of time.Tests are assumed to be carried out according to a Poisson point process (see, e.g. (Klenke, 2014; Anderson and Kurtz,2011) for a detailed introduction) with (possibly time varying) intensity function λ ( t ) . The intensity function λ ( t ) can beinterpreted as the rate at which members of the population are tested for the diseased which is assumed to be known. Forexample, an increase in test availability would be reflected in this model with an increase in λ ( t ) . Additionally, the incidence ofpositive tests being carried out is itself a Poisson point process with intensity λ + ( t ) , and likewise the incidence of negative testsis a Poisson point process with intensity λ − ( t ) , subject to the relation λ ( t ) = λ + ( t ) + λ − ( t ) . In this manuscript, λ ( t ) is takento be constant unless otherwise noted.Each time a test is performed, we determine the status of the person tested. Under the assumption of unbiased testing,we categorize a person into one of three pools. The probability of a test result depends on the respective proportion of the opulation which belongs to each pool. That is, the probability the unbiased test is performed on an asymptomatic infectedperson is P ( this testee is asymptomatic infected ) = I ( t ) I ( t ) + I ( t ) + H ( t ) . (1)Likewise, the probability that an unbiased test is performed on a symptomatic infected person is P ( this testee is symptomatic infected ) = I ( t ) I ( t ) + I ( t ) + H ( t ) (2)and the probability that an unbiased test is performed on a non-infected person is P ( this testee is non-infected ) = H ( t ) I ( t ) + I ( t ) + H ( t ) . (3)We also account for the possibility that a test may give a false-positive or false-negative with some constant probability. Let ε ∈ [ , ] be the false-negative probability of the test and ε ∈ [ , ] be the false-postive probability. Then, for any given testwe can combine eqs. (1) to (3) to see the following: P ( positive test ) = ( − ε ) P ( person is infected ) + ε P ( person is not infected ) (4) P ( negative test ) = ε P ( person is infected ) + ( − ε ) P ( person is not infected ) (5)Testing for disease in a population is not done uniformly at random. Instead, an individual displaying symptoms of thedisease is much more likely to be tested for it than one who is not. We may model a bias in testing by adjusting eqs. (1) to (3).Let B ( t ) be some function of time, with B ( t ) ≥ t ∈ [ , T ] . Then, we can reflect the bias of the testing procedure byre-weighting the population for each test performed. We replace eqs. (1) to (3) with the following: P ( this testee is asymptomatic infected ) = I ( t ) I ( t ) + B ( t ) I ( t ) + H ( t ) (6) P ( this testee is symptomatic infected ) = B ( t ) I ( t ) I ( t ) + B ( t ) I ( t ) + H ( t ) (7) P ( this testee is non-infected ) = H ( t ) I ( t ) + B ( t ) I ( t ) + H ( t ) (8)and combine eqs. (6) to (8) with eqs. (4) and (5) to determine the probability of a single test result. The result is that the numberof positive and negative tests that have been carried out up to time t are each non-homogeneous Poisson point processes withintensity functions λ + ( t ) = λ ( t ) (cid:20) ( − ε ) (cid:18) I ( t ) I ( t ) + B ( t ) I ( t ) + H ( t ) + B ( t ) I ( t ) I ( t ) + I ( t ) + H ( t ) (cid:19) + ε (cid:18) H ( t ) I ( t ) + B ( t ) I ( t ) + H ( t ) (cid:19)(cid:21) (9) λ − ( t ) = λ ( t ) (cid:20) ε (cid:18) I ( t ) I ( t ) + B ( t ) I ( t ) + H ( t ) + B ( t ) I ( t ) I ( t ) + B ( t ) I ( t ) + H ( t ) (cid:19) +( − ε ) (cid:18) H ( t ) I ( t ) + B ( t ) I ( t ) + H ( t ) (cid:19)(cid:21) . (10)We note that the model described above essentially assumes an infinite total population size. Practically, this means thattesting is done on an insignificant proportion of the population, or equivalently that members of a population are immediatelyeligible to be re-tested after being tested. In appendix E, we describe a modification for this model which accounts for smallpopulation size and non-immediate retesting. Simulation with the Stochastic Simulation Algorithm
Both the initial model described in this manuscript and the model adjusted for small populations can be written as the sumof Poisson point processes with time-varying intensities. They can therefore be simulated using a slight adjustment to theStochastic Simulation Algorithm (Gillespie, 1976, 1977). This adjustment accounts for possible changes in the intensity unctions between points in the Poisson processes, for example from variations in λ ( t ) or B ( t ) . This adjustment is made bychoosing event times according to the maximum values of any time-varying functions, and allowing for the possibility of anon-event at each event time, a procedure often called thinning (Asmussen and Glynn, 2007; Anderson and Yuan, 2019).Simulation of the models as described allow us to perform Monte-Carlo estimations of the confidence that can be assumedin the calculation of various statistics from data. We perform Monte-Carlo estimation by repeatedly simulating sampling overthe course of an epidemic with given dynamics and determining the success rate or average error of in determining a metricfrom simulated sampling when compared to determining the same metric from the known underlying dynamics. Estimating trends in disease spread.
To assess trends in data simulated according to our model, we discretize the time interval [ , T ] into evenly spaced intervals(e.g. into single day increments) ending at times t , t , ..., t N = T . We then compute positive-test proportions for these intervals,simply defined as the proportion of tests carried out within the interval that were positive. This allows us to make sense of thesimulated data even as testing capacity λ ( t ) varies in time.We estimate N -day trends in disease spread using a linear least-squares fit to N consecutive days of simulated positive testproportions. We define the linear trend of the simulated data to be the slope of this fitted line. This can then be compared to alinear fit to the infected proportion I ( t ) + I ( t ) I ( t ) + I ( t ) + H ( t ) over the same time interval, computed using time-discretized dynamics. Finding peaks in disease spread.
We attempt to find peaks in the data by estimating the time derivative of the positive-test proportions simulated. We then identifypeaks as points at which the derivative crosses from positive to negative. That is, we estimate the change in true positive-testproportion from day to day, and identify when this proportion stopped increasing and began to decrease.To estimate the time derivative of positive proportions, we first compute a numerical derivative over each time interval. Wethen blur this discretized derivative to reduce noise (and therefore false peaks) using a one-dimensional Gaussian filter.
Underlying dynamics
Our model is designed to simulate sampling of an epidemic with any non-negative underlying compartmental dynamics. Thismeans that the number of healthy, symptomatic infected, and asymptomatic infected members of a population can be anynon-negative known functions of time. As a consequence of this feature, some set of known dynamics must be either chosendirectly or generated by another dynamic model. Confidence in metric sampling is then measured against the known dynamics.To demonstrate our method, we use the SIR model (Edelstein-Keshet, 2005; Hethcote, 2000; Kermack and McKendrick,1927) with a time-variable rate of disease spread to generate underlying dynamics, as well as a similar compartmental modelthat allows for asymptomatic individuals, which we refer to as the SAIR model. See appendix D for details. These modelsrepresent a popular, simple choice of dynamic epidemic model with parameters often reported by the lay news media.
RESULTS
Sampling all infected
One obviously crucial role of testing of a disease outbreak is to identify patients for treatment and possible quarantine. Fromthat perspective, it is crucial to identify all or most infected members of a community. We therefore assess the number of teststhis would require, noting that these tests must be performed in a small enough time frame so that epidemic dynamics and biasin testing can be assumed constant. Using the cumulative distribution function of the hypergeometric distribution, we see thatin a community of 300 million with 5% infection, we need approximately 27 , ,
300 unbiased tests in a short time interval tohave 90% confidence that we have found 90% of the cases. For higher bias, fewer tests are needed, as seen in fig. 1. This islikely an intractable number of tests to be performed in a short enough time interval (perhaps 1 day) so that epidemic dynamicsand bias in testing can be assumed constant. In fact, during the COVID-19 epidemic of 2020, about 500 , − ,
000 testswere performed each day across the United States by the end of June (Lipton et al., 2020).
Simulated sampling
Rather than attempt to find all COVID-19 patients, we may wish to simply have an accurate estimate of the number of infected.In order to asses our ability to do this, we simulated sampling as described above with underlying dynamics generated by ODEmodels of outbreaks. We use the positive test proportion per day as our sampled variable, in order to account for variations in igure 1.
Number of samples needed to identify 90% of infected individuals with 90% confidence, computed using ahypergeometric distribution. In orange, we show the limiting case in which every person tested is infected, which can beinterpreted as infinite bias. In this case, 1 , ,
000 tests are necessary. igure 2.
Positive proportion of tests and scaled positive tests per day. Testing capacity is initially 1350 tests per day but risesto 2700 midway through the simulation.testing capacity. Figure 2 demonstrates that day-to-day variations in testing capacity can mask the real dynamics when onlypositive test count is considered. In contrast, the positive proportion of tests does captures the dynamics. In this simulation, weused a time-varying intensity λ ( t ) which was taken to be λ ( t ) = . + .
25 tanh (cid:18) t − (cid:19) . (11) Biases in testing
Biased sampling
With false positive and false negative rates equal to 0, biased sampling causes an overestimation in the proportion of a largepopulation that is infected. If tests are taken at random in the population (and so not biased), the proportion of positive tests willon average be the same as the proportion of infected individuals r ( t ) , which is the sum of the proportions of symptomatic andasymptomatic infected individuals: r ( t ) = I ( t ) + I ( t ) I ( t ) + I ( t ) + H ( t ) = I ( t ) I ( t ) + I ( t ) + H ( t ) + I ( t ) I ( t ) + I ( t ) + H ( t ) = r I ( t ) + r I ( t ) . (12)Biasing the tests towards symptomatic individuals is analogous to sampling a population with extra symptomatic individualsadded: r B ( t ) = Br I + r I + ( B − ) r I > r ( t ) (13) ampling per day900 2700 8100BiasedSAIR 0.1541 0.2198 0.2959BiasedSIR 0.1865 0.3144 0.5210ExactSAIR 0.1155 0.1924 0.3045ExactSIR 0.1228 0.2080 0.3651PerfectSAIR 0.1871 0.2566 0.3386PerfectSIR 0.2400 0.3924 0.6132UnbiasedSAIR 0.0827 0.1034 0.1568UnbiasedSIR 0.0818 0.1101 0.1783 Table 1.
Precision of peak finding for various samplingand dynamics. Sampling per day900 2700 8100BiasedSAIR 0.8955 0.9000 0.8475BiasedSIR 0.8870 0.9370 0.9740ExactSAIR 0.8565 0.9135 0.9370ExactSIR 0.8580 0.9045 0.9575PerfectSAIR 0.9060 0.8955 0.7990PerfectSIR 0.8985 0.9570 0.9890UnbiasedSAIR 0.8195 0.8505 0.8970UnbiasedSIR 0.8005 0.8545 0.9180
Table 2.
Recall of peak finding for various samplingand dynamics.and we note that this overestimation will not lessen with a higher number of samples. Including the rate of false positive andnegative tests, this becomes r B ( t ) = Br I ( − ε ) + r I ( − ε ) + ε ( − r I − r I ) + ( B − ) r I > r ( t ) (14)In order to estimate the spread of the disease in a community (i.e., estimate the percentage of the population which is infected),we can estimate the bias B if we have unbiased sample data as well. To do this with only positive/negative test data, we mustassume that the ratio of symptomatic to asymptomatic infected members of the community is constant (i.e., r = cr ). Here, weare making an estimation analogous to a Monte-Carlo method, and so for better accuracy we need to reduce the variance in ourestimate of B . We simulated an estimate of B with various biased and unbiased sample capacities. Variances for these estimatesare shown in fig. 3. We see there that with 1000 unbiased and 1000 biased tests, we can estimate B with a variance (and so errorin the estimate of B ) of less than 2. Depending on the magnitude of B , this is likely reasonable error. In simulation, this wastested with a true bias of B =
10, meaning that a variance of 2 represents a 20% error. In fig. 4, we see the effect of this bias inthe overestimation of the infected proportion of the population.
Linear Trends
Dependence on dynamic slope
The ability to correctly characterize a linear trend in the dynamics from sampled data depends on the strength of that trend aswell as the nature of the population sampling. Confidence in trends is reported as the proportion of five-day intervals in 1000simulations which correctly identified the sign (positive/negative) of the linear fit to the dynamics. Sampling identifies the signof the trend robustly for large enough absolute slope of the dynamics. See appendix B for details. We also see that increasingsampling improves identification of trends in data.In fig. 5, we show how this dependence on the underlying dynamics effects confidence in trend identification over thecourse of a simulated outbreak. Here, we see that trends can be identified with good confidence with 8100 samples per day formost time-intervals. Those intervals in which trends could not be confidently identified were those that included local maxima(peaks) or minima (valleys) in the epidemic dynamics.
Confidence with bias and errors
In fig. 6, we see that biased sampling actually improves the confidence of an estimate of the sign of a trend (i.e., is the trendpositive or negative). This is because biased sampling magnifies a trend in infections which are relatively rare, meaning thetrend appears stronger in the biased data. Biased testing allows for higher trend confidence in error free (no false positive/falsenegative) testing as well as testing with 10% error rate.
Peak Finding
Tables 1 and 2 give precision and recall for peak finding with two sets of dynamics (SIR generated and SAIR generated)with various sampling assumptions (with and without bias and errors). We observe that identification of peaks in data usingthe smoothing method described above has a high chance of finding the peaks in the dynamics, but has very poor precision,providing many false peaks. See appendix C for further details. igure 3.
Variance of bias estimate for various sampling rates. This variance represents the error in estimation from a singleday of biased and unbiased tests. In this experiment, we do not have asymptomatic infected individuals. igure 4.
Simulation shows how biasing testing towards symptomatic individuals overestimates infected proportion of thepopulation, even with no false-positive or false-negative tests, and a high rate of testing. Here, the bias parameter is set to B ( t ) =
10 and the testing rate is λ ( t ) = t . igure 5. Trend sign confidence changes with the underlying dynamics. igure 6.
Confidence in linear trends identified for sampling of SIR and SAIR dynamics with and without bias and errors. rends in COVID-19 data
Overall, our model suggests that five-day trends can be used with confidence if bias was constant for testing period. Forexample, we have confidence in five-day trends of the outbreak in the state of Minnesota using data from
The COVID TrackingProject (Lipton et al., 2020) to compute, shown in fig. 7, with data from March 6 to July 14, 2020.Our approach demonstrates that epidemic sampling data is more difficult to interpret accurately when the bias in testingvaries with time. Unfortunately, such a variation is suggested by a significant negative correlation between positive testpercentage and number of tests performed in many states. This can be explained by a reduction in bias as more tests becomeavailable (i.e., an increased willingness to test asymptomatic members of the population). In fact, a strong negative correlationcould indicate that testing may have been initially used in a more restrictive, and therefore more heavily biased, manner. Wehypothesize that as testing increases, testing bias will approach some limit that represents the preferred policies of healthcareand government organizations. It may be that even with high testing capacity, some bias will still exist due to testing practicesand patient self-selection. Changes in policies will result in future changes in testing bias. We note that our model can simulatea change in bias with a time-dependent bias parameter B in eqs. (6) to (8). Our model is built with this problem in mind,allowing a time-varying total intensity function λ ( t ) . However, determining λ ( t ) remains a challenge. This may require otherthan strictly testing data, such as test production data or self-reported testing bias from healthcare providers.It is worth noting that some day-to-day variation may be the result of irregularities in negative test reporting, as evidenced bydays with 100% positive rate. Pearson correlations are shown in fig. 8, with significance computed as p-value of the correlationcoefficient. DISCUSSION
Confidence in any data analysis technique must be carefully assessed in light of the numerous confounding variables in epidemicsampling data, including bias in testing and limits to testing capacity. Our model provides realistic simulated data that canbe used to assess confidence in conclusions based on sampled data, and even to calibrate and engineer novel data analysistechniques. As a relevant test of our approach and an exploration of real data from COVID-19, we use data from the
COVIDTracking Project (Lipton et al., 2020) for the state of Minnesota to compute five-day trends for that data, shown in fig. 7. Fordata collected after mid July 2020, testing capacity was generally over 2000 samples per day, and so these estimates can beseen as somewhat reliable with bias assumed to be approximately constant.We also show the correlation between positive test percentage and number of tests performed for each state in fig. 8. We seefor example that in Minnesota, this correlation is approximately − .
25. In most states, there is a significant anti-correlationbetween between positive test percentage and number of tests performed. This may be due to changes in the policies ofhealthcare providers and government organizations as tests become available. We must therefore account for this change inbias when discussing trends in the spread of the disease. Additionally, some of this correlation may be due to the occurrenceof days on which negative tests are not reported or under-reported. We may model changes in bias simply by choosing sometime-dependent bias functions B ( t ) in eqs. (6) to (8).While the multiple purposes of infectious disease testing would be satisfied if all or almost all infected individuals could beidentified, the amount of testing needing to have a high level confidence that almost all infected individuals have been identifiedis prohibitively high. For example, the hypergeometric distribution suggest that if we have a population of 300 million with 5%infection, then we need about 27 million unbiased tests per day for 90% confidence that we have found 90% of the cases. ForCOVID-19, it remains very unlikely that case numbers reported represent an accurate estimate of the extent of disease spread.Furthermore, these numbers cannot be compared from place to place or time to time because of changes in testing bias (Liptonet al., 2020). As an example of how testing bias can affect perception of a trend, we simulate of an artificial scenario wheretesting capacity (i.e., λ ( t ) ) increases drastically part-way through the course of an infection in fig. 2, and demonstrate thatconsidering only the number of positive tests per day completely obscures a peak in the dynamics. On the other hand, simplyconsidering the proportion of tests in a day which are positive reveals the true dynamics. This emphasizes the importance ofproportion of positive tests over the number of positive tests.Testing for COVID-19 is clearly biased toward finding infected individuals. While reduced testing has drawbacks foraddressing particular scenarios, such as screening healthcare workers, it also has important benefits for tracking the populationlevel trends that inform policy decisions. As an intuitive example, consider a population with a very small proportion ofCOVID-19 cases, as would be expected in the very early or very late stages of an outbreak. Heavily biased testing helps betterdetect the infection by focusing on where the cases are rather than spending the vast majority of tests on negative results. Inthis sense, biased testing is a form of importance sampling. Furthermore, biased testing reduces the number of tests neededto identify all or most infected individuals. In fig. 1, we show the number of samples needed for bias parameters rangingfrom B = B =
50. Bias in testing is the natural result of the role of testing in the healthcare setting, and this confirms theadvantages bias has for detecting population trends. However, using bias in testing as shown in eqs. (6) to (8), we see in fig. 4that the spread of the infection will be overestimated significantly by biased testing. In other words, to accurately estimate the igure 7. five-day trend of positive test proportion in the state of Minnesota, using data from
The COVID tracking project (Lipton et al., 2020) igure 8.
Correlation between positive test percentage and tests performed in each state. Significance indicates p < . p < . p < . p < . pread of disease, we must estimate the bias parameter B . This can be done by conducting a separate set of unbiased tests andusing the relationship given by eq. (14). If we assume that the testing bias is constant (which is reasonable for a single day), thisis a Monte-Carlo estimator where the error in this estimation is determined by the variance in the estimate. We simulated with abias B =
10, and show the variance of single-day estimates for B with various biased and unbiased testing capacities in fig. 3.From that simulation, we conclude that it is not unreasonable to estimate bias with 1000 unbiased samples per day, in additionto a larger capacity of biased testing.Policy changes during an outbreak, such as the recent activation and deactivation of stay-at-home orders in the USA, appearto be based on trends in the disease dynamics, i.e., whether disease spread is accelerating or decelerating, or if there havebeen changes to the rate of acceleration. Our work shows that we can account for testing bias and successfully determine theunderlying trend in disease dynamics. Moreover, we show that the overall positive or negative trend is not overly sensitiveto the bias, meaning that assuming an approximately constant bias may work for most estimates. It is important to note thatdetermining the sign of trends in the disease is easier when the trends are larger in magnitude, as shown in appendix B. The lesschange there is in infection rate, such as those through smaller policy changes, the more testing is needed to identify an effect.As an example, we see in fig. 5 that 8100 samples per day is enough to give good confidence in the estimated sign of a five-daytrend in disease dynamics for most of the course of an outbreak. This confidence is low when the trend is very weak, meaningthe true dynamics are at a local maximum (peak) or minimum (valley). Finally, we see in fig. 6 that a constant bias in testingactually improves our ability to detect the sign of a trend in the dynamics. This is because biased testing magnifies trends in thedata, as can be seen in fig. 4. We see again that 8100 tests per day gives high confidence in the sign of five day trends in thedata as long as that data is done with a constant high bias.Policy may also be based on other metrics in sampling data, such as the occurrence of peaks or more complicated modelfitting. Our model of sampling provides a method for testing the confidence of these metrics. As an example we show that theexact peaks in an outbreak can be found, as seen in table 2, but there will likely be a large number of false peaks, as seen intable 1. However, with the right smoothing, critical points in the dynamics can be identified with some confidence.As written, our model assumes that the dynamics of an epidemic can be characterized by tracking three compartmentswithin a society which we refer to as “symptomatic”, “asymptomatic” and “healthy”. Thus, any dynamics must be recast ascounts of individuals who have a disease and show symptoms, those who have a disease and do not show symptoms, and thosewho do not have a disease. For finer grained models of epidemic sampling, these compartments can still be determined and ourmodel used, but information may be lost. It may be beneficial then to tailor a model analogous to ours to simulate sampling ona more detailed model of epidemic sampling. This can be done simply by increasing the number of compartments in the modeland calculating equations similar to eqs. (6) to (8). CONCLUSION
We provide a model of sampling in a disease outbreak in order to simulate data analysis in different outbreak situations and toassess infection testing strategies. Clearly, we should account for the confidence we have in the measurements of metrics usedto set policy. This confidence is affected by testing capacity, errors, and bias. Our model provides a method to assess confidencewith time-varying testing capacity and bias by simulating sampling over the course of an epidemic. This model demonstratesthe importance of tracking testing capacity, estimating possible changes in bias, and tracking positive test percentage ratherthan raw number of positive tests. Our model provides an essential tool in designing an effective response to the outbreak of aninfectious disease.
REFERENCES
Alvarez, F. E., Argente, D., and Lippi, F. (2020). A simple planning problem for covid-19 lockdown. Technical report, NationalBureau of Economic Research.Anderson, D. F. and Kurtz, T. G. (2011). Continuous time markov chain models for chemical reaction networks. In
Design andanalysis of Biomolecular Circuits , pages 3–42. Springer.Anderson, D. F. and Yuan, C. (2019). Low variance couplings for stochastic models of intracellular processes with time-dependent rate functions.
Bulletin of Mathematical Biology , 81(8):2902–2930.Asmussen, S. and Glynn, P. W. (2007).
Stochastic Simulation: Algorithms and Analysis , volume 57. Springer Science &Business Media.Bai, Y., Yao, L., Wei, T., Tian, F., Jin, D.-Y., Chen, L., and Wang, M. (2020). Presumed asymptomatic carrier transmission ofCOVID-19.
Jama , 323(14):1406–1407.Centers for Disease Control and Prevention (2020). Coronavirus (COVID-19). .Accessed: 2020-07-08. howdhury, R., Kevin Heng, M., Shawon, S. R., Goh, G., Okonofua, D., Ochoa-Rosales, C., Gonzalez-Jaramillo, V., Bhuiya,A., Reidpath, D., Prathapan, S., Shahzad, S., Althaus, C. L., Gonzalez-Jaramillo, N., Franco, O. H., and The Global DynamicInterventions Strategies for COVID-19 Collaborative Group (2020). Dynamic interventions to control COVID-19 pandemic:a multivariate prediction modelling study comparing 16 worldwide countries.
European Journal of Epidemiology , pages1–11.Edelstein-Keshet, L. (2005).
Mathematical Models in Biology . SIAM.Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemicalreactions.
Journal of Computational Physics , 22(4):403–434.Gillespie, D. T. (1977). Exact stochastic simulation of coupled chemical reactions.
The Journal of Physical Chemistry ,81(25):2340–2361.Hethcote, H. W. (2000). The mathematics of infectious diseases.
SIAM Review , 42(4):599–653.Kermack, W. O. and McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics.
Proceedings of theRoyal Society of London. Series A, Containing Papers of a Mathematical and Physical Character , 115(772):700–721.Klenke, A. (2014).
The Poisson Point Process , pages 543–561. Springer London, London.Lee, J. C., Mervosh, S., Avila, Y., Harvey, B., and Matthews, A. L. (2020). See how all 50 states are reopening (and closing again). . Ac-cessed: 2020-07-08.Lipton, Z., Ellington, J., and Riley, K. (2020). The COVID tracking project. https://covidtracking.com .Mizumoto, K., Kagaya, K., Zarebski, A., and Chowell, G. (2020). Estimating the asymptomatic proportion of coronavirusdisease 2019 (COVID-19) cases on board the diamond princess cruise ship, yokohama, japan, 2020.
Eurosurveillance ,25(10):2000180.Occupational Safety and Health Administration (2020). Guidance on preparing workplaces for COVID-19. . Accessed: 2020-07-08.Piguillem, F. and Shi, L. (2020). Optimal COVID-19 quarantine and testing policies.Poll´an, M., P´erez-G´omez, B., Pastor-Barriuso, R., Oteo, J., Hern´an, M. A., P´erez-Olmeda, M., Sanmart´ın, J. L., Fern´andez-Garc´ıa, A., Cruz, I., Fern´andez de Larrea, N., Molina, M., Rodr´ıguez-Cabrera, F., Mart´ın, M., Merino-Amador, P., Le´onPaniagua, J., Mu˜noz-Montalvo, J. F., Blanco, F., Yotti, R., Blanco, F., Guti´errez Fern´andez, R., Mart´ın, M., Mezcua Navarro,S., Molina, M., Mu˜noz-Montalvo, J. F., Salinero Hern´andez, M., Sanmart´ın, J. L., Cuenca-Estrella, M., Yotti, R., Le´onPaniagua, J., Fern´andez de Larrea, N., Fern´andez-Navarro, P., Pastor-Barriuso, R., P´erez-G´omez, B., Poll´an, M., Avell´on, A.,Fedele, G., Fern´andez-Garc´ıa, A., Oteo Iglesias, J., P´erez Olmeda, M. T., Cruz, I., Fernandez Martinez, M. E., Rodr´ıguez-Cabrera, F. D., Hern´an, M. A., Padrones Fern´andez, S., Rumbao Aguirre, J. M., Navarro Mar´ı, J. M., Palop Borr´as, B., P´erezJim´enez, A. B., Rodr´ıguez-Iglesias, M., Calvo Gasc´on, A. M., Lou Alcaine, M. L., Donate Su´arez, I., Su´arez ´Alvarez, O.,Rodr´ıguez P´erez, M., Cases Sanch´ıs, M., Villaf´afila Gomila, C. J., Carbo Saladrigas, L., Hurtado Fern´andez, A., Oliver, A.,Castro Feliciano, E., Gonz´alez Quintana, M. N., Barrasa Fern´andez, J. M., Hern´andez Betancor, M. A., Hern´andez Febles,M., Mart´ın Mart´ın, L., L´opez L´opez, L.-M., Ugarte Miota, T., De Benito Poblaci´on, I., Celada P´erez, M. S., Vall´es Fern´andez,M. N., Mat´e Enr´ıquez, T., Villa Arranz, M., Dom´ınguez-Gil Gonz´alez, M., Fern´andez-Natal, I., Meg´ıas Lob´on, G., Mu˜nozBellido, J. L., Ciruela, P., Mas i Casals, A., Dolad´e Bot´ıas, M., Marcos Maeso, M. A., P´erez del Campo, D., F´elix de Castro,A., Lim´on Ram´ırez, R., El´ıas Retamosa, M. F., Rubio Gonz´alez, M., Blanco Lobeiras, M. S., Fuentes Losada, A., Aguilera,A., Bou, G., Caro, Y., Marauri, N., Soria Blanco, L. M., del Cura Gonz´alez, I., Hern´andez Pascual, M., Alonso Fern´andez, R.,Merino-Amador, P., Cabrera Castro, N., Tom´as Lizcano, A., Ram´ırez Almagro, C., Segovia Hern´andez, M., Ascunce Elizaga,N., Ederra Sanz, M., Ezpeleta Baquedano, C., Bustinduy Bascaran, A., Iglesias Tamayo, S., Elorduy Otazua, L., BenarrochBenarroch, R., Lopera Flores, J., and V´azquez de la Villa, A. (2020). Prevalence of SARS-CoV-2 in spain (ENE-COVID): anationwide, population-based seroepidemiological study.
The Lancet .Reich, N., Wattanachit, N., Ray, E., Niemi, J., Le, K., House, K., Cornell, M., Brennen, A., and Bracher, J. (2020). COVID-19forecast hub. https://github.com/reichlab/covid19-forecast-hub .World Health Organization (2020). Global surveillance for COVID-19 caused by human infection with COVID-19 virus. file:///Users/m197894/Documents/COVID_Model/WHO-2019-nCoV-SurveillanceGuidance-2020.6-eng.pdf . A DETAILED DESCRIPTION OF THE BIAS PARAMETER
Below, we include a calculation to more provide better intuition about the nature of the bias parameter B . To explain thisparameter, we consider the rate at which compartments of the population are tested for a disease. In an infinitesimal time-interval [ t , t + h ) , there is some probability p h that an asymptomatic or healthy individual will be tested, and some probability p h thata symptomatic infected individual will be tested. The nature of Poisson point processes is that (assuming for simplicity no rrors in the tests) λ + ( t ) = p I ( t ) + p I ( t ) (15)and λ − ( t ) = p H ( t ) (16)and that the total rate of testing is λ ( t ) = p ( I ( t ) + H ( t )) + p I ( t ) . We then have from eq. (9) that p I ( t ) + p I ( t ) = ( p ( I ( t ) + H ( t )) + p I ( t )) I ( t ) + BI ( t ) I ( t ) + BI ( t ) + H ( t ) . (17)We can rewrite this as p I ( t ) + p I ( t )( p ( I ( t ) + H ( t )) + p I ( t )) = I ( t ) + BI ( t ) I ( t ) + BI ( t ) + H ( t ) (18)and see that B = p p . (19)Thus the bias B can be interpreted as the increased rate of testing of symptomatic individuals over asymptomatic individuals, aspresumably p > p . The effect of this difference is that the apparent population sampled is an adjusted version of the truepopulation, with apparent total I ( t ) + BI ( t ) + H ( t ) and apparent number of infected I ( t ) + BI ( t ) . B SAMPLING FROM DIFFERENT UNDERLYING DYNAMICS
In fig. A1, we test sampling’s ability to identify a constant trend (i.e., linear increase or decrease) in the infected proportion. Wesee that sampling identifies the sign of the trend robustly for large enough absolute slope.
C RESULTS OF PEAK FINDING
Tables 1 and 2 were generated with a smoothing parameter of 5. Here, we demonstrate that this can be improved with anoptimal choice of smoothing. Figure A2 uses smoothing from 1 to 10, with 10 giving the highest precision and recall.
D MODELS USED TO GENERATE DYNAMICS
SIR model
We generate SIR dynamics by considering three pools of individuals: those that are susceptible, those that are infected, andthose that have recovered. Individuals transition between these pools according to mass action dynamics given in the ODEmodel: dx S dt = − β x I x S dx I dt = β x I x S − γ x i dx R dt = γ x I where x S , x I , x R represent the proportion of the population that is susceptible, infected, or recovered from the disease. We allow β = β ( t ) to be a function of time, and choose a dynamic parameterization which allows us to generate an infection with morethan one peak time, as can be seen in fig. 4. We do this simply by varying β (intuitively varying the virulence of the disease) sothat it decreased until t = : β ( t ) = (cid:32) − exp (cid:32) − (cid:18) t − (cid:19) (cid:33)(cid:33) (20)and γ = . This dynamic parameterization allows us to generate an infection with more than one peak time, as can be seen infig. 4. igure A1. Trend fitting for five-day interval of dynamics with constant slope. Top: Confidence in the sign of the estimatedslope as actual slope varies. Bottom: Error in estimated slope as slope varies. igure A2.
Precision and recall for peak finding with various smoothing parameters. These simulations sought peaks in acourse of infection dynamics generated by the SIR model with a variable β ( t ) . The underlying dynamics were the same asshown in fig. 4. Sampling was done with bias of 10 and false positive/false negative rates of 10%. his model can be interpreted as stating that individuals transition from susceptible to infected at the rate β x I x S , andtransition from infected to recovered at the rate γ x i . The well known “ R ( t ) ” parameter is defined as R ( t ) = N βγ (21)where N is the total population size (Edelstein-Keshet, 2005).Using these dynamics, we take I ( t ) = x I ( t ) , I ( t ) =
0, and H ( t ) = x S ( t ) + x R ( t ) . SAIR model
We also test dynamics that include an asymptomatic infected population. We generate SAIR dynamics by considering threepools of individuals: those that are susceptible, those that are infected, and those that have recovered. Individuals transitionbetween these pools according to mass action dynamics given in the ODE model: dx S dt = − ( β + β ) x I x S − ( β + β ) x I x S dx I dt = β x I x S + β x I x S − γ x I − δ x I dx I dt = β x I x S + β x I x S − γ x I + δ x I dx R dt = γ ( x I + x I ) where x S , x I , x I , x R represent the proportion of the population that is susceptible, asymptomatic infected, symptomatic infected,or recovered from the disease.This model can be interpreted as stating that individuals transition from susceptible to asymptomatic infected at the rate β x I x S + β x I x S , from susceptible to symptomatic at the rate β x I x S + β x I x S , from asymptomatic to symptomatic at therate δ x I , and recover at the rate γ x I if asymptomatic and γ x I if symptomatic. Note that if we take β = β = β = β wemay again define the intrinsic reproduction rate R ( t ) as in the SIR model.Using these dynamics, we have I ( t ) = x I ( t ) , I ( t ) = x I ( t ) , and H ( t ) = x S ( t ) + x R ( t ) . E ACCOUNTING FOR SMALL POPULATIONS.