Counting the uncounted: how many victims did every list miss?

import os

for v_ in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v_] = "1"  # nutpie runs chains as threads: a BLAS pool per chain would fight over the cores

import itertools as it
import logging
import sys
import time
import warnings

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import nutpie
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
from IPython.display import Markdown, display
from matplotlib.patches import Rectangle
from scipy.optimize import minimize
from scipy.special import gammaln, logsumexp

from pymc_challenges import data
from pymc_challenges import public as pub

RANDOM_SEED = 109
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")
logging.getLogger("pymc").setLevel(logging.WARNING)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning, module="arviz")
warnings.filterwarnings("ignore", message="There are no gridspecs with layoutgrids")
sys.modules.pop("jax", None)  # arviz leaves a lazy jax stub that makes matplotlib import all of JAX
pd.set_option("display.width", 160)
T_START = time.time()
TIMES = {}
INK, MUTED, LIGHT, HI, LO = pub.INK, pub.MUTED, pub.LIGHT, pub.HI, pub.LO
GREY = "#555555"
print(f"PyMC {pm.__version__}, ArviZ {az.__version__}, nutpie {nutpie.__version__}")

PyMC 6.3.2, ArviZ 1.3.0, nutpie 0.16.11

#In plain words

The question. After the killings in Kosovo in 1999, four groups each made a list of the dead, and together they named 4,400 people. How many of the dead are on no list at all?

What we did. The overlaps carry the answer: if every list found almost the same people, few were missed; if each found mostly people the others did not, many were. But some people were easier for everyone to find, so the answer depends on how the lists lean on each other. We fitted two very different models of that dependence and weighted them by how well each predicts overlaps it has not seen.

What we found. For every 10 people named on at least one list, about 13 more were on none (somewhere between 10 and 18). That puts the total near 10,000, close to the published estimate of 10,356.

How sure are we. It is very likely (95%) that at least 4,300 of the people killed were on no list. We have medium confidence in that figure: the two models put the total at about 9,300 and 10,400, within about a tenth of each other, but the answer rests on how the lists depend on each other, which the data can only partly check.

HEADLINE = display(Markdown("*The headline figure is drawn in part G from the fits below and "
                            "appears here once the notebook has run.*"), display_id=True)

#The setting

Between March and June 1999 Yugoslav forces and the Kosovo Liberation Army fought in Kosovo, NATO bombed, and most of the Kosovar Albanian population fled. Afterwards several organisations recorded the dead: forensic teams exhuming graves (EXH), the American Bar Association's Central and East European Law Initiative interviewing refugees (ABA), the OSCE's human-rights mission (OSCE) and Human Rights Watch (HRW). Each list is a record of what one organisation could reach. Ball, Betts, Scheuren, Dudukovich & Asher (2002) matched the four lists name by name - about 10,000 reports came down to 4,400 distinct people - and asked the question this notebook asks: how many people were killed whom nobody recorded? Their answer, 10,356 (95% interval 9,002-12,122), was presented as evidence at the International Criminal Tribunal for the former Yugoslavia in the trial of Slobodan Milosevic, together with the timing of the killings and the refugee flows, which did not follow the pattern of NATO bombing or KLA activity.

The method is multiple-systems estimation (MSE), the human-rights name for what ecologists call capture-recapture (E28 counted hares with it). Every person has a capture history, a 0/1 string saying which lists name them. With $K$ lists there are $2^K - 1$ observable histories and one cell, "on no list", that is empty by construction. Fienberg (1972) showed how a log-linear model of the observed cells predicts that missing cell. The prediction rests entirely on how the lists depend on each other, and Lum, Price & Banks (2013) explain why that is the central difficulty in human-rights work: lists are made by organisations that talk to each other, reach the same villages and miss the same people. Two Bayesian answers stand out. Madigan & York (1997) average over the graphical models that say which lists depend on which; Manrique-Vallier (2016) lets people fall into hidden classes with their own visibility, which is how the Human Rights Data Analysis Group now works (the LCMCR package). Silverman (2020) compared both on the UK's modern-slavery lists and got answers a factor of two apart; Chan, Silverman & Vincent (2021) studied what goes wrong when some lists never overlap at all.

We work through both datasets, Kosovo first, because it is the cleaner case: four lists, every combination seen, and a published answer to compare with.

part	question	tool
A	What do the lists look like?	the 15 observed histories, list sizes and overlaps
B	Why does overlap measure the unseen, and how does it fail?	Lincoln-Petersen; a simulation with dependent lists; the six two-list estimates
C	How do log-linear models fill the empty cell?	independence, pairwise and three-way models; a regularised horseshoe; cell-by-cell checks; which prior family
D	What if people differ in how visible they are?	the latent-class model, summed over classes; label switching
E	Which assumption should we believe?	averaging over 113 hierarchical log-linear models; stacking with exact leave-one-cell-out; the headline mixture; "what you would have to believe"
F	Does it work with six sparse lists?	UK 2013: the Home Office model, horseshoe against latent class, the general-public list
G	How do we explain it to everyone?	the iceberg, a HOP GIF, a dotplot column, a 100-world swarm; the statistician's plot for contrast

#A · Four lists and the people on none

The file holds one row per observed capture history for five datasets, with the list names in the order of the history string.

data.describe("e109_mse_lists")
MSE = data.load("e109_mse_lists")
print(MSE.groupby("dataset").agg(lists=("lists", "first"), histories=("count", "size"),
                                 people=("count", "sum")).to_string())


def table(ds, drop=(), merge=None):
    """The capture-history table of one dataset: list names, all 2^K - 1 non-empty histories (rows of 0/1),
    their counts (zero where nobody had that history) and the history strings. `drop` removes lists,
    `merge=(a, b, new)` merges two lists into one."""
    g = MSE[MSE.dataset == ds]
    lists = g.lists.iloc[0].split("|")
    d = pd.DataFrame([list(map(int, h)) for h in g.history], columns=lists)
    d["n"] = g["count"].to_numpy()
    keep = [l for l in lists if l not in drop]
    if merge:
        a, b, new = merge
        d[new] = np.maximum(d[a], d[b])
        keep = [new if l == a else l for l in keep if l != b]
    cnt = d.groupby(keep)["n"].sum()
    K = len(keep)
    H = np.array(list(it.product([0, 1], repeat=K)))[1:]
    n = cnt.reindex([tuple(h) for h in H]).fillna(0).astype(int).to_numpy()
    return keep, H, n, ["".join(map(str, h)) for h in H]


LISTS, H, n, KEY = table("kosovo")
NOBS = int(n.sum())
K = len(LISTS)
sizes = {l: int(n[H[:, k] == 1].sum()) for k, l in enumerate(LISTS)}
print(f"\nKosovo: {NOBS:,} people on {K} lists; list sizes {sizes}")
print("people on exactly 1 / 2 / 3 / 4 lists:", [int(n[H.sum(1) == j].sum()) for j in range(1, 5)])

e109_mse_lists: One row per OBSERVED capture history (the all-zero history is never observed): dataset (kosovo, uk2013, netherlands, new_orleans, western_us), lists (list names joined by '|', in history order), history (0/1 string, one character per list; read as a string), count (people with that history). Kosovo 15 rows / 4,400 people; UK 25 / 2,744; Netherlands 24 / 8,234; New Orleans 19 / 185; Western 13 / 345. Histories with zero people are omitted.
  source: Capture-history counts from two CRAN packages (GPL >= 2), the data from published reports. Kosovo: LCMCR 0.4.14 (Manrique-Vallier) `kosovo_aggregate`, 4,400 killings in Kosovo, 20 March - 22 June 1999, on four lists - EXH (exhumations), ABA (American Bar Association / CEELI interviews), OSCE, HRW (Human Rights Watch) - analysed by Ball, Betts, Scheuren, Dudukovich & Asher (2002, AAAS/ABA-CEELI report to the ICTY). UK 2013: SparseMSE 2.0.1 (Chan, Silverman & Vincent) `UKdat`, 2,744 potential victims of trafficking in the National Crime Agency's 2013 strategic assessment on six lists (LA local authorities, NG NGOs, PF police forces, GO government organisations, GP general public, NCA National Crime Agency), the basis of the Home Office's 2014 estimate (Silverman 2014). Also from SparseMSE: the Netherlands (`Ned`), Greater New Orleans (`NewOrl`, anonymised lists) and a US western site (`Western`, anonymised). DERIVED by tools/build_e109_mse.py from the two tarballs (the url is the first; SparseMSE_2.0.1.tar.gz is the second) - keep the cached file in data/.
  url:    https://cran.r-project.org/src/contrib/LCMCR_0.4.14.tar.gz
                          lists  histories  people
dataset                                           
kosovo         EXH|ABA|OSCE|HRW         15    4400
netherlands         I|K|O|P|R|Z         24    8234
new_orleans     A|B|C|D|E|F|G|H         19     185
uk2013       LA|NG|PF|GO|GP|NCA         25    2744
western_us            A|B|C|D|E         13     345

Kosovo: 4,400 people on 4 lists; list sizes {'EXH': 1910, 'ABA': 1528, 'OSCE': 1786, 'HRW': 685}
people on exactly 1 / 2 / 3 / 4 lists: [3218, 882, 273, 27]

