Skip to content

Biostatistical-review response: parity/two-period detector, covariate model retired + oblast SMR attribution, config single-source, dead-code audit#3

Open
sgrabovyi wants to merge 44 commits into
mainfrom
biostat-review-fixes
Open

Biostatistical-review response: parity/two-period detector, covariate model retired + oblast SMR attribution, config single-source, dead-code audit#3
sgrabovyi wants to merge 44 commits into
mainfrom
biostat-review-fixes

Conversation

@sgrabovyi

@sgrabovyi sgrabovyi commented Jul 10, 2026

Copy link
Copy Markdown
Collaborator

Response to the biostatistical review, plus the methodology rework and cleanup it led to. One branch, reviewed as one PR.

Detector (methodology)

  • Parity, two-period detector. Both axes (level vs national, trend vs own past) are tested at parity — credibly above the reference, no arbitrary multiplier. The old 2.0 / 1.5 multipliers are gone (a simulation showed any multiplier above parity detects nothing on this sparse data).
  • Calling confidence 0.80, the Richardson et al. (2004) disease-mapping rule D(0.8,1), confirmed as the sensitivity/false-positive knee by validation/threshold_simulation.py.
  • Testing-intensity normalization in the detector, so site turnover (wartime closures/openings) does not leak into the trend.
  • Cold-cell gate: a unit with no recent events in either window has no trend to evaluate → it is Normal (grey), not a spurious "decrease".
  • Reliability = posterior CV behind a hard convergence gate (R-hat < 1.01, ESS ≥ 400 — Vehtari et al. 2021; > 5% divergences fatal). The ad-hoc 40/30/30 blend is gone.

Config & correctness

  • Single source of truth: DEFAULT_CONFIG is derived from config.json, so the two can never drift (that drift was the root cause of a stale-resolution run). Standardised on res3.
  • Fixed the broken scripted runner, aligned README / QUICK_START / HOW_TO_RUN with the code, tightened the map caption.

Covariate model retired → oblast SMR attribution

  • On this data the risk-group composition does not explain recency (the groups have near-equal recent-infection rates), so the covariate model's adjustment was a near-no-op. It is removed entirely (model, map, wiring).
  • Its explanatory intent is now met by an oblast-level indirect standardization written automatically to summary/Group_Attribution_Report.xlsx every run: the tested case-mix shift over time (high-risk share fell ~31% → ~18%, p≈1e-98), per-group recency rates with a significance test, and an SMR (observed / expected-given-mix, with CI) that cross-checks the detector — e.g. Zaporizka SMR = 4.25 [3.29–5.41].

Cleanup

Removed dead code and stale docs left by the retired models: the covariate/network/outbreak helpers, the dormant use_loo_ic / LOO-IC machinery (unreachable prompt), vestigial wizard prompts, and retired-model "tombstone" comments. Every step verified with a full green pipeline run.

🤖 Generated with Claude Code

sgrabovyi and others added 22 commits July 2, 2026 00:32
Territories with zero (or a single) recent events in the current window were
being misclassified as "Emerging hotspot": the SIR axis compares the current
rate to a territory's own history, and when that history is a zero on a large
test denominator the Empirical-Bayes baseline collapses to its floor, inflating
SIR by construction. Odesa/Kharkiv/Mykolaiv hexes with 0 recent events, and a
Poltava hex with a single event, surfaced as hotspots.

Adopt the two-part (presence x intensity) decomposition of Fu et al.
(Euro Surveill. 2023;28(14):2200581), implemented in our PyMC Beta-Binomial
framework:

- BayesianAnalyzer.run_two_part_model: zero-inflated Beta-Binomial via
  pm.CustomDist (marginalised over a latent presence indicator). Presence pi is
  shared across the current and baseline windows (separate intensities p, p_hist)
  so it is identified from pooled evidence; the presence-weighted current rate
  theta = pi * p drives SMR/SIR. Reports presence_prob. Enabled by the
  two_part_model config flag; routed in bayesian_dispatch.

- taxonomy: presence gate. A territory needs at least
  MIN_HOTSPOT_CURRENT_EVENTS (=2) current recent events, and presence_prob >= 0.5
  when available, to receive any increase-hotspot label; otherwise it is Normal.
  Epidemiologically, a recency hotspot must rest on observed current signal.

- excel_report: surface the presence_prob column.

- territory_info: deterministic tie-break in the hex -> community/district/oblast
  labelling (sort by admin id) instead of order-dependent iloc[0], so a hex whose
  centroid sits on a boundary keeps a stable label across runs. Counts were never
  affected; the case->hex aggregation is deterministic.

- config: enable two_part_model; set analysis_period to 2026-03-01..2026-05-31.

Validated on real data (2026-03..05, res4+res3, both models): 0 divergences,
zero phantom hotspots, presence_prob discriminates structural zeros (0.69) from
present territories (0.85); standard --test passes.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Ship a fully synthetic dataset in the input folder
(data/synthetic_input_data.xlsx) that stands in for the confidential case data,
and a public notebook (synthetic_data/synthetic_pipeline_public.ipynb) that
turns it into the file the pipeline runs on
(data/synthetic_input_data_double.xlsx) and demonstrates the geographic-privacy
tools (reversible coordinate encryption and donut geomasking) - all with no
real data at any point. Add a light data-exploration notebook, a step-by-step
getting-started guide (synthetic_data/README.md), and links from README.md and
data/README.md. The generation scripts are documented in English.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… tests, reliability caveat)

