Tracking an Outbreak with Mark and Recapture

You can order print and ebook versions of Think Bayes 2e from Bookshop.org and Amazon.

Tracking an Outbreak with Mark and Recapture#

In epidemiology, the true number of people infected in an outbreak is often larger than the number of cases that appear on any official list. If we have several incomplete lists, and we can match individuals across them, the pattern of overlaps is informative: few repeats suggest many missing cases; many repeats suggest we have seen most of the cases.

The key to this approach is mark and recapture modeling, a method originally developed for ecology, but applicable in many other domains. In this excellent paper, Anne Chao reviews the approach and presents an example: a hepatitis A outbreak among students at a college in northern Taiwan in 1995.

Three incomplete lists of cases were available:

  1. P-list (135 cases): serum tests from the Institute of Preventive Medicine of Taiwan.

  2. Q-list (122 cases): cases reported by local hospitals to the National Quarantine Service.

  3. E-list (126 cases): questionnaires collected by epidemiologists.

Across the three lists, 271 distinct students were identified. The question is how many infected students were never listed.

In Chapter 15 of Think Bayes, I used this problem as an exercise. The solution there uses a discrete grid over two parameters: the total number of infections, N, and a shared probability, p, that any case appears on any given list. But this simple model is limited and, as we’ll see, it systematically underestimates the number of cases. Then, we’ll try a more realistic model and we’ll see that the results are better.

Click here to run this notebook on Colab.

The Data#

We can summarize the data by listing the number of cases that appeared on all three lists, denoted k111, the number that appear on each combination of lists, and the number that appear on none of the lists, k000.

P  Q  E   count
1  1  1   28
1  1  0   21
1  0  1   17
1  0  0   69
0  1  1   18
0  1  0   55
0  0  1   63
0  0  0   ??

Of course k000 is unknown. That’s what we’ll estimate!

Here’s the data in array form.

# k000 is unknown; the remaining seven cells are observed
observed = np.array([63, 55, 18, 69, 17, 21, 28])
num_seen = observed.sum()
num_seen

271

So 271 distinct cases appear on at least one list. Once we estimate k000, we can compute N = num_seen + k000.

Independence model#

We’ll start with what I’ll call the independence model, because it assumes that each infected person appears on list P, Q, or E with probabilities p_P, p_Q, and p_E, independently across people and lists.

With that assumption, the eight possible combinations have probabilities equal to various products of p_j and 1 - p_j (see the model below).

For the prior distribution of k000 we’ll use a Poisson distribution with a weakly informative hyperprior on the Poisson mean. The probabilities are drawn from a uniform distribution (Beta with parameters 1, 1).

Here’s the model.

coords = {"list": ["P", "Q", "E"]}

with pm.Model(coords=coords) as model:
    p = pm.Beta("p", alpha=1, beta=1, dims="list")
    p_P, p_Q, p_E = p[0], p[1], p[2]

    mu = pm.HalfNormal("mu", sigma=100)
    k000 = pm.Poisson("k000", mu=mu)
    N = pm.Deterministic("N", num_seen + k000)

    q_P, q_Q, q_E = 1 - p_P, 1 - p_Q, 1 - p_E
    probs = pm.math.stack(
        [
            q_P * q_Q * q_E,  # 000
            q_P * q_Q * p_E,  # 001
            q_P * p_Q * q_E,  # 010
            q_P * p_Q * p_E,  # 011
            p_P * q_Q * q_E,  # 100
            p_P * q_Q * p_E,  # 101
            p_P * p_Q * q_E,  # 110
            p_P * p_Q * p_E,  # 111
        ]
    )
    counts = pm.math.stack([k000, *observed])

    pm.Potential(
        "likelihood",
        pm.logp(pm.Multinomial.dist(n=N, p=probs), counts),
    )

As a technical note, you might wonder why we can’t write the likelihood like this:

likelihood = pm.Multinomial("likelihood", n=N, p=probs, observed=counts)

The problem is that the value we pass as observed has to be constant – it can’t be a random variable or contain any random variables – and counts contains the random variable k000. Instead we use Potential to compute the log likelihood of the data directly.

Here’s a graphical representation of the model.

pm.model_to_graphviz(model)
_images/1c393e81224fb0b081f71a24c63fe5ec73a1a842e7a332536df59ef4cca04f6f.svg

Notice that k000 affects the likelihood of the data in two ways: directly as one of the counts and indirectly as a term of N = num_seen + k000.

