Estimating SARS-CoV-2 Seroprevalence

During the COVID-19 pandemic, governments and public health authorities have used seroprevalence studies to estimate the proportion of persons within a given population who have antibodies to SARS-CoV-2. Seroprevalence is crucial for estimating quantities such as the infection fatality ratio, proportion of asymptomatic cases, and differences in infection rates across population subgroups. However, serologic assays are prone to false positives and negatives, and non-probability sampling methods may induce selection bias. In this paper, we consider nonparametric and parametric seroprevalence estimators that address both challenges by leveraging validation data and assuming equal probabilities of sample inclusion within covariate-defined strata. Both estimators are shown to be consistent and asymptotically normal, and consistent variance estimators are derived. Simulation studies are presented comparing the finite sample performance of the estimators over a range of assay characteristics and sampling scenarios. The methods are used to estimate SARS-CoV-2 seroprevalence in asymptomatic individuals in Belgium and North Carolina.



There are no comments yet.


page 1

page 2

page 3

page 4


Large sample properties of the Midzuno sampling scheme

Midzuno sampling enables to estimate ratios unbiasedly. We prove the asy...

Simple estimators for network sampling

Some conceptually simple estimators for network sampling are introduced....

Estimating the reciprocal of a binomial proportion

As a classic parameter from the binomial distribution, the binomial prop...

Semi-Parametric Estimation of Incubation and Generation Times by Means of Laguerre Polynomials

In epidemics many interesting quantities, like the reproduction number, ...

Statistical Considerations for Cross-Sectional HIV Incidence Estimation Based on Recency Test

Longitudinal cohorts to determine the incidence of HIV infection are log...

Estimating the potential to prevent locally acquired HIV infections in a UNAIDS Fast-Track City, Amsterdam

Amsterdam and other UNAIDS Fast-Track cities aim for zero new HIV infect...

Checking validity of monotone domain mean estimators

Estimates of population characteristics such as domain means are often e...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

SARS-CoV-2 seroprevalence has been estimated in hundreds of studies in various populations worldwide during the COVID-19 pandemic. Seroprevalence studies, which estimate the proportion of persons in a population with detectable antibodies to a pathogen, are crucial for understanding the pandemic. Knowledge of this proportion enables the quantification of, for instance, the total number of infected people; differences in SARS-CoV-2 infection by demographic, geographic, or occupational characteristics; and the infection fatality rate (Arora et al., 2020). These factors have guided governmental and public health policy decisions about when to lock down and reopen sectors of society during the COVID-19 pandemic. Thus, accurate seroprevalence estimates are needed to inform these crucial decisions. However, seroprevalence studies usually suffer from at least two sources of bias: measurement error due to false positives and negatives, and selection bias because convenience sampling designs are often used due to time and cost constraints.

Typically, blood tests for antibodies result in a continuous measure of antibody response which is dichotomized into a positive or negative classification based on a pre-specified cutoff value. This dichotomization almost always produces measurement error in the form of false positives and false negatives (Bouman et al., 2021). The following example from Sempos and Tian (2021) demonstrates how this measurement bias can often be large. Suppose that true seroprevalence is 1% and antibody tests are performed using an assay which perfectly identifies true positives as positive, so with 100% sensitivity, and nearly perfectly identifies true negatives as negative, with 99% specificity. Despite this assay’s high sensitivity and specificity, naively using the sample proportion of positive test results as a seroprevalence estimator would, in expectation, lead to an estimate of nearly 2% rather than 1%.

To appropriately account for measurement error, sensitivity and specificity should be estimated and incorporated into a seroprevalence estimator. The sensitivity and specificity of a diagnostic test are commonly estimated by performing the assay on samples of known (or gold standard) positives and negatives, respectively. For example, samples from patients who had a case of COVID-19 confirmed with reverse transcription polymerase chain reaction (RT-PCR) testing are often assumed to be true positives. Remnant blood samples that were drawn in 2019 or earlier, before major spread of SARS-CoV-2, are often assumed to be true negatives. Sensitivity and specificity are estimated from these validation data as the proportion of samples appropriately classified as positive and negative, respectively. These estimates of sensitivity and specificity can then be incorporated into the seroprevalence estimator, often with a method popularized by Rogan and Gladen (1978) (see also Marchevsky, 1979), to appropriately adjust the seroprevalence estimate for assay characteristics.