- taxonomy: gate hotspot labels on sir_informative. When local history is too
  thin to inform the trend axis, a still-elevated territory becomes the new
  "Elevated vs national (trend uncertain)" label (not counted as a hotspot),
  and a rise-only signal collapses to Normal. New label wired through the
  colour map, UA labels, summary, dashboard.
- stratified outbreak/composition tests: apply a family-wide Benjamini-Hochberg
  FDR correction across all territories (apply_fdr_correction) instead of a raw
  per-territory p < 0.05, and annotate the explanation as descriptive.
- reporting: add a standing low-reliability caveat to the Excel disclaimer sheet
  and metadata, and a visible footnote on the main map. Colours and labels are
  unchanged (all territories stay coloured).

Verified on synthetic data (exit 0): pipeline completes, BH correction fires
(18 outbreak + 9 composition tests, 1 retained), caveat renders in Excel and on
the map, all output files non-empty.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…e aggregation, LOO national rate)

- covariates fallback: fix the BayesianAnalyzer constructor call, which passed
  gdf_cases/analysis_start/analysis_end positionally and did not match
  BaseHotspotAnalyzer.__init__ (it raised TypeError whenever the fallback fired).
  Now uses (config, mode_suffix, orchestrator) and copies gdf_cases.
- retire the Truncated Binomial ("Hurdle") model: it fit a weakly-identified
  per-territory random slope (one Binomial observation per territory) under a
  less robust sampling configuration. Both run_hurdle_model methods are removed;
  the dispatch and the iterative orchestrator route any Hurdle request to the
  two-part zero-inflated Beta-Binomial model, which is the proper replacement.
- covariate territory-level aggregation: replace the unweighted mean across
  strata with exact pooling for proportions (recomputed from summed counts) and
  a denominator-weighted average for the posterior summaries and exceedance
  probabilities. Add the SMR/SIR taxonomy columns to the merge-back so reports
  and the watch-list rate axis see them (previously only the derived label was
  merged).
- add a leave-one-out national rate as a sensitivity option (smr_leave_one_out):
  each territory's SMR denominator excludes its own counts, so a large city can
  no longer partly mask itself. Wired through the two-part and standard paths.

Verified on synthetic data (exit 0): no Hurdle run, covariate model completes
and its report now carries the SMR/SIR columns, all output files non-empty.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…-honest names

- analysis levels: the interactive selection now offers only res3 (larger
  hexagons, more facilities per unit) and adm1 oblasts; res4 and finer
  resolutions are removed because most finer hexes hold zero or one active site.
  Default is res3 + adm1; the standard-mode fallback and the iterative-mode
  H3-id generation move from res4 to res3.
- report labels: the SMR/SIR-style columns keep their internal identifiers but
  are relabelled in the Excel report to "National reference ratio" (current vs
  national) and "Historical trend ratio" (current vs own history), since the
  classic SMR/SIR abbreviations wrongly imply a population-standardised
  mortality/incidence ratio. A new Glossary sheet defines every term, including
  "hotspot" as a testing-network signal (not a residential/population hotspot).

Verified on synthetic data (exit 0): both res3 (26 hexes) and adm1 (13 oblasts)
levels run, reports carry the new column names and the Glossary sheet, no res4.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
- add BayesianAnalyzer.run_two_period_model: fits both windows at once with
  logit(p_it) = mu_t + a_i + delta_i * I(current), reading the trend directly
  from a hierarchical per-territory current-period change delta_i instead of the
  two-step SIR-from-EB-baseline ratio. The "vs national" axis comes from the
  current-period rate, the trend axis from exp(delta_i). Enabled with the config
  flag two_period_model; the dispatch routes to it when set.
- make the crude-model diagnostics robust to model-specific parameter names
  (alpha/beta for the standard/two-part model, a/delta for the two-period model)
  so R-hat/ESS extraction no longer assumes fixed names.
- add validation/two_period_validation.py: a head-to-head simulation of trend
  detection (two-period delta vs legacy SIR) on data with a known injected trend,
  using the pipeline's own FDR cutoff and compute_smr_sir so the rules match
  production.

Finding: on the sparse-event regime the two-period model is properly
conservative -- it controls the false discovery rate and makes no unsupported
trend calls -- while the legacy SIR axis reaches higher nominal sensitivity only
by tolerating a high false-discovery rate. Neither reliably detects trends in
this power-limited regime, which reinforces the study's central power-limit
finding. The two-period model is therefore kept available (flag) but NOT made
the default; the trend axis is reported with explicit uncertainty.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… workbook

Two reporting bugs surfaced by auditing every output file:

- PPC plot for the two-part (crude) model was never produced: az.plot_ppc
  raises on the zero-inflated CustomDist predictive (ragged shapes / degenerate
  KDE). Add a robust manual fallback that overlays the observed recent-count
  distribution with the posterior predictive as plain histograms, so a PPC
  figure is always rendered. The standard-likelihood models keep the arviz plot.
- Diagnostics_<period>.xlsx was written once AFTER the level loop, so with more
  than one analysis level only the LAST level's folder received a workbook. Move
  the write inside the loop and pass each level's own diagnostics, so every level
  (res3 and adm1, crude and covariates) gets its Diagnostics workbook.