Four lists of different sizes: the exhumations name about 1,900 people, the OSCE about 1,800, the ABA interviews about 1,500 and Human Rights Watch about 700. Most people, about 3,200 of the 4,400, appear on one list only, and only 27 appear on all four. The chart shows every history, largest first, with the dots below saying which lists it involves. On the far left is the column every MSE analysis is about: people on no list, whose height is unknown.

order_c = np.argsort(-n)


def overlap_chart(ax_bar, ax_dot, hidden=None, title=""):
    """UpSet-style chart of the observed histories, with the never-observed history in column 0:
    a question mark, or a 20-dot quantile column of the hidden count."""
    xs = np.arange(len(n) + 1)
    ax_bar.bar(xs[1:], n[order_c], color=GREY, width=0.72)
    for x, v in zip(xs[1:], n[order_c]):
        ax_bar.text(x, v + 25, f"{v:,}", ha="center", va="bottom", fontsize=8.5, color=INK)
    if hidden is None:
        ax_bar.bar([0], [n.max()], color="none", edgecolor=HI, ls=(0, (3, 2)), width=0.72)
        ax_bar.text(0, n.max() / 2, "?", ha="center", va="center", fontsize=26, color=HI)
        top = n.max() * 1.15
    else:
        q = np.quantile(hidden, (np.arange(20) + 0.5) / 20)
        ax_bar.scatter(np.zeros(20), q, s=26, color=HI, zorder=3, clip_on=False)
        ax_bar.text(0.55, q.max(), "on no list:\n20 equally likely\nvalues", ha="left", va="top",
                    fontsize=9.5, color=HI)
        top = q.max() * 1.05
    ax_bar.set_ylim(0, top)
    ax_bar.set_ylabel("people")
    ax_bar.set_title(title, loc="left")
    for j, l in enumerate(LISTS):
        for x in xs:
            h = np.zeros(K, int) if x == 0 else H[order_c[x - 1]]
            ax_dot.scatter(x, K - 1 - j, s=42, color=INK if h[j] else LIGHT, zorder=3)
        ax_dot.text(-0.9, K - 1 - j, l, ha="right", va="center", fontsize=10)
    for x in xs[1:]:
        on = np.where(H[order_c[x - 1]] == 1)[0]
        ax_dot.plot([x, x], [K - 1 - on.max(), K - 1 - on.min()], color=INK, lw=1.5, zorder=2)
    ax_dot.set_xlim(-1.5, len(n) + 0.6)
    ax_dot.set_ylim(-0.6, K - 0.4)
    ax_dot.axis("off")


fig, (a1, a2) = plt.subplots(2, 1, figsize=(10, 5.6), sharex=True, gridspec_kw={"height_ratios": [3, 1.2]})
overlap_chart(a1, a2, title="The 15 overlaps of the four Kosovo lists, and the one nobody can see")
a1.set_xlim(-1.5, len(n) + 0.6);

#B · Two lists, and how they mislead

#Lincoln-Petersen

Take two lists, A and B, of sizes $n_A$ and $n_B$, and let $m$ people be on both. If the lists were drawn independently from a population of $N$, the chance of being on A is $n_A/N$ whether or not one is on B, so the share of B's people who are also on A, $m/n_B$, estimates $n_A/N$:

$$\hat N = \frac{n_A\, n_B}{m}.$$

This is the Lincoln-Petersen estimator. Its weak point is the word independently. If some people are easy for every organisation to find - a village near a road, a family that went to the refugee camps early - the lists overlap more than chance, $m$ is too big, and $\hat N$ is too small. If the lists divide the territory between them, they overlap less than chance and $\hat N$ is too big. The simulation below keeps the two list sizes the same in all three worlds (30% and 40% of 10,000 people) and changes only how they depend on each other.

N_TOY, reps = 10_000, 2000
worlds = {  # (share of people, P(on A), P(on B)) for two groups; same list sizes in every world
    "independent": [(1.0, 0.30, 0.40)],
    "same people easy to find\n(positive dependence)": [(0.5, 0.50, 0.60), (0.5, 0.10, 0.20)],
    "lists split the territory\n(negative dependence)": [(0.5, 0.50, 0.20), (0.5, 0.10, 0.60)],
}
toy = {}
for w, groups in worlds.items():
    est = np.zeros(reps)
    for r in range(reps):
        nA = nB = m = 0
        for share, pa, pb in groups:
            size = int(share * N_TOY)
            a = rng.random(size) < pa
            b = rng.random(size) < pb
            nA, nB, m = nA + a.sum(), nB + b.sum(), m + (a & b).sum()
        est[r] = nA * nB / m
    toy[w] = est
    print(f"{w.splitlines()[0]:28s} median estimate {np.median(est):7,.0f}  (89%: "
          f"{np.quantile(est, 0.055):,.0f}-{np.quantile(est, 0.945):,.0f})  truth {N_TOY:,}")

fig, ax = plt.subplots(figsize=(10, 3.4))
for (w, est), col in zip(toy.items(), [MUTED, HI, LO]):
    ax.hist(est, bins=np.arange(6000, 17000, 100), color=col, alpha=0.75, label=w)
    peak = np.histogram(est, bins=np.arange(6000, 17000, 100))[0].max()
    ax.text(np.median(est), peak * 1.03, w, ha="center", va="bottom", fontsize=10, color=col)
ax.axvline(N_TOY, color=INK, lw=1)
ax.set_xlabel("Lincoln-Petersen estimate of a population of 10,000")
ax.set_yticks([])
ax.set_ylim(0, max(np.histogram(e, bins=np.arange(6000, 17000, 100))[0].max() for e in toy.values()) * 1.35)
ax.set_title("Same list sizes, three kinds of dependence: only independent lists give the truth",
             loc="left");

independent                  median estimate  10,001  (89%: 9,710-10,310)  truth 10,000
same people easy to find     median estimate   7,501  (89%: 7,334-7,670)  truth 10,000

lists split the territory    median estimate  15,008  (89%: 14,370-15,618)  truth 10,000

With independent lists the estimate centres on the truth. When the same people are easy for both lists to find, it centres near 7,500, a quarter too low, and no amount of data fixes it: the bias is in the assumption. When the lists split the territory it centres near 15,000. The lists' sizes are identical in all three worlds; only their dependence differs, and two lists alone cannot tell us which world we are in.

#The six two-list estimates for Kosovo

Ball et al. (2002, appendix 2) printed the Lincoln-Petersen estimate for each of the six pairs of Kosovo lists. Recomputing them is also a check that the table we have is the one they used.

BALL_LP = {("EXH", "ABA"): 7245, ("ABA", "HRW"): 9689, ("ABA", "OSCE"): 5970, ("EXH", "HRW"): 6777,
           ("EXH", "OSCE"): 7135, ("OSCE", "HRW"): 5461}
rows = []
for a, b in it.combinations(range(K), 2):
    nA, nB = n[H[:, a] == 1].sum(), n[H[:, b] == 1].sum()
    m = n[(H[:, a] == 1) & (H[:, b] == 1)].sum()
    rows.append(dict(pair=f"{LISTS[a]}-{LISTS[b]}", n_A=nA, n_B=nB, both=m, estimate=round(nA * nB / m),
                     ball_2002=BALL_LP[(LISTS[a], LISTS[b])]))
print(pd.DataFrame(rows).to_string(index=False))

    pair  n_A  n_B  both  estimate  ball_2002
 EXH-ABA 1910 1528   403      7242       7245
EXH-OSCE 1910 1786   478      7137       7135
 EXH-HRW 1910  685   193      6779       6777
ABA-OSCE 1528 1786   457      5972       5970
 ABA-HRW 1528  685   108      9691       9689
OSCE-HRW 1786  685   224      5462       5461

Every estimate agrees with the published table to within 3 people, so the CRAN table is the one the report used. The six estimates run from about 5,500 to 9,700, and all of them are below the 10,356 that Ball et al. finally reported - most of them far below. That pattern is what positive dependence looks like: each pair of lists overlaps more than independent lists would, so each pair concludes that few people were missed. Ball et al. read it the same way ("a great deal of positive dependence between the lists"). With four lists we can do better than any pair, because some of that dependence shows up in the other lists.

#C · Log-linear models of the incomplete table

#The model

Write $n_h$ for the number of people with history $h \in \{0,1\}^4$, and let the counts be Poisson. A log-linear model writes the log of each cell's expected count as a sum of terms:

$$\log \mu_h = \lambda_0 + \sum_k \lambda_k h_k + \sum_{k<l} \lambda_{kl} h_k h_l + \sum_{k<l<m} \lambda_{klm} h_k h_l h_m .$$

The empty history $h = 0000$ switches every term off but $\lambda_0$, so the expected number of people on no list is $e^{\lambda_0}$, and the total is $N = n_{\text{obs}} + n_{0000}$ with $n_{0000} \sim \text{Poisson}(e^{\lambda_0})$. The main effects $\lambda_k$ are the log odds of being on list $k$. The pair terms $\lambda_{kl}$ are log odds ratios: $e^{\lambda_{kl}} = 2$ means being on list $k$ doubles the odds of being on list $l$. Without interactions the model is independence, and every interaction is an assumed dependence. The four-way term cannot be in the model at all: it would need the missing cell to estimate it, which is the untestable assumption at the heart of MSE. Everything in this part assumes it is zero.