Many seroprevalence studies are conducted by non-probability sampling methods such as convenience sampling, which may lead to selection bias when characteristics that drive participation in the study are also risk factors for SARS-CoV-2 infection. Probability-based sampling studies are ideal because they are usually more representative and lead to less biased estimates than convenience samples analyzed with post-hoc statistical adjustments (Shook-Sa, Boyce, & Aiello, 2020; Accorsi et al., 2021). However, probability-based sampling may not always be feasible. Convenience sample based seroprevalence estimators often rely on the assumption that each person in a covariate-defined stratum has an equal probability of being in the sample (Shook-Sa et al., 2020). Under this assumption, population prevalence can be estimated with direct standardization, though complex survey methods like calibration have been used as well (e.g., Bajema et al., 2020).

In this paper, methods are considered which combine standardization and the Rogan-Gladen adjustment to account for both measurement error and selection bias. The article is organized as follows. Section 2 reviews prevalence estimation under measurement error. Nonparametric and parametric standardized prevalence estimators and their large-sample properties are described in Section 3. Section 4

presents simulation studies to evaluate the empirical bias and 95% confidence interval (CI) coverage of the standardized estimators across a range of assay characteristics and bias scenarios. The methods are then applied in Section 

5 to estimate seroprevalence of SARS-CoV-2 among all residents of Belgium in 2020 and, separately, in asymptomatic residents of North Carolina during the spring of 2020. Section 6 concludes with a discussion. Proofs are in the Appendix.

2 Seroprevalence estimation under measurement error

2.1 Problem setup

Let the true serology status for an individual in the target population be denoted by , with if the individual is seropositive and otherwise. Our goal is to draw inference about the population seroprevalence . Because of error in the serology assay, is not observed directly. Let the result of the serology assay be denoted by , with if the individual tests positive and otherwise. Three key quantities are sensitivity, the probability that a true positive tests positive, denoted by ; specificity, the probability that a true negative tests negative, denoted by ; and the population expectation of the serology assay outcome, denoted by . Unless the assay has perfect sensitivity and specificity with , typically will not equal and will be a mismeasured version of .

Suppose validation studies are conducted to estimate the sensitivity and specificity by applying the assay to known positives and negatives. Specifically, the serology assay is conducted on independent and identically distributed (iid) units from strata of the population where and on iid units from strata where . Thus copies of are observed to estimate sensitivity and copies of are observed to estimate specificity. To estimate seroprevalence in a target population, a ‘main’ study with iid copies of is then conducted, among which true infection status is unknown.

Assume, as is realistic in many SARS-CoV-2 studies, that there is no overlap between the units in each of the three studies. Let be an indicator of which study the th individual’s sample is from, with for the sensitivity study, for the specificity study, and for the main study. Note that for , where and here and throughout summations are taken from to unless otherwise specified. Assume as .

2.2 Estimators and statistical properties

Let . Consider the estimator , where , , , and . The prevalence estimator is motivated by rearranging the identity that and is sometimes referred to as the Rogan-Gladen (1978) estimator.

The estimator can be expressed as the solution (for

) to the estimating equation vector

where here and below 0 denotes a column vector of zeros. Since the samples were selected from three different populations, the data are not identically distributed and care must be taken to derive the large sample properties of . In Appendix A.1, the estimator is shown to be consistent and asymptotically normal. Specifically, as , and assuming (as discussed below), where is a covariance matrix with bottom right element


The proof is similar to standard estimating equation theory (e.g., Boos and Stefanski, 2013, Equation 7.10), but because the data are not identically distributed the Lindeberg-Feller Central Limit Theorem (CLT) is used in place of the classical Lindeberg-Lévy CLT.

Let denote the plug-in estimator defined by replacing , and in (1) with , and for , and note that is the variance estimator proposed by Rogan and Gladen (1978). By the continuous mapping theorem, is consistent for the asymptotic variance assuming and can be used to construct Wald-type CIs that asymptotically attain nominal coverage probabilities.

2.3 Truncation into

In finite samples sometimes yields estimates outside of when (i) , (ii) , or (iii) . Indeed, (ii) occurred in the ScreenNC study discussed in Section 5.2. Estimates are typically truncated to be inside because the true population prevalence must exist in (Hilden, 1979). In this article, all point estimates and bounds of interval estimates are so truncated. Note, though, that as the three sample sizes grow large the estimator yields estimates inside almost surely unless . In practice, settings where may be very unlikely; in such scenarios, the probability of a positive test result is higher for seronegative persons than for seropositive persons, meaning such an assay performs worse in expectation than random guessing. Throughout this manuscript, we assume .