- Strip the heavy non-serialisable objects (model, trace, ppc, y_obs) from the
  diagnostics before writing the workbook, keeping it a clean scalar table.

Verified on synthetic data (res3 + adm1): both levels of both models now produce
a PPC figure and their own Diagnostics workbook; no failed-plot warnings.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…trong hotspots are detected

The formal classifier missed even extreme hotspots: an engineered hex with 45
recent infections out of 250 tests (18%, vs a ~2% national rate) was estimated
at ~5% and classified Normal (reference ratio ~1.9). Root cause: the kappa prior
Gamma(3, 0.2) (mean 15) forced the Beta-Binomial to be overdispersed even on
Binomial-like data, which weakened the likelihood so the hierarchical model
shrank every signal toward the national rate.

- Replace kappa ~ Gamma(3, 0.2) with a weakly-informative Gamma(2, 0.01) in the
  crude, two-part, two-period and covariate models, so the data pick the
  concentration (near-Binomial when there is no real overdispersion).
- Add a config-driven fast_sampling preset (default off) for validation runs.

Validated end-to-end: the same 45/250 hex is now classified "Stable high-burden"
(a hotspot) with reference ratio 8.24 and an estimated rate of 13%, while a
normal background hex stays Normal (no new false positives).

Adds validation/taxonomy_validation.py, which drives one engineered target per
taxonomy category through the pipeline and tabulates how each model classifies it.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…presence-weighted theta

Even after loosening the concentration prior, a clear hotspot could still read as
Normal: the presence probability pi is shrunk below 1 (about 0.82) even for a hex
with 45 recent events, so the presence-weighted rate theta = pi * p deflated the
SMR and widened its posterior, dropping the exceedance probability just under the
FDR cutoff (e.g. a 45/250 = 18% hex gave exc_prob_smr 0.947 vs a 0.950 cutoff).

Feed the intensity p (the recency rate given presence) into the SMR/SIR
computation, the exceedance probability and the reported rate; keep presence_prob
(pi) only for the taxonomy presence gate. Structural zeros are already excluded
from the hotspot set by that gate (>=2 recent events / presence_prob > 0.5), so
the pi weighting on the rate was redundant as well as harmful.

Validated on the taxonomy replica (faithful two-part fit + the pipeline's own
compute_smr_sir / classify_with_smr_sir): the 18% hex now clears as a hotspot
(exc_prob_smr 0.96) while normal and low-count hexes stay Normal.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…sing / stable / falling

The trend (SIR) axis is power-limited: with a handful of recent events per window,
the posterior of the current-vs-history change is too wide to clear the
FDR-controlled SIR exceedance (P(rise) sat at 0.6-0.7 even for a clear 2x move).
As a result every confirmed high-burden hex collapsed to "Stable high-burden" --
a genuinely rising Established hotspot and a genuinely falling Declining hex were
indistinguishable.

For a hex ALREADY confirmed elevated vs national (smr_state 'high') and with an
informative local history, describe the trend from the observed current-vs-history
rate ratio (a lean, not an FDR discovery): >=1.3x -> rising, <=1/1.3x -> falling,
else stable. The strict exceedance still governs everything else, so no trend is
invented on a territory that is not a confirmed hotspot, and the
sir-informativeness gate still withholds the trend on thin history
("Elevated vs national (trend uncertain)").

Validated on the taxonomy replica: an 8->45 hex reads Established hotspot, 45->45
reads Stable high-burden, 60->30 reads Declining from high-burden -- the three were
previously all "Stable high-burden".

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…national rate

Replace the arbitrary SMR>2 / SIR>1.5 exceedance thresholds with a parity standard: a cell is elevated when its posterior rate credibly exceeds the national rate (ratio > 1), and rising/falling when the two-period change delta credibly differs from zero -- no invented multipliers on either axis. Drive the detector from the joint two-period model (activated; two-part retired), and compute the national comparison with leave-one-out so a dominant hotspot no longer masks itself in its own denominator on sparse recency counts.

config.json: detection thresholds set to 1.0 (parity), including the low side; two_period_model on, two_part_model off; smr_leave_one_out on; sampling cores 1. bayesian.py: the two-period SMR call now reads smr_low_threshold from config. taxonomy.py: remove the raw-ratio soft trend 'lean' (the last 1.5x); the trend axis is now read solely from the FDR-controlled change exceedance.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ance)

The posterior-probability confidence for a hotspot call was hard-coded at 0.95. Make it a run parameter (detection.confidence_level) so a sparse-data surveillance run can declare 0.90 -- trading a higher false-discovery share for the power to confirm signals that 0.95 cannot on a handful of recent events. The presence gate (>=2 recent events) still holds, so 0.90 never flags a single-event cell.

base.py: thread confidence_level into the FDR-controlled cut-off computation (anchors both start and floor). interactive.py + interactive_setup.py: add a 90/95 dialog prompt (default 90) to both standard and iterative setup, alongside the existing analysis-window prompt.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…e, period tag

Three fixes in the iterative sliding-window path. First, the iterative dispatcher never knew about the joint two-period model and silently fell back to the standard model; it now runs the two-period model when configured. Second, a for...else on the level loop always logged 'No hotspots found in any iteration' after every run, contradicting the report that was actually written; the spurious clause is removed. Third, period_str is now derived from the analysis window when the interactive setup did not set it, so config-driven report filenames read _YYYYMM instead of _None.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…1) rule)

