header logo

Measuring the prevalence of disease with imperfect diagnostic tests (part II)

Photo credit to Johnathan Campos

Goals of this article

We ended part I with this statement

$$ P\left( \text{observation}|\text{hypothesis}\right) \neq P\left( \text{hypothesis}|\text{observation}\right) $$

Part I was a discussion of \(P\left( \text{observation}|\text{hypothesis}\right)\). In this article we are going to explore \(P\left( \text{hypothesis}|\text{observation}\right) \).

Author note

This article is heavily inspired by Keith Richards-Dinger's article. My intellectual contribution here is minimal and superficial but I think it adds some clarity. If nothing else, I improved the color choices (see bottom).

Sensitivity and specificity are the most common terms in medical diagnositic testing. They are the terms that physicians learn, the quantities the FDA wants reported, and so on. For this reason, I view using other metrics which are linear transforms of sensitivity and specificity as adding unnecessary obfuscation to any work. In reproducing Keith's results, I used sensitivity and specificity rather than false positive/negative rates. This is not an error on Keith's part. His results are still true but I feel strongly that this restatement adds clarity, particularly for physicians.

Consider the weighted coin problem backward

The probability distribution function of the true weight of a coin given \(h\) heads in \(n\) flips is given by a beta distribution. This is stated both by this guy and Keith. However, I do not have an authoritative source (textbook, journal article, etc) to support this claim. I think it is likely true but I am still looking for a first principles derivation to explain why it is true. That said, I will proceed assuming this is correct.

$$ P(x, h, n) = \frac{x^h(1-x)^{n-h}}{B (h+1, n-h+1)} = x^h(1-x)^{n-h} \frac{ \Gamma (n+2)}{\Gamma (h+1) \Gamma (n-h+1)} $$

Notice that the gamma functions are not a function of \(x\). They are only there as a normalization factor for a given \(n\) and \(h\). Therefore we can rewrite \(P(x, h, n)\) with a normalization integral.

$$ P(x, h, n) = \frac{ x^h(1-x)^{n-h} }{ \int_0^1 x^h(1-x)^{n-h} dx} $$

Measuring the prevalence

Let's say we did a series of three experiments to measure the disease prevalence of a population. We run the test on \(x_\text{pos}\) positive patients to measure the sensitivity, \(x_\text{neg}\) healthy people to measure the specificity, and \(x_\text{pop}\) random members of the public to measure the prevalence. Our goal is to produce a probability distribution of the population prevalence, \(P(P_\text{prev})\).

Measurement Variable Number of tests Number positive
Sensitivity \(S_n\) \(x_\text{pos}\) \( \lambda_\text{pos}\)
Specificity \(S_p\) \(x_\text{neg}\) \( \lambda_\text{neg}\)
Prevalence \(P_\text{prev}\) \(x_\text{pop}\) \( \lambda_\text{pop}\)

Each one of these is a weighted-coin flip following the equations in the "Consider the weighted coin problem backward" section at the top. Below, we will show them without the normalizing integrals or gamma functions.

$$ \begin{align} P(S_n, x_\text{pos}, \lambda_\text{pos}) & \propto S_n^{\lambda_\text{pos}} (1-S_n)^{x_\text{pos}-\lambda_\text{pos}} \\ P(S_p, x_\text{neg}, \lambda_\text{neg}) & \propto (1-S_p)^{\lambda_\text{neg} } S_p^{x_\text{neg}-\lambda_\text{neg}} \\ P(P_\text{prev}, x_\text{pop}, \lambda_\text{pop}) & \propto \left(P_\text{prev}S_n+(1-P_\text{prev})(1-S_p)\right)^{\lambda_\text{pop}} \left(1-P_\text{prev}S_n-(1-P_\text{prev})(1-S_p) \right)^{x_\text{pop}-\lambda_\text{pop}} \end{align} $$

For any given set of experiments (ie a fixed \(x_\text{pos}\), \(x_\text{neg}\), \(x_\text{pop}\), \( \lambda_\text{pos}\), \( \lambda_\text{neg}\), and \( \lambda_\text{pop}\)) there will be a probability distribution function \(P(P_\text{prev}, S_n, S_p)\) which is the product of \(P(S_n) \times P(S_p) \times P(P_\text{prev})\).

$$ P(P_\text{prev}, S_n, S_p) \propto \begin{align} & \left(P_\text{prev}S_n+(1-P_\text{prev})(1-S_p)\right)^{\lambda_\text{pop}} \\ \times & \left(1-P_\text{prev}S_n-(1-P_\text{prev})(1-S_p) \right)^{x_\text{pop}-\lambda_\text{pop}} \\ \times & S_p^{\lambda_\text{neg} } (1-S_p)^{x_\text{neg}-\lambda_\text{neg}} S_n^{\lambda_\text{pos}} (1-S_n)^{x_\text{pos}-\lambda_\text{pos}} \end{align} $$

To reduce this to a function of only \(P_\text{prev}\), take the double integral over \(S_n\) and \(S_p\) from 0 to 1. To normalize it, divide by the triple integral over \(P_\text{prev}\), \(S_n\) and \(S_p\) from 0 to 1.