3 Standardized seroprevalence estimation

In some settings it may not be reasonable to assume the copies of from the main study constitute a random sample from the target population. Suppose instead that for each copy of a vector of discrete covariates is observed, with taking on possible values . The covariates are of interest because seroprevalence may differ between the strata; for instance, might include demographic variables such as age group, race, or gender. Denote the mean of in the th stratum as and the sample size for the th stratum as , so .

The distribution of strata in the target population, if known, can be used to standardize estimates so they are reflective of the target population (for a review of direct standardization, see van Belle et al. (2004, Ch. 15)). Denote the proportion of the target population comprised by the th stratum as and suppose that these stratum proportions are known with each and . For instance, if the target population were all the adults in a US state, then could be obtained from US census data. Assume all persons in a covariate stratum defined by have the same probability of inclusion in the sample. Then the covariates in the main study sample have a multinomial distribution with categories, sample size , and an unknown sampling probability vector where . Note that if the main study were a simple random sample from the target population, then the sampling probabilities would be equal to the stratum proportions, i.e., for .

3.1 Nonparametric standardization

First, consider a seroprevalence estimator which combines nonparametric standardization and the Rogan-Gladen adjustment to account for both selection bias and measurement error. Notice that is a weighted average of the stratum-conditional means where each weight is a known stratum proportion , i.e., . A nonparametric standardization estimator for using the sample stratum-conditional prevalences for is . A standardized prevalence estimator accounting for measurement error is , which has been used in recent SARS-CoV-2 seroprevalence studies (Havers et al., 2020; Barzin et al., 2020; Cai et al., 2020).

Letting , the estimator solves the estimating equation vector where , and are defined in Section 2; is a -vector with th element ; and . It follows that is consistent and asymptotically normal and that where


The asymptotic variance can be consistently estimated by the plug-in estimator defined by replacing , and in (2) with , and for and . Consistency of holds by continuous mapping, and a proof of asymptotic normality and justification of (2) are in Appendix B.

Standardization requires estimating the stratum-conditional mean of , . However, when for some strata , the corresponding estimator is undefined, and since is then undefined as well alternative estimation strategies are needed. Values of may equal zero for two reasons. First, the study design may exclude these strata (), a situation referred to as deterministic or structural nonpositivity (Westreich and Cole, 2010). Second, even if , random nonpositivity can occur if no individuals with are sampled, which may occur if is small or if is relatively small. When nonpositivity arises, an analytical approach often employed entails “restriction” (Westreich and Cole, 2010), where the target population is redefined to consist only of strata for which . However, this redefined target population may be less relevant from a public health or policy perspective.

3.2 Parametric standardization

Now, consider a parametric approach for modeling the stratum-conditional means , an alternative strategy which allows inference to the original target population and may also perform better when some strata have small sample sizes. Assume the binary regression model holds where

is an appropriate link function for a binary outcome like the logit or probit function;

is a row vector of regression coefficients with intercept ; and is a user-specified -vector function of the th stratum’s covariate values that may include main effects and interaction terms, with th element denoted and set equal to one to correspond to an intercept. Let be the covariate support in the sample, i.e., with dimension , and assume . (Note that only when there is positivity, and in that case can be used with no restriction needed.)

Under the assumed binary regression model, each is a function of the parameters and the covariates that define the th stratum, denoted . A model-based standardized Rogan-Gladen estimator of is , where and is the maximum likelihood estimator of . Estimating equation theory can again be used to derive large-sample properties by replacing the equations for from Section 3.1 with equations for corresponding to the score equations from the binary regression.

Let . The estimator is the solution to the vector of estimating equations where , and are as in Section 2; is a -vector with th element ; and . It follows that is consistent and asymptotically normal and . The asymptotic variance can be consistently estimated by , the lower right element of the empirical sandwich variance estimator of the asymptotic variance of . A proof of asymptotic normality and the empirical sandwich variance estimator are given in Appendix C.

4 Simulation study

A simulation study was conducted to compare , , and . Four data generating processes (DGPs) were considered, within which different scenarios were defined through full factorial designs that varied simulation parameters , and . These DGPs featured no selection bias (DGP 1), selection bias with two strata (DGP 2), and more realistic selection bias settings with 40 strata and 80 strata (DGPs 3 and 4).