Two lists give three observed cells, enough for $\lambda_0$, $\lambda_A$ and $\lambda_B$ and nothing more, which is why Lincoln-Petersen has to assume independence. Four lists give 15 cells: room for all 6 pair terms and all 4 three-way terms (15 parameters, a saturated fit). The question is how many of them to believe.

#Classical fits, for reference

Maximum likelihood with no prior is three lines of NumPy, and gives three quite different answers depending on which terms go in.

def design(Hm, terms):
    cols = [np.ones(len(Hm))] + [Hm[:, k] for k in range(Hm.shape[1])]
    return np.column_stack(cols + [np.prod(Hm[:, list(t)], axis=1) for t in terms]).astype(float)


def fit_mle(Hm, nm, terms):
    """Poisson log-linear maximum likelihood: returns the coefficients, the hidden count exp(l0) and the deviance."""
    X = design(Hm, terms)
    b0 = np.zeros(X.shape[1])
    b0[0] = np.log(nm.mean())
    r = minimize(lambda b: np.sum(np.exp(X @ b) - nm * (X @ b)), b0, jac=lambda b: X.T @ (np.exp(X @ b) - nm),
                 method="L-BFGS-B", options={"maxiter": 5000})
    mu = np.exp(X @ r.x)
    dev = 2 * np.sum(np.where(nm > 0, nm * np.log(np.maximum(nm, 1) / mu), 0) - (nm - mu))
    return r.x, np.exp(r.x[0]), dev, len(nm) - X.shape[1]


PAIRS = list(it.combinations(range(K), 2))
TRIPLES = list(it.combinations(range(K), 3))
for label, terms in [("independence", []), ("all pairs", PAIRS), ("all pairs + all triples", PAIRS + TRIPLES)]:
    b, N0, dev, df = fit_mle(H, n, terms)
    print(f"MLE {label:24s}: total {NOBS + N0:7,.0f}   deviance {dev:6.1f} on {df:2d} df")

MLE independence            : total   7,395   deviance  245.9 on 10 df
MLE all pairs               : total  14,410   deviance   39.8 on  4 df
MLE all pairs + all triples : total  16,951   deviance    0.0 on  0 df

Independence says about 7,400 and fits terribly (a deviance of about 250 on 10 degrees of freedom, where a model that fits would give about 10). All pairs says about 14,400 and still fits badly. Adding the three-way terms fits the 15 cells exactly and says 17,000. The saturated model is not wrong to fit perfectly - it has as many parameters as cells - but that leaves nothing over to check its estimate of the hidden cell against. We will come back to it.

#Priors

The Bayesian version needs a prior for each kind of term. For $\lambda_0$ we centre the hidden count on the number listed with a wide spread, $\lambda_0 \sim N(\log 4400, 1.5)$, which puts 89% of the prior between about 400 and 50,000 hidden people. Centring on the observed total uses the data once more, in the empirical-Bayes way. For the main effects, $\lambda_k \sim N(-1, 1.5)$: a list catches somewhere between a few percent and most of the people. Part C's Failure 2 widens both of these priors and shows when they matter. The interactions get a horseshoe (Carvalho, Polson & Scott 2010): each term has its own scale $c_j \sim \text{HalfCauchy}(1)$ times a global $\tau \sim \text{HalfNormal}(0.5)$ for the pairs and $\text{HalfNormal}(0.25)$ for the triples. The horseshoe says that most dependencies are probably small and a few may be large, which is how analysts of these lists describe them, and it replaces the stepwise search over models that classical MSE uses with shrinkage.

The code builds all the log-linear models from one function. orders=() is independence; slab turns the horseshoe into the regularised version we will need in a moment; ridge swaps the horseshoe for a plain normal prior, for a comparison; fix_pair and masked are for the sensitivity and the leave-one-out refits later on.

def loglinear(lists, Hm, nm, *, orders=(), tau_sd=(), slab=2.0, ridge=None, l0_sd=1.5, lk_sd=1.5, fix_pair=None,
              masked=False):
    """Poisson log-linear MSE. exp(l0) is the expected number on no list. Interactions of each order get a
    (regularised, when `slab` is set) horseshoe written non-centred, or N(0, ridge^2) when `ridge` is set."""
    Kl = len(lists)
    with pm.Model(coords={"cell": ["".join(map(str, h)) for h in Hm], "list": lists}) as m:
        l0 = pm.Normal("l0", np.log(nm.sum()), l0_sd)
        lk = pm.Normal("lk", -1.0, lk_sd, dims="list")
        eta = l0 + pt.dot(Hm.astype(float), lk)
        for order, tsd in zip(orders, tau_sd):
            terms = list(it.combinations(range(Kl), order))
            X = np.column_stack([np.prod(Hm[:, list(t)], axis=1) for t in terms]).astype(float)
            dim = f"term{order}"
            m.add_coord(dim, ["-".join(lists[i] for i in t) for t in terms])
            if ridge is not None:
                scale = ridge
            else:
                tau = pm.HalfNormal(f"tau{order}", tsd)
                scale = tau * pm.HalfCauchy(f"c{order}", 1.0, dims=dim)
                if slab is not None:  # Piironen & Vehtari (2017): large terms see a Normal(0, slab) instead
                    scale = slab * scale / pt.sqrt(slab**2 + scale**2)
            lam = scale * pm.Normal(f"z{order}", 0.0, 1.0, dims=dim)
            if fix_pair is not None and order == 2:  # one dependence held at an assumed value
                lam = pt.set_subtensor(lam[terms.index(fix_pair)], pm.Data("fixed", 0.0))
            lam = pm.Deterministic(f"lam{order}", lam, dims=dim)
            eta = eta + pt.dot(X, lam)
        pm.Deterministic("N0_mean", pt.exp(l0))
        likelihood(pm.Deterministic("mu", pt.exp(eta), dims="cell"), nm, masked)
    return m


def likelihood(mu, nm, masked):
    if masked:  # exact leave-one-cell-out: a 0/1 mask swapped into one compiled model
        mask = pm.Data("mask", np.ones(len(nm)), dims="cell")
        pm.Potential("ll", pt.sum(mask * pm.logp(pm.Poisson.dist(mu), nm)))
    else:
        pm.Poisson("n", mu, observed=nm, dims="cell")


m_prior = loglinear(LISTS, H, n, orders=(2, 3), tau_sd=(0.5, 0.25))
with m_prior:
    prior = pm.sample_prior_predictive(2000, random_seed=RANDOM_SEED)
hid_prior = prior.prior["N0_mean"].to_numpy().ravel()
obs_prior = prior.prior_predictive["n"].to_numpy().reshape(-1, len(n)).sum(1)
print("prior: hidden count 5.5/50/94.5%:", np.round(np.quantile(hid_prior, [0.055, 0.5, 0.945]), -1))
print("prior: total listed 5.5/50/94.5%:", np.round(np.quantile(obs_prior, [0.055, 0.5, 0.945]), -1),
      f"(observed {NOBS:,})")
fig, axs = plt.subplots(1, 2, figsize=(10, 3))
for ax, v, lab in [(axs[0], hid_prior, "expected number on no list"), (axs[1], obs_prior, "people on any list")]:
    ax.hist(np.log10(np.maximum(v, 1)), bins=60, color=MUTED)
    ax.set_xlabel(f"{lab} (log10)")
    ax.set_yticks([])
axs[1].axvline(np.log10(NOBS), color=HI, lw=2)
axs[1].text(np.log10(NOBS), axs[1].get_ylim()[1] * 0.9, " observed 4,400", color=HI)
axs[0].set_title("Prior predictive: the model allows anything from hundreds to tens of thousands",
                 loc="left");

prior: hidden count 5.5/50/94.5%: [  380.  4410. 49470.]
prior: total listed 5.5/50/94.5%: [   1480.   30090. 1168620.] (observed 4,400)

The prior allows hidden counts from about 400 to 50,000 (89%), and listed totals from about 1,500 to well over a million. The observed 4,400 sits in the lower part of that range (the prior median is about 30,000, because the main effects allow lists that catch most people) but well inside it. The prior is wide and rules out nothing a reader could believe before seeing the overlaps; the 15 cells will do the work.

#Failure 1: the horseshoe's funnel

We fit the pairwise model first with the textbook horseshoe and nutpie's defaults, then change one thing at a time: a higher target_accept, and the regularised horseshoe described below.

def fit(m, ta=0.8, **kw):
    global LAST_MODEL
    LAST_MODEL = m  # summarise() needs it for the log-likelihood
    t0 = time.time()
    with m:
        idata = pm.sample(random_seed=RANDOM_SEED, progressbar=False, target_accept=ta, **kw)
    return idata, time.time() - t0


def worst_rhat(idata, skip=()):
    r = {v: float(az.rhat(idata, var_names=[v])[v].max()) for v in idata.posterior.data_vars if v not in skip}
    v = max(r, key=r.get)
    return f"{r[v]:.3f} ({v})"


