hondius_metaculus_43461.py

"""
Metaculus #43461 -- https://www.metaculus.com/questions/43461/5-non-ship-cases-linked-to-the-hondius-outbreak-before-aug-2026/
Will at least 5 non-passengers be linked to the MV Hondius hantavirus outbreak before August 2026?

Resolution: YES if >=5 hantavirus cases (lab-confirmed OR formally designated
probable/suspected) in people NOT aboard during the March-April 2026 voyage,
linked to the outbreak by a public health authority, before 2026-08-01.

Today: 2026-05-10. Community: 23%.

Model structure
---------------
1. Decompose secondary transmission into independent exposure pools generated
   by the symptomatic widow's post-disembark chain (Airlink StH-JNB flight +
   St. Helena island + JNB hospital + OR Tambo + brief KL592 boarding event),
   plus a household pool for other infectious passengers.
2. Mixture in a low-probability Andes super-spreader event (Epuyen analogue).
3. Apply a public-linkage filter (true infection -> formally designated case).
4. Bayesian update: condition on the observation that 0 non-passenger cases
   have been publicly linked as of today (day ~15 since dominant exposure).
   
NOTE: Model built in squigglepy: https://github.com/peterhurford/squigglepy
"""

import numpy as np
import squigglepy as sq
from squigglepy import bayes


# =============================================================================
# Exposure pools (widow's symptomatic chain dominates)
# =============================================================================

# Airlink St. Helena -> Johannesburg flight, April 25. ~82 passengers + 6
# crew, multi-hour flight with widow symptomatic and deteriorating in-flight.
# This is the dominant long-cabin exposure event.
# Close = rows adjacent within ~1m for hours.
# SAR anchored to Epuyen 2018-19 Patient 1 ("fever and malaise" at the party
# infected 5 of ~100 close attendees over 90 min => ~5%). BUT the widow on
# Airlink was earlier in her viral curve than Epuyen Patient 1 -- GI/systemic
# prodrome, not yet severe respiratory; HPS onset was ~24h later at KL592.
# The Alonso et al. viral-shedding curve has infectiousness ramping in the
# 1-2 days BEFORE HPS onset, so transmissibility on Airlink was probably
# 30-70% of peak. Modest downweight (~2x) from raw Epuyen SAR.
# Five intensive close contacts already identified by GGD; KLM flight
# attendant (one of them) tested NEGATIVE (NL Times, May 8). Calendar-quiet
# observation absorbs this.
n_airlink_close   = sq.lognorm(3, 12)
sar_airlink_close = sq.lognorm(0.003, 0.035)

# Rest of Airlink cabin + crew. No documented Andes transmission outside
# close range, so SAR floor is essentially 0 but non-zero given duration.
n_airlink_casual  = sq.lognorm(55, 90)
sar_airlink_casual = sq.lognorm(0.0001, 0.005)

# KL592 ground boarding event (Johannesburg, April 26). Widow boarded
# briefly, was identified as ill by crew, and was removed BEFORE departure
# (per KLM Newsroom statement). The flight then departed without her.
# So the real risk is to the few ground staff / crew who physically assisted
# her removal -- NOT the AMS-bound passengers.
n_kl592_assist   = sq.lognorm(2, 6)
sar_kl592_assist = sq.lognorm(0.005, 0.04)

# St. Helena community. Widow walked the island symptomatic for ~24h before
# being flown out. Population ~4000 but actual encounters likely 15-80.
# Mostly brief contacts; SAR very low for non-prolonged interaction.
n_sth_local   = sq.lognorm(15, 80)
sar_sth_local = sq.lognorm(0.001, 0.02)

# St. Helena body handlers. Husband died on the ship April 11 and his body
# remained aboard for 13 days before St. Helena offload (April 24). By that
# point the body had been prepared/sealed on board; the high-risk handling
# (peri-mortem care, body bag prep) happened with ship CREW (= passengers,
# excluded from this question) not St. Helena personnel. St. Helena handlers
# faced post-preparation transit handling -- much lower transmission risk
# than acute peri-mortem exposure. SAR halved from raw bodily-fluid prior.
n_sth_body   = sq.lognorm(2, 6)
sar_sth_body = sq.lognorm(0.01, 0.10)

