Binary disease tests are used to estimate disease prevalence in a population. When test kits are scarce, researchers and epidemiologists may choose to test samples pooled from multiple people. Further, with observational studies and meta-analyses, you may end up with results from both indiviudal and pooled samples. To combine these results into a single estimate of the population prevalence, some prior work leapt to using compute-heavy and slow MCMC methods. We (Charles Copley, my colleague from Ancera, and I) found a closed-form posterior distribution for the population prevalence that can be computed in a fraction of a second. This post is an overview of that work. The full derivation + analysis is on arXiv.
Let’s say we are trying to measure a population prevalence \(p\). If you were to sample \(n\) individuals, then you would see about \(np\) positives. Now, what if \(p=.01\)? Well you would need to sample 100 individuals to see on average 1 positive result. We can boost that signal by pooling.
The probability of a pooled sample testing positive when \(q\) individuals are pooled, \(\pi _q\), is equivalent to the probability that at least one of the samples which are combined into the pooled sample is positive. We’ll assume that each indiviudal’s positivity is i.i.d. Bernoulli.
\[\pi_q = 1 - \Pr(\text{all }q\text{ samples negative}) = 1 - (1-p)^q\]Or, if you want to see it visually:
Let’s quantify the difference in information a single individual test yields vs. a single pooled test.
Say \(X\) is a Bernoulli random variable representing the result of an individual test. Its likelihood function over possible outcomes \(k\) is
\[f(k \mid p) = p^k (1-p)^{1-k}\]Say \(Y\) is a Bernoulli random variable representing the result of a pooled test. Its likelihood function over possible outcomes \(k\) is
\[g(k \mid p) = \bigg(1 - (1-p)^q\bigg)^k \bigg((1-p)^q\bigg)^{1-k}\]Fisher information \(I(p)\) is calculated as the negative expected value of the second derivative of the log-likelihood. It is a measurement of how much information reads of a RV give us about an unknown parameter of that RV’s distribution.
For the individual test,
\[I_X(p) = \frac{1}{p} + \frac{1}{1-p} = \frac{1}{p(1-p)}\]For the pooled test (thank you Wolfram Alpha),
\[\frac{\partial d^2}{\partial d p^2} \ln g(k \mid p) = -\frac{q\left(k(1-p)^q + kq(1-p)^q - k - 2(1-p)^q + (1-p)^{2q} + 1\right)}{(p-1)^2\left((1-p)^q - 1\right)^2}\]so
\[I_Y(p) = \frac{q^2 (1-p)^{q-2}}{1-(1-p)^q}\]Notice that our \(q=1\) case matches \(I_X (p)\). So, the information ratio between these is
\[\frac{I_Y(p)}{I_X(p)} = \frac{q^2 p (1-p)^{q-1}}{1-(1-p)^q}\]The ratio approaches \(q\) as \(p\rightarrow0\), meaning that at low prevalences a pooled test of \(q\) samples is up to \(q\) times as informative as an individual sample, meaning that you can get the same info for about \(q\) times as few samples.
Define \(P\) as a random variable for population prevalence. We will use \(p\) to denote outcome values of that random variable.
Define \(m\) as the number of individual tests conducted.
Define \(Y_1,Y_2,...Y_m\) as random variables for the binary results of these individual tests.
Define \(y\) as the number of observed positive individual tests.
Define \(n\) as the number of pooled tests conducted.
Define \(Z_1,Z_2,...Z_n\) as random variables for the binary results of these pooled tests.
Define \(z\) as the number of observed positive pooled tests.
As before, we’ll use \(\pi_q\) for the probability of a pooled sample testing positive when \(q\) individuals are pooled.
In general, we assume a beta prior distribution for \(P\) \(P \sim Beta(\alpha, \beta)\) so
\[Pr(P=p) = \frac{ p^{\alpha - 1} (1-p)^{\beta - 1}}{B(\alpha, \beta)}\]where \(B(\alpha, \beta)\) is the beta function of \(\alpha\) and \(\beta\).
Bayes’ Theorem can be stated for this problem as
\[Pr(P=p|Y=y, Z=z) = \frac{Pr(Y=y, Z=z|P=p)Pr(P=p)}{Pr(Y=y, Z=z)}\] \[Y \sim \text{Binomial}(m, p), \qquad Z \sim \text{Binomial}(n, \pi_q), \qquad P \sim \text{Beta}(\alpha,\beta)\]At this point, the prior work plugged these equations into an MCMC solver. We will solve \(\Pr(P=p \mid Y=y, Z=z)\) analytically instead.
For the full derivation, you should check out our arxiv post, but the gist is that you can use a binomial expansion of \(\pi_q^z\)
\[\label{pibin} \pi _q ^{z} = [1 - (1-p)^{q}]^{z} = \sum _{i=0} ^{z} \binom{z}{i} (-1)^{i} (1-p)^{qi}\]and linearity of integration to solve a closed form for the “evidence”, \(Pr(Y=y, Z=z)\). After a few pages of derivations, and substituting \( \gamma = y+\alpha \) and \( \delta = m-y+\beta+qn-qz \), we can derive the posterior distribution.
\[\Pr(P=p \mid y, z) = \frac{\displaystyle\sum_{i=0}^{z}\binom{z}{i}(-1)^{i}\,p^{\gamma-1}(1-p)^{\delta+qi-1}} {\displaystyle\sum_{i=0}^{z}\binom{z}{i}(-1)^{i}\,B(\gamma,\delta+qi)}\]This function is exact, integrable, and evaluates way faster than a MCMC chain. Further, if you need to include sensitivity or specificity, those generalize in the same way, even if you allow them to vary along another beta distribution.
One cool thing about this equation is we can rewrite it to reveal a bit of structure. The posterior is a affine combination of beta distributions!
Let \(f(p, \gamma, \delta + qi)\) be the PDF of the beta distribution with parameters \(\gamma\), and \(\delta + qi\), evaluated at \(p\).
So we can also write the posterior probability distribution for \(P\) in this form:
\[\Pr(P=p \mid y, z) = \frac{\displaystyle\sum_{i=0}^{z}\binom{z}{i}(-1)^{i}\,B(\gamma,\delta+qi)\,f(p;\gamma,\delta+qi)} {\displaystyle\sum_{i=0}^{z}\binom{z}{i}(-1)^{i}\,B(\gamma,\delta+qi)}\]The widget below shows a given posterior’s weighted components and their sum. It also shows the closest-fitting single beta distribution.
Explore the posterior! Drag the sliders to set the experimental setting and see the resulting posterior. Hit resample to change the drawn data and resulting fit. You can see that pooled tests are more useful at low \(p\) and that degraded sensitivity and specificity result in a wider posterior distribution.
So, practically, if you have a fixed test budget (who doesn’t) and you expect low prevalence, you should consider pooling. Our paper shows extensive simulation of the tradeoffs here, so check it out if you’re curious.