for lab, slab, ta in [("plain horseshoe", None, 0.8), ("plain horseshoe", None, 0.99),
                      ("regularised horseshoe", 2.0, 0.8)]:
    idata_hs, sec = fit(loglinear(LISTS, H, n, orders=(2,), tau_sd=(0.5,), slab=slab), ta)
    print(f"{lab:22s} target_accept {ta}: {sec:.0f} s, sampled parameters "
          f"{sum(idata_hs.posterior[v][0, 0].size for v in ('l0', 'lk', 'tau2', 'c2', 'z2'))}, divergences "
          f"{int(idata_hs.sample_stats['diverging'].sum())}, worst r_hat {worst_rhat(idata_hs)}")
del idata_hs

plain horseshoe        target_accept 0.8: 4 s, sampled parameters 18, divergences 92, worst r_hat 1.013 (c2)

plain horseshoe        target_accept 0.99: 3 s, sampled parameters 18, divergences 9, worst r_hat 1.003 (z2)

regularised horseshoe  target_accept 0.8: 2 s, sampled parameters 18, divergences 0, worst r_hat 1.004 (c2)

About ninety divergences, from a model with eighteen sampled parameters (the intercept, four main effects, the global scale, and six local scales and six standardised terms for the pairs). The horseshoe's half-Cauchy scales let a term's prior sd reach almost zero, where its coefficient must be almost exactly zero, and also very large values; NUTS cannot find one step size for both ends of that funnel. Raising target_accept to 0.99 makes the steps smaller and cuts the count to about ten, but does not remove them. What does is the regularised horseshoe of Piironen & Vehtari (2017): the scale is passed through $s \mapsto \text{slab}\cdot s/\sqrt{\text{slab}^2 + s^2}$, so very large terms see a $N(0, 2^2)$ "slab" rather than a Cauchy tail. With it the default target_accept already gives no divergences. A log odds ratio of 2 is a factor of 7 in the odds, already a strong dependence, so the slab costs nothing we believe. We keep target_accept=0.99 for the models with interaction terms below anyway, as cheap insurance on these small models.

#Three models, fitted

We fit independence (M0), pairs (M1) and pairs plus triples (M1b), all with the regularised horseshoe, keep the posterior of the hidden count and the expected counts of every cell, and compute PSIS leave-one-out for later.

def summarise(name, idata, nm, sec, skip=()):
    """Diagnostics, draws of the total N (observed + Poisson hidden count), cell means and PSIS-LOO."""
    nobs = int(nm.sum())
    N = nobs + rng.poisson(idata.posterior["N0_mean"].to_numpy().ravel())
    mu = idata.posterior["mu"].to_numpy().reshape(-1, len(nm))
    yrep = rng.poisson(mu)
    lo, hi = np.quantile(yrep, [0.055, 0.945], axis=0)
    miss = [(KEYS[i], int(nm[i]), int(lo[i]), int(hi[i])) for i in range(len(nm)) if not lo[i] <= nm[i] <= hi[i]]
    q = np.quantile(N, [0.055, 0.5, 0.945])
    print(f"{name}: {sec:.0f} s, tune {idata.posterior.attrs.get('tuning_steps')}, divergences "
          f"{int(idata.sample_stats['diverging'].sum())}, r_hat(N0) "
          f"{float(az.rhat(idata, var_names=['N0_mean'])['N0_mean']):.3f}, worst r_hat {worst_rhat(idata, skip)}, "
          f"ESS(N0) {float(az.ess(idata, var_names=['N0_mean'])['N0_mean']):.0f}")
    print(f"   total N: {q[1]:,.0f}  (89%: {q[0]:,.0f}-{q[2]:,.0f});  cells outside their 89% predictive "
          f"range: {len(miss)} of {len(nm)}  {miss if len(miss) <= 5 else ''}")
    out = dict(N=N, yrep=yrep)
    if "n" in idata.observed_data:
        pm.compute_log_likelihood(idata, model=LAST_MODEL, progressbar=False)
        loo = az.loo(idata, pointwise=True)
        out["k_bad"] = int((loo.pareto_k > 0.7).sum())
        print(f"   PSIS-LOO: {out['k_bad']} of {len(nm)} cells with Pareto k > 0.7")
    return out


KEYS = KEY
KOS = {}
SPEC = {  # name: (model builder, target_accept)
    "M0 independence": (lambda **kw: loglinear(LISTS, H, n, **kw), 0.8),
    "M1 pairs": (lambda **kw: loglinear(LISTS, H, n, orders=(2,), tau_sd=(0.5,), **kw), 0.99),
    "M1b pairs + triples": (lambda **kw: loglinear(LISTS, H, n, orders=(2, 3), tau_sd=(0.5, 0.25), **kw), 0.99),
}
t0 = time.time()
for name, (mk, ta) in SPEC.items():
    idata, sec = fit(mk(), ta)
    KOS[name] = summarise(name, idata, n, sec)
    if name.startswith("M1b"):
        lam2 = idata.posterior["lam2"]
        lam3 = idata.posterior["lam3"]
        ABA_HRW_FREE = np.exp(lam2.sel(term2="ABA-HRW").to_numpy().ravel())
        print("   pair terms as odds multipliers (posterior median):",
              dict(zip(lam2.term2.values, np.round(np.exp(lam2.median(("chain", "draw")).to_numpy()), 2))))
        print("   three-way terms (posterior median):",
              dict(zip(lam3.term3.values, np.round(lam3.median(("chain", "draw")).to_numpy(), 2))))
    del idata
TIMES["log-linear fits"] = time.time() - t0

M0 independence: 1 s, tune 400, divergences 0, r_hat(N0) 1.001, worst r_hat 1.002 (mu), ESS(N0) 1346
   total N: 7,396  (89%: 7,193-7,604);  cells outside their 89% predictive range: 10 of 15  
   PSIS-LOO: 7 of 15 cells with Pareto k > 0.7

M1 pairs: 3 s, tune 400, divergences 0, r_hat(N0) 1.000, worst r_hat 1.003 (z2), ESS(N0) 2517
   total N: 13,658  (89%: 12,264-15,341);  cells outside their 89% predictive range: 5 of 15  [('0110', 217, 219, 289), ('1001', 106, 66, 105), ('1011', 42, 43, 73), ('1100', 177, 180, 243), ('1110', 181, 118, 167)]
   PSIS-LOO: 11 of 15 cells with Pareto k > 0.7

M1b pairs + triples: 4 s, tune 400, divergences 0, r_hat(N0) 1.004, worst r_hat 1.005 (lk), ESS(N0) 757
   total N: 9,282  (89%: 8,281-10,537);  cells outside their 89% predictive range: 0 of 15  []
   PSIS-LOO: 14 of 15 cells with Pareto k > 0.7
   pair terms as odds multipliers (posterior median): {'EXH-ABA': np.float64(0.97), 'EXH-OSCE': np.float64(1.03), 'EXH-HRW': np.float64(1.48), 'ABA-OSCE': np.float64(1.37), 'ABA-HRW': np.float64(0.81), 'OSCE-HRW': np.float64(2.13)}
   three-way terms (posterior median): {'EXH-ABA-OSCE': np.float64(1.25), 'EXH-ABA-HRW': np.float64(0.01), 'EXH-OSCE-HRW': np.float64(-0.03), 'ABA-OSCE-HRW': np.float64(0.06)}

All three sample cleanly now (no divergences, r_hat at most 1.01). Their answers do not agree. Independence (M0) says about 7,400 people were killed in all and misses 10 of the 15 observed cells: it is the Lincoln-Petersen assumption applied to four lists, and the lists plainly are not independent. The pairwise model (M1) says about 13,600 and misses 5 cells. The model with three-way terms (M1b) reproduces all 15 cells and says about 9,300 (89%: 8,300-10,500); with 15 coefficients for 15 cells it is saturated, so fitting every cell is expected of it and is no evidence in its favour. PSIS leave-one-out is unreliable for about half the cells of M0 and for most of M1 and M1b (the Pareto k counts are printed); part E comes back to that.

Where M1 fails is informative. The figure shows each cell's observed count against the 89% range each model predicts for it.

def ppc_panel(ax, out, nm, keys, title, log=True):
    order = np.argsort(-nm)
    lo, mid, hi = np.quantile(out["yrep"], [0.055, 0.5, 0.945], axis=0)
    x = np.arange(len(nm))
    inside = (lo <= nm) & (nm <= hi)
    ax.vlines(x, lo[order], hi[order], color=LIGHT, lw=6)
    ax.scatter(x, nm[order], s=22, color=np.where(inside[order], INK, LO), zorder=3)
    ax.set_xticks(x, [keys[i] for i in order], rotation=90, fontsize=7.5)
    if log:
        ax.set_yscale("log")
    ax.set_title(f"{title}: {int((~inside).sum())} of {len(nm)} cells outside", loc="left", fontsize=11)


fig, axs = plt.subplots(1, 3, figsize=(12, 3.6), sharey=True)
for ax, (name, out) in zip(axs, KOS.items()):
    ppc_panel(ax, out, n, KEY, name)
axs[0].set_ylabel("people (log scale)")
fig.suptitle("Grey: the range each model predicts for a cell (89%). Dots: what was observed "
             "(orange when outside). Histories are EXH ABA OSCE HRW", x=0.01, ha="left", fontsize=11);

The pairwise model's clearest miss is history 1110, the 181 people found by the exhumations, the ABA interviews and the OSCE but not by Human Rights Watch: it predicts 118 to 167. Those three lists share more people than their pairwise dependencies explain, a three-way dependence, and M1b's posterior puts a large positive term on exactly that triple (its median is printed above, a log odds ratio above 1). Once it is in, every cell fits.