# Johannesburg hospital pool, split protected vs unprotected.
# Per Daily Maverick reporting: widow collapsed at Kempton Park hospital
# (brief ER first-response, pre-diagnosis -> minimal PPE). Separate Andes-
# positive patient at Sandton ICU (already isolated with PPE).
# Post-COVID HCW PPE compliance is generally high, so most hospital exposure
# is well-controlled. Main risk: 2-5 unprotected first responders.
n_jnb_unprotected   = sq.lognorm(2, 5)
sar_jnb_unprotected = sq.lognorm(0.03, 0.12)
n_jnb_protected     = sq.lognorm(5, 15)
sar_jnb_protected   = sq.lognorm(0.003, 0.02)

# OR Tambo airport contacts (brief transit encounters).
n_or_tambo   = sq.lognorm(20, 100)
sar_or_tambo = sq.lognorm(0.0001, 0.003)

# Other infectious-but-undetected disembarked passengers -> household contacts.
# 6 confirmed cases are now hospitalized/isolated. Risk comes from undetected
# subclinical/mild infections among the other ~140 ship occupants going home
# with normal household contact. Many countries imposed quarantine (NL: 6
# weeks), which suppresses this pool meaningfully.
# SAR ~Epuyen secondary household rate (~5-15%).
n_other_infectious_pax = sq.lognorm(0.5, 4)
n_household_per_pax    = sq.lognorm(2, 5)
sar_household          = sq.lognorm(0.02, 0.15)

# =============================================================================
# Super-spreader event (Epuyen analogue)
# =============================================================================
# Empirical base rate from ~30 years of ANDV surveillance:
#   - 1996 El Bolson: 16 epi-linked cases (first H2H documented).
#   - 1993-2005 southern Argentina: 9 small clusters, 20 total cases
#     (Padula et al, PubMed 17370522) -- mean ~2.2/cluster, mostly 1
#     secondary, NOT super-spreaders.
#   - 2018-19 Epuyen: 34 cases driven by 3 super-spreaders attending
#     crowded events (NEJM 2020). Characterized as the FIRST documented
#     true super-spreading in Andes history.
# => ~2 clear >=5-secondary events in ~30 years against an unknown but
# substantial number of "candidate" index cases with prodromal social
# exposure. Even accounting for detection bias FAVORING super-spreader
# observation (single H2H goes unnoticed), empirical rate looks 3-8%.
# The widow's exposure profile is more diverse than typical Patagonian
# index cases (multiple venues, dense settings), nudging up within range.
P_SUPERSPREAD = 0.05
# Burst size given super-spreader event. Epuyen 2018-19 had ~5-10 immediate
# secondaries per super-spreader event (Patient 1's birthday party generated
# 5 close-contact secondaries; the two other ss events similar scale).
# Hondius has aggressive contact tracing from ~day 14 vs Epuyen's months-long
# recognition lag, so the runaway-chain regime is capped earlier. Tightened
# upper tail accordingly.
ss_burst = sq.lognorm(3, 10)

