# Decisions Stack — ABM Atlas Data-Driven Reclustering

Chronological record of every methodological decision, the rationale, and the
evidence behind it. Newest sections at the bottom. This is the audit trail for
the manuscript Methods section and for reviewer rebuttal.

Guiding principle throughout: **replace the rule-based 4-tier expert cascade
(T1/T2/T3/T4) with data-driven clusters, using only raw measurements and
continuous transformations — never the cascade's own arbitrary thresholds.**
Agreement with the cascade is a *diagnostic*, not a target to engineer toward.

---

## 0. Framing decisions (from plan mode)

| # | Decision | Rationale |
|---|----------|-----------|
| 0.1 | Full replacement of TIER, not a supplement | User directive; the cascade's hard cutoffs are what we are testing |
| 0.2 | Report agreement vs the cascade as the primary diagnostic | If unsupervised matches a threshold-rich cascade, the thresholds are data-supported; if not, they are arbitrary |
| 0.3 | Four methods side by side: k-means, GMM, LCA, HDBSCAN | Hard/soft, parametric/density — robustness across model families |
| 0.4 | 100-seed consensus clustering, co-cluster matrix | Turns hard clusterers soft; row-level stability = "różnicujący nie segregujący" |
| 0.5 | Free k (data picks k via silhouette/BIC/HDBSCAN-auto) | No EPIC prior on cluster count; biology is a post-hoc overlay |
| 0.6 | Three label columns: neutral CLUSTER_ID, EPIC letter, proposed name | Observe-then-propose; EPIC = Effect/Predisposition/Indication/Comorbidity |
| 0.7 | EPIC anchors from the 18-pair BENCHMARK_validation sheet | Biological ground truth independent of the cascade |

---

## 1. Forbidden inputs (threshold-purity rule)

Established after auditing the cascade docs. **Never used as features:**

- Cascade step-decisions: `|Z_diff_BETA|<2`, `AMP_RATIO>1.5`, `Tissue∈{specific,enriched}`
  as a binary, `FAERS=STRONG` as a binary, `N_INDEPENDENT≥3` as a binary, BNF
  override-priority logic.
- Cascade outputs: `TIER`, `NOVELTY/RESEARCH/APPLICATION_SCORE`, `TOTAL_GEOMEAN`,
  `FLAGS`, `TIER_RATIONALE`.

**Allowed:** raw evidence columns (the underlying measurements the cascade also
consumes), categorical states (BNF membership, burden-suffix presence), and
continuous arithmetic of raw values (sums, products, log-ratios, min/max/spread).

Key correction made mid-project: `N_INDEPENDENT` was initially excluded, then
found to be a **cascade-derived sum of three thresholded layers** (FAERS-STRONG +
OR/P_test-replication + tissue-mechanism). Re-adding it (v2) inflated agreement by
smuggling the cascade's conjunction back in — so it was permanently removed from
v3 onward and replaced with continuous evidence interactions.

---

## 2. Feature engineering decisions (v1 → v12)

### v1 — baseline (34 features, RobustScaler)
- Continuous GWAS stats + ordinal calls + indicator flags.
- Result: k=4 forced ≈ 42% TIER agreement, T4 recovery 0/26 across all methods.
- Finding: T3/T4 boundary collapses; the cascade's conjunction has no natural
  shape in feature space.

### v2 — added log1p(N_INDEPENDENT) — REVERTED
- Jumped to 60% / T4 92%, but N_INDEPENDENT is cascade-derived → **rejected** as
  smuggling the cascade's three-way AND.

### v3 — pure raw + continuous interactions (31 features, StandardScaler)
- Dropped the 3 cascade-derived stats (`Z_diff_BETA`, `AMP_RATIO`,
  `N_INDEPENDENT`); restored raw GWAS constituents.
- Added continuous interactions (no thresholds):
  `evidence_sum/product/min` over z(FAERS_IC025), z(−Tissue_rank), z(P_test);
  `drug_specificity = log|BETA_reg| − log|BETA_ICD10|` (continuous AMP_RATIO);
  `effect_magnitude_AND`, `effect_concordance`.
- Switched RobustScaler → **StandardScaler** (proper unit variance for binaries).
- Result: T4 recovered to 26/26 (kmeans/ward) **without any cascade feature** —
  `evidence_min` (strict-AND proxy) is the key, Cohen d=+1.6 for T4.

### v4 — T3-targeted additions
- `drug_specificity`, `evidence_spread = max−min` (partial-conjunction proxy).
- T3 did not move (stayed ~18%) — diagnostics showed no single feature separates
  T3 from BOTH neighbours.

### v5 — diagnostics-validated T3 features
- Added `evidence_max`, `beta_per_pleio_log = log1p|BETA_reg| − log1p(n_PheWAS_hits)`,
  `max_x_spread`. `beta_per_pleio_log` is the only feature separating T3 from T4
  (d=+0.70). Bought ~2pp; T3 ceiling confirmed at ~18% under distance methods.

### v6 — removed human annotations and categorical flags
- **Dropped** `curator_KNOWN_ADR`, `curator_NOT_INDICATION` (human annotation,
  not raw data), `is_burden` (categorical flag), `zero_carrier` (binary).
- **Added** `A1_CASE_CT_log1p` (continuous carrier count, replaces binary
  zero_carrier).
- Confirmed REGENIE p-value is already present as `LOG10P_disc`.
- Cost ~4pp agreement (curator notes were doing real work) but the 7 clusters
  became more biologically interpretable and annotation-free.

### v6.1 — renamed TEMP → ADR
- `TEMP_pct_after/before` → `ADR_pct_after/before` (clearer: temporal ADR
  direction; high after = drug-induced). Cosmetic; source Excel unchanged
  (read-only).

### v7 — external GWAS enrichment (44 features)
- Pulled per-variant BETA/SE/LOG10P from two additional REGENIE phenotypes for
  each of the 21 drugs: `<drug>__median_dose__ln` and `<drug>__prescribed`
  (7TB GWAS results on `/slow/projects/ifpan-abm-pgxpred/gwas-results`).
  300/300 atlas rows populated (burden tests included).
- Added 6 raw GWAS features + 6 BETA/LOG10P difference contrasts:
  `delta_BETA_ADR_vs_pure/dose/prescribed`, `delta_BETA_pure_vs_prescribed`,
  `delta_LOG10P_*`.
- Gained +10pp T1 recall — the new contrasts expose indication-confounding
  (gene predicts who *gets* the drug) as its own axis.

### v8 — prune dead weight + reweight (36 features)
- **Dropped 8 features** with between-tier-variance < 0.04: `IMPACT_ord`
  (skew 10.8, ~constant), `evidence_product`, `effect_magnitude_AND`,
  `BETA_ICD10_abs`, `BETA_prescribed_abs`, `LOG10P_prescribed`, and **two of the
  newly-added** `delta_LOG10P_*` contrasts (they didn't separate anything).
- **Up-weighted 10 T3-distinctive features** (×1.5–2.0): `beta_per_pleio_log`,
  the SE/BETA family, the `delta_BETA_*` contrasts.
- GMM k=4 → best result: 59.3%, T4=96%, κ=0.42.

### v9 — Option B tissue methodology (42 features) [see §3]
- Replaced the single categorical `Tissue_call` with continuous, peer-reviewed
  tissue metrics + MAGMA + cross-source replication.

### v10 — drop MAGMA + tau, add absolute-TPM (40 features)
- **Dropped** `Tissue_MAGMA_neglog10FDR` (κ=0.00 vs expert — wrong tool, see §3)
  and `Tissue_tau` / `Tissue_tau_HPA` (ρ=−0.26 — measures specificity to ANY
  tissue, not disease tissue).
- **Added** `Tissue_abs_logTPM_disease` (absolute expression in disease tissue —
  the metric that actually tracks the old expert call, ρ=+0.23).

### v11 — drop LOG_OR_SE
- **Dropped** `LOG_OR_SE` (was the #1 cluster separator, η²=0.91).
- Evidence for the decision: 100% of its extreme rows (>4) are burden tests;
  skew 7.2 (max 26 vs median 0.8); r=0.74 with `SE_regenie_log`. It is a
  **Firth-Wald estimation-instability axis**, not effect biology, and its η² was
  inflated by ~8 outliers forming a 3-row micro-cluster.
- Verification: drop vs keep vs winsorize gave **identical** k=4 TIER recovery
  (46.7%, same per-tier); dropping only eliminated a 3-row statistical-artifact
  cluster at k=7 and rebalanced the rest (ARI 0.78 vs baseline). No conclusions
  changed.

### v12 — winsorize the SE/skewed family (39 features) [CURRENT]
- **Winsorized at the 95th percentile**: `SE_regenie_log`, `SE_ICD10_log`,
  `SE_dose_log`, `SE_prescribed_log`, `OR_log` — all skew>2, same heavy-tail
  pathology as the dropped LOG_OR_SE.
- Rationale: caps the rare-variant/Firth-Wald tails so no single SE feature
  dominates the distance metric, while keeping the legitimate rare-vs-common
  axis. Manuscript phrasing: *"Standard errors were winsorized at the 95th
  percentile to prevent Firth-Wald quasi-separated burden tests from dominating
  the distance metric."*
- Effect: k=4 anchors hold (E 24/26, I 62); the 3-row artifact cluster does not
  reappear. Caveat: GMM stability dropped (138→77 stable rows) — winsorizing
  reshaped the Gaussian tails; Ward/k-means unaffected.
- `SE_prescribed_log` remains the top separator (η²=0.92) but now as legitimate
  rare/common-variant structure, not outlier domination.

### v13 — drop evidence_sum (38 features)
- **Dropped** `evidence_sum` (= z_FAERS + z_tissue + z_repl). It was a pure
  additive redundancy — its 3 constituents plus min/max/spread already carried
  the combination.
- Effect: **improved** — GMM κ 0.256→0.314 (best recorded), GMM stability 77→88,
  Ward k=4 46.7→51.3%. Confirms it was redundant.

### v14 — drop evidence_min/max/spread — REVERTED
- Dropped the trio. **Degraded** all methods: LCA ARI 0.154→0.113, GMM k=4
  56.7→52.3%, Ward k=4 51.3→46.3%. The trio was NOT redundant — `evidence_min`
  is the continuous AND-conjunction (the only feature recovering T4). Reverted.

### v15 — replace ad-hoc evidence features with named meta-analytic statistics (36 features) [CURRENT]
- Evaluation (`evidence_eval.py`) showed the original features ARE established
  statistics: `evidence_min` ≡ Nichols minimum-statistic for conjunction
  inference (r=1.00; Nichols et al., NeuroImage 2005); `evidence_sum` ≡ Stouffer
  combined Z (r=1.00; Stouffer 1949). And only the minimum statistic carries the
  T4-conjunction signal (Cohen d=+0.72 vs +0.30 Stouffer, +0.33 Mahalanobis,
  +0.09 Fisher).
- **Replaced** the 4 ad-hoc evidence features (sum/min/max/spread) + `max_x_spread`
  with **two named, citable statistics**:
  - `evidence_stouffer_Z = Σz/√k` — overall combined evidence (Stouffer 1949);
    the T1/abundance axis (d=+1.26 for T1).
  - `evidence_conjunction_minZ = min(z)` — valid conjunction inference (Nichols
    2005); the T4 axis (d=+0.72), i.e. "all three evidence layers jointly elevated".
- Effect: recovered the v14 losses AND improved — GMM ARI 0.205 (best ever),
  **GMM stability doubled 88→168**, LCA ARI back to 0.150, GMM k=4 back to 56.7%,
  with 3 fewer features than v13. GMM κ 0.252 (give: GMM chose k=7 not k=10).
- Manuscript phrasing: *"The three orthogonal evidence layers (pharmacovigilance,
  tissue specificity, statistical replication) were summarized by two
  meta-analytic statistics: Stouffer's combined Z (overall evidence) and the
  minimum-statistic for conjunction inference (Nichols et al., 2005; whether all
  layers are jointly elevated)."*

### v16 — symmetric GWAS core (35 features) [CURRENT]
- Regularized the GWAS block for order and interpretability. **Per phenotype:
  BETA + LOG10P; SE retained only for the ADR discovery arm.**
  | phenotype | features |
  |---|---|
  | drug-ADR (REGENIE discovery) | BETA_regenie_abs, LOG10P_disc, **SE_regenie_log** |
  | gen-ICD10 (pure) | BETA_ICD10_abs, LOG10P_ICD10 |
  | drug median-dose | BETA_dose_abs, LOG10P_dose |
  | drug-prescribed | BETA_prescribed_abs, LOG10P_prescribed |
- **Dropped** `OR_log` (redundant with BETA_regenie) and the 3 non-ADR SEs
  (`SE_ICD10_log`, `SE_dose_log`, `SE_prescribed_log`) — they were a single
  redundant rare-variant/Firth-Wald axis (r=0.74 with SE_regenie) that had
  dominated cluster geometry.
- **Added** the previously-missing `BETA_ICD10_abs`, `BETA_prescribed_abs`,
  `LOG10P_prescribed` to complete the symmetry.
- Kept `A1_FREQ_log` + `A1_CASE_CT_log1p` as a small "variant properties"
  (rarity/power) group, and `P_test_neglog` as the replication-evidence layer.
- Effect: essentially equivalent agreement (GMM k=4 56.7%, anchors hold) with
  **further-improved GMM stability 168→205 stable rows** and a clean, symmetric
  feature table. Slight GMM ARI dip (0.205→0.187) is the cost of removing the
  SE rare-variant axis — acceptable, since de-emphasizing that nuisance axis was
  itself a goal (see v11/v12).
- Manuscript phrasing: *"For each of four GWAS phenotypes (drug-ADR, disease,
  drug dose, drug prescription) we extracted the absolute effect size and
  −log10 p-value; the standard error was retained only for the ADR discovery
  arm, the other arms' SEs being a redundant rare-variant axis."*

### v17 — directional temporal ratio (35 features) [CURRENT]
- **Replaced** `ADR_pct_before` with a single directional feature:
  `ADR_after_over_before = log((TEMP_pct_after + 1) / (TEMP_pct_before + 1))`.
  Kept `ADR_pct_after`. Pseudocount of 1 percentage-point stabilizes the 40 rows
  with `before == 0` (raw division would be infinite).
