Figure 5 — Reparam (nokappa-v3) version¶

Mirror of Figure4_Predictive_Performance.ipynb with reparam (non-centered) predictions instead of centered (nolr).

Panels¶

Panel Status Source
A. Performance bars (multi-disease, multi-method) ✓ already rendered talk_figs/performance_comparison_nokappa_mar4.pdf (sourced from claudefile/results_feb18_complete/)
B. ROC curves (ASCVD, 1-yr at offsets 0–9) rendered below Dropbox/age_offset_nokappa_v3/pi_nokappa_v3_age_offset_*.pt
C. Washout (1mo/3mo/6mo, held-out LOO) rendered below Dropbox/washout_comparison_reparam_loo_10k/
D. Calibration (predicted vs observed at-risk event rates) rendered below Dropbox/enrollment_predictions_nokappa_v3_loo_all40/pi_enroll_fixedphi_sex_FULL.pt

All outputs go to results/paper_figs/fig5_reparam/ so the original fig5/ is untouched.

Panel A — Performance bars (multi-disease)¶

Already rendered from feb18 reparam metrics:

  • PDF: /Users/sarahurbut/aladynoulli2/talk_figs/performance_comparison_nokappa_mar4.pdf
  • CSV: /Users/sarahurbut/aladynoulli2/talk_figs/performance_summary_table_nokappa_mar4_detailed.csv

Three Aladynoulli columns: 1-yr (Baseline), 1-yr (Median), 10-yr (Best of static/dynamic). Numbers verified against claudefile/results_feb18_complete/summary_all_metrics.csv.

Panel B — ROC curves for ASCVD (1-yr at age-offsets 0–9, reparam)¶

Calls the parallel script generate_age_offset_predictions_reparam.py which loads reparam pi tensors from Dropbox/age_offset_nokappa_v3/.

In [ ]:
%run /Users/sarahurbut/aladynoulli2/pyScripts/dec_6_revision/new_notebooks/pythonscripts/generate_age_offset_predictions_reparam.py --max_offset 9 --start_idx 0 --end_idx 10000

Panel C — Washout analysis (held-out LOO reparam, batch 0 / 10K patients)¶

Methodological improvement over the published centered Figure 5C. The published version trained and tested on the same 10K patients (master phi/psi was pooled from all 400K, which includes the 10K). Here we use proper LOO: phi/psi/gamma/kappa are pooled from batches 1–39 only, then delta is refit on batch 0 with each washout E matrix. Train and test are held out.

Pipeline (already executed; pi tensors saved):

  1. claudefile/run_loo_washout_reparam.py — produces pi_washout_{no_washout,1month,3month,6month}_0_10000.pt in Dropbox/washout_comparison_reparam_loo_10k/
  2. pyScripts/dec_6_revision/new_notebooks/pythonscripts/evaluate_washout_predictions.py — produces 1-yr and 10-yr AUC CSVs + washout_comparison_1yr_10yr.csv in results/washout_evaluation_reparam/
  3. plot_washout_results_reparam.py — renders the bar plot below to results/paper_figs/fig5_reparam/

Result vs centered (head-to-head, batch 0): same bar-shape pattern (small 1-yr drop, ~zero 10-yr drop). Reparam absolute AUCs are ~5–10 points higher across the board. Breast Cancer is the biggest 1-yr drop in both (centered +0.058, reparam +0.088). Heart Failure, Colorectal, and All Cancers essentially zero drop in both.

In [ ]:
%run /Users/sarahurbut/aladynoulli2/pyScripts/dec_6_revision/new_notebooks/pythonscripts/plot_washout_results_reparam.py

Panel D — Calibration (full 400k, reparam)¶

Same code as Figure4_Predictive_Performance.ipynb, only the pi-tensor path is swapped from enrollment_predictions_fixedphi_fixedgk_nolr_vectorized/ (centered) to enrollment_predictions_nokappa_v3_loo_all40/ (reparam).

In [ ]:
# ============================================================================
# Reparam Calibration: 400k patients (subsetted to 50k for plot tractability,
# matching the centered-version cell)
# ============================================================================
import torch
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path

# REPARAM pi tensor (replaces the centered nolr path)
REPARAM_PI = '/Users/sarahurbut/Library/CloudStorage/Dropbox/enrollment_predictions_nokappa_v3_loo_all40/pi_enroll_fixedphi_sex_FULL.pt'

pi_full = torch.load(REPARAM_PI)[:50000]
print(f"Loaded reparam pi: {pi_full.shape}")