For each DGP and set of simulation parameters, sensitivity and specificity validation samples of size and were generated with distributed Bernoulli with a mean of or , respectively. In DGPs 1 and 2 a main study of size was then generated where was Bernoulli with mean and was Bernoulli with mean ; in DGPs 3 and 4 was generated directly from the distribution of , as will be explained. Simulation parameter values were selected based on the seroprevalence studies described in Section 5. Sensitivity was varied in , specificity in , and prevalence in . Sample sizes were , , and . The full factorial design led to 120 scenarios per DGP, and within each scenario 1,000 simulations were conducted. Performance was measured across the 1,000 simulations by (a) mean bias, computed as for each estimator ; and (b) empirical coverage, i.e., whether the 95% Wald-type CIs based on each variance estimator contained the true prevalence.

4.1 No selection bias

Simulations for DGP 1 were conducted to assess the performance of when no selection bias was present. The estimator was generally unbiased, as seen in Figure A.1. Performance improved as and tended toward 1, with specificity being a stronger determinant of bias. An exception to these results was for low prevalence (0.05 or lower), when overestimated the true prevalence. Wald CIs based on generally attained nominal coverage, as seen in Figure A.2. However, when prevalence was near zero, coverage of the 95% CIs did not equal the nominal level. Coverage also decreased as increased toward 1. Coverage properties of Wald CIs for a single binomial parameter are known to be variable when the parameter is near 0 or 1 (Brown, Cai & DasGupta, 2001). This issue could be exacerbated for Wald CIs based on , as they estimate three binomial parameters in , , and . Indeed, these results concord with previous simulation studies evaluating (Lang and Reiczigel, 2014). Alternative methods for constructing CIs are mentioned in Section 6.

4.2 Low-dimensional selection bias

In DGP 2, the target population was comprised of two strata defined by a covariate with proportions . Within the main study, was generated from a multinomial distribution of sample size and sampling probability vector . Individuals’ serostatuses were generated from the conditional distribution , which was such that and for each value of . In each simulation and and their corresponding 95% CIs were computed.

The nonparametric standardized estimator was empirically unbiased for true prevalences , as seen in Figure A.3, and 95% CIs based on attained approximately nominal coverage, as seen in Figure A.4. As with in DGP 1, CI coverage fell slightly below the nominal level for very low prevalences, and more noticeably for . Figures A.3 and A.4 show that performed poorly under selection bias, with large negative bias and CI coverage far below the nominal level in most cases.

4.3 More realistic selection bias

DGPs 3 and 4 compared and in scenarios with larger numbers of strata.

4.3.1 Dgp 3

Three covariates were defined as , , and , leading to strata with proportions . Within the main study, was again generated as multinomial with size and known sampling probability vector. Figure 1(a) shows the structure of selection bias in DGP 3 by comparing the stratum proportions and sampling probabilities. Some low-prevalence strata that frequently occur in the population were oversampled, while most remaining strata were undersampled. Individuals’ test results were generated from the conditional distribution , where

The parameters and were set to reflect differential prevalences by stratum, while a “balancing intercept” (Rudolph et al., 2021) was set to different values so that prevalence equaled (approximately) . The nonparametric estimator and corresponding CI were computed using a restricted target population when random nonpositivity arose; the values of prevalence used to compute bias and coverage were based on the total (unrestricted) population, which is the parameter of interest. The parametric estimator

was computed with a correctly-specified logistic regression model.

Figure 1: Panels A and B represent selection bias in the simulation studies of DGPs 3 and 4, described in Sections 4.3.1 and 4.3.2, respectively. Circle size is proportional to prevalence. Points are jittered slightly for legibility, and the diagonal lines denote equality between (stratum proportion) and (sampling probability).

Both and performed well in this scenario. Figure 2 shows that the estimators were generally empirically unbiased, though modest bias occurred when and prevalence was low. Figure A.5 shows 95% CIs based on either or attained nominal coverage, with slight under-coverage for . On average across all 120 scenarios, 89% (range of 86%-92%) of simulated datasets had positivity. As in DGP 2, was seriously biased and the corresponding CIs did not attain nominal coverage.

Figure 2: Empirical bias of the Rogan-Gladen (), nonparametric standardized (), and logistic regression standardized () estimators from simulation study for DGP 3, described in Section 4.3.1. The six facets correspond to a given combination of sensitivity (‘Sens’) and specificity (‘Spec’).

4.3.2 Dgp 4