A simulation on the real site volumes (validation/threshold_simulation.py) traced the false-positive / sensitivity trade-off across confidence levels at parity: the false-positive rate stays ~0 down to 0.85 and is still ~1% at 0.80, while sensitivity to a 3x outbreak roughly doubles (0.26 -> 0.50) between 0.90 and 0.80. 0.80 is the knee, and it coincides with the established disease-mapping decision rule D(0.8, 1) of Richardson, Thomson, Best & Elliott (Interpreting posterior relative risk estimates in disease-mapping studies, Environ Health Perspect 2004;112(9):1016-1025), which flags an area when Pr(relative rate > reference) >= 0.80. Set the config default to 0.80, cite the rule in base.py, and offer 80/85/90/95 in the setup dialog (default 80).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…e vestigial 40/30/30 weighting

The per-territory reliability that drives the map is the posterior CV (interval width) plus a hard convergence gate; the 40/30/30 data/sample/model blend fed nothing (calculate_overall_score was dead code and comparison_report only printed the weights as cosmetic labels). Remove the blend entirely -- the reliability_weights config block, the dead calculate_overall_score method, and the weight fields in the report -- so there is no arbitrary component weighting left to defend.

Re-point validation/reliability_calibration.py at the parity threshold (rate > national, not SMR > 2) and the project national rate. It confirms the 80/60 cut-offs separate actionable from non-actionable estimates cleanly: HIGH cells give a decisive interval 100% of the time, MODERATE 61%, LOW only 3%.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
validation/threshold_simulation.py generates recent counts under a known ground truth on per-hex testing volumes (synthetic input by default; SIM_DATA points at the real file), fits the joint two-period model, and traces the false-positive / sensitivity trade-off across the SMR ratio and confidence grids. It is the evidence behind two decisions: the rate axis stays at parity (any multiplier above 1 catches nothing on this sparse data) and the confidence level is 0.80 (the sensitivity/false-positive knee, matching Richardson et al. 2004 D(0.8,1)). Portable paths; runs from the repo without arguments.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The two axes are now described honestly (level vs country at parity, trend from the joint two-period delta) instead of SMR>2 / SIR>1.5; the config example and the detection notes carry the parity thresholds, confidence_level 0.80 (Richardson D(0.8,1)), smr_leave_one_out and two_period_model; the Models table marks the joint two-period model as the primary detector and records that the Hurdle and two-part branches are retired. .gitignore excludes the simulation run artifacts. (The Classification narrative section still uses the SMR/SIR wording and is queued for a follow-up pass.)

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…odology)

Rewrite the remaining narrative to the honest parity framing: the seven-label cross is trend x level, new units are classified on the level axis only, the report/watch-list columns are the National reference ratio (internal smr_*), and the Methodology section states both axes are tested at parity with confidence 0.80. No functional change.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…or-SMR mention

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The levels prompt now offers numbered options -- 1 = res3, 2 = adm1, 3 = both -- so a user can type an index instead of the level names. The comma-separated name input (res3,adm1) still works, and the default (Enter) stays 'both'.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…variates / Both (1/2/3)

Auto used to pick between the retired GLM and the Bayesian model, which no longer applies -- the Bayesian model is always the detector. The prompt now offers only 1 = Bayesian only, 2 = Bayesian with Covariates, 3 = Both (default). 'both' leaves neither force-flag set, so the orchestrator runs the detector and the explanatory covariate layer side by side, exactly as the old default did.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…e detector

The explanatory covariate layer passed only smr_threshold / sir_threshold to compute_smr_sir, so its below-national cut defaulted to 0.5 and its national reference was not leave-one-out -- inconsistent with the crude detector now that both are parity + LOO. Pass smr_low_threshold and smr_leave_one_out from config so the covariate classification uses the same reference as the detector it explains.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@sgrabovyi sgrabovyi changed the title Biostatistical-review improvements: honesty gates, bug fixes, res3+ADM1, two-period model, output audit Biostatistical-review improvements: parity/two-period detector, calibrated 0.80 confidence, reliability + terminology cleanup Jul 13, 2026
sgrabovyi and others added 2 commits July 13, 2026 21:47
…dispatch branches

The two-part zero-inflated Beta-Binomial model is retired (the joint two-period model is the detector; zeros are handled by the denominator filter + presence gate), so its ~270-line run_two_part_model method was dead code. Delete it, drop the two-part branch from the standard and iterative dispatchers, and remove the now-unused level_use_hurdle / level_hurdle_threshold parameters from run_bayesian_dispatch. The remaining Hurdle wizard/CLI plumbing (inert) is removed next.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
With the two-part model gone, the Hurdle flag it routed to was dead threading: the --use-hurdle / --hurdle-threshold CLI args, the orchestrator constructor params, the wizard's structural-zero use_hurdle computation, the audit-trail 'Hurdle Model Configuration' decision, and the DEFAULT_TRUNCATED_BINOMIAL_STRUCTURAL_ZEROS_PCT constant. All removed; the wizard/audit now carry only use_loo_ic. Verified end-to-end (--test, both models, exit 0).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@sgrabovyi
sgrabovyi changed the base branch from two-part-model-zero-handling to main July 13, 2026 19:12
sgrabovyi and others added 3 commits July 13, 2026 23:04
…wo-period, LOO)