The fitted pair terms are also printed as odds multipliers, each given the other lists and the three-way term. OSCE-HRW overlap about twice as often as chance, EXH-HRW and ABA-OSCE about one and a half times, and two pairs sit at chance. ABA-HRW comes out at about 0.8, slightly less than chance. We will lean on that pair in part E.

#Failure 2: the prior family decides the answer

The saturated maximum-likelihood fit said 17,000; M1b, the same terms with priors, says about 9,300. Both reproduce the 15 cells, so the difference must come from the priors. To find which part of the prior matters, we refit the same ten interaction terms under plain normal priors of three widths, and the horseshoe with three of its other dials turned: the triples' global scale raised from 0.25 to 1, the sd of $\lambda_0$ doubled to 3, and the sd of the main effects $\lambda_k$ raised from 1.5 to 5. The last dial is also tried under the two flattest normal priors. The nearly flat fits mix slowly along the ridge where $\lambda_0$ and the main effects trade off, so every fit here runs 1,500 warm-up steps and 2,000 draws per chain (a few seconds each), and we print the r_hat and ESS of $\lambda_0$.

t0 = time.time()
PRIORS = {"normal, sd 1": dict(ridge=1.0), "normal, sd 3": dict(ridge=3.0), "normal, sd 10": dict(ridge=10.0),
          "normal, sd 3, main sd 5": dict(ridge=3.0, lk_sd=5.0), "normal, sd 10, main sd 5": dict(ridge=10.0, lk_sd=5.0),
          "horseshoe, triples' scale 1": dict(tau_sd=(0.5, 1.0)), "horseshoe, sd of l0 3": dict(l0_sd=3.0),
          "horseshoe, main sd 5": dict(lk_sd=5.0)}
j_ah = PAIRS.index((LISTS.index("ABA"), LISTS.index("HRW")))
b_sat = fit_mle(H, n, PAIRS + TRIPLES)[0]
print("saturated MLE, pair terms as odds multipliers:",
      {f"{LISTS[i]}-{LISTS[j]}": round(float(np.exp(v)), 2) for (i, j), v in zip(PAIRS, b_sat[1 + K:1 + K + len(PAIRS)])})
print(f"saturated MLE, EXH-ABA-OSCE three-way term: {b_sat[1 + K + len(PAIRS)]:.2f}")
for lab, kw in PRIORS.items():
    idata, sec = fit(loglinear(LISTS, H, n, orders=(2, 3), **{"tau_sd": (0.5, 0.25), **kw}), 0.99,
                     tune=1500, draws=2000)
    Nd = NOBS + np.random.default_rng(3).poisson(idata.posterior["N0_mean"].to_numpy().ravel())
    q = np.quantile(Nd, [0.055, 0.5, 0.945])
    ah = idata.posterior["lam2"].to_numpy().reshape(-1, len(PAIRS))[:, j_ah]
    l0d = idata.posterior["l0"].to_numpy().ravel()
    print(f"{lab:27s}: div {int(idata.sample_stats['diverging'].sum())}, r_hat(l0) "
          f"{float(az.rhat(idata, var_names=['l0'])['l0']):.3f}, ESS(l0) {float(az.ess(idata, var_names=['l0'])['l0']):4.0f};"
          f"  N {q[1]:6,.0f} (89%: {q[0]:,.0f}-{q[2]:,.0f});  ABA-HRW odds {np.exp(np.median(ah)):.2f} "
          f"({np.exp(np.quantile(ah, 0.055)):.2f}-{np.exp(np.quantile(ah, 0.945)):.2f}), corr. with l0 "
          f"{np.corrcoef(ah, l0d)[0, 1]:.2f}; main effects {np.round(np.median(idata.posterior['lk'].to_numpy(), (0, 1)), 1)}")
    del idata
TIMES["prior families"] = time.time() - t0

saturated MLE, pair terms as odds multipliers: {'EXH-ABA': 2.32, 'EXH-OSCE': 2.7, 'EXH-HRW': 3.85, 'ABA-OSCE': 3.44, 'ABA-HRW': 1.51, 'OSCE-HRW': 5.39}
saturated MLE, EXH-ABA-OSCE three-way term: 0.39

normal, sd 1               : div 0, r_hat(l0) 1.008, ESS(l0)  795;  N 10,283 (89%: 8,187-13,242);  ABA-HRW odds 0.79 (0.50-1.19), corr. with l0 0.80; main effects [-1.7 -1.9 -1.8 -3. ]

normal, sd 3               : div 0, r_hat(l0) 1.003, ESS(l0)  578;  N 11,230 (89%: 8,050-16,606);  ABA-HRW odds 0.89 (0.49-1.58), corr. with l0 0.88; main effects [-1.8 -2.1 -2.  -3.1]

normal, sd 10              : div 0, r_hat(l0) 1.009, ESS(l0)  574;  N 11,451 (89%: 7,958-17,777);  ABA-HRW odds 0.92 (0.48-1.66), corr. with l0 0.89; main effects [-1.8 -2.1 -2.  -3.1]

normal, sd 3, main sd 5    : div 0, r_hat(l0) 1.004, ESS(l0)  715;  N 13,964 (89%: 9,520-22,095);  ABA-HRW odds 1.18 (0.66-2.15), corr. with l0 0.88; main effects [-2.1 -2.4 -2.3 -3.4]

normal, sd 10, main sd 5   : div 0, r_hat(l0) 1.008, ESS(l0)  523;  N 15,094 (89%: 9,505-26,220);  ABA-HRW odds 1.29 (0.66-2.44), corr. with l0 0.90; main effects [-2.2 -2.5 -2.4 -3.6]

horseshoe, triples' scale 1: div 0, r_hat(l0) 1.002, ESS(l0) 1191;  N  9,215 (89%: 8,152-10,486);  ABA-HRW odds 0.78 (0.54-1.01), corr. with l0 0.61; main effects [-1.5 -1.8 -1.6 -2.8]

horseshoe, sd of l0 3      : div 0, r_hat(l0) 1.002, ESS(l0) 1503;  N  9,260 (89%: 8,229-10,441);  ABA-HRW odds 0.80 (0.57-1.01), corr. with l0 0.56; main effects [-1.5 -1.8 -1.6 -2.8]

horseshoe, main sd 5       : div 0, r_hat(l0) 1.005, ESS(l0) 1281;  N  9,384 (89%: 8,426-10,869);  ABA-HRW odds 0.83 (0.60-1.02), corr. with l0 0.54; main effects [-1.5 -1.8 -1.7 -2.8]

The three horseshoe dials barely matter: raising the triples' scale, doubling the sd of $\lambda_0$ or widening the main effects to sd 5 leaves the answer near 9,200-9,400. With most interactions shrunk, the single-list cells pin the main effects down (their medians run from about -1.5 to -2.8, well inside the prior). The family of the interaction prior is what matters, and once that prior is nearly flat, the main-effect prior matters too. With normal priors of sd 3 or 10 and the main effects at sd 1.5, the median is about 11,200-11,500, with 89% ranges from about 8,000 to 16,600-17,800. Widen the main effects to sd 5 and the median rises to about 14,000-15,000, with ranges reaching 22,000-26,000. The reason is a ridge. In the saturated model each single-list cell fixes $\lambda_0 + \lambda_k$, so a larger hidden count needs smaller main effects, and with the interactions free only the prior on $\lambda_k$ resists. Centred at -1, a catch rate of about 27%, that prior pulls towards fewer hidden people. The saturated MLE of 17,000 lies at or just beyond the top of the sd 1.5 ranges and well inside the sd 5 ones. So "flat" is not one prior, and two reasonable versions of it differ by about a third. A normal prior with sd 1 gives about 10,300. The horseshoe, which expects most dependencies to be small and a few large, gives about 9,300. It gets there by moving the dependence from the pairs into a single triple: the EXH-ABA-OSCE term rises from 0.39 at the MLE to 1.25 (printed above), while every pair falls well below its MLE value (OSCE-HRW from 5.4 to 2.1). With four lists some prior has to decide the barely identified terms, and between sparse, ridge and flat the median runs from about 9,200 to about 15,000. A classical analyst makes the same decision by dropping terms from the model; the Bayesian analyst at least writes it down.

#D · Latent classes: people who differ in visibility

#The model

The log-linear models describe dependence between lists. A different story explains it: people differ. Some victims died in places every organisation reached, others where none did. Manrique-Vallier (2016) models this directly. Each person belongs to one of $C$ hidden classes; within a class the lists are independent, with class-specific capture probabilities $q_{ck}$:

$$P(h) = \sum_{c=1}^{C} \pi_c \prod_{k=1}^{K} q_{ck}^{h_k} (1 - q_{ck})^{1-h_k}.$$

Mixing over classes produces dependence between lists, of any order, without a parameter for each dependence. The weights $\pi$ come from a truncated stick-breaking prior with $C = 10$ (LCMCR's default is 5) and concentration $\alpha \sim \text{Gamma}(0.25, 0.25)$ (LCMCR's default), so the data decide how many classes are used, and $q_{ck} \sim \text{Beta}(1, 1)$. The total enters through a Poisson trick: if $n_h \sim \text{Poisson}(\Lambda P(h))$ for every $h$, the hidden count is $\text{Poisson}(\Lambda P(0))$, and with $p(\Lambda) \propto 1/\Lambda$ this matches the multinomial form with the usual scale-free prior on $N$. The class of each person is summed out in $P(h)$, so NUTS samples only continuous parameters.