# =============================================================================
# Public-linkage filter
# =============================================================================
# A true infection must (a) produce symptoms severe enough to seek care,
# (b) get hantavirus testing ordered, (c) test positive OR get a formal
# probable/suspected designation linked to Hondius.
# Pipeline decomposition under intense active surveillance:
#   P(symptomatic | infected)              ~0.85-0.95 (Andes is high-CFR,
#       asymptomatic seroprevalence 5-15% mostly from mouse exposure; H2H
#       secondaries in Epuyen 2018-19 were all symptomatic).
#   P(seeks care | symptomatic)            ~0.95 (severe respiratory disease).
#   P(hantavirus test ordered | care)      ~0.95 in traced pools, ~0.3-0.6
#       in untraced tail. Vast majority of expected infections are in traced
#       pools (KLM flight, 62 SA contacts, St Helena, JNB HCWs).
#   P(positive + publicly linked | tested) ~0.90+ (any Andes positive right
#       now will be attributed to Hondius by default; sequencing is fast).
# Weighted: ~0.80-0.90. Biological ceiling ~0.90 from asymptomatic floor.
# Beta(10, 2): mean 0.83, P5~0.61, P95~0.97.
# Limitation: applies uniformly across pools. In reality the household pool
# (disembarked pax -> home contacts in random countries) has weaker linkage
# than actively-traced pools (KLM/JNB/St Helena). Effect estimated ~1pp;
# left as single-rate filter for simplicity.
p_link = sq.beta(a=10, b=2)

# =============================================================================
# Pipeline timing: exposure -> incubation -> diagnosis -> public report
# =============================================================================
# Empirical ANDV incubation distribution, H2H-SPECIFIC.
# H2H secondary cases tend toward the upper end of the general range:
#   - CDC general ANDV: 4-42 days from exposure to symptoms.
#   - WHO: 1-6 weeks typical, up to 8 weeks in rare cases.
#   - CDC EID person-to-person Andes paper (Padula et al, 2005): incubation
#     periods of 15-24 days observed in analyzed H2H chains.
#   - Epuyen 2018-19 secondaries (NEJM Martinez 2020): onset 17-40d post-
#     exposure event.
# lognorm(14, 42) gives 90% CI [14, 42], median ~24d, P95~42d. Matches
# H2H-specific literature center; only ~5-10% mass before day 15.
# Note: the audit's pushback that earlier-mass distributions over-update
# the "0 cases at day 15" observation is correct -- this is the fix.
incubation_dist = sq.lognorm(14, 42)

# Onset -> tested -> confirmed -> publicly reported as Hondius-linked.
# Subcomponents: seek care (mild prodrome can delay 1-2d) + clinician orders
# hantavirus PCR (suspicion requires epi link known) + lab turnaround (4-24h) +
# report to public health authority + news pickup (~12h for a high-attention
# outbreak with daily WHO updates).
# lognorm(0.5, 3) gives 90% CI [12h, 3d], median ~1.2d.
report_lag_dist = sq.lognorm(0.5, 3)

# Widow's dominant exposure events: April 23-27.
# Today: May 10 -> 13-17 days since (varies by pool).
# Aug 1 cutoff: May 10 + 83 days.
#
# TODO / NOTE: THIS IS HARDCODED! IT - AND OTHER PARTS OF THE CODE - MUST BE UPDATED TO RE-RUN ON A LATER DATE!
#
DAYS_SINCE_EXPOSURE = sq.uniform(13, 17)
DAYS_TO_AUG1 = 83


# =============================================================================
# Joint model
# =============================================================================

def true_infections():
    """Sum of secondary infections across all pools + super-spreader mix.
    Pre-linkage, pre-observation filter."""
    pool_rates = [
        ~n_airlink_close      * ~sar_airlink_close,
        ~n_airlink_casual     * ~sar_airlink_casual,
        ~n_kl592_assist       * ~sar_kl592_assist,
        ~n_sth_local          * ~sar_sth_local,
        ~n_sth_body           * ~sar_sth_body,
        ~n_jnb_unprotected    * ~sar_jnb_unprotected,
        ~n_jnb_protected      * ~sar_jnb_protected,
        ~n_or_tambo           * ~sar_or_tambo,
        ~n_other_infectious_pax * ~n_household_per_pax * ~sar_household,
    ]
    infections = sum(~sq.poisson(r) for r in pool_rates)
    if sq.event(P_SUPERSPREAD):
        infections += ~sq.poisson(~ss_burst)
    return infections