DEFAULT_CONFIG (used by --test and any config-less run) was stale: hex_resolutions was [4] and detection still carried the old SMR>2 / SIR>1.5 with no confidence_level, no leave-one-out and no two_period_model. So --test silently ran the OLD methodology on res4. Set hex_resolutions to [3]; make detection parity (1.0) with the low thresholds and confidence_level 0.80; add smr_leave_one_out, two_period_model=true, two_part_model=false -- mirroring config.json. Verified: --test now runs res3 with the joint two-period model and parity cut-offs, no res4, exit 0.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Wave 3 removed res4 only from the interactive level menu; it lingered elsewhere and could still run. Fixes: (1) iterative_resolution defaulted to 4 in four places (orchestrator x2, first_window_calibration, iterative_hotspots_report) -> 3, so a config-less iterative run no longer silently uses res4; (2) removed every res4/res5 reference from code and config -- LEVEL_NAMES, the recommendations res4 branch, the analysis-type/iterative dialog text, docstrings/comments, and the unused res4_path/res5_path (both h3_hexagons blocks reduced to res3_path + h3_id_col); (3) removed the dead two_part_model config key (nothing reads it after the two-part model was deleted); (4) removed sampling.cores=1 from config.json -- a debugging leftover that forced single-core and slowed real runs; sampling now adapts. Verified: grep for res4/res5 is empty, and --test runs res3 with the two-period parity methodology (output res3 only, exit 0).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…json

The stale-res4 run was caused by two hand-maintained config copies drifting:
config.json (fixed to res3) and a duplicated DEFAULT_CONFIG literal in
constants.py (still res4 + old thresholds). --test read the second one, so it
silently ran retired methodology.

DEFAULT_CONFIG is now LOADED from config.json at import (_load_repo_config,
path resolved relative to the module, not the CWD) and overlays only the two
fields the wizard would otherwise set for a headless run: hex_resolutions and
admin_levels. Drift between the two is now impossible by construction.

Side effects fixed in the same pass:
- The UA map/oblast label tables moved out of the config literal into code
  constants MAP_STRINGS_UA / OBLAST_NAMES_UA; the three renderers fall back to
  them, so interactive (config.json) runs finally get Ukrainian labels too
  (previously only --test did, via DEFAULT_CONFIG).
- Added the missing "Elevated vs national (trend uncertain)" colour to
  config.json so the single source is complete.
- Removed the config_universal.json phantom that --test probed first, and the
  now-unused os import.

Verified: DEFAULT_CONFIG == config.json on every shared key; full pipeline
(--test methodology, one core) runs Hex_Res3 + Joint Two-Period, res3-only
output, exit 0.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
sgrabovyi and others added 17 commits July 14, 2026 09:34
…ighten map caption

Systematic comb-through after the res4/config-drift incident, second pass.

Confirmed break:
- validation/service_run.py passed use_hurdle/hurdle_threshold kwargs the
  orchestrator no longer accepts -> TypeError on launch. Removed both (and the
  --use-hurdle/--hurdle-threshold flags); the scripted runner works again.

Retired-model dead code / stale references:
- taxonomy.py: removed the presence_prob gate (no model populates presence_prob
  since the two-part model was deleted; the branch was always inert). The
  minimum-current-events presence rule stays.
- excel_report.py: removed the presence_prob report column (never emitted now).
- taxonomy_validation.py: rewritten to run only the production two-period model
  (it previously looped over the deleted two_part model).
- data_quality_audit.py: the audit trail no longer advises "Consider Hurdle
  model"; it describes the actual handling (Beta-Binomial + presence gate).
- wizard_and_audit.py docstring: now describes the two decisions it records
  (spatial structure, model selection), not the removed Hurdle decision.
- Comment cleanups (base/bayesian/covariates/diagnostics/plots_guide): drop
  references to the retired hurdle/two-part fits; keep the honest "retired"
  tombstones and the (behaviour-agnostic) robust-PPC fallback.
- test_convergence_gate.py: Hex_Res4 label -> Hex_Res3.

Consistency:
- covariates.py hardcoded cores=4; now honours config['sampling']['cores'] the
  same way the crude model does (defaults to 4 when unset, so the user's
  config.json behaviour is unchanged). One knob controls both sampling paths.

Map:
- anomaly_map.py: anchor the data-quality caveat directly beneath the map axes
  (axes fraction) instead of the figure bottom. Ukraine is ~2:1 wide, so
  equal-aspect leaves an empty band under the map; the figure-bottom caption
  sat far below it. Now flush under the map, whitespace trimmed by tight bbox.

overall_score: reviewed, NOT changed -- it is the honest mean per-territory
reliability (already documented as un-blended), not a leftover of the retired
40/30/30 composite.

Verified: all changed files compile and import; full pipeline runs green
(Hex_Res3, Joint Two-Period, res3-only output); map caption renders flush under
the map.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…eality

The narrative docs still described machinery the code no longer has, which is
exactly what would confuse a CDC reader diffing docs against code.

- res4/res5 as selectable/primary grids -> res3 is the standardised grid;
  res4/res5 geojson kept in data/ for reference only, not used. Fixed in
  README, docs/QUICK_START, docs/HOW_TO_RUN, data/README.
- Removed the --use-hurdle / --hurdle-threshold rows and examples: those flags
  no longer exist in run_hotspots.py (only --test, --use-loo-ic, logging).