Data were generated as in DGP 3, but the inclusion of a fourth covariate now led to 80 strata. The conditional distribution was such that , where here contains the same terms as in DGP 3 plus a main effect for with corresponding coefficient . Regression parameters were a balancing intercept , , , , , , and . The larger value for , as compared to , led to a stronger relationship between and than was present in DGP 3. Figure 1(b) displays selection bias in DGP 4. Some of the highest-prevalence and most commonly-occurring strata were undersampled to a greater degree than occurred in DGP 3, so in this sense there was more selection bias in DGP 4.

Figure 3 shows that only was empirically unbiased in this scenario, while typically had a moderately negative bias. Nonpositivity almost always occurred (in either all or all but one of the simulations, for each of the 120 scenarios). The worse performance of may be explained by restriction leading to bias under nonpositivity. CIs based on typically attained nominal coverage, unlike those based on or , as seen in Figure A.6. However, the trend of empirical coverage being below the nominal level for low prevalence and was more noticeable compared to other DGPs, perhaps due to greater selection bias. Note that was negative for two simulations in a single scenario, and these “Heywood cases” (Kolenikov and Bollen, 2012) were ignored in the coverage calculation for that scenario.

Figure 3: Bias results from simulation study on DGP 4, described in Section 4.3.2. Figure layout is as in Figure 2.

In summary, both the nonparametric and parametric standardized estimators and had low empirical bias and close to nominal 95% CI coverage when there was little nonpositivity. As the number of covariates, amount of selection bias, and potential for nonpositivity increased, the (correctly-specified) parametric generally maintained its performance while had greater empirical bias and the corresponding CIs did not attain nominal coverage levels.

4.4 Model misspecification

The performance of was also assessed in scenarios similar to DGPs 3 and 4, but under mild model misspecification. Here, the true conditional distributions of were and , where and are the specifications used in the models for in DGPs 3 and 4, respectively. The test results were generated from as in DGPs 1 and 2. The degree of misspecification in each scenario thus depended on the values of , , and . Figures A.7 through A.10 show that, in terms of bias and 95% CI coverage, was robust to the misspecification considered.

5 Applications

5.1 Belgium seroprevalence study

The standardized Rogan-Gladen methods were applied to a nationwide SARS-CoV-2 seroprevalence study in Belgium conducted across seven week-long collection rounds between March and October 2020 (Herzog et al., 2021). Residual sera were collected in a stratified random sample from private laboratories encompassing a wide geographical network, with stratification by age group (10 year groups from 0-9, 10-19, …, 90-plus), sex (male or female), and region (Wallonia, Flanders, or Brussels). The presence of SARS-CoV-2 IgG antibodies was determined using a semi-quantitative EuroImmun ELISA. Based on validation studies of RT-PCR confirmed COVID-19 cases and pre-pandemic negative controls, sensitivity and specificity were estimated to be and (Herzog et al., 2021, Table S1.1). The number of samples for assessing seroprevalence varied between 2,960 and 3,910 across the seven collection rounds.

In our analysis, we estimated nationwide seroprevalence in Belgium during each collection round standardized by age group, sex, and province (11 total), using 2020 stratum proportion data from the Belgian Federal Planning Bureau (2020). Province was used rather than region to match the covariates selected for weighting in Herzog et al. (2021). In six of the seven collection rounds nonpositivity arose, with between 2 to 15 of the 220 strata not sampled, so restricted target populations were used for computation of . For , each logistic regression model had main effects for age group, sex, and province, as well as an interaction term between age group and sex.

Figure 4 displays point estimates and CIs for , , and alongside those for the unadjusted, or naive, sample prevalence for each collection round (with exact Clopper-Pearson 95% CIs). The naive estimates typically are greater than the other three estimates and have shorter CIs. The differences between and are greater than those between and the two standardized estimators, suggesting that the magnitude of measurement error here may have been larger than that of selection bias.

Figure 4: Estimates and corresponding 95% confidence intervals for each of seven collection rounds for the 2020 Belgian seroprevalence study (Herzog et al., 2021), described in Section 5.1.

5.2 North Carolina seroprevalence study

The standardization methods of Section 3 were also applied to ScreenNC, which tested a convenience sample of 2,973 asymptomatic patients age 20 and older in North Carolina (NC) for antibodies to SARS-CoV-2 between April to June 2020 (Barzin et al., 2020). These patients were seeking unrelated medical care at eleven sites in NC associated with the University of North Carolina (UNC) Health Network. The presence of antibodies was determined with the Abbott Architect SARS-CoV-2 IgG assay. Based on validation studies of RT-PCR confirmed positive patients and pre-pandemic serum samples assumed to be negative, sensitivity was estimated as and specificity as .