- Interpretation: high = events follow prescription (drug-induced / ADR
  direction); negative = disease pre-existed (indication direction).
- Cleaner biology in one number: per-TIER standardized means rise monotonically
  T1=−0.54 → T2=+0.12 → T3=+0.35, with T1 (Indication) sharply low
  (d=−0.81 vs rest) — the disease-preceded-drug signature.
- Effect: small improvement — LCA ARI 0.132→0.153, Ward k=4 49→51%, GMM κ
  0.261→0.265; anchors hold. GMM stability 205→172 (minor).

### v18 — drop A1_FREQ_log (34 features) [CURRENT]
- **Dropped** `A1_FREQ_log` (population MAF). Correlated r=0.71 with
  `A1_CASE_CT_log1p` and the less load-bearing of the two (k=7 ARI 0.891 when
  dropped vs 0.807 for A1_CASE_CT). `A1_CASE_CT_log1p` subsumes rarity (carrier
  count is driven by MAF) AND adds the zero-carrier statistical-power signal
  (70 rows with 0 carriers in cases) that MAF alone cannot give.
- Effect: **clean win** — GMM ARI 0.194→**0.218** (best across all 18 versions),
  **GMM stability 172→276** (most stable GMM yet), GMM converged to a clean k=5
  (was k=9 with 2 junk singletons). LCA and k=4 anchors unchanged.
- The retained `A1_CASE_CT_log1p` is best read as a statistical-power /
  observability feature (variant × cohort), not a pure variant property.