- config h3_hexagons example: res4_path -> the real res3_path/h3_id_col.
- resolution_sigma_multiplier example Hex_Res4 -> Hex_Res3.
- posterior SMR -> posterior National reference ratio (internal smr_*).

Same-class res4 leftovers fixed in code:
- multiseed_stability.py --resolution default 4 -> 3 (and its usage example).
- iterative_windows.py / orchestrator.py docstrings said "attach H3 res-4 IDs"
  while the code attaches res3 IDs; corrected.
- orchestrator large-grid warning recommended "Use H3 Res 4 instead" (nonsense,
  the warning fires AT res4+); now recommends the coarser res3.

The Hurdle/two-part "retired" tombstones and the res4+ guard comment are left
as accurate history. All changed .py files compile.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The covariate model adjusts the level for each unit's high-risk testing share,
but that share was only logged, not written per unit. It is now carried from the
stratified frame into the per-territory frame, merged back onto gdf_admin, and
added to the covariate Excel report, so a reader can see the risk-group
composition behind (or, on sparse data where beta_risk is ~0, not behind) each
covariate-adjusted level. Crude reports do not carry the column, so the
conditional in excel_report leaves them unchanged.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Sites enter and leave the network (wartime closures/openings), which moves raw
testing volume for reasons unrelated to incidence. The detector now normalises
for it, so the trend is net of testing effort.

- aggregate_stats (territory.py): compute per-territory test-months (tests per
  active site-month) for both windows, reusing the count spatial joins. Adds
  testing_intensity_curr/hist and n_active_months_curr/hist to the main frame
  (previously these lived only in the covariate/stratified frames).
- run_two_period_model (bayesian.py): standardized log test-months enter the
  LEVEL of each window (pooled standardization, one shared beta_intensity), so
  equal effort maps to an equal shift and the trend delta is net of effort.
  Toggle with config 'intensity_adjustment' (default on); logs beta_intensity.

Gold-standard validation (intensity ON vs OFF, res3, 1 core, real data):
- beta_intensity = -0.350 [-1.012, 0.344] -- weak, CI spans zero on this data.
- The persistent Dnipropetrovska hotspot (2 recent / 54) stays Stable
  high-burden both ways (level 2.92 -> 2.98).
- Only 2/21 labels flip, both zero-event hexes moving Significant decrease ->
  Normal (conservative). No hotspot gained or lost.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…l (grey)

The increase/decrease axes describe a TREND, but a unit with zero recent events
in BOTH the analysis and the baseline window has no trend to evaluate: there is
nothing to have risen or fallen from. Such units were being labelled "Significant
decrease" purely because their current rate is credibly below the national rate
(the ('norm','low') and ('low','low') label-map cells), which reads as a downward
trend that does not exist — e.g. a Mykolaivska hex with 0/164 recent in history
and 0/36 now, ~1.5 years without a single recent infection.

Add a gate mirroring the hotspot presence gate: recent_count_curr == 0 AND
recent_count_hist == 0 -> "Normal". A genuine fall (history had recent events,
current has none) is unaffected because recent_count_hist > 0 and keeps its
decrease label.

Verified by direct classification of crafted rows: cold cell -> Normal; genuine
fall (5 hist, 0 curr) -> Significant decrease; real hotspot -> Established
hotspot; 0 tested -> No Data.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…le detector

On this data the risk-group composition does not explain recency -- the two
groups have near-equal recent-infection rates (2.9% high-risk vs 2.0% general),
so indirect standardization gives a composition ratio of ~1.0 everywhere and the
covariate beta_risk is ~0. The covariate model's adjustment target is thus a
near-no-op, and the one function it carried uniquely (testing-intensity
normalization) already moved into the main detector. So it is retired.

This commit stops it running and drawing its map:
- interactive.py: remove choose_model_selection (no model choice remains).
- interactive_setup.py: fix manual_model_selection = 'bayesian'.
- config.json: manual_model_selection = 'bayesian' (so --test/scripted runs are
  detector-only too).

Verified: pipeline run logs "Skipping Bayesian Covariates (Bayesian only mode)",
output has only the bayesian/ subtree, no bayesian_covariates/, exit 0. Dead
covariate code is deleted in the follow-up commit.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Follow-up to the flow-level disable. Now that nothing runs or imports the
covariate model, remove it.

Deleted:
- pipeline/analyzers/covariates.py, _covariates_runtime.py
- pipeline/orchestration/bayesian_covariates_dispatch.py,
  bayesian_covariates_plots.py, interpret_bayesian_covariates.py

orchestrator.py: dropped the covariate analyzer instantiation, the covariate
dispatch/interpret/plots/maps/finalize block in run_for_mode, the per-level
comparison-report call (it compared two models; there is one now), the covariate
diagnostics-workbook save, the manual model-selection flags
(_force_bayesian_only / _force_bayes_cov_only), and the now-dead bayesian_cov /
bayesian_cov_interpretation variables. Removed the covariate imports.

__init__ exports: dropped BayesianCovariatesAnalyzer (analyzers) and the three
covariate orchestration functions (orchestration).

Verified: package imports; full pipeline runs green (one core), exit 0, zero
tracebacks, only the bayesian/ output subtree, logs "Skipping Bayesian
Covariates". A few harmless covariate-only helpers remain in aggregation
(aggregate_stats_stratified/hard, aggregate_covariates) and base.py and can be
tidied in a follow-up; they are dead but not wired to anything.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Packages the explanatory analysis that replaces the retired covariate model, as
pipeline.standardization.group_attribution + a validation runner.