In our analysis, seroprevalence was estimated in two relevant target populations. First, we standardized to patients accessing the UNC Health Network during a similar timeframe (21,901 patients from February to June of 2020). The main study sample differed from this UNC target population in terms of age group, race, and sex characteristics, as seen in Table 1, and meta-analyses suggested that prevalence of COVID-19 infections differed between levels of these covariates in some populations (Pijls et al., 2021; Mackey et al., 2021), supporting the covariates’ use in standardization. Note that several racial classifications, including patient refused and unknown, were reclassified as ‘Other’. Second, we standardized to the 2019 NC population over the age of 20 (7,873,971 persons) using covariate data from the American Community Survey (U.S. Census Bureau, 2019). The assumption of equal sampling probabilities may be less reasonable for this target population because not all NC residents are in the UNC Health Network and because there were some geographic areas where no patients in the study sample were from. There was no sample data in the main study for two covariate strata that existed in the UNC Health Network, so restriction was used for . Logistic regression models with main effects for sex, race, and age group were used to compute ; interaction effects were not included as the small number of positive test results could have led to model overfit.

ScreenNC UNC Hospitals NC ACS
% % %
2,973 100 31,095 100 7,873,971 100
Sex Female 1,955 66 13,962 64 4,108,603 52
Male 1,018 34 7,975 36 3,765,368 48
Race Asian 67 2 460 2 230,759 3
Black or Af.-Am. 395 13 5,109 23 1,640,311 21
Other 311 10 1,843 6 455,600 6
White or Cauc. 2,200 74 21,074 68 5,547,301 70
Age 20-29 342 12 2060 9 1,400,918 18
30-39 599 20 2,763 13 1,344,647 17
40-49 518 18 3,382 15 1,351,156 17
50-59 602 20 4,200 19 1,360,357 17
60-69 489 17 4,548 21 1,228,123 16
70-79 310 11 3,325 15 806,002 10
80+ 77 3 1,623 7 382,768 5
Table 1: Demographic comparisons of the ScreenNC study sample, UNC Hospitals patient population, and North Carolina population aged 20+. Data on the NC population are from the 2019 American Community Survey (ACS). Several racial classifications including Patient Refused and Unknown were reclassified as Other. Sample size is denoted by . Some column totals do not sum to 100% because of rounding. Abbreviations: Af.-Am for African-American and Cauc. for Caucasian.

The sample proportion of positive tests was . The sample false positive rate was , so the data are, at first appearance, consistent with a population prevalence of 0%. Indeed, the Rogan-Gladen seroprevalence estimate was % (95% CI 0%, 1.00%). Likewise, the UNC target population had nonparametric and parametric standardized estimates of % (0%, 1.11%) and % (0%, 1.13%), and the NC target population had corresponding estimates of 0% (0%, 1.10%) and 0% (0%, 1.11%). All estimates were truncated into . The closeness of the standardized and unstandardized results may be due to the small number of positive test results and similarities between the sample and the target populations.

6 Discussion

In this paper we examined nonparametric and model-based standardized Rogan-Gladen estimators, deriving their large-sample properties and consistent variance estimators. Simulation studies demonstrated that both methods had low empirical bias and nominal CI coverage in the majority of practical settings. The empirical results in Section 4 highlight the tradeoffs inherent in choosing which method to use for a seroprevalence study. The parametric standardized estimator was empirically unbiased even when the number of strata and covariates, and with them the potential for random nonpositivity, increased. A drawback to is the need to correctly specify the form of a regression model. On the other hand, the nonparametric standardized estimator does not require model specification and performed well in scenarios with lower amounts of selection bias and nonpositivity. As the number of strata and covariates grew, however, was empirically biased and its corresponding 95% CIs did not attain nominal coverage.

For practical use of either method, careful choice of covariates is necessary. Including additional covariates may make the assumption of equal probability of sampling within strata more reasonable. However, such inclusion also makes covariate-defined strata smaller and random nonpositivity more likely. An alternative strategy is to collapse smaller strata with few or no persons to create larger strata, which may make random nonpositivity less likely. However, if strata with sufficiently different sampling probabilities were collapsed together, then the assumption of equal probability of sampling within strata would be violated.