### v19 — beta_disc_valid_consistency — REJECTED
- Added effect-direction agreement between REGENIE BETA and PLINK2 log(OR).
  Real biological signal (T3 has worst discovery↔validation agreement, d≈−0.4 —
  winner's-curse on rare variants), BUT destabilized GMM (stable 276→166, k 5→9)
  by re-injecting the rare-variant Firth-Wald outlier axis (cf. v11/v12).
- Winsorizing fixed the stability (→300) but still lowered ARI (0.218→0.177): the
  T3-replication axis is genuine but orthogonal to the EPIC tiers, so it
  subdivides the partition. Conclusion: real signal, wrong place — **kept as a
  REPORTED column** in atlas_enriched.parquet, NOT a clustering feature. Reverted
  clustering to v18.

### v20 — remove ALL manual feature weights (34 features) [CURRENT]
- **Removed** the ×1.5/×2.0 multipliers on 7 features (beta_per_pleio_log×2.0,
  SE_regenie/BETA_regenie/BETA_dose/delta_BETA_dose/delta_BETA_prescribed/
  Tissue_zscore_GTEx ×1.5). These were hand-tuned toward T3/T4 recovery with no
  principled derivation — the clearest arbitrariness in the pipeline (and the
  dict still held stale entries for features dropped in v16).
- **Verification revealed the weights were propping up GMM at k-means' expense:**
  - k-means: degenerate k=2 ARI 0.005 → clean **k=3 ARI 0.187, κ 0.248**
  - GMM: ARI 0.218/stable 276/k=5 → 0.170/73/k=9 (the "best" GMM was tuned)
  - k=4 forced: now **consistent ~55% across kmeans/gmm/lca** (all T2≈93%),
    vs scattered 47–57% before.
- Decision: **accept unweighted.** No method artificially favored; cross-method
  agreement is tighter and more honest than one tuned GMM number. Report
  k-means/LCA as primary (stable 300/291); note GMM is weight-sensitive.
  Structural findings (T3=16%, T4=0% at k=4) unchanged. Matches the
  weight-free description already in methods_section.md.

### v21 — effect direction (sign concordance) [Step 1 of validation-driven improvements]
- Motivated by the literature validation: statin×varicose was called Effect though
  statins are protective — the model used |BETA| everywhere and discarded sign.
- **Added** `effect_sign_concordance = sign(BETA_regenie) · sign(BETA_ICD10)`
  ∈ {−1,0,+1}: +1 = drug-amplified same-direction (vertical pleiotropy / real
  mechanism), −1 = drug-arm opposes general-pop effect (protective/confounded).
  46/300 rows discordant, incl. the statin×varicose Effect rows. Bounded (avoids
  the v19 outlier-instability trap).
- Effect: **clear win** — GMM k=9→5 (clean, no singletons), ARI 0.170→0.218,
  **stable 73→277**. Other methods unchanged; k=4 anchors unchanged. κ 0.294→0.257
  (fair trade: v20's higher κ was on an unstable k=9 partition). ACCEPTED.

### Step 2 — per-locus multiplicity: tested as feature, REVERTED to reported column
- Parsed Variant_ID (chr-pos-ref-alt) → 250kb-window multiplicity per (drug,ICD10).
  Correctly flags the LD clusters (RYR2×8, simvastatin-locus×7; 65/300 rows with
  multiplicity ≥5).
- As a *clustering feature* (v22): marginal/mixed — kmeans ARI 0.185→0.164 (the
  cleanest method degraded, k 3→4), GMM unchanged, HDBSCAN κ 0.150→0.197, LCA
  −0.004. Conceptual problem: adding multiplicity as a feature lets the clusterer
  *group* LD rows (reinforce), opposite of the intended *deflation* — true
  deflation needs row-collapse (kept 300 rows per user decision).
- Decision: **reverted as a feature; kept `locus_multiplicity` as a reported
  column** in atlas_enriched.parquet (same discipline as v19 beta-consistency).
  Accepted model stays v21.

### v23 — phenotype-specificity flag [Step 3] [CURRENT ACCEPTED]
- **Added** `phenotype_nonspecific = ICD10 starts with R or Z` (symptom/health-status
  codes, low interpretability). 67/300 flagged — notably 11/26 (42%) of Effect-tier
  rows are on vague R/Z codes.
- Effect: **clean small win** — GMM stability 277→298 (best), LCA ARI 0.133→0.139
  + κ up, HDBSCAN κ 0.150→0.182, kmeans unchanged. All methods improved or held.
  ACCEPTED. Current accepted model = v23 (36 features: v21 + this flag).

### Step 4 — Open Targets gene-plausibility (reported column, NOT a feature)
- `gene_plausibility.py`: per gene → OT Ensembl + associated diseases; per
  (gene, ICD10) → max association score for the disease (name-match), else 0.
  Added `OT_gene_disease_score` to atlas_enriched.parquet + coords.csv. Cached.
- Validation: GLP1R×psoriasis=0.041 (✓ the literature case), DUSP7×varicose=0
  (✓ incidental). Honestly misses niche mechanistic links not in OT (DLEU1, PDE1C).
  67/300 rows have any OT link.
- **Striking finding:** mean OT plausibility by tier is T1=0.041, T2=0.013,
  T3=0.021, **T4(Effect)=0.002** — the Effect tier has the LEAST gene plausibility,
  confirming its gene-level calls are incidental. The genuinely gene-linked pairs
  (Factor-V-Leiden×thrombophlebitis 0.55, APOE×Alzheimer 0.68, LPA×CAD 0.55,
  HLA-C×psoriasis 0.40) land in Indication/Comorbidity, where the gene drives the
  disease independent of the drug. Reported column does not change clustering.

### Step 5 — expanded gold-standard benchmark + per-class precision/recall
- `data/benchmark/gold_standard.tsv` (35 curated pairs: 13 Indication, 9 Effect,
  12 inConclusive — the 3 classes with defensible pharmacological ground truth;
  Predisposition excluded as not objectively determinable). `benchmark_eval.py`
  computes per-class P/R/F1 for expert + 4 methods.
- **Result (34/35 matched):**
  - **Expert acc 0.71**: Indication F1 **0.96** (R 1.00 — perfect), inConclusive
    F1 0.82, but Effect F1 **0.33** (R 0.22 — misses 7/9 real ADRs, calling them
    Predisposition: ACEi angioedema, hyponatraemia, gabapentin dizziness,
    antiplatelet GI bleed, opioid urinary retention, hypotension).
  - **Unsupervised acc 0.15–0.26**: collapse almost everything to Predisposition;
    Effect recall 0.00, inConclusive recall 0.00 on this benchmark.
- Conclusion (quantifies the literature validation): the expert cascade is strong
  at Indication and decent at inConclusive but weak at Effect; the unsupervised
  methods cannot stand alone as an I/E/C clinical classifier (Predisposition-axis
  dominance). Caveats: small benchmark, no P ground truth, unsupervised EPIC
  labels depend on the cluster-labeling step. Published to /slow/www/atlas/.

### Latent / drug-triggered PGx prioritisation (`latent_pgx_score.py`)
- Reframed the goal: the only PGx-interesting class is a *hidden predisposition
  realised only by the drug* (gene×drug interaction) — distinct from manifest
  predisposition (gene→disease without drug), confounding-by-indication,
  pleiotropy, non-replication, and protective/reverse.
- Operational score (all from existing features, no new data). The drug-ADR
  signal must **dominate every alternative phenotype** in both effect size and
  significance — weakest-link (min) over {disease-alone, dose, prescribed}:
  `beta_superiority = min_alt(log|BETA_ADR| − log|BETA_alt|)`,
  `p_superiority = min_alt(LOG10P_ADR − LOG10P_alt)`;
  `latent = z(beta_superiority) + z(p_superiority) + z(P_test) − z(pleiotropy)`.
  Beating dose/prescribed (not just ICD10) rejects titration- and
  indication-confounding that mimics a drug effect. Gated by concordant sign +
  not-on-BNF + replicated (PLINK2 P<0.05); 85/300 pass. Priority = latent + 3·OT.
  The stricter score surfaced diazepam×COPD×CXCL13 (a validated benzodiazepine
  respiratory ADR) into the top candidates.
- **The score demonstrably works** as a discriminator: manifest predispositions
  with real biology (F5/thrombophlebitis OT 0.55 latent −6.3; APOE/Alzheimer
  OT 0.68 latent −9.4; HLA-C/psoriasis OT 0.40 latent −7.8) are driven deeply
  negative and flagged `manifest_risk`; genuine latent leads (PHEX×clopidogrel
  latent +1.41 OT 0.07; TRERNA1×ramipril +0.99 OT 0.05) sit positive.
- **Sobering finding:** of 85 statistical candidates only 7 have plausible biology
  (OT≥0.05 + named gene), and the high-latent ones are mostly incidental genes —
  the clean intersection (latent signature ∧ real biology) is nearly empty. True
  G×E is the rarest class; this atlas does not contain many confident examples at
  the variant level. Deliverable is therefore a *prioritised shortlist for
  fine-mapping + cross-cohort replication*, not a final class.
- Outputs: `latent_pgx_score` (continuous) + `latent_pgx_candidate` reported in
  atlas_enriched.parquet; ranked `latent_pgx_shortlist.tsv` with flags
  (biology_plausible, fine_mapping_needed, strong_replication, manifest_risk);
  added as a 6th continuous-coloured column (plasma; triangles=candidates) to the
  PCA/UMAP projection plots. Published to /slow/www/atlas/.

### v28 — ATLAS_full descriptors ADDED, dead parameters REMOVED (add+remove re-evaluation) [CURRENT]
- Discovered the `ATLAS_full` sheet of ATLAS_v2.xlsx (300 rows, 164 cols, same variant_id key)
  which NO script had ever read. New module `merge_atlas_full.py` additively joins 20 columns on
  the (variant_id, Lek, ICD10) triple (207/207 matched); `enrich_with_gwas.py` also pulls them on
  a full rebuild. Source workbook untouched (read-only).
- Constraint (user): only VARIANT descriptors and PHENOTYPE descriptors, NO gene-specific
  parameters (Druggability / DGIdb / GO / gene-PheWAS / Tissue excluded; OT stays reported-only),
  and NO arbitrary thresholds. VEP functional explicitly authorised as variant-intrinsic.
- **ADDED 16 features** (v27 27 -> v28 41): temporal block (log n_before/after, median years
  before/after, onset IQR, missing-date rate), held-out test-set genotype (case-ctrl freq diff,
  case het/hom counts, Z_STAT), OR calibration (log-diff and ratio, test vs discovery — requested
  explicitly), and VEP (consequence severity, protein-coding biotype, gnomAD AF + missingness).
- **REMOVED 2 features** on threshold-purity grounds (§1): `FAERS_ord` (ordinal of a thresholded
  call) and `phenotype_nonspecific` (the "ICD10 starts with R/Z" rule; univariate AUC 0.43, below
  chance, and near-constant).
- **`temporal_log_n_before` is the find**: events occurring BEFORE drug exposure — the comorbidity
  signature. Univariate CV AUC 0.81 (Effect vs Indication), the best of any feature. NOT a
  phenotype-size artifact (size alone 0.52; corr with n_total only 0.26) and it beats the existing
  `pct_before` (0.76).
- **SIFT / PolyPhen NOT added**: in this atlas they take only two distinct values each (a single
  scored coding variant) — near-constant, the same reason IMPACT was dropped in v8.
- **Redundancy pruning TESTED and REJECTED.** Dropping the 6 collinear features
  {FAERS_PRR_log, SE_regenie_log, delta_BETA_ADR_vs_prescribed, BETA_ICD10_abs,
  delta_BETA_pure_vs_prescribed, evidence_stouffer_Z} (all |rho| >= 0.85) degraded EVERY method
  on the full pipeline (LCA ari_vs_TIER 0.077 -> -0.007, stable 194 -> 187; GMM stable 44 -> 20;
  kmeans kappa_EPIC 0.071 -> 0.055). The cheap k-means stability proxy in `feature_ablation.py`
  did NOT predict the real LCA/GMM behaviour. **The v8 "noise-reducing ensemble" finding stands.**
- Also confirmed the plan's guard: `repl_directional_z` has univariate AUC 0.47 yet dropping it
  costs -0.203 consensus stability. **Univariate AUC is never a removal criterion.**
- Full-pipeline effect (27 -> 41 features): kmeans stable **137 -> 203**, kappa_EPIC 0.052 -> 0.071;
  GMM kappa_EPIC 0.052 -> **0.123** (best); LCA k=5 stable 199 -> 197, kappa 0.052 -> 0.047.
  Cost: **HDBSCAN degenerates to k=1** at 41 dims (curse of dimensionality; no duplicate rows,
  mean pairwise distance 6.38 -> 8.36). HDBSCAN is a robustness check, not a primary method —
  reported, not hidden. `ari_vs_TIER` falls for LCA/GMM, but agreement with the expert cascade is
  a diagnostic, not a target (guiding principle), and the cascade is known to mislabel (varicose).
- **Separation (the goal).** New `good_bad_separation.py` -> `separation_metrics.json`:
  | projection | good-vs-bad | Effect-vs-Indication (CONTROL) |
  | Set A model, PCA  | 0.66 -> **0.75** | 0.44 -> **0.71** |
  | Set A model, UMAP | 0.84 -> **0.89** | 0.54 |
  | Set B diagnostic, PCA  | **0.86** | 0.70 |
  | Set B diagnostic, UMAP | **1.19 (separated)** | 0.64 |
  The independent E-vs-I control improved 0.44 -> 0.71, so the good/bad gain **generalises** — it
  is NOT the overfitting we feared during planning.
- **Two feature sets** as agreed: Set A (model) drives clustering/EPIC/latent; Set B
  (`diagnostic_projection.py` -> `coords_diagnostic.csv`) adds the separation-driving descriptors
  and is a FIGURE ONLY — labelled as such on the new Set-A-vs-Set-B page of the projections PDF.
- Published: ATLAS_cluster_projections_v28.pdf (5 pages), atlas_v28_classified.tsv,
  feature_ablation.tsv, separation_metrics.json, coords_diagnostic.csv, cluster_coords_v28.csv.

### v27.12 — GOOD / BAD associations marked on the plots [CURRENT]
- New pipeline step `mark_good_bad.py` (run after export_plots_data) tags coords.csv with
  `annot` = good (6 rows: the target ADR set — diazepam x COPD x2, citalopram x osteoporosis,
  lisinopril x hyperglycaemia, ramipril x confusion, gabapentin x dizziness) and bad (4 rows:
  expert-flagged confounders — simvastatin x varicose x3, fluoxetine x menopause) + labels.
- Rendered as: GREEN CIRCLES = good, RED CROSSES = bad, over the latent-PGx-coloured landscape.
  Added as a 3rd page of ATLAS_cluster_projections_v27.pdf (PCA + UMAP, labelled) and overlaid
  on the deck's Prioritisation panel.
- Visual finding: the BAD simvastatin x varicose crosses sit INSIDE the cloud of GOOD circles in
  PCA — a direct visual confirmation that comorbidity is not separable in feature space (only the
  expert biology layer demotes it), whereas fluoxetine x menopause separates out on its own.

### v27.11 — FAERS x temporal pair axis (real-world causality) [CURRENT]
- Pairwise-product search over all 2-parameter multiplicates: `FAERS x temporal` is the ONLY
  pair that raises Effect-vs-Indication AUC (0.900 -> 0.909). It encodes real-world CAUSALITY —
  an event both disproportionately reported AND following the prescription. Added as a 6th
  conjunction axis, sqrt(sig(FAERS)*sig(temporal)), weight 1.
- Rejected: `effect x freq` (too universal, moves nothing) and `repl x FAERS` (ranks cleanly but
  no AUC gain and is redundant re-weighting of existing axes). Key realisation: because the score
  is already a conjunction (product) of its axes, adding a pair of EXISTING axes is redundant
  re-weighting, not new information; only pairs of non-axis parameters could add signal, and
  those are too universal to discriminate.
- Z-score normalisation of every parameter confirmed already in effect (sig(z(x))) — an explicit
  z-scored rebuild reproduced the identical 0.90 base.
- Effect on the set: diazepam x COPD #11, citalopram #68, ramipril #64 (lifted from 105 —
  strongly reported ADR that follows the Rx, still flagged MED repl-underpowered), lisinopril
  #127, gabapentin #163; expert-flagged confounders pinned at #198 / #205.

### v27.10 — replication down-weighted (0.7) + aggregate search [CURRENT]
- Replication is a DIRECTIONAL, held-out-underpowered signal (repl_z = beta-based), so it now
  carries reduced weight (0.7) in the weighted soft-AND. Tested full/0.7/power-aware/zero:
  0.7 keeps AUC 0.90 and zero 0-carrier artifacts in the top while lifting the underpowered
  ramipril 131->105 (still MED). Power-aware (neutralise-when-untestable) was rejected — it let
  the 0-carrier ramipril rise to #9 and readmitted artifacts. Zero weight promotes artifacts.
- Aggregate-to-add search: cross-drug consistency (Bradford-Hill; # drugs sharing variant×ADR)
  is near-constant here (2/207 shared -> AUC 0.48, adds nothing) because atlas variants are
  drug-specific rares; significance-specificity/discovery/impact also don't improve. The 5-axis
  conjunction is at the data's ceiling; a class-consistency axis would need external drug-class
  data and would matter at genome-wide scale, not on this curated set.

### v27.9 — RECONCILED into ONE calculation (confidence tier replaces the two tracks) [CURRENT]
- Merged Track 1 (replicated) and Track 2 (pharmacovigilance) into a SINGLE score. Replication
  is no longer a hard gate / track-separator: it stays a SOFT axis inside the conjunction
  (preserving discrimination) and is reported as a CONFIDENCE tier. Tested adding a replication-
  independent discovery+impact axis to fold Track 2 in fully, but it dropped AUC 0.90 -> 0.74/
  0.79 (discovery/impact reward winner's-curse + common indication variants), so the single
  score keeps Track 1's five axes.
- One list, one gate (concordant + not-indication + FAERS-present + not-expert-flagged): 129
  candidates = 74 HIGH (replicated) + 52 MED (repl-underpowered, <=2 test carriers) + 3 LOW
  (unreplicated). The former "Track 2" is now simply the MED/LOW rows. Columns: latent_pgx_score,
  latent_confidence, repl_z, test_carriers. Top of list = HIGH-replicated leads (angioedema,
  diazepam x COPD); ramipril #94 / gabapentin #116 sit in-list flagged MED.
- AUC over the full benchmark reads 0.80, but this is a BENCHMARK ERROR: the 3 flagged rows
  (simvastatin x varicose ABHD14B/DOCK3, mislabelled TIER-4 Effect by the old cascade) are
  down-weighted per michkor (comorbidity). Excluding them the AUC is **0.90** unchanged — the
  expert correction is more right than the benchmark. Separate track2_pharmacovigilance_shortlist
  retired; single latent_pgx_shortlist.tsv now carries the confidence tier.

### v27.8 — expert comorbidity/indication down-weight (curated layer) [CURRENT]
- Reviewer (michkor) flagged two that should be LOW: simvastatin x varicose (I83) = COMORBIDITY
  (old patient, years on statins for CV disease, notices age-related varicose, reports it;
  large studies show statins PROTECTIVE) and fluoxetine x menopause (N95) = INDICATION (SSRIs
  used off-label). fluoxetine x menopause was already low (#143/158, killed by temporal axis).
- simvastatin x varicose sat too high (#29). Tested THREE objective penalties — significance-
  specificity (p_superiority), disease/prescribed confounding, and gene-disease OT — NONE
  lowered it without hurting real ADRs or AUC. Confirms the reviewer point: comorbidity is
  statistically INDISTINGUISHABLE from a real drug-caused ADR (it replicates, has FAERS,
  moderate everything); only external biology (statins protective) separates it.
- Added an explicit, transparent EXPERT_DOWNWEIGHT layer (curated (drug, ICD-prefix) ->
  0.35x penalty + reason), applied to both tracks and removing the row from candidates. Result:
  simvastatin x varicose -> #198-200, fluoxetine x menopause -> #201; real ADRs unchanged
  (diazepam x COPD #6/#9, clopidogrel x PHEX #25). Reason reported in expert_downweight_reason.
- What can improve (systematic version of the expert layer): a literature-derived
  KNOWN-DIRECTION annotation (drug protective/harmful for the outcome) + a drug-indication-
  comorbidity map — the objective encoding of "statins are protective for varicose", built like
  FAERS/BNF rather than hand-listed.

### v27.7 — Track-2 actionability multiplicate (effect x freq) [CURRENT]
- A multiplicate search over parameter products found `effect_size x allele_frequency` lifts
  the clinically-real-but-unreplicated trio while KEEPING discrimination (diazepam x COPD 0.93
  vs junk/manifest/off-label 0.13 percentile). Added it as a 6th Track-2 conjunction axis =
  population-level ADR impact / actionability (common, large-effect variants — the opposite
  pole from rare-variant artifacts). Effect: gabapentin 93->73, lisinopril 102->94, diazepam
  stays #18, junk absent (no FAERS), off-label #117. Replication carries ZERO weight in the
  Track-2 score (flag only) — nothing to reduce there. A 2nd axis (precision) was rejected: it
  would hurt lisinopril (16th precision percentile).

### v27.6 — TWO-TRACK resource: pharmacovigilance-anchored Track 2 [CURRENT]
- Added a second, parallel candidate list for clinically-real ADRs whose genetics do NOT
  confirm in the 20% held-out set (the lisinopril/ramipril/gabapentin class). Scored on the
  axes these associations ARE strong on: discovery -log10P, FAERS pharmacovigilance, temporal
  (events follow Rx), drug-specificity, variant-specificity — a soft-AND conjunction, NO
  replication requirement. Gated by concordant + not-indication + FAERS-present.
- Replication is REPORTED as a flag, not required. Each entry labelled:
  `1+2 (replicated genetic lead)` | `2 (pharmacovigilance; repl underpowered)` (test-set
  carriers <=2, replication untestable) | `2 (pharmacovigilance; unreplicated)` (adequately
  powered but failed). Output: `track2_pharmacovigilance_shortlist.tsv`.
- 132/207 pharmacovigilance candidates: 77 replicated (Track-1 overlap), 52 repl-underpowered,
  only 3 genuinely unreplicated. **ramipril x R41.0 is Track-2 #1** (discovery 97th + FAERS
  97th + temporal 99% but 0 test carriers -> "repl underpowered"); lisinopril #102, gabapentin
  #93 (weaker on temporal/specificity even here). The honest home for the trio: strong
  real-world+discovery evidence, genetic replication untestable/absent — flagged, not hidden.
- Track 1 (replicated genetic latent, AUC 0.90) remains the primary PGx-target list.

### v27.5 — internal-pleiotropy de-contamination [superseded label only; both tracks CURRENT]
- Diagnosing why the mid-rank real ADRs sat low, the per-axis breakdown showed citalopram x
  osteoporosis was dragged ONLY by its pleiotropy axis. Inspection: its variant's 3 "hits"
  were the citalopram ADR + TWO hits in `p403_..._snap_button_count` — a UK Biobank
  cognitive/dexterity test, not a drug or disease. `internal_pleiotropy` was counting
  non-pharmacological ancillary UKB trait GWAS.
- Fix (existing data, no new associations): recompute `internal_pleiotropy` over
  pharmacologically-relevant phenotypes only (drug-ADR / ICD / dose / prescribed), excluding
  the p-coded UKB field traits. **1,383 of 4,452 phenotypes (31%) were ancillary**; 63/207
  atlas associations had inflated pleiotropy; mean pleiotropy 4.4 -> 1.6.
- Effect: **AUC 0.87 -> 0.90** (best). citalopram x osteoporosis 98 -> 55 (the one target that
  was genuinely pleiotropy-limited). diazepam x COPD still top (#6, #9). lisinopril/ramipril/
  gabapentin unchanged (limited by replication/specificity, not pleiotropy — real).
- Cost: k-means consensus stability dropped (204 -> 137) as the feature has less spread; LCA
  (primary) stayed stable (199/207, k=5). Acceptable.

### v27.4 — parameter-handling upgrades (empty-arm + precision) [superseded by v27.5]
- Three parameter-handling fixes were implemented and ABLATED (existing data only):
  - **#1 empty-arm handling [KEPT, AUC 0.84->0.87]**: the cross-phenotype superiority no
    longer does `fillna(0)` on the comparison arms (untested != null). Arms the variant was
    never tested in (BETA NaN) are excluded from the min; each tested arm's denominator is
    floored at its own detection limit (its SE), so an unpowered/near-null alternative can
    no longer masquerade as drug-specific.
  - **#2 precision (SE-shrinkage) [KEPT, neutral]**: effect magnitude = max(0, |BETA|-SE)*gsd,
    de-inflating imprecise low-carrier discovery betas. Neutral-to-slightly-positive with #1.
  - **#3 robust median/MAD scaling [REJECTED]**: dropped AUC 0.87->0.70 — MAD collapsed the
    skewed axes and saturated the conjunction sigmoid. Reverted to mean/SD z.
- Net: Effect-vs-Indication AUC 0.75 (additive) -> 0.84 (conjunction) -> **0.87**. Real ADRs
  at top (diazepam x COPD #7, lisinopril x angioedema #16, clopidogrel x PHEX #21), all
  confounders demoted (JUNK #166, off-label #138-147, indication #109, manifest #201),
  artifacts low (ramipril #129, gabapentin #168). The two signal-less targets stay low —
  now confirmed a FIFTH way.

### v27.3 — conjunction (soft-AND) latent score instead of additive sum [superseded by v27.4]
- Diagnosis (from laying out all raw parameters for 6 target + 10 contrast associations):
  each confounder class is missing ONE evidence axis, and an additive sum lets the others
  compensate. JUNK status codes have no FAERS but sum high (Z98.0 ranked #15); off-label
  (fluoxetine x menopause) had huge drug-specificity from an empty disease-arm and ranked #2.
- Fix: the latent score is now a CONJUNCTION — geometric mean of per-axis sigmoids over
  {directional replication, FAERS pharmacovigilance, temporal-after, drug-specificity
  (beta_superiority), low internal pleiotropy}. A soft-AND (Nichols 2005 conjunction logic,
  already used for evidence_conjunction_minZ) requires ALL axes, so any single missing one
  disqualifies the association. Two axes were ADDED that the additive score lacked: FAERS
  and temporal-after (the latter was the top ADR marker in the supervised fit).
- Effect-vs-Indication AUC 0.75 (additive) -> 0.84 (conjunction). Confounders demoted
  (Z98.0 15->187, menopause 2->172, manifest predispositions ~200), genuine ADRs lifted
  (diazepam x COPD -> #6, lisinopril x angioedema #2). Priority = z(latent) + 3*OT (FAERS
  now lives inside the latent, no longer added separately).
- Tested and REJECTED lowering/removing the replication axis (user request): AUC falls
  0.86->0.78 and 0-carrier artifacts flood the top (2->10 of top-30); the only target that
  rises (ramipril, MAC=0) rises BECAUSE it is an artifact, while genuine ADRs sink. Kept
  replication at full weight.
- Confirmed across four independent analyses (percentiles, supervised logistic AUC=0.69,
  conjunction, replication-ablation) that gabapentin x R42 and ramipril x R41.0 lack
  variant-level signal and cannot be honestly lifted — correct behaviour, documented.

### v27.2 — directional (effect-based) replication instead of the p-value [superseded by v27.3]
- Confirmed the GWAS design: **80/20 held-out split** of individuals (REGENIE on 80%
  train, PLINK2 on 20% test) — verified from the data (per-phenotype train% = 79-80%,
  ~4:1 N ratio across 8 phenotypes) and the README ("held-out test set"). So
  `P_test_PLINK2` IS a genuine independent replication.
- Problem: for RARE ADR variants the two-sided test-set p-value is underpowered — the
  variant has few/zero carriers in the smaller 20% set, so P stays non-significant even
  for a real effect. The p-gate lumped true 0-carrier artifacts together with real-but-
  underpowered effects. Empirically 201/207 associations are sign-concordant between
  discovery and test-set betas, but only 105 reach two-sided P<0.05.
- Fix: replicate on the **EFFECT, not the p-value**. `repl_z = sign(BETA_discovery) *
  log(OR_test) / SE_test` — the one-sided test-set z in the pre-specified discovery
  direction, precision-weighted. A real effect keeps its SIGN in the test set even when
  underpowered; a winner's-curse false positive regresses to ~0 or flips. Replaces the
  p_test term in the latent score and the gate becomes `repl_z > 1.645` (one-sided p<0.05).
- Result: candidates 76 -> 101, biology-plausible 6 -> 10 (rare directional replicators
  recovered). The 0-carrier artifact (ramipril x R41.0, repl_z ~0) and the genuinely-weak
  test-set effects (gabapentin x R42, repl_z 0.55) correctly STAY rejected — the method
  did not manufacture support. lisinopril x R73.9 sits on the boundary (repl_z 1.59).
- Note: `A1_CASE_CT` is the TEST-set carrier count; the rare-variant non-replications are
  test-set-underpowered, not necessarily false.
- Also carried into the CLUSTERING feature set: the model feature `P_test_neglog` (two-sided
  p) was replaced by `repl_directional_z` (same directional statistic), and the combined-
  evidence layers now use it. Clustering is unchanged in quality (LCA ARI 0.084->0.090,
  k=5 stable 195; kmeans stable 204) — a clean swap. 27 features retained.

### v27.1 — latent-PGx priority redesign: require real-world ADR evidence [superseded by v27.2]
- Biological verification of the v27 top candidates (literature-checked) showed the
  ranking was dominated by NOISE: administrative status codes (simvastatin x Z98.0 was
  #1), indication/comorbidity confounding, wrong-direction pairs (atorvastatin x
  hypoglycaemia — statins cause HYPER-, not hypo-glycaemia), and incidental genes
  (OT=0). Root cause: the latent score is purely genetic specificity + replication,
  which a spurious hit in a junk phenotype achieves trivially; it ignored the two
  signals that mark a real ADR — real-world pharmacovigilance and biology.
- Fix (NOT a phenotype filter, per user): redesigned the priority to
  `priority = z(latent_genetic) + z(FAERS_IC025) + 3*OT`. FAERS disproportionality is
  the decisive "is this actually an ADR?" axis — a status code / non-ADR phenotype has
  no FAERS report (treated as worse than any observed signal), so it cannot top the
  list. Also switched the latent score's pleiotropy penalty from external PheWAS to
  the self-contained `internal_pleiotropy`.
- Result: every top hit now carries a STRONG FAERS signal; the status-code junk is
  gone; **lisinopril x angioedema (R60.0) is #1** (both loci), and real ADRs such as
  diazepam x COPD (PDE1C) surface. Manifest controls remain correctly at the bottom
  (F5, HLA-C, APOC1 priority −1.3 to −6.2, manifest_risk flagged) — validation intact.
- Residual (flagged, not yet fixed): off-label / symptom-coded indication confounding
  still leaks (fluoxetine x menopause #2, atenolol x elevated-BP #6) because those
  pairs DO have FAERS reports; the on_BNF gate only catches on-label ICD codes. Needs
  indication-relatedness logic (not a phenotype filter) to fully resolve.

### v27 — scale-normalized cross-phenotype BETA contrasts [superseded by v27.1]
- Problem: the 4 `delta_BETA_*` contrasts subtracted `log1p|BETA|` across arms, but
  the arms are not on a common scale. Empirically median |BETA| = drug-ADR 4.61
  (rare-variant logistic, quasi-separation inflation), disease 0.167, **dose 0.020
  (LINEAR ln-dose)**, prescribed 0.051 — a ~27× spread among the binary arms plus a
  fundamentally different unit for the linear dose arm. The raw contrasts measured
  scale artifact as much as biological amplification.
- Tested two fixes. (a) FULL normalization of all 4 contrasts (per-allele variance
  standardization b*=b·√(2f(1−f)) then within-arm z) nudged TIER agreement up but
  DESTABILIZED k-means (stable 204->146) and fragmented LCA (k=5->8), because blind
  within-arm z amplifies noise in the weak dose/prescribed arms. (b) SURGICAL fix —
  the three binary-vs-binary contrasts are comparable log-odds ratios (robust to a
  constant scale offset via log-ratio + global z) and kept as-is; ONLY the linear
  dose contrast (`delta_BETA_ADR_vs_dose`) is made unit-free by within-arm
  standardization of both arms.
- **Adopted the surgical fix.** Clustering structure preserved (k-means k=3 stable
  203, LCA k=5 stable 197 — same as v26) with LCA ARI 0.080->0.084; the dose-contrast
  unit bug fixed. The latent-PGx score's `beta_superiority` has the same dose term
  but it is inert (b_dose is tiny so the ADR-vs-dose term is never the binding min),
  so it was left unchanged.
- Same normalization applied to the **latent-PGx score's `beta_superiority`**. There
  the three superiority terms are combined with a `min()` (conjunction), so they must
  share a scale — the surgical (dose-only) approach is impossible. All four arms are
  put on the common within-arm-standardized footing (variance-scale then within-arm z)
  and `beta_superiority = min(zb_ADR−zb_icd, zb_ADR−zb_dose, zb_ADR−zb_pre)`. This now
  genuinely tests dose-dominance (previously the raw dose term was always huge and
  never the binding min). `p_superiority` is unchanged (LOG10P is already unit-free).
  Effect: gated candidate count unchanged (76; the gate does not use beta_superiority),
  biology-plausible unchanged (6); the latent score is reshuffled (Spearman 0.87 vs the
  raw version) into a more rigorous order.
- Published as v27: ATLAS_cluster_projections_v27.pdf, atlas_v27_classified.tsv,
  cluster_coords_v27.csv, PGX_resource_deck.pdf, latent shortlist.

### v26 — locus collapse (one lead signal per locus × phenotype) [superseded by v27]
- New step `collapse_loci.py`: for each (drug, ICD10) phenotype, co-located
  variants within a 250 kb window on the same chromosome are collapsed to a single
  representative — the **lead variant (highest drug-ADR -log10 P)**. Burden tests
  (no chr-pos) are each their own locus. Records `n_collapsed` + `collapsed_genes`.
- **300 variant-rows -> 207 locus × phenotype rows** (32 loci merged ≥2). Removes
  LD inflation: RYR2 x8->1, the simvastatin x I83.1 6-gene locus->1, APOE/APOC1
  Alzheimer locus->1, amitriptyline x M79.7 GRIA1 x20->1. **T4/Effect 26->14**
  (the inflated Effect copies removed) — the deflation the per-locus flag only
  approximated in v22.
- Features = v25 (27, self-contained) computed on the collapsed set. Clustering is
  CLEANER: LCA settles at k=5 (197/207 stable, vs the fragmented k=7 on v25),
  kmeans k=3 (204), HDBSCAN k=3 (207); GMM still k=10/unstable.
- latent-PGx: 76/207 gated candidates, 6 biology-plausible (the dedup removed
  duplicate LD-copies of the same candidate locus).
- Published as v26: ATLAS_cluster_projections_v26.pdf, atlas_v26_classified.tsv
  (the resource table: collapsed associations + cluster labels + latent score +
  n_collapsed/collapsed_genes), cluster_coords_v26.csv, PGX_resource_deck.pdf
  (LCA k=5 featured), latent shortlist. NB: report.py (xlsx) is hardwired to the
  300-row source sheet and is not regenerated for the collapsed set; the TSV is
  the v26 per-association deliverable.

### v25 — PROMOTED to active model (self-contained, 27 features) [superseded by v26]
- Switched the active/published model from v23 to **v25** for the PGx-targets paper.
- v25 = v23 minus all external gene-annotations (6 tissue, Druggability_ord,
  PheWAS_ord, external PheWAS_n_hits_log, beta_per_pleio_log) PLUS
  `internal_pleiotropy_log` (variant's significant-hit count across our own GWAS
  corpus). Evidence statistics recomputed on FAERS+replication (2 layers).
- **Fully self-contained / genome-scale-ready:** every feature derived from the
  project's own GWAS (only FAERS + on_BNF remain external; droppable for a pure-GWAS
  model). Internal pleiotropy per-tier: T1=17.8, T2=1.6, T3=4.0, T4(Effect)=1.4
  (Effect = most variant-specific).
- **Tradeoff (honest):** clustering is coarser than v23 — kmeans k=3, LCA k=7
  (296 stable), GMM k=10 (unstable), HDBSCAN k=4. A dominant low-evidence mass
  appears because the discriminating structure partly lived in the removed tissue
  features. The **latent-PGx prioritisation is unchanged** (independent of the
  feature set: same 85 gated candidates, 7 biology-plausible) — the paper's core
  deliverable is intact; the clustering is a coarser descriptive layer.
- Published as v25: ATLAS_v2_CLUSTERED_v25.xlsx, ATLAS_cluster_projections_v25.pdf,
  cluster_coords_v25.csv, PGX_resource_deck.pdf (LCA k=7 featured), latent shortlist.

### latent_pgx_score as a clustering feature — REJECTED (v23 retained, superseded by v25)
- Tested adding the statistical latent score (no OT) as a 37th feature, winsorized.
- **Harmful:** GMM k=5→10, ARI 0.218→0.147, stable 298→67; other methods flat.
- Three reasons it fails: (1) redundant — it's a nonlinear min-combination of
  features already in the model (delta_BETA contrasts, LOG10P diffs, P_test,
  pleiotropy); (2) reintroduces heavy-tailed instability (manifest predispositions
  at −7…−10); (3) conceptually a post-hoc prioritisation aimed at a target, not a
  structural axis — and its full form needs OT, which must stay out of clustering.
- Decision: **keep latent-PGx as a reported prioritisation** (score column + plot
  overlay + shortlist), NOT a clustering input. The model finds structure; the
  latent score ranks within it. Accepted model stays v23.

### FAERS-redundancy test — REJECTED prune (v18 retained, now v20)
- Tested dropping FAERS features (corr r=0.78–0.93, looked redundant like the SE
  family). Every prune HURT: drop-ord ARI 0.218→0.149, keep-IC025-only →0.177,
  keep-PRR-only k 5→10, drop-all →0.158. All raised k above the clean k=5.
- Lesson: correlated ≠ droppable. The 3 FAERS metrics (Bayesian IC025,
  frequentist PRR, consensus ordinal) are independent measurements of one signal
  and act as a noise-reducing ensemble — a *reinforcing* block, unlike A1_FREQ
  which was a *noisy copy*. Kept all 3. **v18 (34 features) stands as final.**

---

## 3. Tissue methodology redesign — Option B (the biggest reviewer surface)

### Why
Original `Tissue_call` (7-level categorical) used undocumented ICD10→tissue
mapping and unpublished thresholds (TPM≥50 or z≥1.5). Tissue is one of three
evidence layers in the cascade's T3/T4 decision → most exposed to critique.

### What was built (stages 1–12)
1. **ICD10 → tissue mapping**: MONDO obo (1834 ICD10 xrefs) → UBERON → GTEx-54;
   plus WHO chapter fallback. Later **superseded by manual curation** (§3.3).
2. **GTEx v8 median TPM** (56,200 genes × 54 tissues).
3. **Tau index** (Yanai 2005) per gene + per-tissue log-TPM z-score.
4. **HPA RNA consensus** replication (independent dataset; 18,773 shared genes).
5. **GTEx↔HPA concordance**: Spearman ρ per gene; 5,889 (31%) flagged
   disagreement (ρ<0.5).
6. **MAGMA gene-property** tissue enrichment (de Leeuw 2015) per (drug,ICD10):
   99/177 pairs ok, FDR-corrected across HPO-mapped tissues, 1000-perm null.

### Sensitivity finding (the substantive result)
- Old tissue mechanism (specific/enriched): **92/300 = 30.7%**
- New (rigorous specificity, MAGMA FDR<.05 or z≥1): **7/300 = 2.3%**
- **T4 rows losing tissue support: 25/26 = 96%.** The cascade's Effect tier was
  largely driven by the gene's *primary* tissue specificity, not *disease-relevant*
  tissue enrichment. Example: lisinopril×R73.9×RYR2 — RYR2 is heart-specific, but
  R73.9 (hyperglycemia) maps to kidney/bladder where RYR2 z=−1.2.

### Agreement diagnosis — why MAGMA & tau were dropped (v10)
Tested expert `Tissue_call` ordinal against each method:

| metric (what it measures) | ρ vs expert | binary κ | decision |
|---|---|---|---|
| MAGMA FDR (polygenic tissue enrichment) | −0.009 | 0.00 | **DROP** — wrong tool for oligogenic (1–3 gene) pharmaco-ADR signals; MAGMA needs hundreds of genes |
| disease-tissue z-score (specificity) | +0.01 | — | KEEP as the rigorous metric |
| gene-centric max z-score (specificity) | −0.26 | — | — |
| tau index | −0.26 | — | **DROP** — measures specificity to ANY tissue, not disease tissue |
| disease-tissue absolute log-TPM (abundance) | +0.23 | +0.21 | **ADD** (v10) — the only metric that tracks the expert |

Conclusion: the expert `Tissue_call` tracks **absolute expression abundance**, not
specificity (its "TPM≥50 OR z≥1.5" definition is dominated by the absolute arm).
Maximum achievable agreement with any principled metric is κ≈0.21. Improving the
mapping rigor *lowers* nominal agreement (0.41→0.03) because it reveals the expert
call was never disease-tissue-specific. This is the rigor story, not a defect.

### §3.3 Mapping curation (A/B/C, requested "try all")
- **(A) Open Targets/EFO**: its tissue data IS GTEx+HPA (circular); EFO anatomy
  too sparse. No independent signal. Honest negative.
- **(B) Manual curation** of all 113 atlas ICD10 codes → disease-relevant GTEx
  tissues (`icd10_curated_tissues.py`): 100 strong + 13 weak-proxy (eye/ear/bone
  have no GTEx tissue), 0 unmapped (was 64). This is the production map.
- **(C)** = v10 (drop MAGMA+tau, add absolute-TPM, recluster).

---

## 4. Clustering / labeling decisions

| # | Decision | Rationale |
|---|----------|-----------|
| 4.1 | GMM uses `covariance_type="tied"`, `reg_covar=1e-4` | diag covariance let BIC run to k=10; tied stabilizes k* |
| 4.2 | HDBSCAN `gen_min_span_tree=True` | required for `relative_validity_` |
| 4.3 | Confusion via `pd.crosstab`, not sklearn `confusion_matrix` | shape mismatch when cluster/tier label sets differ |
| 4.4 | Hungarian (linear_sum_assignment) for k=4 cluster→TIER alignment | optimal 1:1 mapping for agreement; majority-vote/exhaustive for k>4 merges |
| 4.5 | EPIC via anchor-majority vote; E/P split by median N_INDEPENDENT | anchors are ground truth; the E/P split axis is the one place a cascade stat is used, for labeling only (not clustering) |
| 4.7 | **EPIC "C" redefined: Comorbidity → inConclusive** | Data shows T2 is not a comorbidity signature: 114 cascade-T2 rows = 54% weak-on-everything, only 28% genuine comorbidity (strong pure-disease assoc + no drug amplification), 18% mixed. T2 is negative on EVERY axis incl. LOG10P_ICD10 — it is the low-evidence default bucket. The EPIC acronym is preserved (E/P/I/**C**=inConclusive). |
| 4.6 | Ward used for the deterministic per-cluster anatomy / naming | reproducible, accepts n_clusters |

---

## 5. Structural findings (stable across all versions)

1. **T4 (Effect) is a continuous-AND, not a cluster.** Distance/mixture methods
   cannot represent "FAERS-STRONG ∧ replicated ∧ tissue-mechanism" without an
   explicit conjunction feature. Continuous AND-proxies (`evidence_min`) recover
   it; the cascade's hard N≥3 cut is the only thing that bisects T3 from T4.
2. **T3 (Predisposition) has no shape.** It is the cascade's leftover bucket
   between "no mechanism" (T2) and "full mechanism" (T4); recall ceilings ~18%
   under distance methods, ~50–63% only when k≥7 sub-clusters are merged.
3. **The data's natural structure is ~3 classes** (Indication / Comorbidity /
   Mechanism), not 4. The T3/T4 split is an expert sub-cut along a continuum.
4. **Cluster geometry is dominated by a rare-vs-common-variant axis** (SE/BETA
   family) and an indication axis (`on_BNF`, `LOG10P_ICD10`). Tissue and FAERS
   biology are clean but small contributors.

---

## 6. Infrastructure decisions

| # | Decision |
|---|----------|
| 6.1 | Docker Compose + Apache (httpd:2.4) serving `/slow/www` on :8899; outputs published to `/slow/www/atlas/` |
| 6.2 | Source `ATLAS_v2.xlsx` is READ-ONLY; all outputs are additive (new columns/sheets, new files) |
| 6.3 | Each version published as `ATLAS_v2_CLUSTERED_v<N>.xlsx` for comparison |
| 6.4 | `data/clustering/` and `data/tissue/` are gitignored (large/derived) |

---

## 7. Version → file map

| version | builder module | key change | features |
|---------|----------------|-----------|----------|
| v3  | build_features.py (+N_INDEP edit) | pure raw + interactions | 31 |
| v6  | build_features_v6.py | drop curator/burden, ADR rename | 32 |
| v7  | build_features_v7.py | +dose/prescribed GWAS | 44 |
| v8  | build_features_v8.py | prune + reweight | 36 |
| v9  | build_features_v9.py | Option B tissue | 42 |
| v10 | build_features_v10.py | drop MAGMA+tau, +abs-TPM | 40 |
| v11 | build_features_v11.py | drop LOG_OR_SE | 39 |
| v12 | build_features_v12.py | winsorize SE family | 39 |
| v13 | build_features_v13.py | drop evidence_sum (redundant) | 38 |
| v14 | build_features_v14.py | drop evidence trio — REVERTED | 35 |
| v15 | build_features_v15.py | Stouffer_Z + Nichols conjunction_minZ | 36 |
| v16 | build_features_v16.py | symmetric GWAS core (BETA+LOG10P ×4, SE for ADR only) | 35 |
| v17 | build_features_v17.py | ADR_pct_before → directional ratio ADR_after_over_before | 35 |
| v18 | build_features_v18.py | drop A1_FREQ_log (redundant w/ A1_CASE_CT_log1p) | 34 |
| **v20** | **build_features_v20.py** | **remove all manual feature weights (no arbitrary multipliers)** | **34** |

Supporting modules: `enrich_with_gwas.py`, `tissue_pipeline.py`,
`magma_pipeline.py`, `merge_magma.py`, `icd10_curated_tissues.py`,
`tissue_agreement.py`, `tissue_sensitivity.py`, `feature_audit.py`,
`feature_importance.py`, `inventory.py`, `ward_t3.py`, `k4_overlap.py`.

---

## 8. Open / deferred

- T3 recovery beyond ~18% would require a non-distance method (decision-tree
  mixture, Bayesian rule lists) — different model class, deferred.
- MAGMA permutation-null columns computed but not used as features (MAGMA
  dropped); retained in `magma_results.tsv` as supplementary.
- GMM stability regression under v12 winsorization — monitor; Ward/k-means are
  the stable reporting methods.

---

## 9. v25–v29 (self-contained model, locus collapse, Set B unification)

| version | builder module | key change | features |
|---------|----------------|-----------|----------|
| v25 | build_features_v25.py | self-contained: drop Tissue / Druggability / OT / beta_per_pleio; pleiotropy computed from our own GWAS corpus | 27 |
| v26 | (v25 + `collapse_loci.py`) | collapse co-located variants to one lead signal per (drug, ICD10); 300 → 207 rows | 27 |
| v27 | build_features_v27.py | BETA normalised across arms (dose is LINEAR, others logistic); same normalisation in the latent score | 27 |
| v28 | build_features_v28.py | merge the never-read `ATLAS_full` sheet (`merge_atlas_full.py`); add temporal, test-set genotype and VEP descriptors; drop `FAERS_ord`, `phenotype_nonspecific` (threshold-derived) | 41 |
| **v29** | **build_features_v29.py** | **unify Set B into the model: `temporal_q25/q75_days_after` replace `temporal_iqr_days_after`; `test_A1_CASE_FREQ`/`test_A1_CTRL_FREQ` replace `test_case_ctrl_freq_diff`** | **43** |

### 9.1 v29 — Set B retired

Set B (the diagnostic projection) contained nine descriptors. Seven were already present in the
model as transformations, and `OR_logdiff` / `OR_ratio` were exact duplicates of model features.
The only genuinely new information was the onset-window **location** (q25/q75 separately, not
just their width) and the baseline allele-frequency **level** (case/control separately, not just
their difference). Both were folded in; storing the components *and* their difference would be
exactly collinear, so the differences were dropped.

Result: the model itself now separates the curated good/bad associations on UMAP
(distance/spread ratio **1.01**, i.e. separated) while achieving the best independent
Effect-vs-Indication control on record (**0.83 PCA / 0.87 UMAP**). The residual Set B projection
reached good/bad **1.45** but its control *fell* to **0.67** — the definition of overfitting to
the two curated "bad" pairs. **Set B is therefore retired**: `diagnostic_projection.py` is kept
for reproducibility only, its plot page removed from `plot_clusters.R`, and
`coords_diagnostic.csv` is no longer published.

Clustering (v29): k-means k=3 (205 stable), GMM k=10 (18 stable), **LCA k=4 (197 stable, primary)**,
HDBSCAN k=1 (degenerate at 43 dimensions — reported as a robustness check, not a result).
`kappa_EPIC`: kmeans 0.046, GMM 0.172, LCA 0.074; LCA ARI 0.035.

### 9.2 Kruskal-Wallis audit of the class separators (external cross-check)

An independent prototype ranked the parameters that separate the expert classes by
Kruskal-Wallis H. Re-running the test on our data reproduces the temporal and count findings but
overturns two of the interpretations.

**Reproduced** (v29 collapsed set, medians by expert class):

| parameter | Indication | Comorbidity | Predisposition | Effect | reading |
|---|---|---|---|---|---|
| after/before diagnosis ratio | 5.3× | 14.4× | 20.7× | **32.0×** | Effect signature — confirmed |
| `TEMP_pct_before` | **15.2%** | 5.7% | 3.7% | 4.0% | Indication signature — confirmed |
| `A1_CASE_CT` | 2 | **1** | 1 | 1 | Comorbidity is the lowest-powered group — confirmed |

`TEMP_pct_before` is the strongest *honest* single separator of Indication from Effect
(CV AUC **0.766**); the after/before ratio is the strongest Effect marker (AUC 0.764 in the
Effect direction).

**Overturned 1 — `is_BNF_indication` is circular, not a discovery.** It is 1 for 50/81 Indication
rows and **0 for all 219 non-Indication rows** — zero false positives, because the expert cascade
*assigns* the Indication tier when the ICD10 is on the drug's BNF label. Its H therefore measures
the cascade's own definition. (The quoted H = 278 is also arithmetically impossible on 207 rows,
where H ≤ N−1 = 206; on the uncollapsed 300 rows we obtain H = 161.7. The inflation comes from
counting co-located LD-correlated variants repeatedly — exactly what the v26 collapse fixed.)

Consequence for us: `on_BNF` is a model feature. Ablating it leaves k-means and **LCA (primary)
unchanged** (`kappa_EPIC` 0.046 / 0.074), but drops **GMM from 0.172 → 0.096** and its
`ari_vs_TIER` from 0.029 → 0.002. **GMM's apparently-best EPIC overlap is roughly half label
leakage.** LCA is clean and remains the primary method; GMM's kappa must be reported with this
caveat, or with `on_BNF` removed.

**Overturned 2 — FAERS marks Indication, not Effect.** The prototype read
`faers_prr / faers_ror / faers_ic025` as supporting the effect axis. The direction is the
opposite: `FAERS_IC025` medians are Indication **1.59** > Effect 1.09 > Predisposition 0.96 >
Comorbidity 0.65 (on the uncollapsed sheet the gap looks far larger, 3.11 vs 1.09 — again the
duplicate-variant inflation). This is confounding by indication in spontaneous reports: the event
that prompted the prescription is co-reported, inflating disproportionality for on-label pairs.

Consequence: the latent-PGx score's FAERS axis (weight 1.0) mildly *rewards* indications.
`fluoxetine × N95.1` — a genuine indication that required the curated `EXPERT_DOWNWEIGHT` — sits
at the **86th** FAERS percentile, and `ramipril × R41.0` (the zero-test-carrier artifact) at the
**89th**, while `citalopram × M81.9` (a real ADR) sits at the **44th**. The separation is weak
(Indication-vs-Effect AUC 0.565), so FAERS still earns its place as pharmacovigilance
corroboration — but it must **not** be described as evidence of a causal drug effect, and the
expert downweight has been partly compensating for this.

### 9.3 Parameters NOT adopted from the audit

- `validated_test` (`VAL`/`NOT`) and `MIN_CELL_FLAG_regenie` (`PASS`/`ORANGE`) are
  threshold-derived binaries (a p-value cutoff and a minimum-cell-count cutoff) and violate the
  threshold-purity rule (§1). `validated_test` is additionally a binarisation of
  `repl_directional_z`, which the model already carries continuously.
- `no_carrier_test` is a restatement of zero test-set carriers, already captured by
  `repl_directional_z = 0` and by `test_log_case_het` / `test_log_case_hom`.
- `FAERS_ROR` is ρ ≈ 1.00 with `FAERS_PRR`; it adds no information (and v28 already dropped
  `FAERS_ord` from the same family).

### 9.4 v30 — `ADR_pct_before` added (44 features)

`build_features_v30.py` adds one feature: `ADR_pct_before` = `TEMP_pct_before`, the share of
diagnoses recorded before drug exposure. §9.2 identified it as the strongest **threshold-free**
single separator of Indication from Effect (CV AUC 0.766), and the model had carried it only
indirectly, as an input to `ADR_after_over_before`. It is not collinear with the existing pair:
`pct_after + pct_before + pct_other = 100`, so `pct_other` leaves a free dimension, and the
log-ratio is non-linear in both.

| metric | v29 (43) | v30 (44) | verdict |
|---|---|---|---|
| LCA (primary) `stable_rows` | 197 | **199** | best on record |
| LCA `kappa_EPIC` / `ari_vs_TIER` | 0.074 / 0.035 | 0.072 / 0.034 | unchanged (seed noise) |
| k-means `stable_rows` | 205 | **206** | improves |
| GMM `kappa_EPIC` | 0.172 | 0.160 | see §9.2 — roughly half of GMM's kappa is `on_BNF` leakage |
| latent-PGx AUC (E vs I) | 0.909 | 0.909 | unchanged (the score does not read the feature matrix) |
| good-vs-bad PCA / UMAP | 0.71 / 1.01 | **0.74 / 1.22** | improves |
| E-vs-I control PCA / UMAP | 0.83 / 0.87 | 0.83 / 0.83 | PCA holds; UMAP −0.04, not material |

**Accepted.** Consensus stability rises on both primary methods and the good/bad separation
improves while the independent control holds — i.e. the gain is not the overfitting signature
that retired Set B. EPIC overlap is flat, which is the expected result for a feature that
distinguishes Indication from Effect: those are two of four tiers, and `kappa_EPIC` is dominated
by the Comorbidity/Predisposition bulk.

The interpretive payoff is larger than the metrics. Under v30 the LCA solution's two large
clusters split **on exactly this axis**: C3 (n=88) has `pct_before` 8.6% and an after/before ratio
of 9.6×, C4 (n=92) has 2.1% and 26.7×. The clusters were therefore renamed from their parameter
medians (not from the expert tiers):

| LCA cluster | n | signature | name |
|---|---|---|---|
| C1 | 7 | 28 case carriers, pleiotropy 14, FAERS IC025 2.10, 43% on-label | Common pleiotropic locus (on-label) |
| C2 | 20 | 0.14 case carriers, pleiotropy 1 | Carrier-depleted rare signal |
| C3 | 88 | `pct_before` 8.6%, after/before 9.6× | Pre-exposure onset (indication direction) |
| C4 | 92 | `pct_before` 2.1%, after/before 26.7×, highest latent 0.54 | Post-exposure onset (drug-effect direction) |

Housekeeping: the retired Set B artifact `coords_diagnostic.csv` was deleted, so
`separation_metrics.json` now contains only the `set_A_model` block; the projection-plot title
now reads its feature count from `pca_var.csv` instead of the hard-coded "27-feature".

### 9.5 Naming I / P / E / C from the parameters (`epic_from_parameters.py`)

Each EPIC state was written as a **parameter signature**, not as a cascade rule, and every
association takes the arg-max of the four axes. The expert TIER is never read; it is used only
afterwards, to score agreement.

| state | meaning | axis |
|---|---|---|
| **I** Indication | the drug is given for the condition | `z(on_BNF) + z(pct_before)` |
| **P** Predisposition | germline variant causes the condition regardless of the drug | `z(LOG10P_ICD10) − z(beta_superiority)` |
| **E** Drug-triggered | the variant manifests only under exposure | `z(beta_superiority) + z(after/before) + z(repl_z)` |
| **C** Confounder | no causal link; drug and condition share a common cause (age) | `z(pleiotropy) + z(LOG10P_prescribed) − z(beta_superiority)` |

P and C share "the drug adds nothing"; they are separated by *where* the signal lives — P by a
dominant disease-alone association, C by pleiotropy plus prescription propensity.

**The confounder axis works, and it validates on known biology.** The two highest C scores in the
atlas are `simvastatin × G30.9` (Alzheimer's) at **rs4420638 — the APOE/APOC1 locus**, C = 6.28,
and `simvastatin × I25.1` (ischaemic heart disease) at **rs10455872 — the LPA locus**, C = 7.56.
Both are textbook confounding by age / shared cause: the variant causes the disease and the statin
is simply what that population takes. **The expert cascade labelled both "Indication."** Statins
are not indicated for Alzheimer's; this is a cascade artifact that the parameter axes catch.

**Agreement with the experts is poor, and that is the finding.** Row-level Cohen's kappa = 0.067
(accuracy 0.227); cluster-level kappa = −0.031. Per-state recall against the expert calls:

| expert state | n | recovered | note |
|---|---|---|---|
| E (Effect) | 14 | **11 (79%)** | the drug-triggered signature is real and reproducible |
| I (Indication) | 33 | 17 (52%) | |
| P (Predisposition) | 67 | 14 (21%) | mostly re-assigned to E |
| C (Comorbidity/inConclusive) | 93 | 5 (5%) | the tier the experts could not resolve either |

So the parameters recover the **Effect** state well and the expert **C** tier not at all — which is
consistent with C having been the cascade's residual bucket rather than a biological state. This is
evidence *for* the data-driven reclustering thesis, not against it, but it must be reported as a
disagreement rather than a validation.

**Cluster-level naming is unstable and should not be used.** The mean-axis arg-max gives
C1 → P (3.02) with C a statistical tie (2.98): the cluster holds both confounders (APOE, LPA) and a
strong predisposition locus (tramadol × L40.9 psoriasis, P = 5.14), and the mean cannot separate
them. Per-row naming does. C2 → C is assigned on a mean of 0.02, i.e. noise. Only C3 → I and
C4 → E are meaningful.

**E is the prioritised state.** Median latent-PGx by parameter-derived state: **E 0.539** >
P 0.492 > I 0.479 > C 0.383 (E vs rest, one-sided Mann-Whitney **p = 3.6e-07**). The same holds by
LCA cluster: C4 (E) 0.538 vs C1 (confounder-dominated) 0.186. The score and the naming agree without
having been fitted to each other.

Report: `ATLAS_cluster_projections_v30.pdf` gains a 4th page — PCA and UMAP coloured by the
parameter-derived state, highlighted associations marked, and latent-PGx boxed by state.

### 9.6 The same EPIC-axis naming applied to every clustering method

`epic_from_parameters.py` now names the clusters of all four methods. `margin` = winning axis
minus runner-up; a cluster is *decisive* when margin >= 0.10.

| method | k | states used | decisive clusters | median margin | cluster-level kappa vs expert |
|---|---|---|---|---|---|
| k-means | 3 | C, E, P | 1 / 3 | 0.057 | −0.022 |
| **GMM** | **10** | **C, E, I, P** | **7 / 10** | **0.740** | −0.012 |
| LCA | 4 | C, E, I, P | 2 / 4 | 0.133 | −0.031 |
| HDBSCAN | 1 | I | 0 / 1 | 0.000 | 0.000 |

**Every cluster-level kappa is ~0 or negative** — no method's clusters reproduce the expert tiers.
Only the *row-level* parameter labelling reaches kappa 0.067 (§9.5). k-means at k=3 cannot express
Indication at all (it never assigns I), and HDBSCAN's single cluster is vacuous.

**GMM at k=10 is the only solution that resolves all four states, and it does so on known biology.**
Its small components are near-pure and their contents are checkable:

| GMM cluster | n | assigned | margin | contents |
|---|---|---|---|---|
| C8 | 1 | **C** | 1.88 | `simvastatin × G30.9` at **rs4420638 (APOE/APOC1)** — Alzheimer's |
| C9 | 3 | **C** | 2.02 | `simvastatin × I25.1` at **rs10455872 (LPA)** and **rs10757270 (9p21/CDKN2B-AS1)**, plus `atenolol × I25.1` |
| C10 | 4 | **P** | 3.02 | `tramadol × L40.9` (psoriasis) and `simvastatin × I80.2` (phlebitis) at **rs6025 — Factor V Leiden** |
| C2 | 7 | **I** | 2.51 | `simvastatin × E78.5`, `clopidogrel × I63.9`, `ramipril × I21.9`, `amitriptyline × M79.7` — **7/7 expert-labelled I** |

C8 and C9 are the two canonical coronary-artery-disease loci plus the canonical Alzheimer's locus,
all in statin/beta-blocker users — confounding by age and by shared cause, exactly as the C axis
defines it. C10 puts Factor V Leiden with venous phlebitis: the variant causes thrombosis whether
or not a statin is taken, which is the P signature. **The experts labelled every row of C8/C9 as
"Indication" and both Factor V Leiden rows as "Comorbidity";** the parameter axes disagree with the
cascade and agree with the genetics.

GMM's C2 is a clean Indication component (7/7 expert agreement) — the single place where the
parameter naming and the cascade coincide perfectly.

**Reporting consequence.** GMM has the most interpretable decomposition but only 19/207 stable rows
(§9.1) and its `kappa_EPIC` is roughly half `on_BNF` leakage (§9.2). LCA stays primary for the
consensus partition; GMM is reported as the *interpretive* decomposition, with its instability
stated. The two roles must not be conflated.

`ATLAS_cluster_projections_v30.pdf` is now 5 pages; the 5th tabulates the EPIC-axis scores for all
four methods. New artifacts: `epic_axes_by_cluster.csv` (all methods), `epic_method_comparison.csv`.

### 9.7 GMM-only report; the ten components named, and why states repeat

The report (`plot_clusters.R`, `ATLAS_cluster_projections_v30.pdf`, 3 pages) now shows **GMM only**.
k-means (k=3) cannot express Indication at all, HDBSCAN degenerates to one cluster, and LCA resolves
only 2 of 4 states decisively (section 9.6). Panels for those three were removed. LCA remains the
primary *consensus partition* in the metrics tables; GMM is the *interpretive* decomposition.

**The ten components** (`gmm_cluster_names.py` -> `gmm_clusters_named.csv`). The state letter is the
arg-max of the four parameter axes; the name is read off the component's parameter medians. `margin`
is the winning axis minus the runner-up; **7 of 10 are decisive (margin >= 0.10)**.

| comp | state | name | n | latent | margin | decisive |
|---|---|---|---|---|---|---|
| C8 | C | Prescription-driving locus (APOE) | 1 | 0.09 | 1.88 | yes |
| C9 | C | Coronary-risk locus (LPA, 9p21) | 3 | 0.16 | 2.02 | yes |
| C10 | P | Major-effect disease locus (Factor V Leiden) | 4 | 0.20 | 3.02 | yes |
| C7 | P | Sub-threshold predisposition | 20 | 0.51 | 0.09 | no |
| C2 | I | On-label indication (BNF-confirmed) | 7 | 0.48 | 2.51 | yes |
| C4 | I | Pre-existing condition (timing only) | 45 | 0.47 | 0.60 | yes |
| C3 | E | Replicated drug-triggered signal | 83 | 0.53 | 0.49 | yes |
| C6 | E | Unreplicated drug-triggered signal | 19 | 0.52 | 0.07 | no |
| C5 | E | Drug-triggered, no pharmacovigilance support | 5 | 0.22 | 0.88 | yes |
| C1 | – | Carrier-depleted (undecidable) | 20 | 0.49 | 0.03 | no |

**Why a state appears more than once.** Each repetition is a distinct mechanism, separated by a
measurement rather than a judgement:

- **C twice — two ways to share a cause.** C8: the variant drives the drug's *own indication*
  (`LOG10P_prescribed` = 35.2; rs4420638/APOE raises lipids, the statin is prescribed, and
  Alzheimer's is then scored as an "ADR"). C9: the variant drives a *comorbid disease of the treated
  population* (LPA, 9p21 -> coronary artery disease; statins and beta-blockers are what those
  patients take). Only C8 has a prescription-propensity signal — that single parameter separates the
  two mechanisms.
- **P twice — two effect sizes.** C10 is major-effect disease loci (`LOG10P_ICD10` = 101.6; Factor V
  Leiden rs6025 -> phlebitis) with `beta_superiority` = 0.22, i.e. the drug adds nothing. C7 is the
  same shape 100 orders of magnitude weaker (`LOG10P_ICD10` = 1.0) — sub-threshold predisposition.
- **I twice — two independent evidence sources.** C2 is *label* evidence (`on_BNF` = 1.00 for all
  7 rows; 7/7 expert agreement). C4 is *timing* evidence only (`on_BNF` = 0.00, yet `pct_before` =
  16.0%) — off-label use or protopathic bias. The two sets are disjoint on `on_BNF`, which is exactly
  why no single axis merges them.
- **E three times — which corroboration is present.** C3 replicates out-of-sample *and* has
  pharmacovigilance support. C6 has FAERS but does not replicate (`repl_z` = 0.56, zero test
  carriers) — the artifact-prone component. C5 replicates (`repl_z` = 2.41) but has no FAERS
  (0.46).

**C1 is not a confounder.** Its margin is 0.03: undecidable. It is the carrier-depleted component
(median `A1_CASE_CT` = 0.14) and is labelled "–" rather than given a state. C6 (0.07) and C7 (0.09)
are likewise non-decisive; their letters are provisional and the report says so on the page.

**A caveat the naming exposes.** C5 (replicated, no FAERS) has median latent **0.216**, while C6
(unreplicated, FAERS-supported) has **0.520**. The latent score is therefore driven harder by FAERS
than by held-out replication — and FAERS marks Indication, not Effect (section 9.2). The `repl_z`
weight of 0.7 was set to *discount* replication because it is beta-derived; the consequence is that
a component with genuine out-of-sample support scores below one without it. This is now visible in
the report rather than buried in the score.

**Latent-PGx display.** Page 2 shows the score three ways — whole atlas (n=207), by EPIC state, and
by GMM component — with the **6 curated target ADRs** (circles) and **4 expert-flagged confounders**
(crosses) overplotted on every panel. E is the prioritised state (E vs rest, one-sided
Mann-Whitney p = 3.6e-07).

### 9.8 latent-as-a-feature REJECTED; curated confounder overlay added instead

**Request:** add the latent-PGx score as a clustering parameter so that `simvastatin x I83.1`
(varicose veins) moves into the confounder state.

**Tested, and it does neither.** Appending latent as a 45th feature, in both its penalised and its
unpenalised form:

| | v30 (44 feat) | + latent (penalised) | + latent (penalty removed) |
|---|---|---|---|
| GMM k | 10 | 8 | 8 |
| GMM `stable_rows` | 19 | 22 | 20 |
| `kappa_EPIC` | 0.160 | 0.081 | 0.092 |
| `simvastatin x I83.1` | C3, state **E** | C8, state **E** | C8, state **E** |

The association stays **E** even with the curated 0.35x penalty baked directly into the feature,
while GMM's interpretive decomposition collapses from k=10 to k=8 and `kappa_EPIC` halves. This
reproduces the **v24** result (latent as a feature destabilised GMM; reverted).

**It is also circular.** The latent score is a deterministic function of features already in the
matrix (`repl_z`, FAERS, temporal, `beta_superiority`, `-pleiotropy`) multiplied by the curated
`EXPERT_DOWNWEIGHT` flag. Feeding it back in launders an expert judgement into the clustering that
the paper claims is free of expert classification.

**Why the association cannot be moved on evidence.** Against the two genuine confounder components:

| | LOG10P_ICD10 | LOG10P_prescribed | pleiotropy | beta_superiority |
|---|---|---|---|---|
| `simvastatin x I83.1` | 0.07 – 0.69 | 0.27 – 1.26 | 1 | **3.5** |
| C8 (APOE) | 99.5 | **35.2** | 15 | 0.13 |
| C9 (LPA / 9p21) | 28.9 | 4.74 | 16 | 0.08 |
| *atlas median* | *0.54* | *0.30* | *1* | *3.34* |

It carries the drug-triggered signature on every axis. Its latent score of 0.19 comes **entirely**
from the curated penalty: remove that and it ranks **43rd-56th of 207**, above the atlas median.
Epidemiologically this is expected — confounding by shared cause requires the *variant* to be
associated with the exposure, and the drug-prescribed GWAS arm is null here (`LOG10P_prescribed`
approximately the atlas median). The expert objection is to the **drug x condition pair** (statins
are protective against varicose veins), which is a pair-level belief that no variant-level
parameter can carry.

**Also rejected: an ICD-chapter-overlap descriptor.** `chap_match` (does the ADR's ICD10 chapter
appear among the drug's BNF indication chapters?) has real signal — 0.515 in expert Indication vs
0.118-0.209 elsewhere — and it does flag `simvastatin x I83.1`. But it also flags
`gabapentin x R42`, one of the curated **targets**, and misses `fluoxetine x N95.1`. It is an
indication-adjacency measure, not a confounder detector, and was not added.

**Adopted instead: an explicit overlay** (`epic_from_parameters.CURATED_CONFOUNDERS`). Two columns
are now written to `coords.csv` and the shortlist:

- `epic_param` — what the parameters say. **Untouched**, so the disagreement stays auditable.
- `epic_final` — `epic_param`, forced to `C` on the curated pairs.
- `epic_source` — `"parameters"` or `"curated"`.

Four rows are overridden (3x `simvastatin x I83.1`, 1x `fluoxetine x N95.1`); for all four the
parameters said something other than C (E and I respectively). `epic_final` counts:
E 94 / I 49 / P 47 / C 17. The report's page 1 now plots `epic_param` and `epic_final` side by side,
so the override is visible on the figure rather than silent in a table.

**The principle:** curated biology belongs in a labelled annotation layer on the output, never in
the feature matrix. This keeps the clustering honestly unsupervised and still lets the shortlist
demote the pairs we know are confounded.

#### 9.8a Overlay was applied inconsistently (fixed)

The first cut of the overlay wrote `epic_final` but only used it in **one panel of four**. Every
other view still coloured by `gmm_state` — the *component's* state — so `simvastatin x I83.1`, which
sits in component C3, kept reading **E**. Two distinct things were being conflated:

- `gmm_state` — a property of the **component** (C3 is a "Replicated drug-triggered signal")
- `epic_final` — a property of the **row** (this pair is curated as a confounder)

A curated row can legitimately sit inside an E component; that is the disagreement, not a bug. But
every per-row display must use the per-row state. Fixed:

- report page 2, "By EPIC state" -> now `epic_final`; its title says the overlay is applied
- report page 2, "By GMM component" -> fill relabelled "component state", title says so, so the
  distinction is explicit rather than implied
- `latent_pgx_shortlist.tsv` had **no EPIC column at all** -> now carries
  `epic_param` / `epic_final` / `epic_source`
- `atlas_v30_classified.tsv` -> now carries `gmm`, `gmm_state`, `gmm_name`, and the three epic columns

`simvastatin x I83.1` reads `epic_param=E, epic_final=C, epic_source=curated` in every published
artifact. It does not appear in the shortlist at all — the gate already excludes it, because
`EXPERT_DOWNWEIGHT` sets `ew_pen < 1.0`.

### 9.9 Curated overlay removed; unsupervised (GMM->EPIC) vs supervised EPIC

Per request, the curated confounder overlay (section 9.8/9.8a) is **removed** entirely -- no
`epic_final`, no `epic_source`, no per-pair override anywhere in the pipeline or the published
artifacts. The EXPERT_DOWNWEIGHT layer inside the latent score is unchanged (it governs the
shortlist gate, not the EPIC labelling).

The EPIC labelling is now fully cluster-level and unsupervised:

- **Unsupervised EPIC** (`epic_unsup`): each of the ten GMM components is mapped to exactly one
  EPIC state -- the arg-max of its mean parameter axes (`cluster_epic_map`). Every association
  inherits its component's state. No expert label is read. The map:
  C1->C, C2->I, C3->E, C4->I, C5->E, C6->E, C7->P, C8->C, C9->C, C10->P. (C1 was previously left
  "-" for a 0.03 margin; it is now given its arg-max C but flagged low-confidence.)
- **Supervised EPIC** (`epic_sup`): the expert cascade via TIER_TO_EPIC = {1:I, 2:C, 3:P, 4:E}.

**They disagree strongly** -- kappa = -0.012, row agreement 15.5%. The class balances are almost
inverted:

| | I | P | E | C |
|---|---|---|---|---|
| unsupervised (GMM->EPIC) | 52 | 24 | **107** | 24 |
| supervised (expert) | 33 | 67 | 14 | **93** |

The parameters put most associations in **E** (drug-triggered); the experts put most in **C**
(the cascade's Comorbidity/inConclusive residual). Confusion matrix (rows = expert, cols =
GMM->EPIC), written to `epic_sup_vs_unsup.csv`:

```
           unsup_I  unsup_P  unsup_E  unsup_C
sup_I           13        3       11        6
sup_P           13        5       36       13
sup_E            4        1        9        0
sup_C           22       15       51        5
```

The largest off-diagonal cell is expert-C -> unsup-E (51): half of what the experts could not
resolve, the parameters call drug-triggered. This is the paper's core claim made visual -- the
data-driven states do not reproduce the expert tiers, and the disagreement concentrates in the
expert residual.

**Report** (`ATLAS_cluster_projections_v30.pdf`, now 4 pages):
1. PCA + UMAP: GMM components, unsupervised EPIC, supervised EPIC, latent -- side by side.
2. latent-PGx: whole atlas, by unsupervised EPIC state, by GMM component.
3. the unsupervised-vs-supervised **confusion-matrix** heatmap (diagonal boxed).
4. the named-component table.

`coords.csv`, `atlas_v30_classified.tsv` and `latent_pgx_shortlist.tsv` now carry `epic_unsup`
and `epic_sup` (the curated columns are gone). Curated markers on the plots (the good/bad circles)
are retained purely as visual reference points; they no longer change any label.

### 9.10 v31 -- LOO removal audit (no feature removed); marker recolour

**Question:** does removing any single parameter better separate the unwanted confounder
`simvastatin x I83.1` from the true target ADRs (diazepam x J44.9, citalopram x M81.9,
lisinopril x R73.9, ramipril x R41.0, gabapentin x R42)?

**Method:** full leave-one-out over the 44 v30 features. Fast proxy = standardized centroid
distance/spread between the 3 simva rows and the 6 target rows; shortlisted removals then run
through the real GMM pipeline, scored on (a) does simva still share a GMM cluster with a target,
(b) `kappa_EPIC` and stability, (c) the E-vs-I separation control.

**Result: no removal helps.**

- The best single-feature raw separation gain was `effect_sign_concordance` (+0.13), but it
  collapses GMM k=10 -> 6 and halves `kappa_EPIC` (0.160 -> 0.052).
- Every other removal moved the simva-vs-target separation by < 0.03 (noise) and **every** removal
  dropped `kappa_EPIC` to 0.05-0.13 -- the v29 noise-reducing-ensemble result: these features earn
  their place jointly, not singly.
- In **every** configuration simva still shared a GMM cluster with at least one target. It cannot
  be isolated by removal because it genuinely carries the drug-triggered signature (rare,
  replicated at OR 26, no drug-free arm; sections 9.8-9.9). The confounding is a pair-level
  population fact, not a variant measurement.

| removed | simva-vs-target UMAP | shares cluster w/ target | kappa_EPIC | E-vs-I control |
|---|---|---|---|---|
| *baseline v30* | 2.26 | yes (C3) | **0.160** | 0.83 |
| effect_sign_concordance | - | yes | 0.052 (k=6) | 0.82 |
| drug_specificity | 2.12 | yes | 0.127 | 0.82 |
| gnomad_AF_missing | 2.88 | yes (C5) | 0.114 | 0.84 |
| ADR_pct_after | 2.09 | yes | 0.100 | 0.84 |
| (all others) | <=2.19 | yes | <=0.13 | ~0.83 |

**Decision:** by the project accept/revert rule, no feature is removed. **v31 is feature-identical
to v30** (44 features); the pipeline reproduces k=10, stable=19, kappa_EPIC=0.160 exactly.
`build_features_v31.py` exists so the LOO decision is versioned and reproducible.

**Also (also age-GWAS, from the prior investigation):** the `age_at_first_prescription` GWAS
(exists for 42 drug x ICD10 pairs incl. `simvastatin__I831`) was tested as a confounder axis --
the variant's age-association vs its ADR-association. It does **not** discriminate: every variant,
confounder and true target alike, shows age LOG10P approximately equal to disease LOG10P
(ratio approximately 1.0), a small-N artifact of 1-2 case carriers driving both signals. Rejected.

**Marker recolour (requested):** on all report pages the WANTED target ADRs are now **cyan**
(`#00BCD4`) and the UNWANTED confounders **orange** (`#FF7F0E`), replacing the black/gray circles.
Report published as `ATLAS_cluster_projections_v31.pdf`.

### 9.11 What actually distinguishes simva x I83.1 (web-search biology + correlation audit)

**Is the association true? No -- it is almost certainly confounded, and external biology says so:**

1. **The effect direction contradicts known pharmacology.** Statins are *protective* against
   varicose veins: atorvastatin/rosuvastatin suppress varicose remodelling by ~80% in mouse vein
   models, and human varicose-vein samples from chronically statin-treated patients show reduced
   smooth-muscle proliferation (Varicose Remodeling..., PubMed 26908399; AHA ATVB Abstract 304).
   Our association has OR ~26 in the *risk* direction -- the opposite of the established protective
   effect. A genuine drug-triggered signal would not reverse the drug's known venous pharmacology.
2. **The genes have no venous biology.** Large varicose-vein GWAS (139 risk loci; ~18% heritable)
   implicate EFEMP1, PIEZO1, SOX9, CASZ1, GJD3, TGFB2, KCNH8, SKAP2, ADAM15 (Nat Commun 2022,
   s41467-022-30765-y; Nat Cardiovasc Res 2022). Our three lead genes -- **ABHD14B, DOCK3, CHRM2**
   (a muscarinic acetylcholine receptor) -- are none of them.

**What distinguishes it is real but lives entirely in EXTERNAL knowledge**, not in any variant or
phenotype measurement we have or could compute internally:
- effect-direction-vs-known-pharmacology needs the fact that statins protect veins;
- gene-disease relevance: `OT_gene_disease_score` is ~0 for simva (0.00 / 0.038 / 0.00) -- but it
  is *also* ~0 for every true target (RYR2, ELOVL5, GJB6, PDE1C, ...), because these are novel PGx
  candidates absent from gene-disease databases. So OT cannot separate them either.
- the age-at-prescription GWAS (section 9.10) does not discriminate.

At the variant level simva x I83.1 is genuinely indistinguishable from a real rare replicated ADR.

**"Too many / correlated parameters dominating the clustering" -- structurally TRUE, but not the
cause, and fixing it hurts.** The 44-feature matrix has effective dimensionality **~10.2**
(participation ratio; top-10 PCs = 78% of variance) and **34 feature pairs at |r| >= 0.7**. A large
correlated block -- `BETA_regenie_abs`, `delta_BETA_ADR_vs_prescribed`, `SE_regenie_log`,
`OR_logdiff_test_minus_disc`, `repl_directional_z`, `test_Z_STAT` (|r| 0.85-0.99) -- all encode one
"big replicated effect" axis and collectively carry ~8x the weight of an independent feature. That
is exactly the axis simva is extreme on, so the geometry does group it with high-effect
associations.

But **decorrelating does not separate simva and degrades the model** (PCA-whitening tested):

| transform | k | kappa_EPIC | simva vs targets | sep_UMAP | E-vs-I control |
|---|---|---|---|---|---|
| v31 raw standardized (44) | 10 | **0.160** | share C3 | 2.26 | **0.83** |
| PCA-whiten 10 | 9 | 0.125 | share C1 | 1.57 | 0.57 |
| PCA-whiten 12 | 10 | 0.106 | share C2 | 0.83 | 0.89 |
| PCA-whiten 15 | 10 | 0.109 | share C1 | 1.52 | 0.55 |
| PCA-whiten 44 (full) | 10 | 0.011 | share C1 | 1.19 | 0.19 |

Whitening down-weights the effect block but simva still shares a cluster with the targets in every
case, and `kappa_EPIC`, the good/bad separation and the E-vs-I control all fall. The redundancy is
the v29 noise-reducing ensemble, not dead weight; the raw standardized space is better. simva's
co-clustering is driven by genuine similarity to the targets, not by a weighting artifact.

**Conclusion.** Confounding here is a pair-level pharmacological fact (statins protect veins) with
no variant-level footprint. No parameter -- added, removed, decorrelated, or age-based -- isolates
it. The only faithful routes are the labelled curated overlay (annotation, not feature) or a new
*external* pair-level descriptor (e.g. a literature-derived expected direction of the drug's effect
on the condition), which would sit outside the "variant and phenotype descriptions only" rule.

### 9.12 simva x I83.1 assigned C via an external known-confounder annotation (age / reverse-pharmacology)

Sections 9.8-9.11 established, by exhaustive test, that NO variant- or phenotype-level parameter
-- added, removed, decorrelated, whitened, or age-GWAS-based -- can separate `simvastatin x I83.1`
from a genuine replicated ADR, because at the variant level it IS one. The confounding is a
population-level biological fact with no variant footprint:
  * varicose veins is a degenerative venous condition of older adults and statins are an
    elderly-population drug, so AGE is the common cause of both the prescription and the diagnosis;
  * statins are moreover KNOWN to be protective against varicose remodelling (PubMed 26908399), so
    the OR~26 *risk* signal contradicts the drug's pharmacology -- it cannot be an honest ADR.

Decision: assign it C through an explicit **external annotation**, not through the features.
`epic_from_parameters.KNOWN_CONFOUNDERS` lists the confounded pairs with a per-pair biological
reason. The pipeline now emits, per association:
  * `epic_unsup` -- the unsupervised model (GMM cluster -> EPIC), LEFT UNTOUCHED;
  * `epic_final` -- `epic_unsup` with the known confounders forced to **C**;
  * `epic_source` -- `model` or `confounder`;
  * `confounder_reason` -- the biological justification.

This keeps the clustering honestly unsupervised (nothing enters the feature matrix), keeps the
disagreement auditable (`epic_unsup` still shows what the data said), and gives the correct final
classification. simva x I83.1: `gmm=C3`, `epic_unsup=E`, **`epic_final=C`**. Four rows across two
pairs are overridden (simva x I83.1 x3, fluoxetine x N95.1 x1); `epic_final` counts
E 104 / I 52 / P 24 / **C 27** (was C 24 unsupervised).

Report `ATLAS_cluster_projections_v31.pdf` page 1 now shows three EPIC columns side by side --
UNSUPERVISED (GMM->EPIC), FINAL (+ known-confounder overlay), and SUPERVISED (expert) -- so the
override is explicit. Page 2's latent-by-state boxplot uses `epic_final`, so simva sits in the C
box. `coords.csv`, `atlas_v31_classified.tsv` and `latent_pgx_shortlist.tsv` all carry
`epic_unsup / epic_final / epic_source`.

Note vs section 9.9: this is the same overlay mechanism removed earlier, now reinstated with an
explicit age/reverse-pharmacology rationale and applied on top of the unsupervised model rather
than the per-row parameter arg-max. It is a labelled annotation layer, never a clustering feature.

### 9.13 v32 -- five-state scheme (E/P/I/C/C2), arbitrariness audits, de-arbitrarised latent

**Five states** (replacing the 4-state EPIC). Each is a parameter signature; each axis is a
coefficient-1 sum of z-scored measurements, standardized to unit variance, arg-max wins.

| state | meaning | signature |
|---|---|---|
| **I** | Indication | on-label AND condition precedes Rx: `z(on_BNF)+z(pct_before)` |
| **E** | drug-triggered Effect | drug-specific, follows exposure, replicates, NO disease-alone: `z(beta_sup)+z(ratio)+z(repl_z)-z(dis)` |
| **P** | Preexisting, enhanced | variant causes condition AND drug amplifies: `z(dis)+z(beta_sup)` |
| **C2** | Confounder | locus drives disease/Rx, drug adds nothing: `z(dis)+z(pleio)+z(rx)-z(beta_sup)` |
| **C** | Comorbidity | weak nonspecific residual: `-z(beta_sup)-z(dis)-z(pleio)` |

Splits: E vs P on disease-alone (dis); P vs C2 on beta_superiority (drug enhances vs adds nothing);
C2 vs C on locus strength. Anchors validate: APOE, LPA/9p21, Factor-V-Leiden -> C2; on-label -> I;
replicated ADR -> E; weak targets (citalopram, lisinopril, ramipril) -> C (honestly nonspecific).

**Arbitrariness audit of the axes.** Under +-25% random perturbation of every term weight
(200 draws), **94.7% of labels are unchanged (min 88%)**. The scheme is robust, not arbitrary.
The earlier `-2*beta_sup` coefficient (which pulled weak ADRs into the confounder state) was
removed in favour of coefficient-1 terms.

**Arbitrariness audit of the latent score -- two hand-tuned elements removed.**
- the `repl_z=0.7` weight was immaterial: E-vs-I AUC is flat (0.866-0.887) across repl_w in [0,1],
  and equal weights scored marginally higher (0.887 vs 0.881, rank-corr 0.99). Dropped -> all
  weights 1.0.
- the FAERS x temporal PAIR axis HURT the metric: removing it raised AUC to 0.896. It had been
  added to prop up ramipril/gabapentin, which have no variant-level signal; that lift was the
  artefact. Dropped.
Result: **equal-weight geometric mean of five axes, no free parameters, AUC 0.881 -> 0.890** (raw);
published score (with gate/downweight) 0.909 -> 0.901. Random +-40% weight perturbation keeps
rank-corr >= 0.94. Honest cost: ramipril (71->131) and gabapentin (75->165) fall, correct given
they are variant-level-underpowered; diazepam rises to rank 6.

**Cluster -> 5-state map (GMM k=10):** C8->C2 (APOE), C9->C2 (LPA/9p21), C10->P (Factor V Leiden),
C2->I / C4->I, C3->E / C5->E, C1/C6/C7->C. Unsupervised balance E 88 / C 59 / I 52 / C2 4 / P 4;
the confounder and predisposition states are minority mechanisms that only small components
isolate at k=10 -- honest, and the per-row arg-max (`epic_param`) resolves more (C2 11 / P 43).
`epic_final` additionally forces the two external known-confounder pairs to C2 (simva x I83.1,
fluoxetine x N95.1), so simva reads **C2** in the final labelling.

**Features unchanged.** v32 == v31's 44-feature matrix (GMM k=10, kappa 0.160); `build_features_v32.py`
is a documented copy -- the v32 changes are labelling + score only. Published as
`ATLAS_cluster_projections_v32.pdf`, `atlas_v32_classified.tsv`, `cluster_coords_v32.csv`.

### 9.14 Naming collision fixed: GMM components renamed C1-C10 -> G1-G10

The v32 five-state scheme introduced **C2 = Confounder**, which collided with the GMM component
ids, which had been `C1..C10` ("C" for cluster). On the same figure `C2` therefore read as
"component 2" in the component panel and as "confounder" in the state panels -- and component C2 is
in fact the *On-label indication* cluster (state **I**), i.e. the two readings were contradictory.

Fix: component ids now use a **G** prefix (`G1..G10`, G = GMM component) in
`export_plots_data.py` and `gmm_cluster_names.py`. The EPIC state letters are unchanged
(I / P / E / C / C2), since those are the user-defined scheme. Verified: the cluster-id set and the
state set are now disjoint.

Also fixed in the same pass: `coords["gmm_state"]` was still being filled from the hardcoded
4-state `GMM_NAMES` tuple and is now taken from the live 5-state `cluster_epic_map`; the
`gmm_cluster_names` docstring references were renumbered to G*.

The report's page-1 subtitle now spells the states out in full
("I=Indication, P=Preexisting enhanced by drug, E=drug-triggered Effect, C=Comorbidity,
C2=Confounder. Components are G1-G10 (not states).") so the figure is self-explanatory.

### 9.15 v33 -- fifteen interpretable dimensions and the Causal Specificity Index (CSI)

The 44 flat model features are good for clustering but hard to reason about: an association is
"interesting" for several distinct REASONS, entangled across the columns. v33 regroups the same
measurements into **15 named dimensions**, each a coefficient-1 mean of z-scored terms standardized
to unit variance (the v32 axis construction), and combines them into one signed composite.

**Scope:** self-contained + variant annotation. Gene-level external knowledge (Tissue, Druggability,
DGIdb, OpenTargets, external PheWAS, GO) remains excluded.

| group | dimensions |
|---|---|
| **Rewarding (11)** | Robustness, Replicability, DrugSpecificity, Temporality, DoseDependence, Pharmacovigilance, EffectMagnitude, FunctionalImpact, DirectionalCoherence, Power, DataQuality |
| **Penalising (3)** | Indication, Predisposition, Confounding |
| **Neutral (1)** | Rarity -- reported, never scored |

**CSI** = `exp( (1/14) * [ sum_rewarding log sigmoid(z(D)) + sum_penalising log sigmoid(-z(D)) ] )`
-- signed soft-AND (Nichols 2005), equal weight, no free parameters. Rarity is excluded from the
product because the whole atlas is rare; scoring it would reward underpowered noise.

**Four new source columns merged** (`ATLAS_FULL_COLS`): `FAERS_a`, `INFO`, `N_regenie`,
`A1_CTRL_CT` -- all numeric, <=5% missing, previously unreachable. Total parameters consumed by the
dimensions: ~44, including six that no builder had ever read (`SE_ICD10`, `SE_dose`,
`P_test_PLINK2`, `TEMP_n_total`, `locus_multiplicity`, `beta_disc_valid_consistency`).

**Acceptance checks -- all pass:**
1. All 15 dimensions finite and unit-variance; smallest distinct count 13 (FunctionalImpact), so
   none is near-constant.
2. **No dimension pair reaches |rho| >= 0.85** -- the dimensions are non-redundant despite being
   defined conceptually rather than orthogonalised.
3. **Sanity anchors reproduce known biology unprompted:** `simvastatin x G30.9` (APOE)
   Confounding **7.84**, CSI rank 205/207; `simvastatin x I25.1` (LPA/9p21) Confounding **7.62**,
   rank 198; `simvastatin x I80.2` (Factor V Leiden) Predisposition **4.65**;
   `simvastatin x E78.5` Indication **2.81**. None of this was told to the model.
4. **Robustness:** +-25% perturbation of every dimension weight over 200 draws leaves CSI
   rank-correlation at **0.991** (min 0.983) vs the 0.9 target.
5. Clustering untouched: GMM k=10, stable 19, kappa_EPIC 0.160.

**The honest negative result.** CSI does **not** separate the curated wanted from unwanted rows.
`simvastatin x I83.1` ranks **25 / 207** -- above every target except one -- while the wanted median
rank is 94 vs an unwanted median of 78. The dimension profile shows exactly why, and it is the same
finding as sections 9.8-9.11, now legible for the first time:

| dimension (mean) | unwanted | wanted |
|---|---|---|
| Replicability | **+0.71** | -0.18 |
| DrugSpecificity | **+1.14** | -0.04 |
| EffectMagnitude | **+0.64** | -0.17 |
| DirectionalCoherence | -1.24 | **+0.29** |

The curated targets are genuinely weak on the evidence dimensions -- ramipril x R41.0 scores
Replicability **-2.4**, gabapentin x R42 **-1.4** -- while simva x I83.1 scores **+1.2 to +1.6**.
CSI is measuring the data correctly; the target list was chosen on grounds the variant data does not
support. This is reported, not tuned away. `fluoxetine x N95.1` IS caught (rank 202), but by
Robustness (-4.2), not by the Indication dimension.

Consequence: the `KNOWN_CONFOUNDERS` external annotation remains the only mechanism that demotes
`simvastatin x I83.1`, exactly as predicted in the v33 plan's risk section.

**Artifacts:** `dimensions.py`, `plot_dimensions.R`, `plots/dimensions.csv`,
`plots/dimension_correlation.csv`, `ATLAS_dimension_report_v33.pdf` (3 pages: CSI scale with the 10
curated rows; a 10 x 15 dimension heatmap explaining each score; the correlation matrix and CSI vs
latent-PGx). Also fixed: `TEMP_pct_missing_rx_date` is **constant 0** across all 207 rows, so the
existing `temporal_pct_missing_dates` feature was only ever the dx-date term -- the DataQuality
dimension uses the dx term alone. Dead `DRUGGABILITY_ORD` / `PHEWAS_ORD` imports and a stale
"Option-B tissue features" log line removed from `build_features_v32.py`.

### 9.16 Literature review of the curated targets: only one survives

Web-literature assessment of the 6 curated targets and the 15 highest-CSI associations, judged on
(1) does the drug pharmacologically cause the event, (2) is the direction right, (3) is the gene
plausible. Full write-up: `target_assessment.md`.

**The five curated pairs:**

| association | verdict |
|---|---|
| gabapentin x R42 (ELOVL5) | **TRUE TARGET** -- dizziness/ataxia is gabapentin's hallmark ADR and ELOVL5 mutations cause spinocerebellar ataxia 38. Drug mechanism and gene mechanism converge on the cerebellum. Best hypothesis in the atlas. |
| citalopram x M81.9 (DLEU1) | **TRUE ADR, gene implausible** -- SSRIs: 1.61x fracture, 1.6x bone loss, serotonin transporters on osteoblasts/osteoclasts. But depression itself lowers BMD (partial confounding) and DLEU1 has no bone biology. |
| ramipril x R41.0 (GJB6) | **real ADR, artifact variant** -- ACE-inhibitor SIADH hyponatremia -> confusion is documented, but this variant has 0 test carriers, OR 0.99, Replicability -2.4, and GJB6 is a deafness gene. |
| diazepam x J44.9 (CXCL13/PDE1C) | **LIKELY COMORBID** -- benzodiazepines worsen COPD (+45% exacerbations, higher mortality) but do not *cause* COPD, which is a smoking disease. An incident J44.9 code after diazepam is exacerbation-coding or confounding by indication. |
| lisinopril x R73.9 (RYR2) | **CONFOUNDER** -- see below. |

**A second reverse-pharmacology confounder found.** ACE inhibitors *improve* insulin sensitivity,
*lower* glucose and *reduce* new-onset diabetes (they can cause hypoglycaemia). Hyperglycaemia as a
lisinopril adverse effect is therefore **directionally backwards** -- the identical signature to
`simvastatin x I83.1`. Lisinopril is prescribed for diabetic nephropathy, so its users are enriched
for diabetes; 11.5% of the hyperglycaemia diagnoses precede the prescription.
**`lisinopril x R73` is a candidate for `KNOWN_CONFOUNDERS`.**

**CSI's own top hits are more defensible than the curated list.** `simvastatin x G62.9`
polyneuropathy (rank 1) is a contested but genuine pharmacovigilance question (4-14x in
case-control studies, null in a 2021 meta-analysis, ~12/100k person-years attributable);
`fluoxetine x R55` syncope (rank 7) is mechanistically clean (SSRI SIADH hyponatremia ->
unsteadiness, falls, syncope, worst in the first 1-3 weeks in the elderly). Neither was on the
curated list. Ranks 2-5, 8-10 and 13-15 are non-specific R/A codes with no documented mechanism --
the expected background of a rare-variant scan.

**A methodological suspicion, tested and REFUTED.** I suspected CSI's top was inflated by tiny event
counts (0% pre-exposure is trivial with 2 events). It is not: the top-20 rests on *more* evidence
than the atlas median (176 vs 126 temporal events; 2.0 vs 1.0 case carriers; none with <=5 events),
and CSI correlates only weakly with event count (rho +0.08) and carriers (rho +0.21). Recorded
because it was tested, not assumed.

**Consequence.** The curated target list does not survive literature scrutiny -- which independently
corroborates §9.15. CSI ranked those pairs mid-table because their evidence really is mid-table, not
because the index is miscalibrated. The generalisable pattern for confounders in this atlas is
**reverse pharmacology**: the drug is known to protect against, or is prescribed because of, the
event. That is a pair-level literature fact, not a variant statistic, which is precisely why no
model parameter can capture it (§9.8-9.11).