# Estimating SARS-CoV-2 Seroprevalence

Samuel Rosin
Department of Biostatistics
University of North Carolina at Chapel Hill
Chapel Hill, NC
<EMAIL_ADDRESS>
&Bonnie E. Shook-Sa
Department of Biostatistics
University of North Carolina at Chapel Hill
Chapel Hill, NC
&Stephen R. Cole
Department of Epidemiology
University of North Carolina at Chapel Hill
Chapel Hill, NC
&Michael G. Hudgens
Department of Biostatistics
University of North Carolina at Chapel Hill
Chapel Hill, NC

###### Abstract

Governments and public health authorities use seroprevalence studies to guide
their responses to the COVID-19 pandemic. These seroprevalence surveys
estimate the proportion of persons within a given population who have
detectable antibodies to SARS-CoV-2. However, serologic assays are prone to
misclassification error due to false positives and negatives, and non-
probability sampling methods may induce selection bias. In this paper, we
consider nonparametric and parametric prevalence 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.

_Keywords_ Covid-19 $\cdot$ diagnostic tests $\cdot$ estimating equations
$\cdot$ seroepidemiologic studies $\cdot$ standardization

##  1 Introduction

Estimating the proportion of people who have antibodies to SARS-CoV-2 is
useful for tracking the pandemic’s severity and informing public health
decisions (Arora et al.,, 2021). Individuals may have detectable antibodies
for different reasons, including prior exposure, infection, or vaccination.
Antibody levels within a person are dynamic, typically increasing after an
exposure, infection, or vaccination, and then eventually decreasing (or
waning) over time. Thus individuals may not have detectable antibodies if
never (or very recently) exposed, infected, or vaccinated, or if their
antibody levels have waned to a level below the limit of detection (of the
assay being employed). To the extent that antibody levels are associated with
protection from infection with SARS-CoV-2 or COVID-19 disease (Earle et al.,,
2021; Khoury et al.,, 2021), seroprevalence estimates may be helpful in
modeling the fraction of a population which may be immune or less susceptible
to Covid-19. Likewise, cross-sectional seroprevalence estimates, combined with
certain modeling assumptions and other data, may permit inference about other
parameters such as the cumulative incidence of SARS-CoV-2 infection, infection
fatality rate, or attack rate (Takahashi et al.,, 2021; Shioda et al.,, 2021;
Buss et al.,, 2021; Perez-Saez et al.,, 2021; Brazeau et al.,, 2020).

Unfortunately, seroprevalence studies often 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 a particular antibody response, such as that of immunoglobulin G,
M, or A (IgG, IgM, or IgA), which is dichotomized into a positive or negative
classification based on a cutoff value. This dichotomization almost always
produces misclassification bias 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 error can lead to biased
seroprevalence estimates. 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, it is straightforward to show that naively using
the sample proportion of positive test results as a seroprevalence estimator
would, in expectation, lead to a seroprevalence estimate of nearly 2% rather
than 1%.

To appropriately account for measurement error, sensitivity and specificity
should be estimated and incorporated into seroprevalence estimators. 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 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 prevalence 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 representative by design and lead to less biased estimates than
convenience samples analyzed with post-hoc statistical adjustments (Shook-Sa
et al.,, 2020; Accorsi et al.,, 2021). However, probability-based sampling may
not always be feasible. Convenience sample based estimators often rely on the
assumption that each person in a covariate-defined stratum has an equal
probability of being in the sample (Elliott and Valliant,, 2017). Under this
assumption, population prevalence can be estimated with direct
standardization, though complex survey methods like calibration have been used
as well in the SARS-CoV-2 setting (e.g., Bajema et al.,, 2021).

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 Appendices.

##  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 $Y$, with $Y=1$ if the individual has antibodies against SARS-CoV-2
and $Y=0$ otherwise. Our goal is to draw inference about the population
seroprevalence $\pi=P(Y=1)$. Because of error in the serology assay, $Y$ is
not observed directly. Let the result of the serology assay be denoted by $X$,
with $X=1$ if the individual tests positive (according to the antibody assay
used) and $X=0$ otherwise. Three key quantities are sensitivity, the
probability that a true positive tests positive, denoted by
$\sigma_{e}=P(X=1\mid Y=1)$; specificity, the probability that a true negative
tests negative, denoted by $\sigma_{p}=P(X=0\mid Y=0)$; and the population
expectation of the serology assay outcome, denoted by
$\rho=\mathbb{E}(X)=P(X=1)$. Unless the assay has perfect sensitivity and
specificity with $\sigma_{e}=\sigma_{p}=1$, $\rho$ typically will not equal
$\pi$ and $X$ will be a misclassified version of $Y$.

Suppose validation studies are conducted to estimate the sensitivity and
specificity by applying the assay to known positives and negatives.
Specifically, measurements are taken on $n_{1}$ independent and identically
distributed (iid) units from strata of the population where $Y=1$ and on
$n_{2}$ iid units from strata where $Y=0$. Thus $n_{1}$ copies of $X$ are
observed to estimate sensitivity and $n_{2}$ copies of $X$ are observed to
estimate specificity. To estimate seroprevalence in a target population, a
‘main’ study with $n_{3}$ iid copies of $X$ 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 $\delta_{i}$ be an
indicator of which study the $i$th individual’s sample $X_{i}$ is from, with
$\delta_{i}=1$ for the sensitivity study, $\delta_{i}=2$ for the specificity
study, and $\delta_{i}=3$ for the main study. Note that $\sum
I(\delta_{i}=j)=n_{j}$ for $j=1,2,3$, where $n=n_{1}+n_{2}+n_{3}$ and here and
throughout summations are taken from $i=1$ to $n$ unless otherwise specified.
Assume $n_{j}/n\to c_{j}\in(0,1)$ as $n\to\infty$.

###  2.2 Estimators and statistical properties

Let $\theta=(\sigma_{e},\sigma_{p},\rho,\pi)^{T}$. Consider the estimator
$\hat{\theta}=(\hat{\sigma}_{e},\hat{\sigma}_{p},\hat{\rho},\hat{\pi}_{RG})^{T}$,
where $\hat{\sigma}_{e}=n_{1}^{-1}\sum I(\delta_{i}=1)X_{i}$,
$\hat{\sigma}_{p}=n_{2}^{-1}\sum I(\delta_{i}=2)(1-X_{i})$,
$\hat{\rho}=n_{3}^{-1}\sum I(\delta_{i}=3)X_{i}$, and
$\hat{\pi}_{RG}=(\hat{\rho}+\hat{\sigma}_{p}-1)/(\hat{\sigma}_{e}+\hat{\sigma}_{p}-1)$.
The prevalence estimator $\hat{\pi}_{RG}$ is motivated by rearranging the
identity that $\rho=\pi\sigma_{e}+(1-\pi)(1-\sigma_{p})$ and is sometimes
referred to as the Rogan-Gladen (1978) estimator.

The estimator $\hat{\theta}$ can be expressed as the solution (for $\theta$)
to the estimating equation vector

$\sum\psi(X_{i};\delta_{i},\theta)=\left(\begin{array}[]{l}\sum\psi_{e}(X_{i};\delta_{i},\theta)\\\
\sum\psi_{p}(X_{i};\delta_{i},\theta)\\\
\sum\psi_{\rho}(X_{i};\delta_{i},\theta)\\\
\psi_{\pi}(X_{i};\delta_{i},\theta)\end{array}\right)=\left(\begin{array}[]{l}\sum
I(\delta_{i}=1)(X_{i}-\sigma_{e})\\\ \sum
I(\delta_{i}=2)\\{(1-X_{i})-\sigma_{p}\\}\\\ \sum
I(\delta_{i}=3)(X_{i}-\rho)\\\
(\rho+\sigma_{p}-1)-\pi(\sigma_{e}+\sigma_{p}-1)\end{array}\right)=0$

where here and below 0 denotes a column vector of zeros. Since the samples
were selected from three different populations, the data $X_{1},\dots,X_{n}$
are not identically distributed and care must be taken to derive the large
sample properties of $\hat{\theta}$. In Appendix A, the estimator
$\hat{\theta}$ is shown to be consistent and asymptotically normal.
Specifically, as $n\to\infty$,
$\sqrt{n}(\hat{\theta}-\theta)\to_{d}\mathcal{N}\left(0,\mathbb{A}(\theta)^{-1}\mathbb{B}(\theta)\mathbb{A}(\theta)^{-T}\right)$
and $\sqrt{n}(\hat{\pi}-\pi)\to_{d}\mathcal{N}\left(0,V_{\pi,RG}\right)$
assuming $\sigma_{e}>1-\sigma_{p}$ (as discussed below), where
$\mathbb{A}(\theta)^{-1}\mathbb{B}(\theta)\mathbb{A}(\theta)^{-T}$ is a
covariance matrix with bottom right element

$V_{\pi,RG}=\left\\{\frac{\pi^{2}\sigma_{e}(1-\sigma_{e})}{c_{1}}+\frac{(1-\pi)^{2}\sigma_{p}(1-\sigma_{p})}{c_{2}}+\frac{\rho(1-{\rho})}{c_{3}}\right\\}(\sigma_{e}+\sigma_{p}-1)^{-2}.$
(1)

The proof of consistency and asymptotic normality is similar to those from
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. Note that the asymptotic variance (1) consists of three
components corresponding to the sensitivity, specificity, and main studies. In
some circumstances, investigators may be able to decrease the variance of
$\hat{\pi}_{RG}$ by increasing the sample sizes of the sensitivity or
specificity studies compared to the main study (Larremore et al.,, 2020).

Let $\hat{V}_{\pi,RG}$ denote the plug-in estimator defined by replacing
$\sigma_{e},\sigma_{p},\rho,\pi$, and $c_{j}$ in (1) with
$\hat{\sigma}_{e},\hat{\sigma}_{p},\hat{\rho},\hat{\pi}_{RG}$, and $n_{j}/n$
for $j=1,2,3$, and note that $\hat{V}_{\pi,RG}/n$ is the variance estimator
proposed by Rogan and Gladen, (1978). By the continuous mapping theorem,
$\hat{V}_{\pi,RG}$ is consistent for the asymptotic variance assuming
$\sigma_{e}>1-\sigma_{p}$ and can be used to construct Wald-type CIs that
asymptotically attain nominal coverage probabilities.

###  2.3 Truncation into $[0,1]$