Topics for future research and broader issues in SARS-CoV-2 seroprevalence studies merit mention. As an alternative to standardization, inverse probability of sampling weights (Lesko et al., 2017) or inverse odds of sampling weights (Westreich et al., 2017) can be considered. Standardization and weighting methods could then be combined to create a doubly robust Rogan-Gladen estimator. Bootstrap (Cai et al., 2020) or Bayesian (Gelman and Carpenter, 2020) CIs are alternatives that may address the limitations of Wald-type CIs; likewise, several alternative approaches to CI construction based on test inversion that account for measurement error (DiCiccio et al., 2021) could be extended to account for selection bias. Methods could also be developed to adjust for “spectrum bias”, which can arise when validation studies for assay sensitivity are sampled from patients that are sicker than the average infected person (Takahashi, Greenhouse, and Rodríguez-Barraquer, 2020).


  • [1] Accorsi, E. K., Qiu, X., Rumpler, E., Kennedy-Shaffer, L., Kahn, R., Joshi, K., et al. (2021). How to detect and reduce potential sources of biases in studies of SARS-CoV-2 and COVID-19. European Journal of Epidemiology, 36, 179-196.
  • [2] Agresti, A. & Coull, B. A. (1998). Approximate is better than “exact” for interval estimation of binomial proportions. The American Statistician, 52(2), 119-126.
  • [3] Arora, R. K., Joseph, A., Van Wyk, J., Rocco, S., Atmaja, A., May, E., et al. (2021). SeroTracker: A global SARS-CoV-2 seroprevalence dashboard. The Lancet Infectious Diseases, 21(4), e75-e76.
  • [4] Bajema, K. L., Wiegand, R. E., Cuffe, K., Patel, S. V., Iachan, R., Lim, et al. (2021). Estimated SARS-CoV-2 seroprevalence in the US as of September 2020. JAMA Internal Medicine, 181(4), 450-460.
  • [5] Barzin, A., Schmitz, J. L., Rosin, S., Sirpal, R., Almond, M., Robinette, C., et al. (2020). SARS-CoV-2 seroprevalence among a southern U.S. population indicates limited asymptomatic spread under physical distancing measures. MBio, 11(5), e02426-20.
  • [6] Boos, D. D. & Stefanski, L. A. (2013). Essential Statistical Inference. Springer New York, Ch. 7, pp. 297-337.
  • [7] Bouman, J. A., Riou, J., Bonhoeffer, S., & Regoes, R. R. (2020). Estimating cumulative incidence of SARS-CoV-2 with imperfect serological tests: Exploiting cutoff-free approaches. PLoS Computational Biology, 17(2), e1008728.
  • [8] Cai, B., Ioannidis, J. P. A., Bendavid, E., & Tian, L. (2020). Exact inference for disease prevalence based on a test with unknown specificity and sensitivity.
  • [9] DiCiccio, T.J., Ritzwoller, D.M., Romano, J.P. & Shaikh, A.M. (2021). Confidence intervals for seroprevalence.
  • [10] Federal Planning Bureau (2020). Population par province et âge, au 1er janvier. Electronic dataset, accessed August 18, 2021.
  • [11] Gelman, A. & Carpenter, B. (2020). Bayesian analysis of tests with unknown specificity and sensitivity. Journal of the Royal Statistical Society: Series C (Applied Statistics), 69(5), 1269–1283.
  • [12] Havers, F. P., Reed, C., Lim, T., Montgomery, J.M., Klena, J.D., Hall, A.J., et al. (2020). Seroprevalence of antibodies to SARS-CoV-2 in 10 sites in the United States, March 23-May 12, 2020. JAMA Internal Medicine 180(12), 1576-1586.
  • [13] Herzog, S., De Bie, J., Abrams, S., Wouters, I., Ekinci, E., Patteet, L., et al. (2021). Seroprevalence of IgG antibodies against SARS coronavirus 2 in Belgium – a serial prospective cross-sectional nationwide study of residual samples (March - October 2020). MedRxiv.
  • [14] Hilden, J. (1979). A further comment on “Estimating prevalence from the results of a screening test.” American Journal of Epidemiology, 109(6), 721–722.
  • [15] Kolenikov, S. & Bollen, K.A. (2012). Testing negative error variances: Is a Heywood case a symptom of misspecification?. Sociological Methods & Research, 41(1), 124-167.
  • [16] Lang, Z. & Reiczigel, J. (2014). Confidence limits for prevalence of disease adjusted for estimated sensitivity and specificity. Preventive Veterinary Medicine, 113(1), 13–22.
  • [17] Lesko, C.R., Buchanan, A.L., Westreich, D., Edwards, J.K., Hudgens, M.G., & Cole, S.R. (2017). Generalizing study results: a potential outcomes perspective. Epidemiology, 28(4), 553.
  • [18] Mackey, K., Ayers, C.K., Kondo, K.K., Saha, S., Advani, S.M., Young, S., et al. (2021). Racial and ethnic disparities in COVID-19–related infections, hospitalizations, and deaths: a systematic review. Annals of Internal Medicine, 174(3), 362-373.
  • [19] Marchevsky, N. (1979). Re: “Estimating prevalence from the results of a screening test.” American Journal of Epidemiology, 109(6), 720–721.
  • [20] Pijls, B.G., Jolani, S., Atherley, A., Derckx, R.T., Dijkstra, J.I., Franssen, G.H., et al. (2021). Demographic risk factors for COVID-19 infection, severity, ICU admission and death: a meta-analysis of 59 studies. BMJ Open, 11(1), e044640.
  • [21] Rogan, W. J. & Gladen, B. (1978). Estimating prevalence from the results of a screening test. American Journal of Epidemiology, 107(1), 71-76.
  • [22] Rudolph, J. E., Edwards, J. K., Naimi, A. I., & Westreich, D. J (2021). Simulation in practice: The balancing intercept. American Journal of Epidemiology, 190(8), 1696-1698.
  • [23] Sempos, C. T. & Tian, L. (2021). Adjusting coronavirus prevalence estimates for laboratory test kit error. American Journal of Epidemiology, 190(1), 109–115.
  • [24] Shook-Sa, B. E., Boyce, R. M., & Aiello, A. E. (2020). Estimation without representation: Early severe acute respiratory syndrome coronavirus 2 seroprevalence studies and the path forward. The Journal of Infectious Diseases, 222(7), 1086-1089.
  • [25] Takahashi, S., Greenhouse, B., & Rodriguez-Barraquer, I. (2020). Are SARS-CoV-2 seroprevalence estimates biased? The Journal of Infectious Diseases, 222(11), 1772-1775.
  • [26] U.S. Census Bureau (2019). American Community Survey 1-Year Estimates, Public Use Microdata Sample. Retrieved from, 2 February 2021.
  • [27] van Belle, G., Fisher, L. D., Heagerty, P. J., & Lumley, T. (2004). Biostatistics: a Methodology for the Health Sciences (2nd ed.). John Wiley & Sons, pp. 640-648.
  • [28] Westreich, D. & Cole, S. R. (2010). Invited commentary: Positivity in practice. American Journal of Epidemiology, 171(6), 674–677.
  • [29] Westreich, D., Edwards, J.K., Lesko, C.R., Stuart, E., & Cole, S.R. (2017). Transportability of trial results using inverse odds of sampling weights. American Journal of Epidemiology, 186(8), 1010-1014.