$$ P(P_\text{prev}) =\frac{\left[ \begin{align}\int_0^1 \int_0^1 &\left(P_\text{prev}S_n+(1-P_\text{prev})(1-S_p)\right)^{\lambda_\text{pop}}\\ &\left(1-P_\text{prev}S_n-(1-P_\text{prev})(1-S_p) \right)^{x_\text{pop}-\lambda_\text{pop}}\\ &(1-S_p)^{\lambda_\text{neg} } S_p^{x_\text{neg}-\lambda_\text{neg}} \\ &S_n^{\lambda_\text{pos}} (1-S_n)^{x_\text{pos}-\lambda_\text{pos}} dS_n dS_p\end{align} \right]} {\left[ \begin{align}\int_0^1 \int_0^1 \int_0^1 &\left(P_\text{prev}S_n+(1-P_\text{prev})(1-S_p)\right)^{\lambda_\text{pop}}\\ &\left(1-P_\text{prev}S_n-(1-P_\text{prev})(1-S_p) \right)^{x_\text{pop}-\lambda_\text{pop}}\\ &(1-S_p)^{\lambda_\text{neg} } S_p^{x_\text{neg}-\lambda_\text{neg}} \\ &S_n^{\lambda_\text{pos}} (1-S_n)^{x_\text{pos}-\lambda_\text{pos}} dS_n dS_p dP_\text{prev} \end{align}\right]} $$

It turns out implementing this numerically or otherwise in javascript to make a cool interactive graph is not possible. It just crashes the browser. I'm also not aware of a good way to reduce the integrals in the numerator or denominator.

Reproducing Keith's results

As stated above, this article is heavily built on Keith Richards-Dinger's work with only some slight improvements. He did the data analysis in R. I believe the Python solution is more elegant. This code reproduces the marginal 2D and 1D plots near the bottom of Keith's article under "Including uncertainties in estimates of sensitivity." He was reanalyzing statistics from [1] using the following experimental results.

Measurement Variable Number of tests Number positive
Sensitivity \(S_n\) \(x_\text{pos} = 122\) \( \lambda_\text{pos} = 103\)
Specificity \(S_p\) \(x_\text{neg} = 401\) \( \lambda_\text{neg} = 2\)
Prevalence \(P_\text{prev}\) \(x_\text{pop} = 3300\) \( \lambda_\text{pop}=50\)


from matplotlib import cm

import numpy as np

import matplotlib.pyplot as plt

# measured results from population

xpop = 3300; lpop = 50

# measured results from known positive samples (sensitivity estimate)

xpos = 122; lpos = 103

# measured results from known negative samples (specificity estimate)

xneg = 401; lneg = 2

res = 300; #resolution of the computation (how many points?)

def joint_likely(pop, sens, spec):

p = (pop*sens + (1-pop)*(1-spec))**lpop * (1- pop*sens - (1-pop)*(1-spec))**(xpop-lpop)

pp = (1-spec)**lneg * (spec)**(xneg-lneg) * sens**lpos * (1-sens)**(xpos-lpos)

return p*pp

def confidence(x,y):

ycum = np.cumsum(y)/np.sum(y)

up = np.where(ycum > 0.9999999999)[0][0]

func = interp1d(ycum[0:up],x[0:up],'cubic');

return [func(0.025), func(0.975)]

def oneD(x,y,ax):

upper,lower = confidence(x,y)

print("95 confidence interval (%.5f,%.5f)"%(upper,lower))

ax.plot(x, y/y.sum()/(x[1]-x[0]), linewidth=2.5)

h = np.max(y/y.sum()/(x[1]-x[0]))

ax.plot([upper,lower], [h/10,h/10], linewidth=3.5, color='k')

ax.set_ylabel('Marginal PDF')

pop = np.linspace(0.00001, .04, res)

sens = np.linspace(0.601, 1, res)

spec = np.linspace(.96, 1, res)

popv, sensv, specv = np.meshgrid(pop, sens, spec)

full = joint_likely(popv, sensv, specv)

fig,axe = plt.subplots(2,3,figsize=(20,11))

ax = axe[0,0]

ax.pcolormesh(popv[1,:,:], specv[1,:,:], full.sum(axis=0), cmap=cm.viridis)

ax.set_xlabel('Prevalence'); ax.set_ylabel('Specificity')

ax = axe[0,1]

ax.pcolormesh( popv[1,:,:], sensv[:,1,:].T, full.sum(axis=2).T, cmap=cm.viridis)

ax.set_xlabel('Prevalence'); ax.set_ylabel('Sensitivity')

ax = axe[1,0]

ax.pcolormesh(specv[1,:,:],sensv[:,:,1], full.sum(axis=1), cmap=cm.viridis)

ax.set_xlabel('Specificity'); ax.set_ylabel('Sensitivity')

oneD(pop, full.sum(axis=0).sum(axis=1), axe[1,1] )


oneD(sens, full.sum(axis=2).sum(axis=1), axe[1,2] )


oneD(spec, full.sum(axis=0).sum(axis=0), axe[0,2] )




>>> exec(open('keith.py').read())

95 confidence interval (0.00126,0.01903) # prevalence

95 confidence interval (0.76685,0.89622) # sensitivity

95 confidence interval (0.98572,0.99842) # specificity

plots reproduced from Keith's work

[Caption] Recreation of Keith's 1D and 2D marginal PDF plots accounting for uncertainty in both sensitivity and specificity. These plots are mirror images of Keith's because he plotted false negative and postive rates while I plotted sensitivity and specificity.

I used the viridis color map because it is the most perceptually uniform color map I am aware of. Jet, the colormap used in Keith's images, is the default colormap in many applications. For more information, see my article on colormaps.


[1] E. Bendavid, B. Mulaney, N. Sood, S. Shah, E. Ling, R. Bromley-dulfano, C. Lai, Z. Weissberg, R. Saavedra, and J. Tedrow, "COVID-19 Antibody Seroprevalence in Santa Clara County, California," MedRxiv, 2020.

Follow @domesticengine7

© MC Byington