with model:
    idata = pm.sample(2000, tune=2000, random_seed=42, target_accept=0.9)

Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [p, mu]
>Metropolis: [k000]



Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 6 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details













az.summary(idata, var_names=["N", "p", "mu", "k000"])













  
    

      
      mean
      sd
      hdi_3%
      hdi_97%
      mcse_mean
      mcse_sd
      ess_bulk
      ess_tail
      r_hat
    

  
  
    

      N
      385.172
      20.810
      351.000
      425.000
      1.256
      0.913
      279.0
      431.0
      1.01
    

    

      p[P]
      0.352
      0.030
      0.296
      0.410
      0.001
      0.001
      642.0
      1371.0
      1.00
    

    

      p[Q]
      0.318
      0.029
      0.266
      0.376
      0.001
      0.000
      740.0
      1733.0
      1.00
    

    

      p[E]
      0.329
      0.029
      0.271
      0.382
      0.001
      0.001
      714.0
      1735.0
      1.00
    

    

      mu
      113.782
      22.871
      73.256
      158.240
      1.263
      0.821
      330.0
      605.0
      1.01
    

    

      k000
      114.172
      20.810
      80.000
      154.000
      1.256
      0.913
      279.0
      431.0
      1.01
    

  







The model samples reasonably well. Here are the posterior distributions.






az.plot_posterior(idata, var_names=["N", "p"])
decorate()








_images/b1dfc5b206cf8cfdbb9a8e295fa4f11ea0ba7492da32332b71f1e8d56e07c1a4.png





And here’s the estimated value of N with a 90% CI.






az.summary(idata, var_names=["N"], hdi_prob=0.9)













  
    

      
      mean
      sd
      hdi_5%
      hdi_95%
      mcse_mean
      mcse_sd
      ess_bulk
      ess_tail
      r_hat
    

  
  
    

      N
      385.172
      20.81
      352.0
      418.0
      1.256
      0.913
      279.0
      431.0
      1.01
    

  







The posterior mean, 381, is consistent with Chao’s independence estimator N̂0 and with the solution in Think Bayes.
Unfortunately, it turns out to be too low. A later campus-wide serum screen put the true number of infections around 545.


The independence model treats each infected student as equally likely to appear on a given list.
It also assumes that, after accounting for the different coverage of the three lists, appearing on one list provides no information about appearing on another.
Both assumptions are questionable.


In reality, some students are easier to detect than others — for example, because they had more severe symptoms or greater contact with medical services.
Such students would be more likely to appear on all three lists, while students with mild or asymptomatic infections might appear on none.


If that’s true, the data will contain more repeated observations than the independence model expects.
In that case, the independence model interprets the observed overlap as evidence that the lists have covered most of the population, and it tends to underestimate the number of infections.
And in this example, we can see that it did.


So let’s build a better model.






Individual heterogeneity#


We’ll assume that each person has a level of detectability, \(z_k\), and each list has a baseline detection level \(\alpha_P\), \(\alpha_Q\), \(\alpha_E\).
For a person with detectability \(z_k\), the chance they appear on list j is p_{kj}, where



\[\mathrm{logit}(p_{kj}) = \alpha_j + \sigma z_k\]


Here \(\sigma\) controls how much influence individual detectability has, compared to the baseline detection levels.
A person with high \(z_k\) is more likely to appear on more than one list; people with low detectability are more likely to be missed entirely.


To implement this model, we don’t need a latent variable for every person.
Instead we’ll define a grid of values for z and assign to each value a weight, w_k, that represents the fraction of people in the population with that value.
Then we can compute the cell probabilities using a weighted sum of the cell probabilities conditioned on z_k.



\[\pi_{abc} = \sum_k w_k \, P(abc \mid z_k)\]


We’ll precompute the values of z and the corresponding weights (using the PDF of a normal distribution to make a bell curve).






z = np.linspace(-4, 4, 21)
weights = norm.pdf(z)
weights /= weights.sum()









And precompute the pattern of detection/non-detection for each cell.






from itertools import product

patterns = np.array(list(product([0, 1], repeat=3)), dtype=float)
patterns









array([[0., 0., 0.],
       [0., 0., 1.],
       [0., 1., 0.],
       [0., 1., 1.],
       [1., 0., 0.],
       [1., 0., 1.],
       [1., 1., 0.],
       [1., 1., 1.]])









Now here’s a function that builds the model.
The argument sigma_prior_scale is the scale of the half-normal prior on \(\sigma\); we’ll vary it later to check sensitivity.