def define_event():
    """For each true infection: decide if it will ever be diagnosed-and-linked
    (link_rate), then draw incubation + reporting lag to determine WHEN it
    surfaces publicly. Count how many have surfaced by today and by Aug 1.
    link_rate sampled once per event so today vs Aug 1 share the linkage rate."""
    n_true = true_infections()
    link_rate = ~p_link
    days_elapsed = ~DAYS_SINCE_EXPOSURE
    reported_today = 0
    reported_aug1 = 0
    for _ in range(n_true):
        if not ~sq.bernoulli(link_rate):
            continue
        pipeline = ~incubation_dist + ~report_lag_dist
        if pipeline <= days_elapsed:
            reported_today += 1
        if pipeline <= days_elapsed + DAYS_TO_AUG1:
            reported_aug1 += 1
    return {
        "n_true": n_true,
        "linked_today": reported_today,
        "linked_aug1": reported_aug1,
    }


# =============================================================================
# Run: unconditional + Bayesian update on today's observation
# =============================================================================

N = 50_000

print("Metaculus #43461 -- non-passenger Hondius-linked cases by Aug 1, 2026")
print("=" * 70)
print(f"Samples: {N:,}\n")

print("[1/2] Sampling unconditional model...")
unc_samples = sq.sample(define_event, n=N, verbose=True)
unc_aug1 = np.array([s["linked_aug1"] for s in unc_samples])
p_unc = float(np.mean(unc_aug1 >= 5))
print(f"\nUnconditional P(>= 5 by Aug 1): {p_unc:.1%}")
print(f"  mean={np.mean(unc_aug1):.2f}, median={np.median(unc_aug1):.1f}, "
      f"P90={np.percentile(unc_aug1, 90):.1f}")

# Condition on: 0 non-passenger cases have been publicly linked as of today.
# Note: KLM flight attendant is hospitalized but no positive test yet, so
# treating her as "not yet linked" for the observation filter.
print("\n[2/2] Conditional on observed 0 linked non-passenger cases today...")
p_yes_cond = bayes.bayesnet(
    define_event,
    find=lambda e: e["linked_aug1"] >= 5,
    conditional_on=lambda e: e["linked_today"] == 0,
    n=N,
    verbose=True,
)

# Resample for percentile detail on the conditional distribution.
print("\nResampling conditional distribution for percentile detail...")
cond_aug1 = []
attempts = 0
while len(cond_aug1) < 20_000 and attempts < 500_000:
    e = define_event()
    attempts += 1
    if e["linked_today"] == 0:
        cond_aug1.append(e["linked_aug1"])
cond_aug1 = np.array(cond_aug1)
print(f"  Acceptance: {len(cond_aug1)}/{attempts} ({len(cond_aug1)/attempts:.1%})")

print()
print("=" * 70)
print("RESULTS (as of May 10)")
print("=" * 70)
print(f"P(YES | observed 0 today) = {p_yes_cond:.1%}")
print(f"P(YES) unconditional      = {p_unc:.1%}")
print()
print("Conditional distribution of linked_aug1:")
print(f"  mean   = {np.mean(cond_aug1):.2f}")
print(f"  median = {np.median(cond_aug1):.1f}")
for q in (10, 25, 50, 75, 90, 95, 99):
    print(f"  P{q:>2} = {np.percentile(cond_aug1, q):.1f}")
print()
print("P(linked_aug1 >= k | observed 0 today):")
for k in (1, 2, 3, 5, 7, 10):
    print(f"  >= {k:>2}: {np.mean(cond_aug1 >= k):.1%}")

# ======================================================================
# RESULTS (as of May 10)
# ======================================================================
# P(YES | observed 0 today) = 3.8%
# P(YES) unconditional      = 4.7%
#
# P(linked_aug1 >= k | observed 0 today):
#  >=  1: 64.8%
#  >=  2: 31.9%
#  >=  3: 14.3%
#  >=  5: 3.4%
#  >=  7: 1.4%
#  >= 10: 0.4%