def latent_class(lists, Hm, nm, *, C=10, masked=False):
    """Latent-class MSE (Manrique-Vallier 2016) with a truncated stick-breaking prior, in Poisson-total form."""
    H_all = np.vstack([np.zeros(len(lists), int), Hm]).astype(float)
    with pm.Model(coords={"cell": ["".join(map(str, h)) for h in Hm], "list": lists, "cls": np.arange(C)}) as m:
        alpha = pm.Gamma("alpha", 0.25, 0.25)
        w = pm.StickBreakingWeights("w", alpha, C - 1, dims="cls")
        q = pm.Beta("q", 1.0, 1.0, dims=("cls", "list"))
        lp_hc = pt.dot(H_all, pt.log(q).T) + pt.dot(1 - H_all, pt.log1p(-q).T)  # (cells, classes)
        logP = pt.logsumexp(lp_hc + pt.log(w)[None, :], axis=1)
        u = pm.Flat("u")  # log Lambda - log n_obs; flat in log Lambda is p(Lambda) ~ 1/Lambda
        logLam = np.log(nm.sum()) + u
        pm.Deterministic("N0_mean", pt.exp(logLam + logP[0]))
        likelihood(pm.Deterministic("mu", pt.exp(logLam + logP[1:]), dims="cell"), nm, masked)
    return m


SPEC["M2 latent class"] = (lambda **kw: latent_class(LISTS, H, n, C=10, **kw), 0.99)
t0 = time.time()
idata, sec = fit(SPEC["M2 latent class"][0](), 0.99)
KOS["M2 latent class"] = summarise("M2 latent class", idata, n, sec, skip=("q", "w", "mu"))
post = idata.posterior
print("   worst r_hat of q (capture probabilities):", round(float(az.rhat(idata, var_names=["q"])["q"].max()), 3))
print("   per-chain median hidden count:", np.round(np.median(post["N0_mean"].to_numpy(), axis=1), -1))
wdraw = post["w"].to_numpy().reshape(-1, 10)
print("   classes holding more than 5% of people, per draw (mean):", round(float((wdraw > 0.05).sum(1).mean()), 1))
qdraw = post["q"].to_numpy().reshape(-1, 10, K)
top = np.argsort(-wdraw, axis=1)[:, :2]  # the two largest classes in each draw, whatever their labels
r_ = np.arange(len(wdraw))[:, None]
for j in range(2):
    qj = qdraw[r_[:, 0], top[:, j]]
    print(f"   class ranked {j + 1} by size: share of all people {np.median(wdraw[r_[:, 0], top[:, j]]):.2f}, "
          f"capture probabilities {dict(zip(LISTS, np.round(np.median(qj, 0), 2)))}, "
          f"P(on no list) {np.median(np.prod(1 - qj, 1)):.2f}")
del idata, post
TIMES["latent class fit"] = time.time() - t0

M2 latent class: 9 s, tune 400, divergences 1, r_hat(N0) 1.003, worst r_hat 1.005 (alpha), ESS(N0) 873
   total N: 10,420  (89%: 9,282-12,795);  cells outside their 89% predictive range: 0 of 15  []
   PSIS-LOO: 15 of 15 cells with Pareto k > 0.7

   worst r_hat of q (capture probabilities): 1.272
   per-chain median hidden count: [6010. 5950. 6000. 6090.]
   classes holding more than 5% of people, per draw (mean): 1.9
   class ranked 1 by size: share of all people 0.86, capture probabilities {'EXH': np.float64(0.15), 'ABA': np.float64(0.12), 'OSCE': np.float64(0.12), 'HRW': np.float64(0.02)}, P(on no list) 0.65
   class ranked 2 by size: share of all people 0.07, capture probabilities {'EXH': np.float64(0.28), 'ABA': np.float64(0.08), 'OSCE': np.float64(0.33), 'HRW': np.float64(0.27)}, P(on no list) 0.23

The latent-class model (M2) says about 10,400 (89%: 9,300-12,800), almost exactly Ball et al.'s 10,356 (9,002-12,122), which came from log-linear models fitted separately by time and region - a different route to the same number. Silverman (2020) got 10,400 from LCMCR itself on this table. M2 reproduces all 15 cells, which with ten classes of four capture probabilities is again easy. On average about two classes hold more than 5% of the people each. Ranking the classes by size in each draw, so that labels do not matter, the largest holds over four fifths of all the people, each list finds at most about one in seven of them (Human Rights Watch one in fifty), and about two thirds are on no list. The next, under a tenth of the people, is found about twice as often by the exhumations, nearly three times as often by the OSCE and more than ten times as often by Human Rights Watch (about one in four rather than one in fifty), a little less often by the ABA, and only about a quarter of them are on no list.

Two diagnostics need reading together. The capture probabilities $q$ have r_hat up to about 1.3: chains disagree about which class is which, because swapping two classes' labels leaves the likelihood unchanged (label switching, as in E69), and the near-empty classes are pure prior. The hidden count does not care about labels: its r_hat is 1.00, and the four chains' medians agree within about 2%. We judge the model by the quantity we want, not by its unidentified parameters. One divergence remains at target_accept=0.99; one draw in 4,000 does not change the hidden count's summary. The limit of the model is its flexibility: with ten classes of four capture probabilities it can fit almost any table of 15 cells, so fitting this one is no test of it.

We now have four models and four answers: 7,400 from a model that fails, 13,600 from one that misses a dependence, and about 9,300 and 10,400 from the two that fit every cell.

#E · Averaging over assumptions

A Bayesian does not have to choose one model. There are two common ways to combine them, and they answer different questions.

#Bayesian model averaging over 113 hierarchical log-linear models

Madigan & York (1997) put a prior over the models themselves and average with posterior model probabilities, $p(M \mid n) \propto p(n \mid M)\, p(M)$. They used decomposable graphical models with conjugate priors; we work in their spirit with a wider set. For four lists we take every hierarchical log-linear model - any set of pair terms, plus any three-way term whose three pairs are all in - which gives 113 models, each with equal prior weight. The marginal likelihood $p(n \mid M)$ of each is computed by a Laplace approximation with $N(0, 1)$ priors on the interactions (the same priors as above for the other terms). For a Poisson log-linear model the posterior is very close to Gaussian, so this takes a fraction of a second for all 113.

def laplace_models(Hm, nm, sd_int=1.0, triples=True, lk_sd=1.5):
    """Laplace evidence and the hidden-count posterior (log-normal) of every hierarchical log-linear model."""
    Kl, nobs = Hm.shape[1], nm.sum()
    pairs = list(it.combinations(range(Kl), 2))
    trips = list(it.combinations(range(Kl), 3)) if triples else []
    res = []
    for mask in range(2 ** len(pairs)):
        P = [p for j, p in enumerate(pairs) if mask >> j & 1]
        okT = [t for t in trips if all(pp in P for pp in it.combinations(t, 2))]
        for r in range(len(okT) + 1):
            for T in it.combinations(okT, r):
                X = design(Hm, P + list(T))
                d = X.shape[1]
                m0 = np.r_[np.log(nobs), np.full(Kl, -1.0), np.zeros(d - 1 - Kl)]
                s0 = np.r_[1.5, np.full(Kl, lk_sd), np.full(d - 1 - Kl, sd_int)]
                b = m0.copy()
                b[0] = np.log(nm.mean() + 1)
                for _ in range(100):  # Newton steps to the posterior mode
                    mu = np.exp(X @ b)
                    Hs = X.T @ (X * mu[:, None]) + np.diag(1 / s0**2)
                    step = np.linalg.solve(Hs, X.T @ (nm - mu) - (b - m0) / s0**2)
                    b = b + np.clip(step, -2, 2)
                    if np.abs(step).max() < 1e-9:
                        break
                mu = np.exp(X @ b)
                Hs = X.T @ (X * mu[:, None]) + np.diag(1 / s0**2)
                ll = np.sum(nm * np.log(mu) - mu - gammaln(nm + 1))
                lp = np.sum(-0.5 * ((b - m0) / s0) ** 2 - np.log(s0) - 0.5 * np.log(2 * np.pi))
                logev = ll + lp + 0.5 * d * np.log(2 * np.pi) - 0.5 * np.linalg.slogdet(Hs)[1]
                res.append((P, list(T), logev, b[0], np.sqrt(np.linalg.inv(Hs)[0, 0])))
    le = np.array([r[2] for r in res])
    return res, np.exp(le - logsumexp(le))


def bma_draws(res, w, nobs, size=20_000, seed=1):
    g = np.random.default_rng(seed)
    j = g.choice(len(res), size, p=w)
    l0 = np.array([res[i][3] for i in j]) + g.standard_normal(size) * np.array([res[i][4] for i in j])
    return nobs + g.poisson(np.exp(l0))


short = lambda P, T: " ".join(["".join(LISTS[i][:2] for i in p) for p in P] +
                              ["[" + "".join(LISTS[i][:2] for i in t) + "]" for t in T])