def make_heterogeneity_model(sigma_prior_scale=1):
    with pm.Model(coords=coords) as model:
        alpha = pm.Normal("alpha", mu=0, sigma=1.5, dims="list")
        sigma = pm.HalfNormal("sigma", sigma=sigma_prior_scale)

        # Probability someone with z=0 appears on each list
        p_typical = pm.Deterministic(
            "p_typical", pm.math.sigmoid(alpha), dims="list"
        )

        # Probability someone with each level of z appears on each list
        # shape: level × list
        logits = alpha[None, :] + sigma * z[:, None]
        p_z = pm.math.sigmoid(logits)

        # Population average probability of appearing on each list
        p_avg = pm.Deterministic(
            "p_avg",
            pm.math.sum(p_z * weights[:, None], axis=0),
            dims="list",
        )

        # Probability of each cell for each level of z and list
        # shape: cell × level × list
        y = patterns[:, None, :]
        cell_prob_given_z = pm.math.prod(
            p_z[None, :, :] ** y * (1 - p_z[None, :, :]) ** (1 - y),
            axis=2,
        )

        # Mixture over levels of z; shape: cell
        probs = pm.math.sum(cell_prob_given_z * weights[None, :], axis=1)
        probs = probs / pm.math.sum(probs)

        # The rest of this model is the same as the previous one
        mu = pm.HalfNormal("mu", sigma=100)
        k000 = pm.Poisson("k000", mu=mu)
        N = pm.Deterministic("N", num_seen + k000)

        counts = pm.math.stack([k000, *observed])
        pm.Potential(
            "likelihood",
            pm.logp(pm.Multinomial.dist(n=N, p=probs), counts),
        )
    return model













heterogeneity_model = make_heterogeneity_model()
pm.model_to_graphviz(heterogeneity_model)








_images/7577003e4202b6aa337ecb2d5a6f6f5c456c870058f63e4cd22ea6197e6f604e.svg









with heterogeneity_model:
    idata_het = pm.sample(2000, tune=2000, random_seed=42)









Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [alpha, sigma, mu]
>Metropolis: [k000]







Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 8 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details









This model doesn’t sample as well, but the results are stable enough for our purposes.






az.summary(
    idata_het,
    var_names=["N", "alpha", "sigma", "p_typical", "p_avg", "mu", "k000"],
)













  
    

      
      mean
      sd
      hdi_3%
      hdi_97%
      mcse_mean
      mcse_sd
      ess_bulk
      ess_tail
      r_hat
    

  
  
    

      N
      494.968
      54.767
      391.000
      589.000
      5.985
      4.051
      84.0
      184.0
      1.06
    

    

      alpha[P]
      -1.146
      0.258
      -1.659
      -0.702
      0.025
      0.014
      104.0
      227.0
      1.04
    

    

      alpha[Q]
      -1.303
      0.257
      -1.776
      -0.832
      0.024
      0.013
      111.0
      267.0
      1.04
    

    

      alpha[E]
      -1.254
      0.260
      -1.752
      -0.789
      0.025
      0.013
      108.0
      214.0
      1.04
    

    

      sigma
      0.935
      0.213
      0.539
      1.336
      0.018
      0.008
      138.0
      376.0
      1.03
    

    

      p_typical[P]
      0.244
      0.046
      0.159
      0.330
      0.005
      0.002
      104.0
      227.0
      1.04
    

    

      p_typical[Q]
      0.217
      0.043
      0.144
      0.303
      0.004
      0.002
      111.0
      267.0
      1.04
    

    

      p_typical[E]
      0.225
      0.044
      0.143
      0.306
      0.004
      0.002
      108.0
      214.0
      1.04
    

    

      p_avg[P]
      0.278
      0.036
      0.214
      0.346
      0.003
      0.002
      115.0
      270.0
      1.04
    

    

      p_avg[Q]
      0.251
      0.033
      0.191
      0.313
      0.003
      0.001
      129.0
      319.0
      1.03
    

    

      p_avg[E]
      0.259
      0.034
      0.196
      0.321
      0.003
      0.001
      122.0
      272.0
      1.04
    

    

      mu
      220.063
      54.572
      121.053
      318.643
      5.791
      3.848
      89.0
      194.0
      1.05
    

    

      k000
      223.968
      54.767
      120.000
      318.000
      5.985
      4.051
      84.0
      184.0
      1.06
    

  











az.plot_posterior(idata_het, var_names=["N", "sigma", "p_avg"])
decorate()








_images/251ad1bb08ea1f70c8ca7414ac08b2b6c40c4c4f5220fe1302bb46cd63b8e12a.png









