Add the naive per-condition baseline arm and manuscript Figure 3

[jaredgalloway]

Closes #294

EPIC #290 Phase 4. Builds the naive per-condition baseline the manuscript
contrasts against joint fitting, and renders manuscript Figure 3 (six
panels, A–F) from it.

The naive approach is the obvious alternative to joint modelling: fit each
condition alone, then call the difference of the two independently inferred
effects a "shift". This PR makes that baseline reproducible inside the
pipeline so the two arms can be compared on equal footing, and — this turned
out to be the substantive part of the work — measures whether the resulting
comparison actually supports the claim the manuscript makes with it.

What landed

Piece	Where
spike.naive config block (times_seen_threshold: 1, three conditions)	all six config/*_downstream.yaml
fit_single_condition, fit_naive_arm, derive_naive_shifts, assert_wt_agreement	notebooks/_downstream.py
8 unit tests pinning the arm's load-bearing invariants	tests/test_naive_baseline.py
naive_baseline.ipynb — six single-condition fits, shift derivation, α sensitivity check	notebooks/
rule naive_baseline + fig3 registry entry + manuscript_figures inputs	Snakefile
Figure 3, panels A–F, one mosaic / one savefig, plus the confound diagnostics	notebooks/manuscript_figures.ipynb
Docs toctree entry	docs/spike_naive_baseline.nblink, docs/index.rst
S16 marks the non-converged fit instead of dropping it (deviation 7)	notebooks/manuscript_figures.ipynb, Snakefile

Three design choices make the comparison fair, and all three are deliberate:

1. Equal convergence. The naive arm inherits the joint arm's tol and
   maxiter from the same config. The original analysis gave the naive fits a
   tenth of the iteration budget, so part of its poor showing was budget
   rather than method. All 6 fits converge (max final error 9.4e-07 against
   tol 1e-06), and manuscript_figures prints that rather than assuming it.
2. One shared mutation index at times_seen_threshold=1. Every scatter
   panel reads the same 5,428 mutations. This matches evaluate.ipynb's
   joint-arm call; filtering the arms differently rigs the comparison —
   joint shift_Delta reads 0.543 at threshold 3 but 0.405 at 1.
3. The intersection premise is verified, not assumed. 5,428 mutations
   checked against a site map carrying 35 sites whose wildtype letter
   differs by background; none reach the shared index.

No refit. config/config.yaml is untouched (git diff main -- <it> is
empty), so the 1.76 GB fit_collection.pkl and the ~2h20m fit are not
invalidated. The DAG is naive_baseline → evaluate → manuscript_figures;
fit_models is absent. evaluate reruns (~7 min) because it legitimately
declares config_downstream.yaml as an input.

Measured results

Against results-prod-292-spike-fit-tuning, λ = 8e-05, shared index
n = 5,428, times_seen ≥ 1.

Panel	Quantity	Measured	Published
C	joint β (BA.1)	0.833	0.84
D	joint shift Delta / BA.2	0.405 / 0.378	0.48 / 0.54
E	naive β (BA.1)	0.843	0.78
F	naive shift Delta / BA.2	0.349 / 0.169	0.14 / 0.16

Deviations from spec

1. The headline statistic is confounded — for both arms

The spec framed the deliverable as reproducing the published gap (naive 0.14
vs joint 0.48, "3.4× worse") and expected the equal-convergence correction to
narrow it. It does narrow it — but the narrowed number rests on the same
broken statistic as the original.

A naive shift is β_cond − β_ref. Both replicates subtract the same
well-replicating reference beta vector, and that shared term manufactures
correlation on its own. Permuting the condition betas identically in both
replicates — destroying every condition-specific signal while leaving β_ref
paired — gives:

	observed R²	permutation null
naive Delta	0.349	0.697
naive BA.2	0.169	0.811

The null is higher than what is observed. Replicate R² of a
difference-of-betas is not a clean measure of shift reproducibility. So do
not report "3.4× collapsed to 1.2×" as the finding.

Residualizing each arm's shift on β_ref isolates the real signal:

	joint	naive	ratio
Delta	0.341	0.155	2.20×
BA.2	0.391	0.186	2.10×

Joint beats naive by a consistent ~2.1–2.2× in both conditions — unlike the
published asymmetry. This is the defensible statement, and the notebook
computes it next to the headline numbers rather than asserting it.

2. The advantage lives entirely in low-coverage mutations

Stratified by min times_seen across both replicates and all conditions:

                    Delta                  BA.2
 threshold    joint  naive  ratio    joint  naive  ratio      n
 ts>=1        0.405  0.349  1.16x    0.378  0.169  2.23x   5,428
 ts>=3        0.529  0.559  0.95x    0.478  0.378  1.27x   4,058
 ts>=5        0.641  0.673  0.95x    0.564  0.502  1.12x   3,059
 ts>=10       0.752  0.793  0.95x    0.684  0.673  1.02x   1,736
 ts>=20       0.843  0.848  0.99x    0.800  0.801  1.00x     642

The joint advantage vanishes by ts >= 10; for Delta it inverts at
ts >= 3. The defensible story is "joint fitting helps where data is thin"
— which is the manuscript's own thesis, but a narrower claim than "3.4×
better replication", and it should be stated as the narrower one.

Supporting evidence: naive β_Delta replicate R² is 0.461 against joint's
0.718, while BA.1 and BA.2 naive betas are fine (0.843, 0.784). Delta has
the least data (28k variants vs 60–69k; median times_seen 8 vs 16), and
independent fitting cannot borrow strength.

3. ⚠ The spec's alpha lower-bound argument is refuted

The issue body argued that per-condition α mismatch is replicate-stable,
therefore adds a correlated component to both replicates' naive shifts,
therefore inflates naive R² — making the measured gap a lower bound on
the true gap.

Measured, this runs backwards. naive_baseline.ipynb now rescales each
condition's betas by its fitted α before subtracting and prints the result:

  Delta          uncorrected 0.340   rescaled 0.416   -> alpha mismatch is suppressing
  Omicron_BA2    uncorrected 0.161   rescaled 0.179   -> alpha mismatch is suppressing

α mismatch is suppressing the naive correlation, not inflating it, and the
gap is not a lower bound. α is confirmed replicate-stable (Delta
11.72/11.41, BA.1 13.06/13.14, BA.2 12.47/11.79) — the premise holds; the
inference from it does not.

Issue #294's body carries a ⚠ CORRECTION block, and the notebook's prose has
been rewritten to state the measured result rather than the refuted argument
— it publishes as method documentation via the docs toctree, so shipping the
false version was not an option. (These two numbers are on the naive arm's
own paired index, 0.340/0.161, not the joint-intersected 5,428 used by the
figure; the direction is identical either way.)

This needs a call on whether the manuscript text carries the same
argument.

4. The joint number is materially λ-dependent

fusionreg   shift_Delta  shift_BA2   both-zero pairs
0.0            0.431       0.190        0 /    0
2e-05          0.390       0.263       60 /  369
8e-05          0.405       0.378      854 / 1671   <- chosen
1.6e-04        0.391       0.439     2855 / 2405
3.2e-04        0.248       0.530     4927 / 4489

Delta's joint R² peaks at λ = 0 (0.431), above the chosen 8e-05. BA.2
climbs monotonically to 3.2e-04, but by then 4,489 of 5,428 shifts are zero
in both replicates. The two conditions want opposite λ. 8e-05 was not
chosen to maximize shift replicate R², so reporting one joint number without
the curve overstates its stability. (Excluding exactly-(0,0) pairs moves
little: Delta 0.405 → 0.392, BA.2 0.378 → 0.376.)

5. Figure 3 reproduces the legacy layout verbatim

Figure 3 was first drawn with derived axis limits and a bespoke panel
layout. It is now a verbatim port of the legacy figure
(matsengrp/SARS-CoV-2_spike_multidms@6c98b7b, spike-analysis.ipynb
cells 84–95), matching the convention Figure 4's heatmap already follows:
legacy structure as written, only the data plumbing adapted to this
pipeline's CSV schema. The mosaic, the three-panel correlation rows, the
hard-coded limits and ticks, the panel letters and titles are the legacy
code.

Three departures remain, each commented at its site:

• One shared mutation index across both arms. Legacy called
  combine_replicate_muts(how="inner") per arm, leaving the joint and
  naive panels on different mutation sets. Every panel here reads the
  same 5,428 mutations — filtering the arms differently rigs the
  comparison the figure exists to make.

• The clip is assigned, not passed inplace. Legacy's
  mut_df_replicates["avg_beta"].clip(lower=-5, inplace=True) is a no-op
  under modern pandas copy-on-write, so the clip silently never happened.

• The distribution block is 2 rows, not the legacy 3. The legacy
  middle row was in-frame codon deletion muts (muts == '-'). This
  pipeline deletes every deletion-bearing variant during data prep
  (prepare_data.ipynb, "Remove deletions and invalid variants" — on
  main since Improve spike data prep: post_count threshold, nonsense truncation, log2 scores #226), so no deletion mutation can reach the model and
  that row rendered as an empty labelled panel. It is dropped.

  Two classes remain because two are all the model has: β is indexed by
  observed amino-acid substitution, and synonymous variants carry no
  substitution — prep folds them into the wildtype row used for
  func-score subtraction, so they have no β to plot. Showing synonymous
  mutations in this block would require both a prep change (retain the
  class) and a full refit; that is a modelling decision, not a figure
  one, and is left to the epic.

  Re-rendered against results-prod-294-naive-baseline-arm: all six
  replicate R² unchanged, confirming the change touches only the
  distribution block.

• Axis limits and bins refit to this pipeline's data. The legacy
  ranges were sized for the legacy fit's much wider parameter spread.
  Carried over verbatim they left every subpanel mostly empty and
  rendered E and F as small off-centre blobs:

  panel	legacy limits	data range	wasted
  C β	[-6, 3]	−1.97 … 1.11	~66%
  D shifts	[-2.8, 2.8]	−1.48 … 1.83	~41%
  E naive β	[-11, 6]	−2.25 … 1.25	~79%
  F naive shifts	[-8.2, 12.2]	−2.74 … 3.59	~69%

  Limits are now the observed range across both replicates padded ~5% and
  rounded to clean ticks, kept square so the identity line stays a
  true diagonal, and shared across the two columns of D and of F so
  those pairs remain directly comparable.

  Distribution bins likewise: the legacy ranges gave only ~5 populated
  0.5-wide bins per panel. Halving the width to 0.25 over the observed
  range resolves the distributions. Verified no mutation falls outside
  the new bins, and the tallest bin is 0.97 — so the legacy [0, 1]
  probability axis still holds.

  Noted at its site: the naive clip(lower=-10) does not bind on this
  fit (β bottoms out at −2.25). Retained as a guard rather than removed,
  since it would matter on a fit where β does run away.

  All six replicate R² unchanged — presentation only.

6. The prod fit was re-run remotely rather than reused from cache

The plan anticipated shipping a cached checkpoint to a remote host, and an
earlier draft of this section reported that no refit was needed. That is
no longer accurate.

The local results-prod-292-spike-fit-tuning cache was destroyed during
this work by a snakemake --forcerun manuscript_figures invocation run
without a preceding --dry-run: config.yaml is declared as a rule
input:, the Phase 2 PR (32d7ffb) had added a comment block to it, and
Snakemake hashes whole files — so it planned a 7-rule rerun and deleted
fit_collection.pkl and cv_fit_collection.pkl at job start.

The remote copies on orca04 were unaffected, and the arm was re-run there
in a fresh worktree (results-prod-294-naive-baseline-arm, commit
147e057): 7/7 rules, 51 minutes, fit_models alone taking 41. The
regenerated training_functional_scores.csv is byte-identical to the Aug
12 version (17,210,919 B) and both pickles match their originals to within
3 bytes of pickle metadata, confirming the config change was inert and the
numbers in this PR are unchanged.

The measured results below were computed against the original cache and
have been re-verified against the re-run.

7. Figure S16 no longer hides the one fit that did not converge

Scope beyond #294, folded in at the user's request after the dashboard and
S16 appeared to contradict each other: fit_convergence.csv reports 19/20 fits
converged, but S16 showed every line dropping below tol.

Both artifacts were correct. S16 called drop_excluded() on the trajectory
frame, removing the 1.28e-3 rung — which contains the study's only
non-converged fit (rep_1, 500-sweep cap, drift_frac 0.10). A convergence
figure that deletes the only non-converged trace asserts the opposite of the
truth.

config_downstream.yaml already stated the correct intent — figures that plot
"all rungs" should "mark it rather than silently showing it as a normal point"
— so this is the code catching up to the documented contract, not new policy.

S16 now iterates the full ladder and marks the excluded rung grey/dashed, with
the non-converged trace annotated in place. Two judgement calls:

• The mark is driven by fit_convergence.csv, not by fusionreg. The
  instability is per-fit, not per-rung: rep_2 converges at 1.28e-3 in 109
  sweeps. Greying its trace as "did not converge" would swap one dishonesty for
  another, so the legend reads 0.00128 (excluded) and only rep_1 is
  annotated did not converge.
• Colours are indexed off the 9-rung analysed ladder, so the existing
  palette is unchanged and the marked trace takes no colour slot. The other
  nine lines are pixel-identical to the previous version.

fit_convergence.csv is added to the rule's input:. It is an evaluate
output in the same downstream tier as convergence_trajectory.csv, which the
rule already consumed, so this adds no new dependency edge and cannot reach the
fit tier.

S9 keeps dropping the rung via drop_excluded(), which is correct there:
its panels are the λ-selection criteria, only defined for converged fits, on a
categorical axis over the analysed ladder. Marking S9's slot needs rung_x()
to gain a tenth position and is deliberately out of scope. drop_excluded()'s
docstring now records why the two figures differ, since the config's
"(S9, S16)" wording elides it.

Everything else: no deviation

D1–D7 landed as written, including validations V1–V3.

Checked, and cleared as non-artifacts

• Index subsetting — joint Delta 0.4097 on its own 5,809 → 0.4049 on the
  shared index. Negligible.
• Outliers — Spearman 0.307; trimming |shift| > 2 gives 0.370.
• Stop codons — R² 0.001–0.011 over n = 137. Dropping them changes nothing.
• Naive shift arithmetic — reproduces β_cond − β_ref to 7e-16; the
  reference condition's shift is identically 0; no clipping at beta bounds.

Split out of this PR

In-frame codon deletions are dropped at the prep stage and have been since
baaa0dc, the first commit of the spike pipeline. The legacy analysis kept
them. Restoring them would add 881 mutations (354 with ≥5-variant
support, led by NTD antigenic-supersite residues P25-, N164-, F140-),
and would require re-deriving λ, since shift_sparsity uses stop codons as
its false-positive denominator and cannot see deletions.

That is out of scope here and is filed as #313. A working restoration edit
is parked as git stash 8d911b5 and referenced from that issue.

Open question for Phase 7

Under-convergence alone does not explain the published naive 0.14. An
under-converged Delta fit shrinks β_Delta toward its init, pushing
β_Delta − β_BA1 toward −β_BA1, which per the decomposition pushes naive
R² up toward the 0.70 null — not down to 0.14. The legacy 0.14 more
likely came from a different index, a different times_seen threshold, or
per-replicate rather than shared-index pairing. Resolving it needs the legacy
naive betas, not just its reported correlation.

Verification

• pixi run pytest --doctest-modules multidms tests → 142 passed (8 in
  test_naive_baseline.py)
• pixi run ruff check . → clean; explicit run on _downstream.py → clean
  (ruff's config excludes experiments/)
• pixi run black --check → clean, 13 files
• pixi run docs → build succeeded, spike_naive_baseline present
• Test profile end-to-end → snakemake --config profile=test --forceall,
  7/7 rules, Figure 3 renders. The naive fits report 0/6 converged there
  by design (2 outer sweeps) — reported, not asserted, per spec.
• Prod run → naive_baseline, evaluate, manuscript_figures;
  fit_models absent from the DAG
• git diff main -- experiments/scv2-spike/config/config.yaml → empty
• S16 re-render (deviation 7) — rendered on orca01 at a4a7536.
  --dry-run gate first: DAG was exactly manuscript_figures + all, reason
  "Updated input files: manuscript_figures.ipynb"; no fit_models,
  cross_validation, evaluate, or naive_baseline. Run took 50 s, 30/30
  cells, zero errors. Post-run parity: all 7 fit-tier CSVs and all 20 figures
  md5-identical local↔remote; both pickles retain their Aug 16 16:39/16:40
  mtimes, proving the fit was never re-run.

A fresh-context reviewer audited the completed diff and found five real
defects, all fixed in 6e71a4c: the shared-index violation (panels C/D and
E/F were on different mutation sets despite the prose claiming otherwise),
the per-replicate union inflating every reported n from 5,428 to 6,386,
V3 shipped as dead code, naive_convergence.csv declared as an input but
never read, and panel B silently dropping the stop class.

Yolo trail: spec gate ✓, plan gate ✓ — full log in _ignore/yolo/294.log

🤖 Generated with Claude Code