Four products, each with its honest power caveat:
- period_composition: the tested case-mix shift over time (well powered) --
  high-risk share fell from ~31% pre-2024 to ~18% in 2024+ (chi2 p~1e-126,
  monotonic, z=21).
- escalating_group_rates: national group recency rates at yearly -> 2-year ->
  whole windows; the 2025-2026 tail needs a 2-year pool to clear >=10 recent/group.
- oblast_group_rates: per-oblast high-risk vs general recency rate with a
  two-proportion test. Where powered, the high-risk group is significantly hotter
  -- Zaporizka 16.6% vs 6.6% (2.5x, p<0.001), Kyiv 3.4% vs 1.8% (1.9x, p=0.007);
  Dnipropetrovska 2.7% vs 2.1% is not significant. (This corrects the earlier
  over-broad "groups are similar" read, which was the national average.)
- oblast_standardization: indirect standardization. comp_ratio (~1 everywhere,
  because national group rates are close) is the mix effect; smr = observed/
  expected tracks the detector's National reference ratio, so it CROSS-CHECKS the
  detector rather than adding a new axis. Whole-pool smr: Zaporizka 4.25
  [3.29-5.41], Mykolaivska 2.08 [1.37-3.03], Odeska 0.48 [0.31-0.71].

Two distinct findings kept separate: the high-risk group IS hotter (per-oblast,
where powered), but the oblast EXCESS is not explained by composition (it sits in
both groups). Methods: Kitagawa (1955); indirect standardization / SMR (Breslow &
Day 1987; Rothman-Greenland-Lash); pooling-for-power after Fay & Herriot (1979).

Runner: validation/group_attribution_report.py prints all oblasts and writes
Group_Attribution_Report.xlsx.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Follow-up to the covariate-model removal. These helpers had no live caller once
the covariate analyzer was deleted:

- aggregation/territory.py: removed aggregate_covariates, aggregate_stats_stratified,
  aggregate_stats_hard_stratified, detect_outbreak_and_artifact, apply_fdr_correction
  (and their FDR helpers) -- ~970 lines. territory.py now holds only aggregate_stats.
  Also removed a REDUNDANT testing-intensity block I had added in step 1: aggregate_stats
  already computes testing_intensity_curr/hist via calculate_testing_intensity further
  down, which overwrote it. The two-period model reads that existing column; re-validated
  (beta_intensity = -0.342 [-0.99, 0.38], unchanged from before).
- diagnostics/bayesian.py: removed calculate_covariates_diagnostics and its stratified
  variant (~358 lines); kept calculate_bayesian_diagnostics.
- Deleted diagnostics/multicollinearity.py.
- base.py: removed the covariate output dir, the MODEL_TYPE=="bayesian_covariates"
  branch, and the aggregate_covariates / aggregate_stats_stratified / hard /
  detect_outbreak_and_artifact wrappers + imports.
- Pruned the now-orphaned exports and imports; fixed the module docstrings.

Verified: full package imports; pipeline runs green (one core), exit 0, zero
tracebacks, only the bayesian/ output. A few harmless network/outbreak wrappers
(analyze_site_profile, soft_fallback_result, etc.) remain wired through base.py
and can go in the upcoming full audit.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
- Removed choose_parametrization: with bayesian.auto_select_parametrization (on
  by default) the pipeline picks non-centered/centered from the data and ignores
  the manual answer, so the prompt was vestigial. The wizard now seeds
  'non_centered' and the auto-selector decides.
- Deleted three dead methods no longer called by any flow: choose_hex_resolutions,
  choose_iterative_resolution (both superseded by choose_levels), and choose_period
  (superseded by choose_analysis_period / choose_iterative_date_range).

The model-selection prompt was already removed with the covariate model. The
wizard now asks only what still has a real choice: analysis type, levels, window/
period, and calling confidence.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The site-profile / network-change / outbreak-fallback helpers were only reached
through the deleted covariate outbreak-detection path; no live code calls them.

- testing_network.py: removed analyze_site_profile, analyze_network_change,
  generate_network_explanation (340 -> 152 lines). Kept the two live functions the
  detector uses: calculate_testing_intensity and classify_network_stability.
- Deleted aggregation/outbreak_defaults.py (soft_fallback_result).
- base.py: removed the analyze_site_profile / analyze_network_change /
  generate_network_explanation / _get_soft_fallback_result wrappers and their imports.
- Pruned the aggregation __init__ exports; fixed the module docstrings.

Verified: full import; pipeline runs green (one core), exit 0, zero tracebacks,
only the bayesian/ output.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
An Explore-agent audit found covariate-only columns (high_outbreak, low_outbreak,
testing_artifact, outbreak_type, proportion_high_risk, imputed_proportion_high_risk)
that are read but never written now that the covariate model is gone, plus docs
that describe the removed model as live.

Dead read-branches removed:
- recommendations.py: the bayesian_covariates/ report + diagnostics globs, the
  high_outbreak/low_outbreak hotspot branch, and the outbreak_type reporting block.
- excel_report.py: the proportion_high_risk column append (I had added it for the
  covariate report) and the high_outbreak/testing_artifact column block.