az.summary(idata_het, var_names=["N"], hdi_prob=0.9)













  
    

      
      mean
      sd
      hdi_5%
      hdi_95%
      mcse_mean
      mcse_sd
      ess_bulk
      ess_tail
      r_hat
    

  
  
    

      N
      494.968
      54.767
      410.0
      584.0
      5.985
      4.051
      84.0
      184.0
      1.06
    

  







The mean is much closer to the true value, 545.
And the 90% CI contains the true value.


But the price we pay for a more flexible model is a wider CI.
Here are the posterior distributions of the two models (the dashed line is the true value).






az.plot_density(
    [idata, idata_het],
    var_names=["N"],
    data_labels=["Independence", "Heterogeneity"],
)
plt.axvline(545, ls='--', color='gray')
decorate(title="Posterior of N")








_images/afb5f3718ad57605afac99c01e6c4b01f4ae9cc100483baee96afd1b722b3ad0.png









summary_ind = az.summary(idata, var_names=["N"], hdi_prob=0.9)
summary_het = az.summary(idata_het, var_names=["N"], hdi_prob=0.9)
pd.DataFrame(
    {
        "mean": [summary_ind.loc["N", "mean"], summary_het.loc["N", "mean"]],
        "hdi_5%": [summary_ind.loc["N", "hdi_5%"], summary_het.loc["N", "hdi_5%"]],
        "hdi_95%": [summary_ind.loc["N", "hdi_95%"], summary_het.loc["N", "hdi_95%"]],
    },
    index=["Independence", "Heterogeneity"],
)













  
    

      
      mean
      hdi_5%
      hdi_95%
    

  
  
    

      Independence
      385.172
      352.0
      418.0
    

    

      Heterogeneity
      494.968
      410.0
      584.0
    

  











Sensitivity to Priors#


With only seven observed cells the baselines \(\alpha_j\), and \(\sigma\) can trade off.
Larger \(\sigma\) creates a subgroup with low probability of detection, which can support a larger unseen population.
So the upper tail of N may be prior-sensitive.


To check, we’ll test a few values of \(\sigma_prior_scale\), which controls the width of the prior distribution of \sigma.






def fit_heterogeneity(sigma_prior_scale, random_seed=42):
    model = make_heterogeneity_model(sigma_prior_scale)
    with model:
        idata_s = pm.sample(
            1000, tune=1000, random_seed=random_seed, progressbar=False
        )
    return idata_s













sensitivity = []
for scale in [0.5, 1.0, 2.0]:
    idata_s = fit_heterogeneity(scale)
    s = az.summary(idata_s, var_names=["N", "sigma"], hdi_prob=0.9)
    sensitivity.append(
        {
            "HalfNormal scale": scale,
            "N mean": s.loc["N", "mean"],
            "N 5%": s.loc["N", "hdi_5%"],
            "N 95%": s.loc["N", "hdi_95%"],
            "sigma mean": s.loc["sigma", "mean"],
        }
    )

pd.DataFrame(sensitivity).round(1)









Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [alpha, sigma, mu]
>Metropolis: [k000]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [alpha, sigma, mu]
>Metropolis: [k000]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [alpha, sigma, mu]
>Metropolis: [k000]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details









  
    

      
      HalfNormal scale
      N mean
      N 5%
      N 95%
      sigma mean
    

  
  
    

      0
      0.5
      476.9
      382.0
      558.0
      0.8
    

    

      1
      1.0
      494.4
      403.0
      572.0
      0.9
    

    

      2
      2.0
      493.3
      383.0
      571.0
      0.9
    

  







If the prior distribution of sigma is too narrow, it restrains the effect of individual variation, which biases the estimate of N downward.


But once it’s large enough, it doesn’t have much effect on the results.






Discussion#


This example demonstrates a feature of Bayesian modeling with PyMC: the ability to develop and test a model incrementally, adding features as needed and comparing successive iterations.


It also demonstrates the impact of AI tools on the development process for Bayesian models. The models in this example are moderately complex. Both of them use a Potential to compute a custom likelihood function, and the second involves array operations that might take some time to write, test, and debug. Using various AI tools, I was able to get both models working quickly – which allowed time to inspect, test, and refine them.
And then I was able to test a few variations to make the models simpler and easier to explain.


I’m not sure I would say it was easy – this kind of work still requires a fundamental understanding of Bayesian statistics and PyMC. But it is a lot easier than it used to be!


Copyright 2025 Allen B. Downey


License: Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)