t0 = time.time()
BMA = {}
for sd in (0.5, 1.0, 2.0):
    res, w = laplace_models(H, n, sd)
    BMA[sd] = bma_draws(res, w, NOBS)
    q = np.quantile(BMA[sd], [0.055, 0.5, 0.945])
    print(f"interaction prior sd {sd}: {len(res)} models, effective number {1 / np.sum(w**2):.1f}, "
          f"N {q[1]:,.0f} (89%: {q[0]:,.0f}-{q[2]:,.0f})")
    if sd == 1.0:
        RES1, W1 = res, w
        for j in np.argsort(-w)[:4]:
            print(f"     weight {w[j]:.2f}  hidden ~{np.exp(res[j][3]):6,.0f}  terms: {short(res[j][0], res[j][1])}")
        inc = lambda test: sum(w[i] for i, r in enumerate(res) if test(r))
        print(f"     P(ABA-HRW pair in) = {inc(lambda r: (1, 3) in r[0]):.2f}, "
              f"P(EXH-ABA-OSCE triple in) = {inc(lambda r: (0, 1, 2) in r[1]):.2f}")
res, w = laplace_models(H, n, 1.0, lk_sd=5.0)
q = np.quantile(bma_draws(res, w, NOBS), [0.055, 0.5, 0.945])
print(f"interaction prior sd 1.0, main effects sd 5: N {q[1]:,.0f} (89%: {q[0]:,.0f}-{q[2]:,.0f}), "
      f"top model's weight {w.max():.2f}")
TIMES["BMA"] = time.time() - t0

interaction prior sd 0.5: 113 models, effective number 6.7, N 9,978 (89%: 8,675-11,542)
interaction prior sd 1.0: 113 models, effective number 2.6, N 10,868 (89%: 9,057-13,173)
     weight 0.58  hidden ~ 7,022  terms: EXAB EXOS EXHR ABOS OSHR [EXABOS] [EXOSHR]
     weight 0.17  hidden ~ 5,728  terms: EXAB EXOS EXHR ABOS OSHR [EXABOS]
     weight 0.09  hidden ~ 6,587  terms: EXAB EXOS EXHR ABOS ABHR OSHR [EXABOS] [EXOSHR]
     weight 0.08  hidden ~ 5,230  terms: EXAB EXOS EXHR ABOS ABHR OSHR [EXABOS]
     P(ABA-HRW pair in) = 0.25, P(EXH-ABA-OSCE triple in) = 1.00

interaction prior sd 2.0: 113 models, effective number 2.1, N 11,232 (89%: 9,178-13,988)
interaction prior sd 1.0, main effects sd 5: N 11,400 (89%: 9,444-13,903), top model's weight 0.73

With $N(0, 1)$ interaction priors, the average over 113 models puts the total at about 10,900 (89%: about 9,100-13,200). Most of the weight sits on a few models, all of which contain the EXH-ABA-OSCE triple (its inclusion probability is essentially 1), and the top one leaves out the ABA-HRW pair - the same two facts the horseshoe fits found. In the codes, EXAB is the EXH-ABA pair and [EXABOS] the EXH-ABA-OSCE triple.

The catch is visible in the three rows. Model averaging is sensitive to the prior on the interactions: a narrower prior (sd 0.5) penalises big models less and spreads the weight over more of them, and the median moves from about 10,000 to 11,200 as the sd goes from 0.5 to 2. The marginal likelihood contains an Occam factor that depends on how wide the prior is, so the prior here does not just regularise, it votes. That is well known (Madigan & York discuss it), and it is one reason to also look at an approach that does not use marginal likelihoods. The main-effect prior votes too: widened to sd 5 it moves the average to about 11,400 (last row).

Silverman (2020) ran Madigan & York's own method (the dga package) on this table and got a median of 12,600 (95%: 9,800-16,500), with one model holding almost all the weight. Their version differs from ours in three ways: decomposable graphs only, conjugate priors on cell probabilities, and no informative prior on the main effects, which the last row suggests is worth about 500 people.

#Stacking needs exact leave-one-cell-out

Stacking (Yao, Vehtari, Simpson & Gelman 2018) chooses weights that maximise the leave-one-out predictive density of the combined model. The unit of data here is a cell: leaving out one history and predicting its count from the other 14. PSIS-LOO approximates that from one fit, and the fits above printed how often it is reliable.

for name, out in KOS.items():
    print(f"{name:20s}: PSIS-LOO Pareto k > 0.7 in {out['k_bad']} of {len(n)} cells")

M0 independence     : PSIS-LOO Pareto k > 0.7 in 7 of 15 cells
M1 pairs            : PSIS-LOO Pareto k > 0.7 in 11 of 15 cells
M1b pairs + triples : PSIS-LOO Pareto k > 0.7 in 14 of 15 cells
M2 latent class     : PSIS-LOO Pareto k > 0.7 in 15 of 15 cells

Unreliable almost everywhere. Leaving out a cell of several hundred people changes the posterior far too much for importance sampling to bridge. So we refit each model 15 times, once without each cell, and score the left-out count exactly. Each model is compiled once by nutpie with the likelihood in a pm.Potential multiplied by a 0/1 mask held in pm.Data; with_data swaps the mask between refits without recompiling (E64 and E103 use the same trick). Sixty refits at 300 draws per chain take a little over a minute.

def exact_loo(spec, nm, draws=300, seed0=100):
    """Leave-one-cell-out log predictive density of every cell, for each model in `spec`."""
    L = np.zeros((len(nm), len(spec)))
    for j, (name, (mk, ta)) in enumerate(spec.items()):
        t0 = time.time()
        comp = nutpie.compile_pymc_model(mk(masked=True))
        divs = 0
        for i in range(len(nm)):
            mask = np.ones(len(nm))
            mask[i] = 0
            tr = nutpie.sample(comp.with_data(mask=mask), chains=4, draws=draws, tune=400, seed=seed0 + i,
                               progress_bar=False, target_accept=ta)
            mu = tr.posterior["mu"].to_numpy()[..., i].ravel()
            L[i, j] = logsumexp(nm[i] * np.log(mu) - mu - gammaln(nm[i] + 1)) - np.log(mu.size)
            divs += int(tr.sample_stats["diverging"].sum())
        print(f"{name:20s}: {len(nm)} refits in {time.time() - t0:4.0f} s, divergences {divs:3d}, "
              f"elpd {L[:, j].sum():7.1f}")
    return L


def stack_weights(L):
    def neg(z):
        return -np.sum(logsumexp(L + np.log(np.exp(z - logsumexp(z)) + 1e-300), axis=1))
    r = minimize(neg, np.zeros(L.shape[1]), method="Nelder-Mead",
                 options={"maxiter": 20000, "xatol": 1e-8, "fatol": 1e-10})
    return np.exp(r.x - logsumexp(r.x))


t0 = time.time()
LOO = exact_loo(SPEC, n)
TIMES["exact LOO"] = time.time() - t0
names = list(SPEC)
w_st = stack_weights(LOO)
elpd = LOO.sum(0)
w_pb = np.exp(elpd - logsumexp(elpd))
print("\n" + pd.DataFrame({"elpd (exact)": elpd.round(1), "stacking": w_st.round(2), "pseudo-BMA": w_pb.round(2)},
                          index=names).to_string())
worst = LOO.argmin(axis=0)
print("each model's worst-predicted cell:", {nm: f"{KEY[i]} (n={n[i]}, lpd {LOO[i, j]:.0f})"
                                             for j, (nm, i) in enumerate(zip(names, worst))})

M0 independence     : 15 refits in    2 s, divergences   0, elpd  -250.9

M1 pairs            : 15 refits in   11 s, divergences   0, elpd  -379.5

M1b pairs + triples : 15 refits in   14 s, divergences   0, elpd   -88.7

M2 latent class     : 15 refits in   48 s, divergences   1, elpd   -84.2

                     elpd (exact)  stacking  pseudo-BMA
M0 independence            -250.9      0.05        0.00
M1 pairs                   -379.5      0.01        0.00
M1b pairs + triples         -88.7      0.22        0.01
M2 latent class             -84.2      0.73        0.99
each model's worst-predicted cell: {'M0 independence': '1110 (n=181, lpd -63)', 'M1 pairs': '0100 (n=845, lpd -161)', 'M1b pairs + triples': '0100 (n=845, lpd -11)', 'M2 latent class': '0100 (n=845, lpd -9)'}

The exact scores are revealing. The pairwise model M1 has by far the worst score, and over two fifths of it comes from one cell: left without the 845 people on the ABA list alone, its pair terms run away and it predicts that cell absurdly. The two models that fit every cell, M1b and M2, score about 4.5 nats apart and far ahead of the other two, by more than 160 nats.

Exact leave-one-out at 300 draws per chain has a Monte Carlo error of a few nats, so the weights built on these scores are themselves uncertain, and the two weighting schemes react to that differently. Pseudo-BMA weights each model by $e^{\text{elpd}}$, so a difference of about 4.5 nats becomes a factor of about 90 and M2 takes almost all the weight; a few nats of noise could reverse it. Stacking asks which mixture predicts best and splits the weight between M1b and M2, since each predicts some cells better than the other: about a fifth on M1b and the rest on M2 here, a split that a few nats of noise could also shift noticeably. Stacking can also give a small weight to a model that is plainly wrong about the quantity we care about - independence gets a few percent - because it judges only the observed cells, and a model can predict a cell of 300 people tolerably while being thousands of people wrong about the cell no one sees.