- reliability.py: the imputed_proportion_high_risk penalties (three sites).
- base.py: the high_outbreak save_report debug logs.

Stale docs fixed:
- base.py / bayesian.py class docstrings (no longer describe a parallel covariate
  analyzer or a bayesian_covariates MODEL_TYPE).
- diagnostics/__init__ (dropped "multicollinearity checks").
- config/interactive.py (docstring now lists the prompts that actually exist:
  analysis type, levels, window/period, confidence -- not hex/parametrization/
  model-selection).
- excel_report.py comment; reliability.py imputation-penalty note.

Kept legitimate "retired Hurdle/two-part" tombstones. Verified: full import;
pipeline runs green (one core), exit 0, zero tracebacks.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
ModelConfigurationWizard asked two questions; the data-type one (site vs
residence) was discarded -- _determine_config ignored it and data_type is set to
'facility_based' in interactive_setup regardless. Removed it. Kept the analysis-
mode question: it toggles use_loo_ic, which still drives real research
diagnostics (the Binomial-vs-Beta-Binomial overdispersion test and sigma
calibration), so it is not dead. Fixed the docstring accordingly.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The "Standard vs Research" prompt that set use_loo_ic was UNREACHABLE: run_for_mode
always built a truthy cli_args ({'use_loo_ic': ...}), and ModelConfigurationWizard
skips its interactive block whenever cli_args is truthy -- so the prompt never
showed. use_loo_ic was therefore only ever set by the --use-loo-ic CLI flag (off by
default), and all it did was run a Binomial-vs-Beta-Binomial LOO check that always
recommends Beta-Binomial -- the likelihood the detector already uses. A dormant
diagnostic that confirms the existing choice.

Removed end to end:
- ModelConfigurationWizard (deleted pipeline/config/wizard.py) and its config export.
- _run_loo_ic_model_selection (orchestrator) and its use_loo_ic __init__ param.
- diagnostics/overdispersion.py (compare_binomial_betabinomial), used only there.
- The use_loo_ic branch in specification_analysis (heuristic AutoSpecificationSystem
  is now the sole path); wizard_and_audit no longer runs a wizard, just records the
  standing exchangeable-structure decision and returns pct_structural_zeros.
- The --use-loo-ic CLI flag from run_hotspots.py and validation/service_run.py.

Kept: the always-on single-model LOO diagnostic in diagnostics/bayesian.py (loo_elpd)
and the Cameron-Trivedi overdispersion test in spec/auto.py -- neither is tied to the
removed flag.

Verified: full import; pipeline runs green (one core), exit 0, zero tracebacks.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… code

Per feedback: the source should describe what the code IS, not narrate what was
removed -- git history holds the "why". Removed the "retired Hurdle / two-part /
zero-inflated / covariate model" archaeology from docstrings and comments in
bayesian.py, base.py, bayesian_dispatch.py, orchestrator.py, current_results.py,
interactive_setup.py. Kept comments that describe live behaviour (zero-inflation
thresholds, leave-one-out, the exchangeable-structure rationale).

Dead code that surfaced alongside:
- force_bayes_cov_only param + the "Bayesian Covariates only" branch in
  bayesian_dispatch (the orchestrator always passed False since the covariate
  model was removed).
- manual_model_selection config key: written (config.json, interactive_setup,
  taxonomy_validation) but never read once the model-selection branch was removed.
- comparison_report.py (generate_comparison_report) and
  DiagnosticInterpreter.generate_recommendations_report -- a crude-vs-covariate
  side-by-side report whose only caller was removed with the covariate model.

Verified: full import; pipeline runs green (one core), exit 0, zero tracebacks;
config.json and DEFAULT_CONFIG still identical on every shared key.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The oblast-level indirect standardization (the SMR that shows, with a CI, which
oblasts are hot beyond their testing mix) previously lived only in a standalone
validation script. Now the pipeline writes summary/Group_Attribution_Report.xlsx
automatically after the recommendations, and logs a one-line summary.

- New pipeline/orchestration/group_attribution_report.py: loads cases + ADM1,
  spatial-joins to oblasts, runs the four group_attribution products
  (composition shift over time, escalating-window group rates, per-oblast group
  rates + two-proportion test, per-oblast standardization for the current window
  and the whole 2021-2026 pool), writes the workbook. Skips gracefully (logs why)
  when the data lack a risk_group column, coordinates, or an ADM1 layer.
- Wired into run_full_pipeline, wrapped so a failure never aborts the run.
- README Output section documents the new report.

Verified on the real data: run is green, exit 0, and the log shows
"Composition shift pre/post-2024: 31%->18% high-risk (p=2.9e-98). Oblasts
credibly hot beyond composition (pooled): Zaporizka SMR=4.25 [3.29-5.41],
Mykolaivska SMR=2.08 [1.37-3.03]." The workbook lands in summary/.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The README output tree and HOW_TO_RUN still listed a bayesian_covariates/ subtree,
left over from the retired covariate model. No run produces it -- output is only
bayesian/ + summary/. Removed both references. (No code change; the covariate
model is genuinely gone.)

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@sgrabovyi sgrabovyi changed the title Biostatistical-review improvements: parity/two-period detector, calibrated 0.80 confidence, reliability + terminology cleanup Biostatistical-review response: parity/two-period detector, covariate model retired + oblast SMR attribution, config single-source, dead-code audit Jul 14, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant