Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
70,716 characters · 21 sections · 37 citation commands
Estimation of Covid-19 Prevalence from Serology Tests: A Partial Identification Approach
{\em Keywords}: partial identification; disease prevalence; serology tests; Covid-19.
{\em JEL classification codes}: C12, C14, I10.
\thispagestyle{empty}
\setcounter{page}{1}
Since December 2019 the world has been facing the Covid-19 pandemic, and its disastrous effects in human life and the economy. Responding to the pandemic, most countries have closed off their borders, and imposed unprecedented, universal lockdowns on their entire economies. The key reason for such drastic measures is uncertainty: we do not yet know the actual transmission rate, the lethality, or the prevalence of this new deadly disease. As governments and policy makers were caught by surprise, there is no doubt that these drastic measures were needed as a first line of defense. The data show that we would have to deal with a massive humanitarian disaster otherwise.
At the same time, as the economic pain mounts, especially for the most vulnerable and disadvantaged segments of the population, there is an urgent need to think of careful ways to safely reopen the economy. Estimating the true prevalence of Covid-19 has been identified as a key parameter to this effort alvarez2020simple. In the United States, the number of confirmed Covid-19 cases is 1,193,813 as of May 7 with 70,802 total deaths. This implies a (case) prevalence of 0.36% (assuming 328m as the US population), and a 5.9% mortality rate of Covid-19, which is even higher than the mortality rate reported at times by the World Health Organization.\footnote{The official mortality rate was revised from 2% in late January to 4% in early March; see also an official WHO situation report from early March: \url{https://www.who.int/docs/default-source/coronaviruse/situation-reports/20200306-sitrep-46-covid-19.pdf?sfvrsn=96b04adf_2}.} However, the true prevalence, that is, the number of people who are currently infected or have been infected by Covid-19 over the entire population is likely much higher, and so the mortality rate should be significantly lower than 5.9%. A growing literature is attempting to estimate these numbers through epidemiological models li2020substantial, flaxman2020report, or structural assumptions hortaccsu2020estimating.
A more robust alternative seems to be possible through randomized serological studies that detect marker antibodies indicating exposure to Covid-19. In the US, there is currently a massive coordinated effort to evaluate the widespread application of these tests. The results are expected in late May of 2020.\footnote{CDC page: \url{https://www.cdc.gov/coronavirus/2019-ncov/lab/serology-testing.html}.} The hope is that these tests will determine the true prevalence of the virus, and thus its lethality, and also determine whether someone is immune enough to return to work (the extent of immunity is still uncertain, however). Furthermore, seroprevalence studies can give information on risk factors for the disease, such as a patient's age, location, or underlying health conditions. They may also reveal important medical information on immune responses to the virus, such as how long antibodies last in people’s bodies following infection, and could also identify those able to donate blood plasma, which is a possible treatment to seriously ill Covid-19 patients.\footnote{Food and Drug Administration (FDA) announcement on serology studies (04/07/2020): \url{https://www.fda.gov/news-events/press-announcements/coronavirus-covid-19-update-serological-tests}.} The development of serology tests is therefore essential to designing a careful strategy towards both effective medical treatments and a gradual reopening of the economy.
Until widespread serology testing is possible, however, we have to rely on a limited number of serology studies that have started to emerge in various areas of the globe, including the US. Table (ref) presents a non-exhaustive summary of such studies around the world. For example, in Germany, serology tests in early April showed a 14% prevalence in a sample of 500 people.\footnote{Report in German: \url{https://www.land.nrw/sites/default/files/asset/document/zwischenergebnis_covid19_case_study_gangelt_0.pdf}.} In the Netherlands, a study in mid-April showed a lower prevalence at 3.5% in a small sample of blood donors.\footnote{Presentation slides in Dutch: \url{https://www.tweedekamer.nl/sites/default/files/atoms/files/tb_jaap_van_dissel_1604_1.pdf}. } \footnote{ See also a summary of these projects in the journal “Science": \url{https://www.sciencemag.org/news/2020/04/antibody-surveys-suggesting-vast-undercount-coronavirus-infections-may-be-unreliable}. } In the US, in a recent and relatively large study in Santa Clara, California, bendavid2020covid estimated an in-sample prevalence of 1.5% from 50 positive test results in a sample of 3330 patients. Using a reweighing technique, the authors extrapolated this estimate to 2-4% prevalence in the general population. A follow-up study in LA County found 35 positives out of 846 tests. What is unique about these last two studies is that data from a prior validation study are also available, where, say, 401 “true negatives" were tested with 2 positive results, implying a false positive rate of 0.5%. Upon publication, these studies received intense criticism because the false positive rate appears to be large enough compared to the underlying disease prevalence. For example, the Agresti-Coull and Clopper-Pearson 95% confidence intervals for the false positive rate are $[0.014\%, 1.92\%]$ and $[0.06\%, 1.79\%]$, respectively. These intervals for the false positive rate are not incompatible even with a 0% prevalence, since a 1.5% false positive rate achieves $0.015 \times 3330 \approx 50$ (false) positives on average, same as the observed value in the sample.
Such standard methods, however, are justified based on approximations, asymptotic arguments, prior specifications (for Bayesian methods), or normality assumptions, which are always suspect in small samples. In this paper, we develop a method that can assess finite-sample statistical significance in a robust way. The key idea is to treat all unknown quantities as parameters, and explore the entire parameter space to assess agreement with the observed data. Our method adopts the partial identification framework, where the goal is not to produce point estimates, but to identify sets of plausible parameter values wooldridge2007s, tamer2010partial, chernozhukov2007estimation, manski2003partial, manski2010partial,manski2007partial, romano2008inference, romano2010inference, honore2006bounds, imbens2004confidence, beresteanu2012partial, stoye2009more, kaido2019confidence. Within that literature, our proposed method appears to be unique in the sense that it constructs a procedure that is valid in finite samples given the correct distribution of the test statistic. Importantly, the choice of the test statistic can affect only the power of our method, but not its validity. Such flexibility may be especially valuable in choosing a test statistic that is both powerful and easy to compute. Thus, the main benefit of our approach is that it is valid with {\em enough computation}, whereas more standard methods are only valid with {\em enough samples}.
The rest of this paper is structured as follows. In Section (ref) we describe the problem formally. In Section (ref) we describe the proposed method on a high level. A more detailed analysis along with a modicum of theory is given in Section (ref). In Section (ref) we apply the proposed method on data from the Santa Clara study, the LA County study, and a recent study from New York state.
Here, we formalize the statistical problem of estimating disease prevalence through imperfect medical tests. Every individual $i$ is associated with a binary status $x_i$: it is $x_i= 1$ if the individual has developed antibodies from exposure to the disease, and $x_i=0$ if not. We will also refer to these cases as “positive" and “negative", respectively. Patient status is not directly observed, but can be estimated with a serology (antibody) test.
This medical antibody test can be represented by a function $t: \{0, 1\} \to \{0, 1\}$, and determines whether someone is positive or negative. As usual, the categorization of the test results can be described through the following table:
We will assume that each test result is an independent random outcome, such that the true positive rate and false positive rate, denoted respectively by $q$ and $p$,\footnote{The terms “sensitivity" and “specificity" are frequently used in practice of medical testing. In our setting, sensitivity maps to the true positive rate $(q)$, and specificity maps to one minus the false positive rate $(1-p)$. In this paper, we will only use the terms “true/false positive rate" as they are more precise and self-explanatory. } are constant:
This assumption may be untenable in practice. In general, patient characteristics, or test target and delivery conditions can affect the test results. For example, bendavid2020covid report slightly different test performance characteristics depending on which antibody (either IgM or IgG) was being detected. We note, however, that this assumption is not strictly necessary for the validity of our proposed inference procedure. It is only useful in order to obtain a precise calculation for the distribution of the test statistic (see Theorem (ref) and remarks).
To determine test performance characteristics, and gain information about the true/false positive rates of the antibody test, there is usually a {\em validation study} where the underlying status of participating individuals is known. In the Covid-19 case, for example, such validation study could include pre-Covid-19 blood samples that have been preserved, and are thus “true negatives". To simplify, we assume that in the validation study there is a set $\mathcal{I}_\mathrm{c}^-$ of participating individuals, where it is known that everyone is a true negative, and a set $\mathcal{I}_\mathrm{c}^+$ where everyone is positive; i.e., $$ x_i=0,~\text{for all}~i\in \mathcal{I}_\mathrm{c}^-,\text{and}~x_i=1,~\text{for all}~i\in \mathcal{I}_\mathrm{c}^+. $$ There is also the {\em main study} with a set $\mathcal{I}_\mathrm{m}$ of participating individuals, where the true status is not known. We assume no overlap between sets $\mathcal{I}_\mathrm{c}^-, \mathcal{I}_\mathrm{c}^+$ and $\mathcal{I}_\mathrm{m}$, which is a realistic assumption. We define $N_\mathrm{c}^- = |\mathcal{I}_\mathrm{c}^-|$ and $N_\mathrm{c}^+=|\mathcal{I}_\mathrm{c}^+|$ as the respective number of participants in the validation study, and $N_\mathrm{m} = |\mathcal{I}_\mathrm{m}|$ as the number of participants in the main study. These numbers are observed, but the full patient sets or the patient characteristics, may not be observed.
We also observe the positive test results in both studies:
Thus, $S_\mathrm{c}^-$ is the number of false positives in the validation study since we know that all individuals in $\mathcal{I}_\mathrm{c}^-$ are true negatives. Similarly, $S_\mathrm{c}^+$ is the number of true positives in the validation study since all individuals in $\mathcal{I}_\mathrm{c}^+$ are known to be positive. These numbers offer some simple estimates of the false positive rate and true positive rate of the medical test, respectively: $\hat p = S_\mathrm{c}^- / N_\mathrm{c}^-$ and $\hat q = S_\mathrm{c}^+ / N_\mathrm{c}^+$. We use $(s_{\mathrm{c}, \mathrm{obs}}^-, s_{\mathrm{c}, \mathrm{obs}}^+, s_{\mathrm{m}, \mathrm{obs}})$ to denote the observed values of test positives $(S_\mathrm{c}^-, S_\mathrm{c}^+, S_\mathrm{m})$, respectively, which are integer-valued random variables.
The statistical task is therefore to use observed data $\{(N_\mathrm{c}^-, N_\mathrm{c}^+, N_\mathrm{m}), (s_{\mathrm{c}, \mathrm{obs}}^-, s_{\mathrm{c}, \mathrm{obs}}^+, s_{\mathrm{m}, \mathrm{obs}})\}$ and do inference on the quantity:
i.e., the unknown disease prevalence in the main study. We emphasize that $\pi$ is a finite-population estimand --- we discuss (briefly) the issue of extrapolation to the general population in Section (ref). The challenge here is that $S_\mathrm{m}$ generally includes both false positives and true positives, which depends on the unknown test parameters, namely the true/false positive rates $q$ and $p$. Since $\pi N_\mathrm{m}$ is the (unknown) number of infected individuals in the main study, we can use Assumptions (A1) and (A2) to write down this decomposition formally:
where $\mathrm{Binom}$ denotes the binomial random variable. For brevity, we define $\boldsymbol{S} = (S_\mathrm{c}^-, S_\mathrm{c}^+, S_\mathrm{m})$ as our joint data statistic, and $\boldsymbol{\theta} = (p, q, \pi)$ as the joint parameter value. The independence of tests implies that the density of $\boldsymbol{S}$ can be computed exactly as follows.
where $\mathrm{d}(k; n, s)$ denotes the probability of $k$ successes in a binomial experiment with $n$ trials and $s$ probability of success. There are several ways to implement Equation (ref) efficiently --- we defer discussions on computational issues to Section (ref).
We begin with an illustrative example to describe the proposed method on a high level. We give more details along with some theoretical guarantees in the section that follows.
Let us consider the Santa Clara study bendavid2020covid with observed data: $$ (N_\mathrm{c}^-, N_\mathrm{c}^+, N_\mathrm{m}) =(401, 197, 3330),~\text{and}~ (s_{\mathrm{c}, \mathrm{obs}}^-, s_{\mathrm{c}, \mathrm{obs}}^+, s_{\mathrm{m}, \mathrm{obs}}) = (2, 178, 50). $$ The unknown quantities in our analysis are $q,p$ and $\pi$: the true positive rate of the test, the false positive rate, and the unknown prevalence in the main study, respectively. Assume zero prevalence ($\pi=0\%$), 90% true positive rate ($q=0.9$), and 1.5% false positive rate ($p=0.015$). We ask the question: “Is the combination $(p, q, \pi) = (0.015, 0.90, 0)$ compatible with the data?". Naturally, this can be framed in statistical terms as a null hypothesis:
To test $H_0$ we have to compare the observed positive test results with the values that {\em could have been observed} if indeed the true model parameter values were $(p, q, \pi) = (0.015, 0.90, 0)$. Our model is simple enough that we can execute this hypothetical analysis exactly based on the density of $f(\boldsymbol{S} | \boldsymbol{\theta})$ in Equation (ref), where $\boldsymbol{S} = (S_\mathrm{c}^-, S_\mathrm{c}^+, S_\mathrm{m})$ is the vector of all positive test results, and $\boldsymbol{\theta}$ is specified as in $H_0$; see Figure (ref).
To simplify visualization, in Figure (ref) we fix the component $S_\mathrm{c}^+$ of $\boldsymbol{S}$ to its observed value $(S_\mathrm{c}^+= 178$), and only plot the density with respect to the other two components, $(S_\mathrm{c}^-, S_\mathrm{m})$; i.e., we plot $f(\boldsymbol{S} \mid H_0, S_\mathrm{c}^+=178)$. One can visualize the full joint distribution $f(\boldsymbol{S} \mid H_0)$ as a collection of such conditional densities for all possible values of $S_\mathrm{c}^+$.
The next step is to decide whether the observed value of $\boldsymbol{S}$, namely $\boldsymbol{s}_\mathrm{obs} = (2,178, 50)$, is compatible with the distribution of Figure (ref). We see that the mode of the distribution is around the point $(S_\mathrm{c}^-, S_\mathrm{m}) = (5, 45)$, whereas the point $(2, 50)$ is at the lower edge of the distribution. If the observed values were even further, say at $(S_\mathrm{c}^-, S_\mathrm{m}) = (2, 80)$, then we could confidently reject $H_0$ since the density at $(2, 80)$ is basically zero. Here, we have to be careful because the actual observed values are still somewhat likely under $H_0$. Our method essentially accepts $H_0$ when the density of this distribution at the observed value $\boldsymbol{s}_\mathrm{obs}$ of statistic $\boldsymbol{S}$ is above some threshold $c_0$, that is, we decide based on the following rule:
The test in Equation (ref) is reminiscent of the likelihood ratio test, the key difference being that our test does not require maximizations of the likelihood function over the parameter space, which is computationally intensive, and frequently unstable numerically --- we make a concrete comparison in the application of Section (ref). Our test essentially uses the density of $\boldsymbol{S}$ as the test statistic for $H_0$, while threshold $c_0$ generally depends on the particular null values being tested. Assuming that the test of Equation (ref) has been defined, we can then test for all possible combinations of our parameter values, $\boldsymbol{\theta} \in \Theta$, in some large enough parameter space $\Theta$, and then invert this procedure in order to construct the confidence set. As usual, we would like this confidence set to cover the true parameters with some minimum probability (e.g., 95%). In the following section, we show that this is possible through an appropriate construction of the test in Equation (ref), which takes into account the level sets of the density function depicted in Figure (ref). The overall procedure is computationally intensive, but is valid in finite samples without the need of asymptotic or normality assumptions. The details of this construction, including the appropriate selection of the test threshold and the proof of validity, are presented in the following section.
Let $\boldsymbol{S} = (S_\mathrm{c}^-, S_\mathrm{c}^+, S_\mathrm{m}) \in \mathbb{S}$ denote the statistic, where $\mathbb{S}=\{0, \ldots, N_\mathrm{c}^-\} \times \{0, \ldots, N_\mathrm{c}^+\} \times \{0, \ldots, N_\mathrm{m}\}$, and let $\boldsymbol{\theta} = (p, q, \pi) \in\Theta$ be the model parameters. We take $\Theta$ to be finite and discrete; e.g., for probabilities we take a grid of values between 0 and 1. Let $\boldsymbol{s}_\mathrm{obs} = (s_{\mathrm{c}, \mathrm{obs}}^-, s_{\mathrm{c}, \mathrm{obs}}^+, s_{\mathrm{m}, \mathrm{obs}})$ denote the observed value of $\boldsymbol{S}$ in the sample. Let $f(\boldsymbol{S} | \boldsymbol{\theta})$ denote the density of the joint statistic conditional on the model parameter value $\boldsymbol{\theta}$, as defined in Equation (ref). Suppose that $\boldsymbol{\theta}_0$ is the true unknown parameter value, and assume that
Assumption (A4) basically posits that our discretization is fine enough to include the true parameter value with probability one. In our application, this assumption is rather mild as we are dealing with parameters that are either probabilities or integers, and so bounded within well-defined ranges. Moreover, this assumption is implicit essentially in all empirical work since computers operate with finite precision. Our goal is to construct a confidence set $\widehat\Theta_{1-\alpha} \subseteq\Theta$ such that $ P(\boldsymbol{\theta}_0 \in\widehat\Theta_{1-\alpha}) \ge 1-\alpha, $ where $\alpha$ is some desired level (e.g., $\alpha=0.05$). Trivially, $\widehat\Theta_{1-\alpha}=\Theta$ satisfies this criterion, so we will aim to make $\widehat\Theta$ as narrow as possible. We will also need the following definition:
Function $\nu$ depends on level sets of $f$, and counts the number of sample data points (over the sample space $\mathbb{S}$) with likelihood at $\boldsymbol{\theta}$ that is smaller than the observed likelihood at $\boldsymbol{\theta}$.
We can now prove the following theorem.
When $z^\ast_{\boldsymbol{\theta}}$ is not a discontinuity point of $g$, for all $\boldsymbol{\theta}$, then our test is exact in the sense that $P(\boldsymbol{\theta}_0\in\widehat\Theta_{1-\alpha}) = 1-\alpha$. In general, however, this condition will not hold for all $\Theta$, and so the confidence set of Equation (ref) may be conservative and lose power. We could potentially achieve more power if instead we define the confidence set as follows:
It is straightforward to see that $\widehat\Theta_{1-\alpha}^{\mathrm{ alt}} \subseteq \widehat\Theta_{1-\alpha}$ almost surely since $$ \sum_{\boldsymbol{s}\in\mathbb{S}} \mathbb{I}\big\{f(\boldsymbol{s} | \boldsymbol{\theta}) \le f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}) \big\} f(\boldsymbol{s} | \boldsymbol{\theta}) \le f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}) \sum_{\boldsymbol{s}\in\mathbb{S}} \mathbb{I}\big\{f(\boldsymbol{s} | \boldsymbol{\theta}) \le f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}) \big\} = f(\boldsymbol{s}_\mathrm{obs}|\boldsymbol{\theta}) \nu(f(\boldsymbol{s}_\mathrm{obs}|\boldsymbol{\theta}), \boldsymbol{\theta}). $$
Since both constructions are valid in finite samples, the choice between $\widehat\Theta_{1-\alpha}$ or $\widehat\Theta_{1-\alpha}^{\mathrm{ alt}}$ should be mainly based on computational feasibility. The construction of $\widehat\Theta_{1-\alpha}$ may be easier to compute in practice as it depends on a summary of the distribution $f(\boldsymbol{s} | \boldsymbol{\theta})$ through the level set function $\nu$, while the construction of $\widehat\Theta_{1-\alpha}^{\mathrm{ alt}}$ requires full knowledge of the entire distribution. If it is computationally feasible, however, $\widehat\Theta_{1-\alpha}^{\mathrm{ alt}}$ should be preferred because it is contained in $\widehat\Theta_{1-\alpha}$ with probability one, as argued above. This leads to sharper inference. See also the applications on serology studies in Section (ref) for more details, where the construction of $\widehat\Theta_{1-\alpha}^{\mathrm{ alt}}$ is feasible.
Theorems (ref) and (ref) imply the following simple procedure to construct a 95% confidence set:
How does our method compare to a more standard frequentist or Bayesian approach? Here, we discuss two key differences. First, as we have repeatedly emphasized in this paper, our method is valid in finite samples under only independence of test results, which is a mild assumption. In contrast, a standard frequentist approach, say based on the bootstrap, is inherently approximate and relies on asymptotics, while a Bayesian method requires the specification of priors and posterior sampling. Of course, our procedure requires more computation, mainly compared to the bootstrap, and can be conservative for marginal inference (see Remark (ref)), but this is arguably a small price to pay in a critical application such as the estimation of Covid-19 prevalence.
A second, more subtle, difference is the way our method performs inference. Specifically, we decide whether any $\boldsymbol{\theta}\in\Theta$ is in the confidence set based on the entire density $f(\boldsymbol{s}|\boldsymbol{\theta})$ over all $\boldsymbol{s}\in\mathbb{S}$, whereas both frequentist and Bayesian methods typically perform inference “around the mode" of the likelihood function $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta})$ with fixed $\boldsymbol{s}_\mathrm{obs}\in\mathbb{S}$ (we ignore how the prior specification affects Bayesian inference to simplify exposition). This can explain, on an intuitive level, how the inferences of the respective methods may differ. Figure (ref) illustrates the difference. On the left panel, we plot the likelihood, $f(\boldsymbol{s}_\mathrm{obs}|\boldsymbol{\theta})$, as a function of $\boldsymbol{\theta}\in\Theta$. Typically, in frequentist or Bayesian methods, the confidence set is around the mode, say $\hat\boldsymbol{\theta}$. We see that a parameter value, say $\boldsymbol{\theta}_1$, with a likelihood value, $f(\boldsymbol{s}_\mathrm{obs} |\boldsymbol{\theta}_1)$, that is low in absolute terms will generally not be included in the confidence set. However, in our approach, the value $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}_1)$ is not important in absolute terms for doing inference, but is only important relative to all other values $\{f(\boldsymbol{s} | \boldsymbol{\theta}_1) : \boldsymbol{s}\in\mathbb{S}\}$ of the test statistic distribution $f(\boldsymbol{s}|\boldsymbol{\theta}_1)$. Such inference will typically include the mode, $\hat\boldsymbol{\theta}$, but will also include parameter values at the tails of the likelihood function, such as $\boldsymbol{\theta}_1$. As such, our method is expected to give more accurate inference in small-sample problems, or in settings with poor identifiability where the likelihood is non-smooth and multimodal. We argue that we actually see these effects in the application on Covid-19 serology studies analyzed in the following section --- see also Section (ref) and Appendix (ref) for concrete numerical examples.
In this section, we apply the inference procedure of Section (ref) to several serology test datasets in the US. Moreover, we present results for combinations of these datasets, assuming that the tests have identical specifications. This is likely an untenable assumption, but it helps to illustrate how we can use our approach to flexibly combine all evidence. Before we present the analysis, we first discuss some data on serology test performance to inform our inference.
An important aspect of serology studies are the test performance characteristics. As of May 2020, there are perhaps more than a hundred commercial serology tests in the US, but they can differ substantially across manufacturers and technologies. In our application, we use data from bendavid2020covid, who applied a serology testing kit distributed by Premier Biotech. bendavid2020covid used validation test results provided by the test manufacturer, and also performed a local validation study in the lab. The combined validation study estimated a true positive rate of 80.3% (95% CI: 72.1%-87%), and a false positive rate of 0.5% (95% CI: 0.1%-1.7%).
To get an idea about how these performance characteristics relate to other available serology tests we use a dataset published by the FDA based on benchmarking 12 other testing kits to grant emergency use authorization (EUA) status. The dataset is summarized in Table (ref). We see that the characteristics of the testing kit used by bendavid2020covid are compatible with the FDA data shown in the table. For example, a true positive rate of 80% is below the mean and median of the point estimates in the FDA dataset. A false positive rate of 0.5% falls between the median and mean of the respective FDA point estimates. A reason for this skewness is likely the existence of one outlier testing kit that performs notably worse than the others (Chembio Diagnostic Systems). Removing this datapoint brings the mean false positive rate down to 0.6%, very close to the estimate provided by bendavid2020covid.
In the Santa Clara study, bendavid2020covid report a validation study and main study, with $(N_\mathrm{c}^-, N_\mathrm{c}^+, N_\mathrm{m}) = (401, 197, 3330)$ participants, respectively. The observed test positives are $\boldsymbol{s}_\mathrm{obs} = (s_{\mathrm{c}, \mathrm{obs}}^-, s_{\mathrm{c}, \mathrm{obs}}^+, s_{\mathrm{m}, \mathrm{obs}}) = (2, 178, 50)$, respectively. Given these data, we produce the 95% confidence sets for $(p, q, \pi)$ following both procedures in (ref) and (ref) described in Section (ref). In Figure (ref) of Appendix (ref), we jointly plot all triples in the 3-dimensional space $\widehat\Theta_{0.95}$ of Equation (ref), with additional coloring based on prevalence values. We see that the confidence set is a convex space tilting to higher prevalence values as the false positive rate of the test decreases. The true positive rate does not affect prevalence, as long as it stays in the range 80%-95%.
To better visualize the pairwise relationships between the model parameters, we also provide Figure (ref) that breaks down Figure (ref) into two subplots, one visualizing the pairs $(\pi, p)$ and another visualizing the pairs $(\pi, q)$. The figure visualizes both $\widehat\Theta_{0.95}$ and $\widehat\Theta_{0.95}^{\mathrm{ alt}}$ to illustrate the differences between the two constructions. From Figure (ref), we see that the Santa Clara study is not conclusive about Covid-19 prevalence. A prevalence of 0% is plausible, given a high enough false positive rate. However, if the true false positive rate is near its empirical value of 0.5%, as estimated by bendavid2020covid, then the identified prevalence rate is estimated in the range 0.4%-1.8% in $\widehat\Theta_{0.95}$. Under this assumption, we see that $\widehat\Theta_{0.95}^{\mathrm{ alt}}$ offers a sharper inference, as expected, with an estimated prevalence in the range 0.7%-1.5%. Even though, strictly speaking, the statistical evidence is not sufficient here for definite inference on prevalence, we tend to favor the latter interval because (i) common sense precludes 0% prevalence in the Santa Clara county (total pop. of about 2 million); (ii) the interval generally agrees with the test performance data presented earlier, and (iii) it is still in the low end compared to prevalence estimates from other serology studies (see Table (ref)). Regardless, pinning down the false positive rate is important for estimating prevalence, especially when prevalence is as low as it appears to be in the Santa Clara study. Roughly speaking, a decrease of 1% in the false positive rate implies an increase of 1.3% in prevalence.
In this section, we aim to discuss how our method practically compares to more standard methods using data from the Santa Clara study. In our comparison we include Bayesian methods, a classical likelihood ratio-based test, and the Monte Carlo-based approach to partial identification proposed by chen2018monte.\footnote{ All these methods are fully implemented in the accompanying code at \url{https://github.com/ptoulis/covid-19}. }
Due to initial criticism, the authors of the original Santa Clara study published a revision of their work, where they use a bootstrap procedure to calculate confidence intervals for prevalence in the range 0.7%-1.8%.\footnote{Link: \url{https://www.medrxiv.org/content/10.1101/2020.04.14.20062463v2.full.pdf}.} Some recent Bayesian analyses report wider prevalence intervals in the range 0.3%-2.1% gelman2020bayesian. In another Bayesian multi-level analysis, levesque2020note report similar findings but mention that posterior summarization here may be subtle, since the posterior density of prevalence in their specification includes 0%. These results are in agreement with our analysis in the previous section only if we assume that the true false positive rate of the serology test was near its empirical estimate ($\sim$0.5%). We discussed intuitively the reason for such discrepancy in Section (ref), where we argued that standard methods typically do inference “around the mode" of the likelihood, and may thus miscalculate the amount of statistical information hidden in the tails.
For a numerical illustration, consider two parameter values, namely $\boldsymbol{\theta}_1 = (0.5\%, 90\%, 1.2\%)$ and $\boldsymbol{\theta}_2 = (1.5\%, 80\%, 0\%)$, where the components denote the false positive rate, true positive rate, and prevalence, respectively. In the Santa Clara study, $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}_1) = 2.2 \times 10^{-3}$ and $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}_2) = 9.58 \times 10^{-8}$, that is, $\boldsymbol{\theta}_2$ (which implies 0% prevalence) maps to a likelihood value that is many orders of magnitude smaller than $\boldsymbol{\theta}_1$. In fact, $\boldsymbol{\theta}_1$ is close to the mode of the likelihood, and so frequentist or Bayesian inference is mostly based around that mode, ignoring the tails of the likelihood function, such as $\boldsymbol{\theta}_2$. For our method, however, the small value of $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}_2)$ is more-or-less irrelevant --- what matters is how this value compares to the entire distribution $f(\boldsymbol{s} | \boldsymbol{\theta}_2)$. It turns out that $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}_2) \nu(f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}_2), \boldsymbol{\theta}_2) = 0.137$, that is, 13.7% of the mass of $f(\boldsymbol{s}|\boldsymbol{\theta}_2)$ is below the observed level $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}_2) = 9.58 \times 10^{-8}$. As such, $\boldsymbol{\theta}_2$ cannot be rejected at the 5% level (see also Appendix (ref)). This highlights the key difference of our procedure compared to frequentist or Bayesian procedures. More generally, we expect to see such important differences between the inference from our method and the inference from other more standard methods in settings with small samples or poor identification (e.g., non-separable, multimodal likelihood).
As briefly described in Section (ref), our test is related to the likelihood ratio test lehmann2006testing. Here, we study the similarities and differences between the two tests, both theoretically and empirically through the Santa Clara study. Specifically, consider testing a null hypothesis that the true parameter is equal to some value $\boldsymbol{\theta}$ using the likelihood ratio statistic,
Since $f$ is known analytically from Equation (ref), the null distribution of $t(\boldsymbol{S} | \boldsymbol{\theta})$ can be fully simulated. An exact p-value can then be obtained by comparing this null distribution with the observed value $t(\boldsymbol{s}_\mathrm{obs}|\boldsymbol{\theta})$. We can see that this method is similar to ours in the sense that both methods use the full density $f(\boldsymbol{S}|\boldsymbol{\theta})$ in the test, and both are exact. The main difference, however, is that our method is using a summary of the density values $f(\boldsymbol{S} | \boldsymbol{\theta})$ that are below the observed value $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta})$, which avoids the expensive (and sometimes numerically unstable) maximization in the denominator of the likelihood ratio test in (ref). Our proposed method turns out to be orders of magnitude faster than the likelihood ratio approach as we get $50$-fold to $200$-fold speedups in our setup --- see Section (ref) for a more detailed comparison in computational efficiency.
To efficiently compare the inference between the two tests, we sampled 5,000 different parameter values from inside $\widehat\Theta_{0.95}$ --- i.e., the 95% confidence set from the basic test in Equation (ref) --- and 5,000 parameter values from $\Theta\setminus \widehat\Theta_{0.95}$, and then calculated the overlap between the test decisions. The likelihood ratio test rejected 3% of the values from the first set, and 98% of the values from the second set, indicating a good amount of overlap between the two tests. The correlation between the p-value from the likelihood ratio test, and the values $f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}) \nu(f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}), \boldsymbol{\theta})$, which our basic test uses to make a decision in Equation (ref), was equal to 0.94. The correlation with the alternative confidence set construction is 0.90, using instead the values $\mathbb{I}\big\{f(\boldsymbol{s} | \boldsymbol{\theta}) \le f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}) \big\} f(\boldsymbol{s} | \boldsymbol{\theta})$ in the above calculation. Since the likelihood ratio test is exact, these results suggest that our test procedures are generally high-powered.
In Figures (ref) and (ref) of Appendix (ref), we plot the 95% confidence sets from the likelihood ratio test described above for the Santa Clara study and the LA county study (of the following section). The estimated prevalence is 0%-1.9% for Santa Clara, which is shorter than $\widehat\Theta_{0.95}$ but wider than $\widehat\Theta_{0.95}^{\mathrm{ alt}}$, as reported earlier; the same holds for LA county. As with our method, prevalence here is estimated through direct projection of the confidence set, which may be conservative. It is also possible that with more samples the likelihood ratio test could achieve the same interval as $\widehat\Theta_{0.95}^{\mathrm{ alt}}$ (we used only 100 samples), but this would come at an increased computational cost. Overall, the likelihood ratio test produces very similar results to our method, but it is not as efficient computationally.
In recent work, chen2018monte proposed a Monte Carlo-based method of inference in partially identified models. The idea is to sample from a quasi-posterior distribution, and then calculate $q_n$, the 95% percentile of $\{ f(\boldsymbol{s}_\mathrm{obs} | \boldsymbol{\theta}^{(j)}),~j=1, \ldots\}$, where $\boldsymbol{\theta}^{(j)}$ denotes the $j$-th sample from the posterior. The 95% confidence set is then defined as:
We implemented this procedure with an MCMC chain that appears to be mixing well --- see Appendix (ref) and Figure (ref) for details. The 95% confidence set, $\widehat\Theta$, is given in Figure (ref) of Appendix (ref). Simple projection, yields a prevalence in the range 0.9%-1.43%. This suggests that our MCMC “spends more time" around the mode of the likelihood, which we back up with numerical evidence in Appendix (ref). Finally, we also tried Procedure 3 of chen2018monte, which does not require MCMC simulations but is generally more conservative. Prevalence was estimated in the range 0.12%-1.65%, which is comparable to our method and the likelihood ratio test.
Next, we analyze the results from a recent serology study in Los Angeles county, which estimated a prevalence of 4.1% over the entire county population.\footnote{\url{http://publichealth.lacounty.gov/phcommon/public/media/mediapubhpdetail.cfm?prid=2328}} We use the same validation study as before since this study was executed by the same team as the Santa Clara one. Here, the main study had $N_\mathrm{m}=846$ participants with $s_{\mathrm{m}, \mathrm{obs}}=35$ positives.\footnote{ This number was not reported in the official study announcement mentioned above. It was reported in a Science article referencing one of the authors of the study: \url{https://www.sciencemag.org/news/2020/04/antibody-surveys-suggesting-vast-undercount-coronavirus-infections-may-be-unreliable}.} For inference, we only use the alternative construction, $\widehat\Theta_{0.95}^{\mathrm{ alt}}$, of Equation (ref) to simplify exposition. The results are shown in Figure (ref).
In contrast to the Santa Clara study, we see that the results from this study are conclusive. The prevalence rate is estimated in the range 1.7%-5.2%. If the false positive rate is, for example, closer to its empirical estimate (0.5%) then the identified prevalence is relatively high, somewhere in the range 3%-5.2%. We also see that the true positive rate is estimated in the range 85%-95%, which is higher than the empirical point estimate of 80% provided by bendavid2020covid. In fact, the empirical point estimate is not even in the 95% confidence set. Finally, as an illustration, we combine the data from the Santa Clara and LA county studies. The assumption is that the characteristics of the tests used in both studies were identical. The results are shown in Figure (ref) of Appendix (ref). We see that 0% prevalence is consistent with the combined study as well. Furthermore, prevalence values higher than 2.5% do not seem plausible in the combined data.
Recently, a quasi-randomized study was conducted in New York state, including NYC, which sampled individuals shopping in grocery stores. Details about this study were not made available. Here, we assume that the medical testing technology used was the same as in the Santa Clara and LA county studies, or at least similar enough that the comparison remains informative.
Under this assumption, we can use the same validation study as before, with $(N_\mathrm{c}^-, N_\mathrm{c}^+) = (401, 197)$ participants in the validation study, and $(s_{\mathrm{c}, \mathrm{obs}}^-, s_{\mathrm{c}, \mathrm{obs}}^+)=(2, 178)$ positives, respectively. The main study in New York had $N_\mathrm{m}=3000$ participants with $s_{\mathrm{m}, \mathrm{obs}}=420$ observed test positives.\footnote{ \url{https://www.nytimes.com/2020/04/23/nyregion/coronavirus-antibodies-test-ny.html}} The $\widehat\Theta_{0.95}^{\mathrm{ alt}}$ confidence set on this dataset is shown in Figure (ref). We see that the evidence in this study is much stronger than the Santa Clara/LA county studies with an estimated prevalence in the range 12.9%-16.6%. The true positive rate is now an important identifying parameter in the sense that knowing its true value could narrow down the confidence set even further.
Finally, in Figure (ref) of Appendix (ref) we present prevalence estimates for a combination of all datasets presented so far, while using both constructions, $\widehat\Theta_{0.95}$ and $\widehat\Theta_{0.95}^{\mathrm{ alt}}$, to illustrate their differences. As mentioned earlier, this requires the assumption that the antibody testing kits used in all three studies had identical specifications, or at least very similar so that the comparison remains informative. This assumption is most likely untenable given the available knowledge. However, we present the results there for illustration and completeness. The general picture in the combined study is a juxtaposition of earlier findings. For example, both false and positive rates are now important for identification. The identified prevalence is in the range 5.2%-8.2% in $\widehat\Theta_{0.95}^{\mathrm{ alt}}$ (and 3.2%-8.9% in $\widehat\Theta_{0.95}$). These numbers are larger than the Santa Clara/LA county studies but smaller than the New York study.
The procedure described in Section (ref) is computationally intensive for two main reasons. First, we need to consider all values of $\boldsymbol{\theta}\in\Theta$, which is a three-dimensional grid. Second, given some $\boldsymbol{\theta}$, we need to calculate $f(\boldsymbol{s} | \boldsymbol{\theta})$ for each $\boldsymbol{s}\in\mathbb{S}$, which is also a three-dimensional grid.
To deal with the first problem we can use parallelization, since the test decisions in step 3 of our procedure are independent of each other. For instance, the results in Section (ref) were obtained in a computing cluster (managed by Slurm) comprised of 500 nodes, each with x86 architecture, 64-bit processors, and 16GB of memory. The total wall clock time to produce all results of the previous section was about 1 hour. The results for, say, the Santa Clara study can be obtained in much shorter time (a few minutes) because they contain few positive test results. To address the second computational bottleneck we can exploit the independence property between $S_\mathrm{c}^-, S_\mathrm{c}^+$, and $S_\mathrm{m}$, as shown in the product of Equation (ref). Since any zero term in this product implies a zero value for $f$, we can ignore all individual term values that are very small. Through numerical experiments, we estimate that this computational trick prunes on average 97% of $\mathbb{S}$ leading to a significant computational speedup. For example, to test one single value $\boldsymbol{\theta}\in\Theta$ takes about 0.25 seconds in a typical high-end laptop, which is a 200-fold speedup compared to 50 seconds required by the likelihood ratio test of Section (ref) --- see Appendix (ref) for more details.
As mentioned earlier, prevalence $\pi$ in Equation (ref) is a finite-population estimand, that is, it is a number that refers to the particular population in the study. Theorem (ref) shows that our procedure is valid for $\pi$ only under Assumption (A4). However, to extrapolate to the general population we generally need to assume that
This is currently an untenable assumption. For example, in the Santa Clara study the population of middle-aged white women was overrepresented, while the population of Asian or Latino communities was underrepresented. The impact from such selection bias on the inferential task is very hard to ascertain in the available studies. Techniques such as post-stratification or reweighing can help, but at this early stage any extrapolation using distributional assumptions would be too speculative. However, selection bias is a well-known issue among researchers, and can be addressed as widespread and carefully designed antibody testing catches on. We leave this for future work.
In this paper, we presented a partial identification method for estimating prevalence of Covid-19 from randomized serology studies. The benefit of our method is that it is valid in finite samples, as it does not rely on asymptotics, approximations or normality assumptions. We show that some recent serology studies in the US are not conclusive (0% prevalence is in the 95% confidence set). However, the New York study gives strong evidence for high prevalence in the range 12.9%-16.6%. A combination of all datasets shifts this range down to 5.2%-8.2%, under a test uniformity assumption. Looking ahead, we hope that the method developed here can contribute to a more robust analysis of future Covid-19 serology tests.
I would like to thank Guanglei Hong, Ali Hortascu, Chuck Manski, Casey Mulligan, Joerg Stoye, and Harald Uhlig for useful suggestions and feedback. Special thanks to Connor Dowd for his suggestion of the alternative construction (ref), and to Elie Tamer for various important suggestions. Finally, I gratefully acknowledge support from the John E. Jeuck Fellowship at Booth School of Business.