In finite samples $\hat{\pi}_{RG}$ sometimes yields estimates outside of
$[0,1]$ when (i) $\hat{\sigma}_{e}<1-\hat{\sigma}_{p}$, (ii)
$\hat{\rho}<1-\hat{\sigma}_{p}$, or (iii) $\hat{\rho}>\hat{\sigma}_{e}$.
Indeed, (ii) occurred in the ScreenNC study discussed in Section 5.2.
Estimates are typically truncated to be inside $[0,1]$ because the true
population prevalence must exist in $[0,1]$ (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
$\hat{\pi}_{RG}$ yields estimates inside $[0,1]$ almost surely unless
$\sigma_{e}<1-\sigma_{p}$. In practice, settings where
$\sigma_{e}<1-\sigma_{p}$ may be very unlikely; in such scenarios, the
probability of a positive test result is higher for seronegative persons than
for seropositive persons, so such a measurement instrument performs worse in
expectation than random guessing. Throughout this manuscript, we assume
$\sigma_{e}>1-\sigma_{p}$.

##  3 Standardized seroprevalence estimation

###  3.1 Problem setup

In some settings it may not be reasonable to assume the $n_{3}$ copies of $X$
from the main study constitute a random sample from the target population.
Suppose instead that for each copy of $X$ a vector of discrete covariates $Z$
is observed, with $Z$ taking on $k$ possible values $z_{1},\ldots,z_{k}$. The
covariates $Z$ are of interest because seroprevalence may differ between the
strata; for instance, $Z$ might include demographic variables such as age
group, race, or gender. Denote the mean of $X$ in the $j$th stratum as
$\rho_{j}=P(X=1\mid Z=z_{j})$ and the sample size for the $j$th stratum as
$n_{z_{j}}=\sum I(\delta_{i}=3,Z_{i}=z_{j})$, so
$\sum_{j=1}^{k}n_{z_{j}}=n_{3}$.

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, Chapter 15)).
Denote the proportion of the target population comprised by the $j$th stratum
as $\gamma_{j}=P(Z=z_{j})$ and suppose that these stratum proportions are
known with each $\gamma_{j}>0$ and $\sum_{j=1}^{k}\gamma_{j}=1$. For instance,
if the target population were all of the adults in a US state, then
$\gamma_{1},\dots,\gamma_{k}$ could be obtained from US census data. Assume
that all persons in a covariate stratum defined by $Z$ have the same
probability of inclusion in the sample. Then the covariates $Z$ in the main
study sample have a multinomial distribution with $k$ categories, sample size
$n_{3}$, and an unknown sampling probability vector $(s_{1},\dots,s_{k})^{T}$
where $\sum_{j=1}^{k}s_{j}=1$. 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, that is, $s_{j}=\gamma_{j}$ for
$j=1,\dots,k$.

###  3.2 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 $\rho$ is a weighted average of the
stratum-conditional means $\rho_{j}$ where each weight is a known stratum
proportion $\gamma_{j}$, i.e., $\rho=\sum_{j=1}^{k}\rho_{j}\gamma_{j}$. A
nonparametric standardization estimator for $\rho$ using the sample stratum-
conditional prevalences $\hat{\rho}_{j}=n_{z_{j}}^{-1}\sum
I(Z_{i}=z_{j},\delta_{i}=3)X_{i}$ for $j=1,\dots,k$ is
$\hat{\rho}_{SRG}=\sum_{j=1}^{k}\hat{\rho}_{j}\gamma_{j}$. A standardized
prevalence estimator accounting for measurement error is
$\hat{\pi}_{SRG}=(\hat{\rho}_{SRG}+\hat{\sigma}_{p}-1)/(\hat{\sigma}_{e}+\hat{\sigma}_{p}-1)$,
which has been used in SARS-CoV-2 seroprevalence studies (Havers et al.,,
2020; Barzin et al.,, 2020; Cai et al.,, 2020).

Letting
$\theta_{s}=(\sigma_{e},\sigma_{p},\rho_{1},\dots,\rho_{k},\rho,\pi)^{T}$, the
estimator
$\hat{\theta}_{s}=(\hat{\sigma}_{e},\hat{\sigma}_{p},\hat{\rho}_{1},\dots,\hat{\rho}_{k},\hat{\rho}_{SRG},\hat{\pi}_{SRG})^{T}$
solves the estimating equation vector
$\sum\psi(X_{i},Z_{i};\delta_{i},\theta_{s})=\left(\sum\psi_{e},\sum\psi_{p},\sum\boldsymbol{\psi}_{\boldsymbol{\rho}},\psi_{\rho},\psi_{\pi}\right)^{T}=0$
where $\sum\psi_{e},\sum\psi_{p}$, and $\psi_{\pi}$ are defined in Section 2;
$\sum\boldsymbol{\psi}_{\boldsymbol{\rho}}$ is a $k$-vector with $j$th element
$\sum\psi_{\rho_{j}}=\sum I(Z_{i}=z_{j},\delta_{i}=3)(X_{i}-\rho_{j})$; and
$\psi_{\rho}=\sum_{j=1}^{k}\rho_{j}\gamma_{j}-\rho$. It follows that
$\hat{\theta}_{s}$ is consistent and asymptotically normal and that
$\sqrt{n}(\hat{\pi}_{SRG}-\pi)\to_{d}\mathcal{N}\left(0,V_{\pi,SRG}\right)$
where

$V_{\pi,SRG}=\left\\{\frac{\pi^{2}\sigma_{e}(1-\sigma_{e})}{c_{1}}+\frac{(1-\pi)^{2}\sigma_{p}(1-\sigma_{p})}{c_{2}}+\sum_{j=1}^{k}\frac{\gamma_{j}^{2}\rho_{j}(1-{\rho_{j}})}{c_{3}s_{j}}\right\\}(\sigma_{e}+\sigma_{p}-1)^{-2}.$
(2)

The asymptotic variance $V_{\pi,SRG}$ can be consistently estimated by the
plug-in estimator $\hat{V}_{\pi,SRG}$ defined by replacing
$\sigma_{e},\sigma_{p},\rho_{j},s_{j},\pi$, and $c_{l}$ in (2) with
$\hat{\sigma}_{e},\hat{\sigma}_{p},\hat{\rho}_{j},n_{z_{j}}/n_{3},\hat{\pi}_{SRG}$,
and $n_{l}/n$ for $j=1,\dots,k$ and $l=1,2,3$. Consistency of
$\hat{V}_{\pi,SRG}$ 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 $X$,
$\rho_{j}=P(X=1\mid Z=z_{j})$. However, when $n_{z_{j}}=0$ for some strata
$j$, the corresponding estimator $\hat{\rho}_{j}$ is undefined, and
$\hat{\rho}_{SRG}$ is then undefined as well. Values of $n_{z_{j}}$ may equal
zero for two reasons. First, the study design may exclude these strata
($s_{j}=0$), a situation referred to as deterministic or structural
nonpositivity (Westreich and Cole, 2010). Second, even if $s_{j}>0$, random
nonpositivity can occur if no individuals with $Z=z_{j}$ are sampled, which
may occur if $s_{j}$ is small or if $n_{3}$ 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 $j$ for which $n_{z_{j}}>0$. However, this
redefined target population may be less relevant from a public health or
policy perspective.

###  3.3 Parametric standardization

We now consider a parametric approach for modeling the stratum-conditional
means $\rho_{j}$, 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 $g(\rho_{j})=\beta
h(z_{j})$ holds where $g$ is an appropriate link function for a binary outcome
like the logit or probit function; $\beta$ is a row vector of $p$ regression
coefficients with intercept $\beta_{1}$; and $h(z_{j})$ is a user-specified
$p$-vector function of the $j$th stratum’s covariate values that may include
main effects and interaction terms, with $l$th element denoted $h_{l}(z_{j})$
and $h_{1}(z_{j})$ set equal to one to correspond to an intercept. Let
$\operatorname{supp}(z)$ be the covariate support in the sample, i.e.,
$\operatorname{supp}(z)=\\{z_{j}:n_{z_{j}}>0\\}$ with dimension
$\dim\\{\operatorname{supp}(z)\\}=\sum_{j=1}^{k}I(n_{z_{j}}>0)$, and assume
$p\leq\dim\\{\operatorname{supp}(z)\\}\leq k$. (Note that
$\dim\\{\operatorname{supp}(z)\\}=k$ only when there is positivity, and in
that case $\hat{\pi}_{SRG}$ can be used with no restriction needed.)

Under the assumed binary regression model, each $\rho_{j}$ is a function of
the parameters $\beta$ and the covariates $z_{j}$ that define the $j$th
stratum, denoted $\rho_{j}(\beta,z_{j})=g^{-1}\\{\beta h(z_{j})\\}$. A model-
based standardized Rogan-Gladen estimator of $\pi$ is
$\hat{\pi}_{SRGM}=(\hat{\rho}_{SRGM}+\hat{\sigma}_{p}-1)/(\hat{\sigma}_{e}+\hat{\sigma}_{p}-1)$,
where
$\hat{\rho}_{SRGM}=\sum_{j=1}^{k}\hat{\rho}_{j}(\hat{\beta},z_{j})\gamma_{j}$
and $\hat{\beta}$ is the maximum likelihood estimator of $\beta$. Estimating
equation theory can again be used to derive large-sample properties by
replacing the $k$ equations for $\rho_{1},\dots,\rho_{k}$ from Section 3.2
with $p$ equations for $\beta_{1},\dots,\beta_{p}$ corresponding to the score
equations from the binary regression.

Let
$\theta_{m}=(\sigma_{e},\sigma_{p},\beta_{1},\dots,\beta_{p},\rho,\pi)^{T}$.
The estimator
$\hat{\theta}_{m}=(\hat{\sigma}_{e},\hat{\sigma}_{p},\hat{\beta}_{1},\dots,\hat{\beta}_{p},\hat{\rho}_{SRGM},\hat{\pi}_{SRGM})^{T}$
is the solution to the vector
$\sum\psi(X_{i},Z_{i};\delta_{i},\theta_{m})=\left(\sum\psi_{e},\sum\psi_{p},\sum\psi_{\beta},\psi_{\rho},\psi_{\pi}\right)^{T}=0$
of estimating equations where $\sum\psi_{e},\sum\psi_{p}$, and $\psi_{\pi}$
are as in Section 2; $\sum\psi_{\beta}$ is a $p$-vector with $j$th element
$\sum\psi_{\beta_{j}}=\sum I(\delta_{i}=3)\left[X_{i}-g^{-1}\left\\{\beta
h(Z_{i})\right\\}\right]h_{j}(Z_{i})$; and
$\psi_{\rho}=\sum_{j=1}^{k}g^{-1}\\{\beta h(Z_{j})\\}\gamma_{j}-\rho$. It
follows that $\hat{\theta}_{m}$ is consistent and asymptotically normal and
$\sqrt{n}(\hat{\pi}_{SRGM}-\pi)\to_{d}\mathcal{N}\left(0,V_{\pi,SRGM}\right)$.
The asymptotic variance $V_{\pi,SRGM}$ can be consistently estimated by
$\hat{V}_{\pi,SRGM}$, the lower right element of the empirical sandwich
variance estimator of the asymptotic variance of $\hat{\theta}_{m}$. 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 $\hat{\pi}_{RG}$,
$\hat{\pi}_{SRG}$, and $\hat{\pi}_{SRGM}$. Four data generating processes
(DGPs) were considered, within which different scenarios were defined through
full factorial designs that varied simulation parameters
$\pi,\sigma_{e},\sigma_{p},n_{1},n_{2}$, and $n_{3}$. 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 $n_{1}$ and $n_{2}$ were generated with $X$
distributed Bernoulli with a mean of $\sigma_{e}$ or $1-\sigma_{p}$,
respectively. In DGPs 1 and 2 a main study of size $n_{3}$ was then generated
where $Y$ was Bernoulli with mean $\pi$ and $X\mid Y$ was Bernoulli with mean
$\sigma_{e}Y+(1-\sigma_{p})(1-Y)$; in DGPs 3 and 4 $X$ was generated directly
from the distribution of $X\mid Z$, as will be explained. Simulation parameter
values were selected based on the seroprevalence studies described in Section
5. Sensitivity was varied in $\sigma_{e}\in\\{.8,.99\\}$, specificity in
$\sigma_{p}\in\\{.8,.95,.99\\}$, and prevalence in
$\pi\in\\{.01,.02,\dots,.20\\}$. Sample sizes were $n_{1}=40$, $n_{2}=250$,
and $n_{3}=2500$. 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
$\hat{\pi}-\pi$ for each estimator $\hat{\pi}$; and (b) empirical coverage,
i.e., whether the 95% Wald-type CIs based on each variance estimator
$\hat{V}_{\pi}$ contained the true prevalence.

###  4.1 No selection bias

Simulations for DGP 1 were conducted to assess the performance of
$\hat{\pi}_{RG}$ when no selection bias was present. The estimator
$\hat{\pi}_{RG}$ was generally unbiased, as seen in Appendix Figure A.1.
Performance improved as $\sigma_{e}$ and $\sigma_{p}$ tended toward 1, with
specificity $\sigma_{p}$ being a stronger determinant of bias. An exception to
these results was for low prevalence $\pi$ (0.05 or lower), when
$\hat{\pi}_{RG}$ overestimated the true prevalence. Wald CIs based on
$\hat{V}_{\pi,RG}$ generally attained nominal coverage, as seen in Appendix
Figure A.2. However, when prevalence $\pi$ was near zero, coverage of the 95%
CIs did not equal the nominal level. Coverage also decreased as $\sigma_{p}$
increased toward 1. These variable CI coverage results concord with previous
simulation studies evaluating $\hat{V}_{\pi,RG}$ (Lang and Reiczigel,, 2014).
Alternative methods for constructing CIs are discussed in Section 6.

###  4.2 Low-dimensional selection bias

In DGP 2, the target population was comprised of two strata defined by a
covariate $Z\in\\{z_{1},z_{2}\\}$ with proportions $\gamma_{1}=\gamma_{2}=.5$.
Within the main study, $Z$ was generated from a multinomial distribution of
sample size $n_{3}$ and sampling probability vector $(.2,.8)$. Individuals’
serostatuses were generated from the conditional distribution $Y\mid Z$, which
was such that $P(Y=1\mid Z=z_{1})=1.5\pi$ and $P(Y=1\mid Z=z_{2})=0.5\pi$ for
each value of $\pi$. In each simulation $\hat{\pi}_{RG}$ and $\hat{\pi}_{SRG}$
and their corresponding 95% CIs were computed.

The nonparametric standardized estimator $\hat{\pi}_{SRG}$ was empirically
unbiased for true prevalences $\pi\geq 0.05$, as seen in Appendix Figure A.3,
and 95% CIs based on $\hat{V}_{\pi,SRG}$ attained approximately nominal
coverage, as seen in Appendix Figure 4 A.4. As with $\hat{\pi}_{RG}$ in DGP 1,
CI coverage fell slightly below the nominal level for very low prevalences,
and more noticeably for $\sigma_{p}=.99$. Appendix Figures A.3 and A.4 show
that $\hat{\pi}_{RG}$ 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 $\hat{\pi}_{SRG}$ and $\hat{\pi}_{SRGM}$ in scenarios
with larger numbers of strata.

####  4.3.1 DGP 3

Three covariates were defined as $Z_{1}\in\\{z_{10},z_{11}\\}$,
$Z_{2}\in\\{z_{20},z_{21},z_{22},z_{23}\\}$, and
$Z_{3}\in\\{z_{30},z_{31},z_{32},z_{33},z_{34}\\}$, leading to $k=40$ strata
with proportions $(\gamma_{1},\dots,\gamma_{40})$. Within the main study, $Z$
was again generated as multinomial with size $n_{3}$ 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 $X\mid Z$, where

$\begin{split}\operatorname{logit}\\{P(X=1\mid
Z)\\}&=\beta_{0}+\beta_{1}I(Z_{1}=z_{11})+\beta_{2}I(Z_{2}=z_{20})+\beta_{3}I(Z_{2}=z_{21})\\\
&+\beta_{4}I(Z_{3}=z_{30})+\beta_{5}I(Z_{3}=z_{31}).\end{split}$

The parameters $\beta_{1}=-1,\ \beta_{2}=-.6,\ \beta_{3}=.8,\ \beta_{4}=.6,$
and $\beta_{5}=.4$ were set to reflect differential prevalences by stratum,
while a “balancing intercept” $\beta_{0}$ (Rudolph et al.,, 2021) was set to
different values so that prevalence $\pi$ equaled (approximately)
$\\{.01,.02,\dots,.20\\}$. The nonparametric estimator $\hat{\pi}_{SRG}$ 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 $\hat{\pi}_{SRGM}$ 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 $\gamma_{j}$ (stratum
proportion) and $s_{j}$ (sampling probability).

Both $\hat{\pi}_{SRG}$ and $\hat{\pi}_{SRGM}$ performed well in this scenario.
Figure 2 shows that the estimators were generally empirically unbiased, though
modest bias occurred when $\sigma_{p}=0.8$ and prevalence was low. Appendix
Figure A.5 shows 95% CIs based on either $\hat{V}_{\pi,SRG}$ or
$\hat{V}_{\pi,SRGM}$ attained nominal coverage, with slight under-coverage for
$\pi<0.05$. On average across all 120 scenarios, 89% (range of 86%-92%) of
simulated datasets had positivity. As in DGP 2, $\hat{\pi}_{RG}$ was seriously
biased and the corresponding CIs did not attain nominal coverage.

Figure 2: Empirical bias of the Rogan-Gladen ($\hat{\pi}_{RG}$), nonparametric
standardized ($\hat{\pi}_{SRG}$), and logistic regression standardized
($\hat{\pi}_{SRGM}$) estimators from simulation study for DGP 3, described in
Section 4.3.1. The six facets correspond to a given combination of sensitivity
$\sigma_{e}$ (‘Sens’) and specificity $\sigma_{p}$ (‘Spec’).

####  4.3.2 DGP 4

Data were generated as in DGP 3, but the inclusion of a fourth covariate
$Z_{4}\in\\{z_{40},z_{41}\\}$ now led to 80 strata. The conditional
distribution $X\mid Z$ was such that $\operatorname{logit}\\{P(X=1\mid
Z)\\}=\nu h(Z)$, where $h(Z)$ here contains the same terms as in DGP 3 plus a
main effect for $I(Z_{4}=z_{41})$ with corresponding coefficient $\nu_{6}$.
Regression parameters were a balancing intercept $\nu_{0}$, $\nu_{1}=-1$,
$\nu_{2}=3.25$, $\nu_{3}=.8$, $\nu_{4}=.6$, $\nu_{5}=.4$, and $\nu_{6}=.1$.
The larger value for $\nu_{2}$, as compared to $\beta_{2}$, led to a stronger
relationship between $X$ and $Z$ 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 $\hat{\pi}_{SRGM}$ was empirically unbiased in this
scenario, while $\hat{\pi}_{SRG}$ 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
$\hat{\pi}_{SRG}$ may be explained by restriction leading to bias under
nonpositivity. CIs based on $\hat{V}_{\pi,SRGM}$ typically attained nominal
coverage, unlike those based on $\hat{V}_{\pi,SRG}$ or $\hat{V}_{\pi,RG}$, as
seen in Appendix Figure A.6. However, the trend of empirical coverage being
below the nominal level for low prevalence and $\sigma_{p}=.99$ was more
noticeable compared to other DGPs, perhaps due to greater selection bias. Note
that $\hat{V}_{\pi,SRGM}$ 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
$\hat{\pi}_{SRG}$ and $\hat{\pi}_{SRGM}$ 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 $\hat{\pi}_{SRGM}$ generally
maintained its performance while $\hat{\pi}_{SRG}$ had greater empirical bias
and the corresponding CIs did not attain nominal coverage levels.

###  4.4 Model misspecification

The performance of $\hat{\pi}_{SRGM}$ was also assessed in scenarios similar
to DGPs 3 and 4, but under mild model misspecification. Here, the true
conditional distributions of $Y\mid Z$ were $\operatorname{logit}\\{P(Y=1\mid
Z)\\}=\beta h(Z)$ and $\operatorname{logit}\\{P(Y=1\mid Z)\\}=\nu h(Z)$, where
$\beta h(Z)$ and $\nu h(Z)$ are the specifications used in the models for
$P(X=1\mid Z)$ in DGPs 3 and 4, respectively. The test results $X$ were
generated from $X\mid Y$ as in DGPs 1 and 2. The degree of misspecification in
each scenario thus depended on the values of $\sigma_{e}$, $\sigma_{p}$, and
$\pi$. Appendix Figures A.7 through A.10 show that, in terms of bias and 95%
CI coverage, $\hat{\pi}_{SRGM}$ 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). The final
collection round took place before the first vaccine authorization in the
European Union in December 2020. 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 $n_{1}=181$ RT-
PCR confirmed COVID-19 cases and $n_{2}=326$ pre-pandemic negative controls,
sensitivity and specificity were estimated to be $\hat{\sigma}_{e}=.851$ and
$\hat{\sigma}_{p}=.988$ (Herzog et al.,, 2021, Table S1.1). The number of
samples for assessing seroprevalence varied between $n_{3}=$ 2,960 and
$n_{3}=$ 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,
(2021). 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 $\hat{\pi}_{SRG}$.
For $\hat{\pi}_{SRGM}$, 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 $\hat{\pi}_{RG}$,
$\hat{\pi}_{SRG}$, and $\hat{\pi}_{SRGM}$ alongside those for the unadjusted,
or naive, sample prevalence $\hat{\rho}$ for each collection round (with exact
Clopper-Pearson 95% CIs). The naive estimates $\hat{\rho}$ were typically
greater than the other three estimates and had narrower CIs. The greatest
differences were between $\hat{\rho}$ and $\hat{\pi}_{RG}$, which can be
attributed to (estimated) measurement error in the assay. Both of the
standardized estimates, $\hat{\pi}_{SRG}$ and $\hat{\pi}_{SRGM},$ were similar
in value to $\hat{\pi}_{RG}$ in most collection periods. These estimates, in
combination with the stratified random sampling design, suggest that the
magnitude of measurement error in this study 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.

