# Trying to Estimate Disease Prevalence from Fragmentary Test Results

In response to the spread of COVID-19 disease, all Californians were ordered on 19 March 2020 to stay at home, except for such necessary errands as trips to grocery stores, pharmacies, etc. On 21 March the San Francisco Chronicle reported that about $$11\,500$$ tests for the virus have been completed so far with $$1206$$ positive results among them.

Can anything useful be said about the prevalence of the virus in California based on these results? If (a) the 11500 tests were from randomly chosen Californians and (b) the tests were of "gold standard" quality, then prevalence would be close to $$1206/11500 = 0.1049 \approx 10.5\%.$$

However, according to the Chronicle article, guidelines have centered on testing people showing symptoms of COVID-19 disease, care workers who have had contact with COVID-19 patients, and people who have symptoms of pneumonia, but test negative for pneumonia. Moreover, tests are "coming from a growing list of private companies, private labs, and public and private universities," so there is no way to know proportions of false positives and false negatives for any one test, let alone the overall rates for whatever mix of tests accounts for the 11,500 reported results.

In my answer below, I will extend my answer to this question to give prevalence values resulting from various assumptions about the sensitivity and specificity of the tests involved. I use the notation of that page. This analysis may help to answer the question for whatever population produced the specimens tested.

Additional answers are welcome, perhaps showing alternative methods or based on additional information that may become available as the pandemic continues

• You may find this post interesting medium.com/@tomaspueyo/…. – Jesper for President Mar 22 at 9:26
• @JesperforPresident. My post is intended to be purely based on probability computations, absent political commentary. – BruceET Jun 13 at 22:34

Notation. Let $$\pi = P(\text{Disease})$$ be the prevalence of the disease in the population and $$\tau = P(\text{Pos Test})$$ be the proportion testing positive. For the test, let $$\eta = P(\text{Pos}|\text{Disease})$$ be the sensitivity and $$\theta = P(\text{Neg}|\text{No Disease})$$ be its specificity. Also, given test results, let $$\gamma = P(\text{Disease}| \text{Pos})$$ and $$\delta = P(\text{No Disease}| \text{Neg})$$ be, respectively, the predictive powers of a positive or negative test.

If a test is of gold standard quality with $$\eta = \theta = 1,$$ then $$\pi = \tau.$$ Tests that accurately sequence the genome of the virus may be gold standard tests. Often the first tests for a virus may have considerably lower values of $$\eta$$ and $$\theta.$$ It is difficult to find accounts of values of $$\eta$$ and $$\theta$$ for any of the tests in current use for COVID-19. (According to one unauthorized report, a test used in China had $$\theta \approx 0.7.)$$

Traditional estimate. First, we look at results for tests with $$\eta = \theta = 0.95.$$ Then for data with $$n= 11\,500; a = 1206,$$ we have $$\hat \tau = t = 0.1049.$$ The associated 95% Wald confidence interval for $$\tau$$ is $$(0.0993, 0.1105)$$ from which one can derive the confidence interval $$(0.0547, 0.0672)$$ for $$\pi.$$ Also, $$\pi = 0.061$$ implies that the predictive power of a positive test is $$\gamma = 0.5523.$$ Only about half of the subjects testing positive are actually infected. Some computations in R follow:

ETA = THETA = .95
n = 11500;  a = 1206;  p0m = -1:1; t = a/n
wald.TAU = t +p0m*1.96*sqrt(t*(1-t)/n);  wald.TAU
[1] 0.09926973 0.10486957 0.11046940
ci.PI = (wald.T + THETA - 1)/(ETA + THETA - 1);  ci.PI
[1] 0.05474415 0.06096618 0.06718822
PI = (t + THETA -1)/(ETA + THETA -1); PI
[1] 0.06096618
GAMMA = PI*ETA/(PI*ETA + (1-PI)*(1-THETA)); GAMMA
[1] 0.5522849


When the traditional estimate is problematic. For a poorer test with $$\eta = \theta = 0.90,$$ this method gives a CI for $$\pi$$ as $$(-0.0009, 0.0131),$$ which has a (nonsensical) negative left endpoint. (We would expect about 1150 false positive tests even with no infected subjects. This is getting close to the the observed number 1206 of positive tests.) In such circumstances, one wonders whether to trust the point estimates $$\hat \pi = 0.0061$$ and $$\hat \gamma = 0.0522.$$

ETA = THETA = .9
n = 11500;  a = 1206;  p0m = -1:1; t = a/n
wald.TAU = t +p0m*1.96*sqrt(t*(1-t)/n);  wald.TAU
[1] 0.09926973 0.10486957 0.11046940
ci.PI = (wald.T + THETA - 1)/(ETA + THETA - 1);  ci.PI
[1] -0.0009128343  0.0060869565  0.0130867473
PI = (t + THETA -1)/(ETA + THETA -1); PI
[1] 0.006086957
GAMMA = PI*ETA/(PI*ETA + (1-PI)*(1-THETA)); GAMMA
[1] 0.05223881


A Gibbs sampler. One useful alternative approach is to assume a beta prior distribution on prevalence $$\pi \sim \mathsf{Beta(\alpha, \beta)}.$$ Even if noninformative with $$\alpha = \beta = 0.5,$$ such a prior distribution excludes value of $$\pi$$ outside $$(0,1).$$ Then we use a Gibbs sampler to find the posterior distribution of $$\pi,$$ given our data $$n = 11\,500, a = 1206.$$ Its steps, within each iteration, are as follows:

• We begin with an arbitrary value of $$\pi_1 \in (0,1)$$ and use it to estimate 'latent' counts of subjects with the disease based on predictive values $$\gamma$$ and $$\delta.$$ We sample counts $$X \sim \mathsf{Binom}(a, \gamma)$$ and $$Y \sim \mathsf{Binom}(n-A, 1 - \delta).$$

• Then with the estimated $$S = X+Y$$ infected subjects, we update the beta prior at this step as $$\pi|S \sim \mathsf{BETA}(\alpha + S, \beta + n - S).$$ Finally, we sample $$\pi_2$$ from this updated distribution. Even with arbitrary $$\pi_1,$$ this new value $$\pi_2$$ is likely to be closer to the truth.

Simulated posterior distribution. Iterating through many such steps we get successive values $$\pi_1, \pi_2, \pi_3, \dots ,$$ in a convergent Markov chain, for which the limiting distribution is the desired posterior distribution of the prevalence. To make sure that the chain has reached steady state, we use only the last half of the many values $$\pi_i$$ generated in this way.

Cutting 2.5% of the probability from each tail of the simulated posterior distribution, we can obtain a 95% Bayesian probability estimate for prevalence $$\pi.$$

set.seed(1020)
m = 10^5                             # iterations
PI = GAMMA = Gamma = numeric(m)      # vectors for results
PI[1] = .5                           # initial value
alpha = .5; beta = .5                # parameters of beta prior
ETA = .9; THETA = .9                 # sensitivity; specificity
n = 11500; A = 1206; B = n - A       # data
for (i in 2:m)
{
num.x = PI[i-1]*ETA; den.x = num.x + (1-PI[i-1])*(1 - THETA)
GAMMA[i] = num.x/den.x
X = rbinom(1, A, num.x/den.x)  # use est of gamma as probability
num.y = PI[i-1]*(1 - ETA); den.y = num.y + (1-PI[i-1])*THETA
Y = rbinom(1, B, num.y/den.y)  # use 1 - est of delta as probability
PI[i] = rbeta(1, X + Y + alpha, n - X - Y + beta)
}
aft.brn = seq(floor(m/2),m)
quantile(PI[aft.brn], c(.025, .975))
2.5%        97.5%
3.329477e-05 1.225794e-02
quantile(PI[aft.brn], .95)
95%
0.01101075
mean(PI[aft.brn])
[1] 0.0049096

quantile(GAMMA[aft.brn], c(.025, .975))
2.5%        97.5%
0.0002995732 0.1004690791
mean(GAMMA[aft.brn])
[1] 0.04176755
quantile(Gamma[aft.brn], c(.025, .975))
2.5%        97.5%
0.0002995732 0.1004690791


Because the two-sided Bayesian probability interval $$(.00003, .0123)$$ has its lower endpoint so near to 0, we also look at the one-sided interval $$(0, .0110)$$ for $$\pi.$$

Because we estimate the predictive power $$\gamma$$ of a positive test at each step of the chain, we capture its values in order to get a 95\% Bayesian probability interval $$(0.0003, 0.1005)$$ for the predictive power $$\gamma$$ of a positive test. If we were to sequester subjects that get a positive result with this test, then only a relatively small proportion of sequestered subjects would actually be infected.

Diagnostic Plots. Not all Gibbs samplers converge as anticipated. Diagnostic plots show that this one does.

• A plot of successive values of $$\pi$$ shows the nature of the convergence of the Markov chain. The history plot of $$\pi$$ shows that the chain 'mixes well'; that is, it moves freely among appropriate values. There are no points of near-absorption.

• The trace of running averages of the $$\pi_i$$ shows smooth convergence to prevalence 0.0049. Vertical blue lines indicate the burn-in period.

• The ACF plot shows that the $$\pi_i$$ are not independent. Among the $$m = 100\,000$$ values, perhaps there are 1000 independent ones. In many Gibbs samplers, Markov dependence 'wears away' more rapidly than here.

• The plot at lower-right is a history plot of the $$\gamma_i.$$

Variations. If we run the same program with $$\eta = \theta = .95,$$ results are nearly the same as for the traditional procedure. If we have useful prior information (or opinions) about the prevalence, we can incorporate that information into the prior distribution on $$\pi.$$

References: (1) Suess, Gardner, & Johnson (2002), "Hierarchical Bayesian model for prevalence inferences and determination of a country’s status for an animal pathogen" Preventive Veterinary Medicine, and its references. (2) Suess & Trumbo (2010) Probability simulation and Gibbs sampling, (Sect. 9.1), Springer.