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.