Y_full = torch.load('/Users/sarahurbut/Library/CloudStorage/Dropbox-Personal/data_for_running/Y_tensor.pt',
                    map_location='cpu', weights_only=False)[:50000]
E_corrected_full = torch.load('/Users/sarahurbut/Library/CloudStorage/Dropbox-Personal/data_for_running/E_matrix_corrected.pt',
                              map_location='cpu', weights_only=False)[:50000]
print(f"Loaded Y: {Y_full.shape}; E_corrected: {E_corrected_full.shape}")

pi_np = pi_full.detach().numpy() if torch.is_tensor(pi_full) else np.asarray(pi_full)
Y_np = Y_full.detach().numpy()
E_np = E_corrected_full.detach().numpy() if torch.is_tensor(E_corrected_full) else E_corrected_full

N, D, T = Y_np.shape
print(f"\nDimensions: N={N}, D={D}, T={T}")

# At-risk mask
print('Building at-risk mask...')
at_risk = (E_np[:, :, None] >= np.arange(T)[None, None, :])
print('At-risk mask done')

# Collect at-risk predictions and observations across all timepoints
all_pred = []
all_obs = []
for t in range(T):
    if t % 10 == 0:
        print(f'  collecting t={t}/{T}')
    mask_t = at_risk[:, :, t]
    if mask_t.sum() > 0:
        all_pred.append(pi_np[:, :, t][mask_t])
        all_obs.append(Y_np[:, :, t][mask_t])
all_pred = np.concatenate(all_pred)
all_obs = np.concatenate(all_obs)
print(f"Collected {len(all_pred):,} (pred, obs) pairs")
print(f"  mean predicted: {all_pred.mean():.2e}")
print(f"  mean observed:  {all_obs.mean():.2e}")

# Log-spaced binning
n_bins = 50
min_bin_count = 10000
bin_edges = np.logspace(np.log10(max(1e-7, all_pred.min())),
                        np.log10(all_pred.max()),
                        n_bins + 1)
bin_means, obs_means, counts = [], [], []
for i in range(n_bins):
    mask = (all_pred >= bin_edges[i]) & (all_pred < bin_edges[i + 1])
    if mask.sum() >= min_bin_count:
        bin_means.append(all_pred[mask].mean())
        obs_means.append(all_obs[mask].mean())
        counts.append(int(mask.sum()))

fig, ax = plt.subplots(figsize=(12, 10), dpi=300)
ax.plot([1e-7, 1], [1e-7, 1], '--', color='gray', alpha=0.5, label='Perfect calibration', linewidth=2)
ax.set_xscale('log')
ax.set_yscale('log')
ax.plot(bin_means, obs_means, 'o-', color='#c8102e', markersize=10, linewidth=2.5,
        label='Reparam observed rates', alpha=0.85)
for x, y, c in zip(bin_means, obs_means, counts):
    ax.annotate(f'n={c:,}', (x, y), xytext=(0, 12), textcoords='offset points',
                ha='center', fontsize=9)

mse = np.mean((np.array(bin_means) - np.array(obs_means)) ** 2)
stats_text = (f'MSE: {mse:.2e}\n'
              f'Mean Predicted: {all_pred.mean():.2e}\n'
              f'Mean Observed:  {all_obs.mean():.2e}\n'
              f'N total: {sum(counts):,}')
ax.text(0.05, 0.95, stats_text, transform=ax.transAxes, va='top',
        bbox=dict(boxstyle='round', facecolor='white', alpha=0.9), fontsize=11)
ax.grid(True, which='both', linestyle='--', alpha=0.3)
ax.set_xlabel('Predicted Event Rate', fontsize=14, fontweight='bold')
ax.set_ylabel('Observed Event Rate', fontsize=14, fontweight='bold')
ax.set_title('Calibration Across All Follow-up (At-Risk Only) — REPARAM\nFull Dataset (50k subset of 400k)',
             fontsize=15, fontweight='bold', pad=18)
ax.legend(loc='lower right', fontsize=12)
plt.tight_layout()

out_dir = Path('/Users/sarahurbut/aladynoulli2/pyScripts/dec_6_revision/new_notebooks/results/paper_figs/fig5_reparam')
out_dir.mkdir(parents=True, exist_ok=True)
out_path = out_dir / 'calibration_plots_full_400k_reparam.pdf'
plt.savefig(out_path, format='pdf', dpi=300, bbox_inches='tight')
print(f'\nSaved: {out_path}')
plt.show()