That last point decides how we build the headline. We stack only M1b and M2: they are more than 160 nats ahead of M0 and M1 on exact leave-one-cell-out, and M0 and M1 also fail the cell-by-cell check of part C. (That M1b and M2 pass it says little, since both have about as many parameters as cells.) The restriction was chosen after seeing those checks; it removes models whose small stacking weight would come from fitting the observed cells while being wrong about the hidden one.

j2 = [names.index("M1b pairs + triples"), names.index("M2 latent class")]
w2 = stack_weights(LOO[:, j2])
print(f"two-model stacking weights: M1b {w2[0]:.2f}, M2 {w2[1]:.2f}")
g = np.random.default_rng(7)
pick = g.random(KOS["M2 latent class"]["N"].size) < w2[0]
N_MIX = np.where(pick, KOS["M1b pairs + triples"]["N"], KOS["M2 latent class"]["N"])
HID = N_MIX - NOBS
qN = np.quantile(N_MIX, [0.025, 0.055, 0.5, 0.945, 0.975])
qH = np.quantile(HID, [0.055, 0.5, 0.945, 0.975])
print(f"total killed: median {qN[2]:,.0f}, 89% {qN[1]:,.0f}-{qN[3]:,.0f}, 95% {qN[0]:,.0f}-{qN[4]:,.0f}")
print(f"on no list:   median {qH[1]:,.0f}, 89% {qH[0]:,.0f}-{qH[2]:,.0f}, 97.5% point {qH[3]:,.0f}")
per10 = 10 * qH[:3] / NOBS
print(f"hidden per 10 listed: about {per10[1]:.1f} (89%: {per10[0]:.1f}-{per10[2]:.1f})")
FLOOR = int(round(qH[0], -2))  # the bottom of the solid block in the iceberg
for lab, p in [("hidden > listed (more on no list than on any)", (HID > NOBS).mean()),
               (f"at least {FLOOR:,} on no list", (HID >= FLOOR).mean()),
               ("total above Ball et al.'s 10,356", (N_MIX > 10356).mean())]:
    print(f"P({lab}) = {p:.3f}  -> {pub.calibrated_word(p)}")
for k in ("M1b pairs + triples", "M2 latent class"):
    p = (KOS[k]["N"] - NOBS >= FLOOR).mean()
    print(f"   {k} alone: P(at least {FLOOR:,} on no list) = {p:.3f}  -> {pub.calibrated_word(p)}")

two-model stacking weights: M1b 0.26, M2 0.74
total killed: median 10,128, 89% 8,720-12,399, 95% 8,431-13,341
on no list:   median 5,728, 89% 4,320-7,999, 97.5% point 8,941
hidden per 10 listed: about 13.0 (89%: 9.8-18.2)
P(hidden > listed (more on no list than on any)) = 0.932  -> very likely
P(at least 4,300 on no list) = 0.948  -> very likely
P(total above Ball et al.'s 10,356) = 0.413  -> about as likely as not
   M1b pairs + triples alone: P(at least 4,300 on no list) = 0.832  -> likely
   M2 latent class alone: P(at least 4,300 on no list) = 0.991  -> virtually certain

With the printed weights, the stacked total is about 10,100 (89%: about 8,700-12,400) and the number of people on no list about 5,700 (89%: about 4,300-8,000). Per 10 people named on a list, that is about 13 more on none, between about 10 and 18 - the sentence in the plain-words summary. It is very likely (95%) that at least 4,300 of the dead are on no list, the bottom of the solid block in the headline figure. That number leans on the mixture: M2 alone makes it virtually certain (99%), M1b alone only likely (83%). That more people were missed than were listed is also very likely (93%), though close to the 90% line between "likely" and "very likely". About four in ten draws have more deaths than the published 10,356.

The bars below show how the weights split the models, all four under each scheme and the two that make the headline.

fig, ax = plt.subplots(figsize=(10, 4.0))
labels = ["M0: lists are independent", "M1: lists depend in pairs", "M1b: pairs and threes",
          "M2: people differ in visibility"]
y = np.arange(4)
w_head = np.zeros(4)
w_head[j2] = w2
for off, w_, col, tcol, lab in [(-0.27, w_head, INK, INK, "stacking, M1b + M2 (the headline)"),
                                (0.0, w_st, HI, HI, "stacking, all four"),
                                (0.27, w_pb, LIGHT, MUTED, "pseudo-BMA, all four")]:
    ax.barh(y + off, w_, height=0.25, color=col, label=lab)
    for yy, v in zip(y, w_):
        ax.text(v + 0.01, yy + off, f"{v:.0%}", va="center", fontsize=9.5, color=tcol)
ax.set_yticks(y, labels)
ax.invert_yaxis()
ax.set_xlim(0, 1.12)
ax.set_xlabel("share of the weight")
ax.legend(loc="center right", frameon=False)
ax.set_title("Three ways to weight the models: both stackings split the answer, pseudo-BMA picks one",
             loc="left");

#What you would have to believe

Averaging handles the dependencies the data can see. A sceptic can still ask: what if two lists depend on each other more than the model suggests? The question has a precise answer. Take the ABA-HRW pair, which M1b puts at an odds multiplier of about 0.8 (printed below), fix it at a chosen value - 1, 1.5, 2, 3 or 4 times chance - and refit M1b with everything else free. One compiled model with the fixed value in pm.Data serves all five fits.

print("ABA-HRW odds multiplier when free (M1b): median {:.2f}, 89% {:.2f}-{:.2f}".format(
    *np.quantile(ABA_HRW_FREE, [0.5, 0.055, 0.945])))
t0 = time.time()
comp = nutpie.compile_pymc_model(loglinear(LISTS, H, n, orders=(2, 3), tau_sd=(0.5, 0.25),
                                           fix_pair=(LISTS.index("ABA"), LISTS.index("HRW"))))
SENS = {}
for mult in (1, 1.5, 2, 3, 4):
    tr = nutpie.sample(comp.with_data(fixed=np.log(mult)), chains=4, draws=500, tune=400, seed=RANDOM_SEED,
                       progress_bar=False, target_accept=0.99)
    SENS[mult] = NOBS + rng.poisson(tr.posterior["N0_mean"].to_numpy().ravel())
    q = np.quantile(SENS[mult], [0.055, 0.5, 0.945])
    print(f"ABA-HRW fixed at x{mult}: divergences {int(tr.sample_stats['diverging'].sum())}, "
          f"N {q[1]:,.0f} (89%: {q[0]:,.0f}-{q[2]:,.0f})")
TIMES["sensitivity"] = time.time() - t0


def ridgeline(ax, rows, ys, xgrid, color=HI, height=0.8):
    """Smoothed densities of each row's draws, drawn as ridges at heights ys (largest peak = height)."""
    for v, y0 in zip(rows, ys):
        h = np.histogram(v, bins=xgrid)[0].astype(float)
        k = np.exp(-0.5 * (np.arange(-12, 13) / 4) ** 2)
        h = np.convolve(h, k / k.sum(), mode="same")
        h = height * h / h.max()
        xc = 0.5 * (xgrid[1:] + xgrid[:-1])
        ax.fill_between(xc, y0, y0 + h, color=color, alpha=0.55, lw=0, zorder=3)
        ax.plot(xc, y0 + h, color=color, lw=1, zorder=4)


fig, ax = plt.subplots(figsize=(10, 4.8))
mults = list(SENS)
ys = np.log2(mults)
lo_f, hi_f = np.quantile(ABA_HRW_FREE, [0.055, 0.945])
ax.axhspan(np.log2(lo_f), np.log2(hi_f), color=LIGHT, alpha=0.7, lw=0, zorder=0)
ax.text(40_000, np.log2(lo_f) + 0.05, f"what M1b, with its horseshoe prior, says about this pair\n"
        f"(89% range: x{lo_f:.1f} to x{hi_f:.1f})", ha="right", va="bottom", fontsize=10, color=MUTED)
ridgeline(ax, [SENS[mm] for mm in mults], ys, np.arange(6000, 42000, 100), height=0.36)
for mm, y0 in zip(mults, ys):
    ax.text(np.quantile(SENS[mm], 0.945) + 400, y0 + 0.04, f"middle: {round(np.median(SENS[mm]), -2):,.0f}",
            ha="left", va="bottom", fontsize=9.5, color=INK)
ax.axvline(np.median(N_MIX), color=INK, lw=1, ls=(0, (3, 2)))
ax.text(np.median(N_MIX), ys[-1] + 0.5, "  headline (stacked) median", ha="left", fontsize=10)
ax.set_yticks(ys, [f"x{mm}" for mm in mults])
ax.set_ylabel("assumed ABA-HRW dependence\n(odds multiplier)")
ax.set_xlabel("total killed")
ax.set_xlim(6000, 40000)
ax.set_ylim(np.log2(lo_f) - 0.25, ys[-1] + 0.7)
ax.set_title(f"To reach about {round(np.median(SENS[2]), -2):,.0f} deaths, two lists would have to overlap "
             "twice as often as chance", loc="left");