We thank Dirk Dittmer, the ScreenNC research team, and participants. We also thank Shaina Alexandria, Bryan Blette, and Kayla Kilpatrick for helpful suggestions and comments. This research was supported by the NIH (Grant R01 AI085073), the UNC Chapel Hill Center for AIDS Research (Grants P30 AI050410 and R01 AI157758), and the NSF (Grant GRFP DGE-1650116). The content is solely the responsibility of the authors and does not represent the official views of the NIH, UNC, or NSF.

Code and data availability statement

R code based on the simulations in Section 4 and data analysis in Section 5, data from the Belgium and North Carolina seroprevalence studies in Section 5, and a Microsoft Excel spreadsheet that computes the estimators and are available at

Appendix A Proofs for Section 2

a.1 Proof of asymptotic normality

Taylor expansion of around the true parameter yields

where and is a remainder term. Rearranging and multiplying by yields


where is a new remainder defined below. It is shown below that , , and , where and are defined below. Therefore, by Slutsky’s theorem, .

First, define

and let

As , by the assumption that for .

Second, for brevity, let denote , and similarly for , and . Define

and let

Note that as and likewise for and . It follows that by the Lindeberg-Feller CLT. In particular, let and let denote the Euclidean norm. Note that as the maximum magnitude of each element of is 1, and for a given observation it is always true that two of the indicators equal zero and the third indicator equals one. Therefore for all ,

implying the Lindeberg condition holds.

Third, it remains to prove . The outline of the proof of Boos and Stefanski (2013) Theorem 7.2 can be followed, but their assumption of identically distributed data must be removed. Consider the second-order Taylor series expansion of the th element of the vector , denoted , around the true value :

where are on the line segment joining and and is a matrix with entry equal to for . Writing these 4 equations in matrix notation yields


where and

is the matrix with th row given by . Note that