Time-specific seropositivity was at its highest level in the third collection
period in May 2020, and thereafter followed a general decreasing trend through
September 2020. This trend could reflect a decrease in seroprevalence in
Belgium during this calendar time period, perhaps due to a decrease in
seroincidence and waning antibodies (Herzog et al.,, 2021). Alternatively, the
apparent trend could reflect changes in care seeking behaviors during the
study timeframe which could have affected the composition of the population
contributing residual sera in each collection period (Herzog et al.,, 2021),
although the standardized estimates should at least partially account for
changes in the composition of individuals contributing sera over time.

###  5.2 North Carolina seroprevalence study

The standardization methods of Section 3 were also applied to ScreenNC, which
tested a convenience sample of $n_{3}=$ 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), before the authorization of vaccines in the
United States. 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 $n_{1}=40$ RT-PCR
confirmed positive patients and $n_{2}=277$ pre-pandemic serum samples assumed
to be negative, sensitivity was estimated as $\hat{\sigma}_{e}=1$ and
specificity as $\hat{\sigma}_{p}=0.989$.

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 (US 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 $\hat{\pi}_{SRG}$. Logistic
regression models with main effects for sex, race, and age group were used to
compute $\hat{\pi}_{SRGM}$; interaction effects were not included as the small
number of positive test results could have led to model overfit.

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 $n$. Some column totals do not sum to 100% because of rounding.  |  | ScreenNC | UNC Hospitals | NC ACS
---|---|---|---|---
|  | $n$ | % | $n$ | % | $n$ | %
|  | 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 African-American | 395 | 13 | 5,109 | 23 | 1,640,311 | 21
Other | 311 | 10 | 1,843 | 6 | 455,600 | 6
White or Caucasian | 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

The sample proportion of positive tests was $\hat{\rho}=24/2973=0.81\%$. The
sample false positive rate was $1-\hat{\sigma}_{p}=1.08\%$, so the data are,
at first appearance, consistent with a population prevalence of 0%. Indeed,
the Rogan-Gladen seroprevalence estimate was $\hat{\pi}_{RG}=0$% (95% CI 0%,
1.00%). Likewise, the UNC target population had nonparametric and parametric
standardized estimates of $\hat{\pi}_{SRG}=0$% (0%, 1.11%) and
$\hat{\pi}_{SRGM}=0$% (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 $[0,1]$. 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

We examined nonparametric and model-based standardized Rogan-Gladen
estimators, deriving their large-sample properties and consistent variance
estimators. While motivated by SARS-CoV-2 seroprevalence studies, the methods
considered here are also applicable to prevalence estimation of any binary
variable for settings where validation data can be used to estimate the
measurement instrument’s sensitivity and specificity and covariate data can be
used for standardization. Simulation studies demonstrated that both
standardized Rogan-Gladen 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
$\hat{\pi}_{SRGM}$ was empirically unbiased even when the number of strata and
covariates, and with them the potential for random nonpositivity, increased. A
drawback to $\hat{\pi}_{SRGM}$ is the need to correctly specify the form of a
regression model. On the other hand, the nonparametric standardized estimator
$\hat{\pi}_{SRG}$ 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, $\hat{\pi}_{SRG}$ 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. In this paper Wald-type confidence intervals are
considered, which have known limitations (Brown et al.,, 2001; Dean and
Pagano,, 2015); alternative types of confidence intervals could be considered
based on the bootstrap (Cai et al.,, 2020), Bayesian posterior intervals
(Gelman and Carpenter,, 2020), or test inversion (DiCiccio et al.,, 2021).
While the approaches here estimate seroprevalence at a fixed point in time,
seroprevalence is a dynamic parameter. For analysis of studies with lengthier
data collection periods, extensions of the estimators in this paper could be
considered which make additional assumptions (e.g., smoothness, monotonicity)
about the longitudinal nature of seroprevalence. Another possible extension
could consider variations in assay sensitivity, which may depend on a variety
of factors such as: the type of assay used; the recency of exposure,
infection, or vaccination of an individual; disease severity in infected
individuals; the type and dose of vaccine for vaccinated individuals; and so
forth. Where additional data are available related to these factors, then
extensions of the standardized Rogan-Gladen estimators which incorporate these
additional data could be developed. As an alternative to standardization,
inverse probability of sampling weights (Lesko et al.,, 2017) or inverse odds
of sampling weights (Westreich and Cole,, 2010) could be considered.
Standardization and weighting methods may possibly be combined to create a
doubly robust Rogan-Gladen estimator.

## References

  * Accorsi et al., (2021) Accorsi, E. K., Qiu, X., Rumpler, E., Kennedy-Shaffer, L., Kahn, R., Joshi, K., Goldstein, E., Stensrud, M. J., Niehus, R., Cevik, M., and Lipsitch, M. (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.
  * Arora et al., (2021) Arora, R. K., Joseph, A., Van Wyk, J., Rocco, S., Atmaja, A., May, E., Yan, T., Bobrovitz, N., Chevrier, J., Cheng, M. P., Williamson, T., and Buckeridge, D. L. (2021).  SeroTracker: a global SARS-CoV-2 seroprevalence dashboard.  The Lancet Infectious Diseases, 21(4):e75–e76.
  * Bajema et al., (2021) Bajema, K. L., Wiegand, R. E., Cuffe, K., Patel, S. V., Iachan, R., Lim, T., Lee, A., Moyse, D., Havers, F. P., Harding, L., et al. (2021).  Estimated SARS-CoV-2 seroprevalence in the US as of September 2020\.  JAMA Internal Medicine, 181(4):450–460.
  * Barzin et al., (2020) Barzin, A., Schmitz, J. L., Rosin, S., Sirpal, R., Almond, M., Robinette, C., Wells, S., Hudgens, M., Olshan, A., Deen, S., et al. (2020).  SARS-CoV-2 seroprevalences among a southern U.S. population indicates limited asymptomatic spread under physical distancing measures.  mBio, 11(5):e02426–20.
  * Boos and Stefanski, (2013) Boos, D. D. and Stefanski, L. A. (2013).  Essential Statistical Inference: Theory and Methods.  Springer.
  * Bouman et al., (2021) Bouman, J. A., Riou, J., Bonhoeffer, S., and Regoes, R. R. (2021).  Estimating the cumulative incidence of SARS-CoV-2 with imperfect serological tests: Exploiting cutoff-free approaches.  PLOS Computational Biology, 17(2):e1008728.
  * Brazeau et al., (2020) Brazeau, N. F., Verity, R., Jenks, S., Fu, H., Whittaker, C., Winskill, P., Walker, P., Riley, S., Schnekenberg, R. P., Hoeltgebaum, H., et al. (2020).  Report 34: COVID-19 infection fatality ratio: Estimates from seroprevalence.  Technical report, Imperial College London.
  * Brown et al., (2001) Brown, L. D., Cai, T. T., and DasGupta, A. (2001).  Interval estimation for a binomial proportion.  Statistical Science, 16(2):101–133.
  * Buss et al., (2021) Buss, L. F., Prete Jr., C. A., Abrahim, C. M. M., Mendrone Jr., A., Salomon, T., Almeida-Neto, C. d., França, R. F. O., Belotti, M. C., Carvalho, M. P. S. S., Costa, A. G., et al. (2021).  Three-quarters attack rate of SARS-CoV-2 in the Brazilian Amazon during a largely unmitigated epidemic.  Science, 371(6526):288–292.
  * Cai et al., (2020) Cai, B., Ioannidis, J. P. A., Bendavid, E., and Tian, L. (2020).  Exact inference for disease prevalence based on a test with unknown specificity and sensitivity.  arXiv. http://arxiv.org/abs/2011.14423.
  * Dean and Pagano, (2015) Dean, N. and Pagano, M. (2015).  Evaluating confidence interval methods for binomial proportions in clustered surveys.  Journal of Survey Statistics and Methodology, 3(4):484–503.
  * DiCiccio et al., (2021) DiCiccio, T. J., Ritzwoller, D. M., Romano, J. P., and Shaikh, A. M. (2021).  Confidence intervals for seroprevalence.  arXiv. http://arxiv.org/abs/2103.15018.
  * Earle et al., (2021) Earle, K. A., Ambrosino, D. M., Fiore-Gartland, A., Goldblatt, D., Gilbert, P. B., Siber, G. R., Dull, P., and Plotkin, S. A. (2021).  Evidence for antibody as a protective correlate for COVID-19 vaccines.  Vaccine, 39:4423–4428.
  * Elliott and Valliant, (2017) Elliott, M. R. and Valliant, R. (2017).  Inference for nonprobability samples.  Statistical Science, 32(2):249–264.
  * Federal Planning Bureau, (2021) Federal Planning Bureau (2021).  Population projections 2020-2070.  Dataset. https://www.plan.be/databases/data-35-en-population_projections_2020_2070.
  * Gelman and Carpenter, (2020) Gelman, A. and 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.
  * Havers et al., (2020) Havers, F. P., Reed, C., Lim, T., Montgomery, J. M., Klena, J. D., Hall, A. J., Fry, A. M., Cannon, D. L., Chiang, C.-F., Gibbons, A., 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.
  * Herzog et al., (2021) Herzog, S., Bie, J. D., Abrams, S., Wouters, I., Ekinci, E., Patteet, L., Coppens, A., Spiegeleer, S. D., Beutels, P., Damme, P. V., 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. https://www.medrxiv.org/content/10.1101/2020.06.08.20125179v6.
  * Hilden, (1979) Hilden, J. (1979).  A further comment on “Estimating prevalence from the results of a screening test”.  American Journal of Epidemiology, 109(6):721–722.
  * Khoury et al., (2021) Khoury, D. S., Cromer, D., Reynaldi, A., Schlub, T. E., Wheatley, A. K., Juno, J. A., Subbarao, K., Kent, S. J., Triccas, J. A., and Davenport, M. P. (2021).  Neutralizing antibody levels are highly predictive of immune protection from symptomatic SARS-CoV-2 infection.  Nature Medicine, 27:1205–1211.
  * Kolenikov and Bollen, (2012) Kolenikov, S. and Bollen, K. A. (2012).  Testing negative error variances: Is a Heywood case a symptom of misspecification?  Sociological Methods & Research, 41(1):124–167.
  * Lang and Reiczigel, (2014) Lang, Z. and Reiczigel, J. (2014).  Confidence limits for prevalence of disease adjusted for estimated sensitivity and specificity.  Preventive Veterinary Medicine, 113(1):13–22.
  * Larremore et al., (2020) Larremore, D. B., Fosdick, B. K., Zhang, S., and Grad, Y. H. (2020).  Jointly modeling prevalence, sensitivity and specificity for optimal sample allocation.  bioRxiv. http://biorxiv.org/lookup/doi/10.1101/2020.05.23.112649.
  * Lesko et al., (2017) Lesko, C. R., Buchanan, A. L., Westreich, D., Edwards, J. K., Hudgens, M. G., and Cole, S. R. (2017).  Generalizing study results: a potential outcomes perspective.  Epidemiology, 28(4):553–561.
  * Mackey et al., (2021) Mackey, K., Ayers, C. K., Kondo, K. K., Saha, S., Advani, S. M., Young, S., Spencer, H., Rusek, M., Anderson, J., Veazie, S., et al. (2021).  Racial and ethnic disparities in COVID-19-Related infections, hospitalizations, and deaths.  Annals of Internal Medicine, 174(3):362–373.
  * Marchevsky, (1979) Marchevsky, N. (1979).  Re: “Estimating prevalence from the results of a screening test".  American Journal of Epidemiology, 109(6):720–721.
  * Perez-Saez et al., (2021) Perez-Saez, J., Zaballa, M.-E., Yerly, S., Andrey, D. O., Meyer, B., Eckerle, I., Balavoine, J.-F., Chappuis, F., Pittet, D., Trono, D., et al. (2021).  Persistence of anti-SARS-CoV-2 antibodies: immunoassay heterogeneity and implications for serosurveillance.  Clinical Microbiology and Infection, 27(11):1695.e7–1695.e12.
  * Pijls et al., (2021) Pijls, B. G., Jolani, S., Atherley, A., Derckx, R. T., Dijkstra, J. I. R., Franssen, G. H. L., Hendriks, S., Richters, A., Venemans-Jellema, A., Zalpuri, S., and Zeegers, M. P. (2021).  Demographic risk factors for COVID-19 infection, severity, ICU admission and death: a meta-analysis of 59 studies.  BMJ Open, 11(1):e044640.
  * Rogan and Gladen, (1978) Rogan, W. J. and Gladen, B. (1978).  Estimating prevalence from the results of a screening test.  American Journal of Epidemiology, 107(1):71–76.
  * Rudolph et al., (2021) Rudolph, J. E., Edwards, J. K., Naimi, A. I., and Westreich, D. J. (2021).  Simulation in practice: The balancing intercept.  American Journal of Epidemiology, 190(8):1696–1698.
  * Shioda et al., (2021) Shioda, K., Lau, M. S., Kraay, A. N., Nelson, K. N., Siegler, A. J., Sullivan, P. S., Collins, M. H., Weitz, J. S., and Lopman, B. A. (2021).  Estimating the cumulative incidence of SARS-CoV-2 infection and the infection fatality ratio in light of waning antibodies.  Epidemiology, 32(4):518–524.
  * Shook-Sa et al., (2020) Shook-Sa, B. E., Boyce, R. M., and 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.
  * Takahashi et al., (2021) Takahashi, S., Peluso, M. J., Hakim, J., Turcios, K., Janson, O., Routledge, I., Busch, M. P., Hoh, R., Tai, V., Kelly, J. D., et al. (2021).  SARS-CoV-2 serology across scales: a framework for unbiased seroprevalence estimation incorporating antibody kinetics and epidemic recency.  medRxiv. http://medrxiv.org/lookup/doi/10.1101/2021.09.09.21263139.
  * US Census Bureau, (2019) US Census Bureau (2019).  American Community Survey 1-year estimates, Public Use Microdata Sample.  Dataset. https://data.census.gov/mdat/#/.
  * van Belle et al., (2004) van Belle, G., Fisher, L. D., Heagerty, P. J., and Lumley, T. (2004).  Biostatistics: A Methodology for the Health Sciences.  John Wiley & Sons, 2nd edition.
  * Westreich and Cole, (2010) Westreich, D. and Cole, S. R. (2010).  Invited Commentary: Positivity in practice.  American Journal of Epidemiology, 171(6):674–677.

## Acknowledgments

We thank Dirk Dittmer and the ScreenNC research team for data access. We also
thank Shaina Alexandria, Bryan Blette, and Kayla Kilpatrick for constructive
suggestions and discussions. 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.

## 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
$\hat{\pi}_{SRG}$ and $\hat{V}_{\pi,SRG}$ are available at
https://github.com/samrosin/rgStandardized.

##  Appendix A Proofs for Section 2

###  A.1 Proof of asymptotic normality

Taylor expansion of $n^{-1}\sum\psi(X_{i};\delta_{i},\hat{\theta})$ around the
true parameter $\theta$ yields

$0=n^{-1}\sum\psi(X_{i};\delta_{i},\hat{\theta})=n^{-1}\sum\psi(X_{i};\delta_{i},\theta)+n^{-1}\sum\dot{\psi}(X_{i};\delta_{i},\theta)(\hat{\theta}-\theta)+R$

where
$\dot{\psi}(X_{i};\delta_{i},\theta)=\partial\psi(X_{i};\delta_{i},\theta)/\partial\theta^{T}$
and $R$ is a remainder term. Rearranging and multiplying by $\sqrt{n}$ yields

$\sqrt{n}(\hat{\theta}-\theta)=\left\\{n^{-1}\sum-\dot{\psi}(X_{i};\delta_{i},\theta)\right\\}^{-1}\sqrt{n}\left\\{n^{-1}\sum\psi(X_{i};\delta_{i},\theta)\right\\}+\sqrt{n}R^{*}$
(A.1)

where $R^{*}$ is a new remainder defined below. It is shown below that
$n^{-1}\sum-\dot{\psi}(X_{i};\delta_{i},\theta)\to_{p}\mathbb{A}(\theta)$,
$\sqrt{n}\\{n^{-1}\sum\psi(X_{i};\delta_{i},\theta)\\}\to_{d}\mathcal{N}\\{0,\mathbb{B}(\theta)\\}$,
and $\sqrt{n}R^{*}\to_{p}0$, where $\mathbb{A}(\theta)$ and
$\mathbb{B}(\theta)$ are defined below. Therefore, by Slutsky’s theorem,
$\sqrt{n}(\hat{\theta}-\theta)\to_{d}\mathcal{N}\left\\{0,\mathbb{A}(\theta)^{-1}\mathbb{B}(\theta)\mathbb{A}(\theta)^{-T}\right\\}$.

First, define

$\displaystyle\mathbb{A}_{n}(X,\delta,\theta)$
$\displaystyle=n^{-1}\sum-\dot{\psi}(X_{i};\delta_{i},\theta)$
$\displaystyle=n^{-1}\sum\left[\begin{array}[]{cccc}I(\delta_{i}=1)&0&0&0\\\
0&I(\delta_{i}=2)&0&0\\\ 0&0&I(\delta_{i}=3)&0\\\
\pi&-1+\pi&-1&\sigma_{e}+\sigma_{p}-1\end{array}\right]$

and let

$\displaystyle\mathbb{A}(\theta)=\left[\begin{array}[]{cccc}c_{1}&0&0&0\\\
0&c_{2}&0&0\\\ 0&0&c_{3}&0\\\
\pi&-1+\pi&-1&\sigma_{e}+\sigma_{p}-1\end{array}\right].$

As $n\to\infty$, $\mathbb{A}_{n}(X,\delta,\theta)\to\mathbb{A}(\theta)$ by the
assumption that $n^{-1}\sum I(\delta_{i}=j)=n_{j}/n\to c_{j}\in(0,1)$ for
$j\in\\{1,2,3\\}$.

Second, for brevity, let $\psi_{e}$ denote
$\psi_{e}(X_{i};\delta_{i},\theta_{i})$, and similarly for
$\psi_{p},\psi_{\rho}$, and $\psi_{\pi}$. Define

$\displaystyle\mathbb{B}_{n}(X,\delta,\theta)$
$\displaystyle=n^{-1}\sum\mathbb{E}\left\\{\psi(X_{i};\delta_{i},\theta)\psi(X_{i};\delta_{i},\theta)^{T}\right\\}$
$\displaystyle=n^{-1}\sum\mathbb{E}\left[\begin{array}[]{cccc}\psi_{e}^{2}&0&0&\psi_{e}\psi_{\pi}\\\
0&\psi_{p}^{2}&0&\psi_{p}\psi_{\pi}\\\
0&0&\psi_{\rho}^{2}&\psi_{\rho}\psi_{\pi}\\\
\psi_{\pi}\psi_{e}&\psi_{\pi}\psi_{p}&\psi_{\pi}\psi_{\rho}&\psi_{\pi}^{2}\end{array}\right]$

and let

$\mathbb{B}(\theta)=\left[\begin{array}[]{cccc}c_{1}\sigma_{e}(1-\sigma_{e})&0&0&0\\\
0&c_{2}\sigma_{p}(1-\sigma_{p})&0&0\\\ 0&0&c_{3}{\rho}(1-{\rho})&0\\\
0&0&0&0\end{array}\right].$

Note that
$\mathbb{B}(\theta)=\lim_{n\to\infty}\mathbb{B}_{n}(X,\delta,\theta)$ as
$\mathbb{E}(\psi_{e}^{2})=\mathbb{E}\\{I(\delta_{i}=1)(X_{i}-\sigma_{e})^{2}\\}=I(\delta_{i}=1)\sigma_{e}(1-\sigma_{e})$
and likewise for $\mathbb{E}(\psi_{p}^{2})$ and $\mathbb{E}(\psi_{\rho}^{2})$.
It follows that
$\sqrt{n}\left\\{n^{-1}\sum\psi(X_{i};\delta_{i},\theta)\right\\}\to_{d}\mathcal{N}\\{0,\mathbb{B}(\theta)\\}$
by the Lindeberg-Feller CLT. In particular, let
$s_{n}^{2}=\sum\\{\operatorname{Var}(\psi_{e})+\operatorname{Var}(\psi_{p})+\operatorname{Var}(\psi_{\rho})+\operatorname{Var}(\psi_{\pi})\\}=n_{1}\sigma_{e}(1-\sigma_{e})+n_{2}\sigma_{p}(1-\sigma_{p})+n_{3}\rho(1-\rho)$
and let $\lVert\cdot\rVert$ denote the Euclidean norm. Note that
$\max_{i}\lVert\psi(X_{i};\delta_{i},\theta)\rVert\leq\sqrt{2}$ as the maximum
magnitude of each element of $\psi(X_{i};\delta_{i},\theta)$ is 1, and for a
given observation $X_{i}$ it is always true that two of the indicators
$\delta_{1},\delta_{2},\delta_{3}$ equal zero and the third indicator equals
one. Therefore for all $\epsilon>0$,

$\lim_{n\to\infty}s_{n}^{-2}\sum\mathbb{E}\left\\{\lVert\psi(X_{i};\delta_{i},\theta)\rVert^{2}I(\lVert\psi(X_{i};\delta_{i},\theta)\rVert\geq\epsilon
s_{n})\right\\}\leq\lim_{n\to\infty}2s_{n}^{-2}\sum P(\sqrt{2}\geq\epsilon
s_{n})=0,$

implying the Lindeberg condition holds.

Third, it remains to prove $\sqrt{n}R^{*}\to_{p}0$. 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 $j$th element of the vector
$n^{-1}\sum\psi(X_{i};\delta_{i},\hat{\theta})$, denoted
$n^{-1}\sum\psi_{j}(X_{i};\delta_{i},\hat{\theta})$, around the true value
$\theta$:

$\displaystyle 0=n^{-1}\sum\psi_{j}(X_{i};\delta_{i},\hat{\theta})$
$\displaystyle=n^{-1}\sum\psi_{j}(X_{i};\delta_{i},\theta)+n^{-1}\sum\dot{\psi}_{j}(X_{i};\delta_{i},\theta)(\hat{\theta}-\theta)$
$\displaystyle+\frac{1}{2}(\hat{\theta}-\theta)^{T}n^{-1}\sum\ddot{\psi}_{j}(X_{i};\delta_{i},\widetilde{\theta}_{j})(\hat{\theta}-\theta)$

where $\widetilde{\theta}_{1},\dots,\widetilde{\theta}_{4}$ are on the line
segment joining $\hat{\theta}$ and $\theta$ and
$\ddot{\psi}_{j}(X_{i};\delta_{i},\theta_{j})$ is a $4\times 4$ matrix with
entry $(j,k)$ equal to
$\partial^{2}\psi_{j}(X_{i};\delta_{i},\theta_{j})/\partial\theta_{j}\partial\theta_{k}$
for $j,k\in\\{1,2,3,4\\}$. Writing these 4 equations in matrix notation yields

$0=n^{-1}\sum\psi(X_{i};\delta_{i},\theta)+\widetilde{R}(\hat{\theta}-\theta)$
(A.2)

where
$\widetilde{R}=\left\\{n^{-1}\sum\dot{\psi}(X_{i};\delta_{i},\theta)+(1/2)\widetilde{Q}\right\\}$
and

$\widetilde{Q}=\left[\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\
\sigma_{e}-\hat{\sigma}_{e}&\sigma_{p}-\hat{\sigma}_{p}&0&0\end{array}\right]$

is the $4\times 4$ matrix with $j$th row given by
$(\hat{\theta}-\theta)^{T}n^{-1}\sum\ddot{\psi}_{j}(X_{i};\delta_{i},\widetilde{\theta}_{j})$.
Note that $\widetilde{Q}\to_{p}0_{4\times 4}$ by the Weak Law of Large
Numbers, where in general $0_{r\times c}$ denotes an $r\times c$ matrix of
zeros. Also note that
$n^{-1}\sum\dot{\psi}(X_{i};\delta_{i},\theta)\to_{p}-\mathbb{A}(\theta)$,
where $-\mathbb{A}(\theta)$ is nonsingular under the assumption that
$\sigma_{e}>1-\sigma_{p}$. It follows that as $n\to\infty$, $\widetilde{R}$ is
invertible with probability one. On the set $S_{n}$ where $\widetilde{R}$ is
invertible, (A.2) can be rearranged to yield

$\hat{\theta}-\theta=\left(-\widetilde{R}\right)^{-1}\left\\{n^{-1}\sum\psi(X_{i};\delta_{i},\theta)\right\\}.$
(A.3)

Define

$\widetilde{R}^{*}=\frac{1}{1+g}\left\\{n^{-1}\sum-\dot{\psi}(X_{i};\delta_{i},\theta)\right\\}^{-1}(1/2)\widetilde{Q}\left\\{n^{-1}\sum-\dot{\psi}(X_{i};\delta_{i},\theta)\right\\}^{-1}$

where
$g=\mathrm{Tr}\left[-(1/2)\widetilde{Q}\left\\{n^{-1}\sum-\dot{\psi}(X_{i};\delta_{i},\theta)\right\\}^{-1}\right]$,
and note that $\widetilde{R}^{*}\to_{p}0_{4\times 4}$ by Slutsky’s theorem
because $\widetilde{Q}\to_{p}0_{4\times 4}$. Since $(1/2)\widetilde{Q}$ has
rank one,
$\left(-\widetilde{R}\right)^{-1}=\left\\{n^{-1}\sum-\dot{\psi}(X_{i};\delta_{i},\theta)\right\\}^{-1}+\widetilde{R}^{*}$
by an application of the Sherman-Morrison-Woodbury formula (Miller, 1981).
Thus, substituting this expression for $(-\widetilde{R})^{-1}$ into (A.3) and
multiplying by $\sqrt{n}$ yields (A.1), where
$R^{*}=\widetilde{R}^{*}\left\\{n^{-1}\sum\psi(X_{i};\delta_{i},\theta)\right\\}$.
Therefore, because $\lim_{n\to\infty}\Pr(S_{n})=1$,
$\widetilde{R}^{*}\to_{p}0_{4\times 4}$, and
$\sqrt{n}\left\\{n^{-1}\sum\psi(X_{i};\delta_{i},\theta)\right\\}\to_{d}\mathcal{N}(0,\mathbb{B}(\theta))$,
by Slutsky’s theorem $\sqrt{n}R^{*}\to_{p}0$.

###  A.2 Computation of asymptotic variance

Since $\mathbb{A}(\theta)$ is lower triangular, it follows that

$\mathbb{A}(\theta)^{-1}=\left[\begin{array}[]{cccc}1/c_{1}&0&0&0\\\
0&1/c_{2}&0&0\\\ 0&0&1/c_{3}&0\\\
-\dfrac{\pi}{c_{1}(\sigma_{e}+\sigma_{p}-1)}&\dfrac{1-\pi}{c_{2}(\sigma_{e}+\sigma_{p}-1)}&\dfrac{1}{c_{3}(\sigma_{e}+\sigma_{p}-1)}&(\sigma_{e}+\sigma_{p}-1)^{-1}\end{array}\right],$

and therefore
$\mathbb{A}(\theta)^{-1}\mathbb{B}(\theta)\mathbb{A}(\theta)^{-T}$ equals

$\left[\begin{array}[]{cccc}\sigma_{e}(1-\sigma_{e})/c_{1}&0&0&*\\\
0&\sigma_{p}(1-\sigma_{p})/c_{2}&0&*\\\ 0&0&{\rho}(1-{\rho})/c_{3}&*\\\
&*&*&V_{\pi,RG}\end{array}\right],$

where * denotes quantities not expressed explicitly and

$V_{\pi,RG}=\left\\{\frac{\pi^{2}\sigma_{e}(1-\sigma_{e})}{c_{1}}+\frac{(1-\pi)^{2}\sigma_{p}(1-\sigma_{p})}{c_{2}}+\frac{\rho(1-{\rho})}{c_{3}}\right\\}(\sigma_{e}+\sigma_{p}-1)^{-2}.$

##  Appendix B Proofs for Section 3.2

###  B.1 Proof of asymptotic normality

The Taylor expansion of
$n^{-1}\sum\psi(X_{i},Z_{i};\delta_{i},\hat{\theta}_{s})$ around the true
parameter $\theta_{s}$ yields

$\sqrt{n}(\hat{\theta}_{s}-\theta_{s})=\left\\{n^{-1}\sum-\dot{\psi}(X_{i},Z_{i};\delta_{i},\theta_{s})\right\\}^{-1}\sqrt{n}\left\\{n^{-1}\sum\psi(X_{i},Z_{i};\delta_{i},\theta_{s})\right\\}+\sqrt{n}R^{*},$

which is similar in form to (A.1), except here the estimating equations are
dependent on covariates $Z$. The remainder $R^{*}$ here is distinct from that
in Appendix A.1, and in general symbols may be reused and notation may not
hold the same meaning across appendices. Below it is established, using an
analogous approach to that of Appendix A.1, that
$n^{-1}\sum-\dot{\psi}(X_{i},Z_{i};\delta_{i},\theta_{s})\to_{p}\mathbb{A}(\theta_{s})$,
$\sqrt{n}\\{n^{-1}\sum\psi(X_{i},Z_{i};\delta_{i},\theta_{s})\\}\to_{d}\mathcal{N}\\{0,\mathbb{B}(\theta_{s})\\}$,
and $\sqrt{n}R^{*}\to_{p}0$. Therefore, by Slutsky’s theorem,
$\sqrt{n}(\hat{\theta}_{s}-\theta_{s})\to_{d}\mathcal{N}\\{0,\mathbb{A}(\theta_{s})^{-1}\mathbb{B}(\theta_{s})\mathbb{A}(\theta_{s})^{-T}\\}$.

First, define

$\displaystyle\mathbb{A}_{n}(X,Z,\delta,\theta_{s})=n^{-1}\sum\left\\{-\dot{\psi}(X_{i},Z_{i};\delta_{i},\theta_{s})\right\\}=\left[\begin{array}[]{cc}A^{\prime}&0_{(k+2)\times
2}\\\ C&D\end{array}\right]$

as a block matrix where
$A^{\prime}=\operatorname{diag}\left(n_{1}/n,n_{2}/n,n_{z_{1}}/n,\dots,n_{z_{k}}/n\right)$
is $(k+2)\times(k+2)$,

$C=\left[\begin{array}[]{ccccc}0&0&-\gamma_{1}&\dots&-\gamma_{k}\\\
\pi&-1+\pi&0&\dots&0\end{array}\right]$

is $2\times(k+2)$, and

$D=\begin{bmatrix}1&0\\\ -1&\sigma_{e}+\sigma_{p}-1\end{bmatrix}.$

Let

$\mathbb{A}(\theta_{s})=\left[\begin{array}[]{cc}A&0_{(k+2)\times 2}\\\
C&D\end{array}\right]$

where $A=\operatorname{diag}\\{c_{1},c_{2},c_{3}s_{1},\dots,c_{3}s_{k}\\}$,
and note that
$\mathbb{A}(\theta_{s})=\lim_{n\to\infty}\mathbb{A}_{n}(X,Z,\delta,\theta_{s})$
since $n_{j}/n\to c_{j}$ for $j\in\\{1,2,3\\}$. Note that $s_{1},\dots,s_{k}$
are all nonzero due to positivity, without which $\hat{\pi}_{SRG}$ is
undefined.

Second, define

$\displaystyle\mathbb{B}_{n}(X,Z,\delta,\theta_{s})=n^{-1}\sum\mathbb{E}\\{\psi(X_{i},Z_{i};\delta_{i},\theta_{s})\psi(X_{i},Z_{i};\delta_{i},\theta_{s})^{T}\\}=\begin{bmatrix}E^{\prime}&0_{(k+2)\times
2}\\\ 0_{2\times(k+2)}&0_{2\times 2}\end{bmatrix},$

where

$E^{\prime}=\operatorname{diag}\\{n_{1}\sigma_{e}(1-\sigma_{e})/n,n_{2}\sigma_{p}(1-\sigma_{p})/n,n_{z_{1}}\rho_{1}(1-\rho_{1})/n,\dots,n_{z_{k}}\rho_{k}(1-\rho_{k})/n\\},$

since
$\mathbb{E}(\psi_{\rho_{j}}^{2})=\mathbb{E}\\{I(Z_{i}=z_{j},\delta_{i}=3)(X_{i}-\rho_{j})^{2}\\}=I(Z_{i}=z_{j},\delta_{i}=3)\rho_{j}(1-\rho_{j})$.
Define $\mathbb{B}(\theta_{s})$ to have the same form as
$\mathbb{B}_{n}(X,Z,\delta,\theta_{s})$ except $E^{\prime}$ is replaced by

$E=\operatorname{diag}\\{c_{1}\sigma_{e}(1-\sigma_{e}),c_{2}\sigma_{p}(1-\sigma_{p}),c_{3}s_{1}\rho_{1}(1-\rho_{1}),\dots,c_{3}s_{k}\rho_{k}(1-\rho_{k})\\}.$

Since
$\mathbb{B}(\theta_{s})=\lim_{n\to\infty}\mathbb{B}_{n}(X,Z,\delta,\theta_{s})$,
$\sqrt{n}\left\\{n^{-1}\sum\psi(X_{i},Z_{i};\delta_{i},\theta_{s})\right\\}\to_{d}\mathcal{N}\\{0,\mathbb{B}(\theta_{s})\\}$
by the Lindeberg-Feller CLT. In particular, let

$\displaystyle v_{n}^{2}$
$\displaystyle=\sum\\{\operatorname{Var}(\psi_{e})+\operatorname{Var}(\psi_{p})+\operatorname{Var}(\psi_{\rho_{1}})+\dots+\operatorname{Var}(\psi_{\rho_{k}})+\operatorname{Var}(\psi_{\rho})+\operatorname{Var}(\psi_{\pi})\\}$
$\displaystyle=n_{1}\sigma_{e}(1-\sigma_{e})+n_{2}\sigma_{p}(1-\sigma_{p})+\sum_{j=1}^{k}n_{z_{j}}\rho_{j}(1-\rho_{j}),$

and note that
$\max_{i}\lVert\psi(X_{i},Z_{i};\delta_{i},\theta_{s})\rVert\leq\sqrt{3}$
because at most one of
$\psi_{e},\psi_{p},\psi_{\rho_{1}},\dots,\psi_{\rho_{k}}$ is nonzero and each
element of $\psi(X_{i},Z_{i};\delta_{i},\theta_{s})$ has a maximum value of
one. Therefore for all $\epsilon>0$,

$\lim_{n\to\infty}v_{n}^{-2}\sum\mathbb{E}\left\\{\lVert\psi(X_{i};\delta_{i},\theta)\rVert^{2}I(\lVert\psi(X_{i};\delta_{i},\theta)\rVert\geq\epsilon
v_{n})\right\\}\leq\lim_{n\to\infty}3v_{n}^{-2}\sum P(\sqrt{3}\geq\epsilon
v_{n})=0,$

implying the Lindeberg condition holds.

Third, let $\widetilde{Q}_{s}$ be analogous to $\widetilde{Q}$ from Appendix
A.1. Denote the $j$th entry of $\theta_{s}$ as $\theta_{s_{j}}$ and note
$\lvert\partial^{2}\psi_{j}(X_{i},Z_{i};\delta_{i},\theta_{s})/\partial\theta_{s_{j}}\theta_{s_{l}}\rvert\leq
1$ for all $j,l\in\\{1,2,\dots,k+4\\}$. Thus there exists a function
$g(X_{i},Z_{i},\delta_{i})$ such that, for each $\theta_{s}^{*}$ in a
neighborhood of $\theta_{s}$,
$\left\lvert\partial^{2}\psi_{j}(X_{i},Z_{i};\delta_{i},\theta_{s})/\partial\theta^{*}_{s_{j}}\partial\theta^{*}_{s_{l}}\right\rvert\leq
g(X_{i},Z_{i},\delta_{i})$ for all $X_{i},Z_{i},\delta_{i}$ where $\int
g(X_{i},Z_{i},\delta_{i})dF(X_{i},Z_{i},\delta_{i})<\infty$; e.g., let
$g(X_{i},Z_{i},\delta_{i})=2$. Therefore, each entry in $\widetilde{Q}_{s}$ is
bounded by $\lVert\hat{\theta}_{s}-\theta_{s}\rVert n^{-1}\sum
g(X_{i},\delta_{i})=o_{p}(1)$, so
$\widetilde{Q}_{s}\to_{p}0_{(k+4)\times(k+4)}$. This fact can be used to prove
$\sqrt{n}R^{*}\to_{p}0$ with the same technique as in Appendix A.1, so we omit
the rest of the proof.

###  B.2 Computation of asymptotic variance

Note that $\mathbb{A}(\theta_{s})$ is block lower triangular, and thus

$\mathbb{A}(\theta_{s})^{-1}=\begin{bmatrix}A^{-1}&0_{(k+2)\times 2}\\\
-D^{-1}CA^{-1}&D^{-1}\end{bmatrix}$

and therefore

$\displaystyle\mathbb{A}(\theta_{s})^{-1}\mathbb{B}(\theta_{s})\mathbb{A}(\theta_{s})^{-T}$
$\displaystyle=\begin{bmatrix}A^{-1}EA^{-1}&-A^{-1}EA^{-1}C^{T}D^{-T}\\\
-D^{-1}CA^{-1}EA^{-1}&D^{-1}CA^{-1}EA^{-1}C^{T}D^{-T}\end{bmatrix}.$

The bottom right submatrix of
$\mathbb{A}(\theta_{s})^{-1}\mathbb{B}(\theta_{s})\mathbb{A}(\theta_{s})^{-T}$,
and specifically the bottom right element of that submatrix, is of primary
interest. Since $A$ and $E$ are both diagonal,

$A^{-1}EA^{-1}=\operatorname{diag}\\{c_{1}^{-1}\sigma_{e}(1-\sigma_{e}),c_{2}^{-1}\sigma_{p}(1-\sigma_{p}),c_{3}^{-1}s_{1}^{-1}\rho_{1}(1-\rho_{1}),\dots,c_{3}^{-1}s_{k}^{-1}\rho_{k}(1-\rho_{k})\\}.$

Next, note

$D^{-1}C=\begin{bmatrix}0&0&-\gamma_{1}&\dots&-\gamma_{k}\\\
\dfrac{\pi}{\sigma_{e}+\sigma_{p}-1}&\dfrac{-1+\pi}{\sigma_{e}+\sigma_{p}-1}&\dfrac{-\gamma_{1}}{\sigma_{e}+\sigma_{p}-1}&\dots&\dfrac{-\gamma_{k}}{\sigma_{e}+\sigma_{p}-1}\end{bmatrix}.$

Therefore, $D^{-1}CA^{-1}EA^{-1}$ equals

$\displaystyle\begin{bmatrix}0&0&-\gamma_{1}c_{3}^{-1}s_{1}^{-1}\rho_{1}(1-\rho_{1})&\dots&-\gamma_{k}c_{3}^{-1}s_{k}^{-1}\rho_{1}(1-\rho_{1}))\\\
\dfrac{\pi\sigma_{e}(1-\sigma_{e})}{c_{1}(\sigma_{e}+\sigma_{p}-1)}&\dfrac{(-1+\pi)\sigma_{p}(1-\sigma_{p})}{c_{2}(\sigma_{e}+\sigma_{p}-1)}&-\dfrac{\gamma_{1}\rho_{1}(1-\rho_{1})}{c_{3}s_{1}(\sigma_{e}+\sigma_{p}-1)}&\dots&-\dfrac{\gamma_{k}\rho_{k}(1-\rho_{k})}{c_{3}s_{k}(\sigma_{e}+\sigma_{p}-1)}\end{bmatrix}$

and thus

$D^{-1}CA^{-1}EA^{-1}C^{T}D^{-T}=\left[\begin{array}[]{cc}*&\ *\\\
&V_{\pi,SRG}\end{array}\right]$

where * denotes quantities not expressed explicitly and

$V_{\pi,SRG}=\left\\{\frac{\pi^{2}\sigma_{e}(1-\sigma_{e})}{c_{1}}+\frac{(1-\pi)^{2}\sigma_{p}(1-\sigma_{p})}{c_{2}}+\sum_{j=1}^{k}\frac{\gamma_{j}^{2}\rho_{j}(1-{\rho_{j}})}{c_{3}s_{j}}\right\\}(\sigma_{e}+\sigma_{p}-1)^{-2}.$

##  Appendix C Proof for Section 3.3

The proof that $\theta_{m}$ is asymptotically normal is given for the case of
logistic regression, but it extends to any link function $g(\cdot)$
appropriate for binary regression. Recall that the $(p+4)$-vector of
estimating equations is

$\sum\psi(X_{i},Z_{i};\delta_{i},\theta_{m})=\left(\sum\psi_{e},\sum\psi_{p},\sum\psi_{\beta},\psi_{\rho},\psi_{\pi}\right)^{T}=0.$

In the above vector $\sum\psi_{e},\sum\psi_{p}$, and $\psi_{\pi}$ are
identical to the equations used in Appendices A and B;
$\psi_{\rho}=\sum_{j=1}^{k}\operatorname{logit}^{-1}\\{\beta
h(Z_{j})\\}\gamma_{j}-\rho$; and $\sum\psi_{\beta}$ is a $p$-vector with $j$th
element $\sum\psi_{\beta_{j}}=\sum
I(\delta_{i}=3)\left[X_{i}-\operatorname{logit}^{-1}\\{\beta
h(Z_{i})\\}\right]h_{j}(Z_{i})$. The Taylor expansion of
$n^{-1}\sum\psi(X_{i},Z_{i};\delta_{i},\hat{\theta}_{m})$ around the true
parameter $\theta_{m}$ yields

$\sqrt{n}(\hat{\theta}_{m}-\theta_{m})=\left\\{n^{-1}\sum-\dot{\psi}(X_{i},Z_{i};\delta_{i},\theta_{m})\right\\}^{-1}\sqrt{n}\left\\{n^{-1}\sum\psi(X_{i},Z_{i};\delta_{i},\theta_{m})\right\\}+\sqrt{n}R^{*}.$

The rest of the proof is similar to Appendices A.1 and B.1. Namely, first
define

$\mathbb{A}_{n}(X,Z,\delta,\theta_{m})=n^{-1}\sum\left\\{-\dot{\psi}(X_{i},Z_{i};\delta_{i},\theta_{m})\right\\}=\left[\begin{array}[]{ccc}A^{\prime}&0_{2\times
p}&0_{2\times 2}\\\ 0_{p\times 2}&B^{\prime}&0_{p\times 2}\\\
C&D&E\end{array}\right]$

as a block matrix where $A^{\prime}=\operatorname{diag}(n_{1}/n,n_{2}/n)$,
$B^{\prime}$ is $p\times p$ with entry $(j,k)$ equal to $n^{-1}\sum
I(\delta_{i}=3)h_{j}(Z_{i})h_{k}(Z_{i})\exp\\{\beta
h(Z_{i})\\}/\left[1+\exp\\{\beta h(Z_{i})\\}\right]^{2}$,

$C=\begin{bmatrix}0&0\\\ \pi&-1+\pi\end{bmatrix},$
$D=\left[\begin{array}[]{ccccc}-\sum_{j=1}^{k}h_{1}(Z_{j})\dfrac{\exp\\{\beta
h(Z_{i})\\}}{\left[1+\exp\\{\beta
h(Z_{i})\\}\right]^{2}}\gamma_{j}&\dots&-\sum_{j=1}^{k}h_{p}(Z_{j})\dfrac{\exp\\{\beta
h(Z_{i})\\}}{\left[1+\exp\\{\beta h(Z_{i})\\}\right]^{2}}\gamma_{j}\\\
0&\dots&0\end{array}\right]$

is $2\times p$, and

$E=\begin{bmatrix}1&0\\\ -1&\sigma_{e}+\sigma_{p}-1\end{bmatrix}.$

Let $\mathbb{A}(\theta_{m})$ have the same form as
$\mathbb{A}_{n}(X,Z,\delta,\theta_{m})$, replacing $A^{\prime}$ with
$A=\operatorname{diag}(c_{1},c_{2})$ and $B^{\prime}$ with $B$, with entry
$(j,k)$ of $B$ equal to $c_{3}\mathbb{E}\left(h_{j}(Z)h_{k}(Z)\exp\\{\beta
h(Z_{i})\\}/\left[1+\exp\\{\beta h(Z_{i})\\}\right]^{2}\right)$. Note that
$\mathbb{A}_{n}(X,Z,\delta,\theta_{m})\to_{p}\mathbb{A}(\theta_{m})$ by the
Weak Law of Large Numbers.

Second, define

$\mathbb{B}_{n}(X,Z,\delta,\theta_{m})=n^{-1}\sum\mathbb{E}(\psi(X_{i},Z_{i};\delta_{i},\theta_{m})\psi(X_{i},Z_{i};\delta_{i},\theta_{m})^{T})=\left[\begin{array}[]{ccc}F^{\prime}&0_{2\times
p}&0_{2\times 2}\\\ 0_{p\times 2}&G^{\prime}&0_{p\times 2}\\\ 0_{2\times
2}&0_{2\times p}&0_{2\times 2}\end{array}\right]$

where
$F^{\prime}=\operatorname{diag}\left\\{n_{1}\sigma_{e}(1-\sigma_{e})/n,n_{2}\sigma_{p}(1-\sigma_{p})/n\right\\}$
and $G^{\prime}$ is a $p\times p$ matrix with entry $(j,k)$ equal to
$(n_{3}/n)\mathbb{E}\left(h_{j}(Z)h_{k}(Z)\left[X-\operatorname{logit}^{-1}\\{\beta
h(Z)\\}\right]^{2}\right)$. Let $\mathbb{B}(\theta_{m})$ be of the same form
as $\mathbb{B}_{n}(X,Z,\delta,\theta_{m})$ except with $F$ and $G$ replacing
$F^{\prime}$ and $G^{\prime}$, where
$F=\operatorname{diag}\\{c_{1}\sigma_{e}(1-\sigma_{e}),c_{2}\sigma_{p}(1-\sigma_{p})\\}$
and $G$ is identical to $G^{\prime}$ except with $c_{3}$ replacing $n_{3}/n$
in each element. Noting that
$\mathbb{B}(\theta_{m})=\lim_{n\to\infty}\mathbb{B}_{n}(X,Z,\delta,\theta_{m})$,
it follows that
$\sqrt{n}\left\\{n^{-1}\sum\psi(X_{i},Z_{i};\delta_{i},\theta_{m})\right\\}\to_{d}\mathcal{N}\left\\{0,\mathbb{B}(\theta_{m})\right\\}$
by the Lindeberg-Feller CLT. Specifically, note that the range of the inverse
logit function is $[0,1]$ and that all elements of $h(Z)$, the user-specified
function of the covariates, must be finite. Then the Lindeberg condition holds
by the same logic as in Appendices A.1 and B.1.

Finally, $\sqrt{n}R^{*}\to_{p}0$ can be proved as in Appendix B.1. Letting
$\theta_{m_{j}}$ denote the $j$th entry of $\theta_{m}$, the proof follows
from the fact that
$\lvert\partial^{2}\psi_{j}(X_{i},Z_{i};\delta_{i},\theta_{m})/\partial\theta_{m_{j}}\partial\theta_{m_{k}}\rvert<\infty$
for all $j,k\in\\{1,2,\dots,p+4\\}$. Therefore, by Slutsky’s Theorem

$\sqrt{n}(\hat{\theta}_{m}-\theta_{m})\to_{d}\mathcal{N}\left\\{0,\mathbb{A}(\theta_{m})^{-1}\mathbb{B}(\theta_{m})\mathbb{A}(\theta_{m})^{-T}\right\\}\
\text{and}\ \sqrt{n}(\hat{\pi}-\pi)\to_{d}\mathcal{N}(0,V_{\pi,SRGM})$

where the bottom right element of
$\mathbb{A}(\theta_{m})^{-1}\mathbb{B}(\theta_{m})\mathbb{A}(\theta_{m})^{-T}$
is denoted as $V_{\pi,SRGM}$. The asymptotic variance of $\hat{\theta}_{m}$
can be consistently estimated by the empirical sandwich variance estimator
$\mathbb{A}(\hat{\theta}_{m})^{-1}\mathbb{B}(\hat{\theta}_{m})\mathbb{A}(\hat{\theta}_{m})^{-T}$;
the bottom right element of the sandwich estimator is $\hat{V}_{\pi,SRGM}$.

## Appendix Figures

Figure A.1:  Empirical bias of the Rogan-Gladen ($\hat{\pi}_{RG}$) estimator
from simulation study for DGP 1, described in Section 4.1. The six facets
correspond to a given combination of sensitivity $\sigma_{e}$ (‘Sens’) and
specificity $\sigma_{p}$ (‘Spec’).  Figure A.2: Confidence interval coverage
of the Rogan-Gladen ($\hat{\pi}_{RG}$) estimator from simulation study for DGP
1, described in Section 4.1. Figure A.3:  Empirical bias of the Rogan-Gladen
($\hat{\pi}_{RG}$) and nonparametric standardized ($\hat{\pi}_{SRG}$)
estimators from simulation study for DGP 2, described in Section 4.2. Figure
A.4: Confidence interval coverage of the Rogan-Gladen ($\hat{\pi}_{RG}$) and
nonparametric standardized ($\hat{\pi}_{SRG}$) estimators from simulation
study for DGP 2, described in Section 4.2. Figure A.5: Confidence interval
coverage of the Rogan-Gladen ($\hat{\pi}_{RG}$), nonparametric standardized
($\hat{\pi}_{SRG}$), and parametric standardized ($\hat{\pi}_{SRGM}$)
estimators from simulation study for DGP 3, described in Section 4.3.1.

Figure A.6: Confidence interval coverage of the Rogan-Gladen
($\hat{\pi}_{RG}$), nonparametric standardized ($\hat{\pi}_{SRG}$), and
parametric standardized ($\hat{\pi}_{SRGM}$) estimators from simulation study
for DGP 4, described in Section 4.3.2.

Figure A.7: Empirical bias of the estimators from simulation study for DGP 3
under model misspecification, described in Section 4.4.

Figure A.8: Confidence interval coverage of the estimators from simulation
study for DGP 3 under model misspecification, described in Section 4.4.

Figure A.9: Empirical bias of the estimators from simulation study for DGP 4
under model misspecification, described in Section 4.4.

Figure A.10: Confidence interval coverage of the estimators from simulation
study for DGP 4 under model misspecification, described in Section 4.4.

## Appendix References

Boos, D. D., & Stefanski, L. A. (2013). Essential Statistical Inference.
Springer New York, Ch. 7, pp. 297-337.

Miller, K. S. (1981). On the inverse of the sum of matrices. Mathematics
Magazine, 54(2), 67-72.