#!/usr/bin/env python3 """Calibration-first modelling protocol for INP-era lab-to-field prediction. Implements the design the pre-analysis plan specifies, rather than a single-model sketch: * Nested models A–D (A: Lighthouse score only; B: + lab metric detail; C: + public field context; D: + repeated-lab-run stability features), each fitted and scored independently so the incremental value of every evidence tier is visible. * Two target families: the binary all-CWV pass probability, and per-metric p75 regressions (LCP/INP/CLS), plus a multiclass model over the divergence taxonomy when the panel carries `divergence_type`. * Calibration-first evaluation: Brier, ECE, calibration slope and intercept (Cox recalibration on the logit scale), AUROC and AUPRC. * Split-conformal prediction intervals for the p75 regressions, with realised coverage reported against the nominal level. * Three holdout regimes: origin-blocked, time-blocked (chronological) and technology-segment. Calibration folds are origin-aware — an integer `cv` on CalibratedClassifierCV silently uses StratifiedKFold, which would put rows from one origin on both sides of the calibration split. * A machine-readable manifest tying every number to the code version, panel file and its digest. Usage: python lab-field-inp-modeling-protocol.py --panel panel.csv --out metrics.json python lab-field-inp-modeling-protocol.py --panel panel.csv --out metrics.json --holdout time """ from __future__ import annotations import argparse import hashlib import json from pathlib import Path import numpy as np import pandas as pd from sklearn.calibration import CalibratedClassifierCV from sklearn.compose import ColumnTransformer from sklearn.ensemble import HistGradientBoostingClassifier, HistGradientBoostingRegressor from sklearn.impute import SimpleImputer from sklearn.linear_model import LogisticRegression from sklearn.metrics import ( average_precision_score, brier_score_loss, f1_score, mean_absolute_error, roc_auc_score, ) from sklearn.model_selection import GroupKFold, GroupShuffleSplit from sklearn.pipeline import Pipeline from sklearn.preprocessing import OneHotEncoder, StandardScaler PROTOCOL_VERSION = '2.0' RANDOM_STATE = 20260622 # --------------------------------------------------------------------------------------------- # Nested model families. Each is a strict superset of the previous tier, so the deltas between # them are the evidence-tier contributions the article claims to measure. # --------------------------------------------------------------------------------------------- MODEL_A_NUM = ['lighthouse_score'] MODEL_B_NUM = MODEL_A_NUM + ['lab_lcp_ms', 'lab_tbt_ms', 'lab_cls', 'lab_fcp_ms', 'lab_speed_index_ms', 'lab_total_bytes'] MODEL_C_NUM = MODEL_B_NUM + ['p75_ttfb', 'rtt_high_fraction', 'navigation_bfcache_fraction', 'navigation_type_navigate', 'navigation_type_reload', 'navigation_type_prerender'] MODEL_D_NUM = MODEL_C_NUM + ['lab_run_count', 'lab_lcp_iqr_ms', 'lab_tbt_iqr_ms', 'lab_cls_iqr'] MODELS: dict[str, list[str]] = { 'A_score_only': MODEL_A_NUM, 'B_lab_detail': MODEL_B_NUM, 'C_plus_field_context': MODEL_C_NUM, 'D_plus_run_stability': MODEL_D_NUM, } CATEGORICAL_FEATURES = ['form_factor'] BINARY_TARGET = 'cwv_good_primary' P75_TARGETS = ['p75_lcp', 'p75_inp', 'p75_cls'] DIVERGENCE_TARGET = 'divergence_type' CONFORMAL_ALPHA = 0.1 # 90% nominal coverage. # --------------------------------------------------------------------------------------------- # Metrics # --------------------------------------------------------------------------------------------- def expected_calibration_error(y_true: np.ndarray, y_prob: np.ndarray, bins: int = 15) -> float: edges = np.linspace(0, 1, bins + 1) ece = 0.0 for lo, hi in zip(edges[:-1], edges[1:]): mask = (y_prob >= lo) & (y_prob < hi if hi < 1 else y_prob <= hi) if mask.sum() == 0: continue ece += (mask.sum() / len(y_prob)) * abs(y_true[mask].mean() - y_prob[mask].mean()) return float(ece) def calibration_slope_intercept(y_true: np.ndarray, y_prob: np.ndarray) -> tuple[float, float]: """Cox recalibration: regress outcomes on the predicted logit. Slope 1 and intercept 0 mean perfectly calibrated. Slope < 1 indicates over-extreme predictions; a non-zero intercept indicates systematic over/under-prediction. """ eps = 1e-6 p = np.clip(y_prob, eps, 1 - eps) logit = np.log(p / (1 - p)).reshape(-1, 1) if len(np.unique(y_true)) < 2: return float('nan'), float('nan') # Effectively unpenalised. `penalty=None` is deprecated from scikit-learn 1.8 and removed in # 1.10; a very large C means the same thing and works across versions. lr = LogisticRegression(C=1e12, solver='lbfgs', max_iter=1000) lr.fit(logit, y_true) return float(lr.coef_[0][0]), float(lr.intercept_[0]) # --------------------------------------------------------------------------------------------- # Pipelines # --------------------------------------------------------------------------------------------- def _preprocessor(numeric: list[str], categorical: list[str]) -> ColumnTransformer: return ColumnTransformer([ ('num', Pipeline([('impute', SimpleImputer(strategy='median')), ('scale', StandardScaler())]), numeric), ('cat', Pipeline([('impute', SimpleImputer(strategy='most_frequent')), ('onehot', OneHotEncoder(handle_unknown='ignore'))]), categorical), ], remainder='drop') def classifier_pipeline(numeric: list[str], categorical: list[str]) -> Pipeline: return Pipeline([('pre', _preprocessor(numeric, categorical)), ('model', HistGradientBoostingClassifier( max_iter=250, learning_rate=0.05, l2_regularization=0.05, random_state=RANDOM_STATE))]) def regressor_pipeline(numeric: list[str], categorical: list[str]) -> Pipeline: return Pipeline([('pre', _preprocessor(numeric, categorical)), ('model', HistGradientBoostingRegressor( max_iter=250, learning_rate=0.05, l2_regularization=0.05, random_state=RANDOM_STATE))]) # --------------------------------------------------------------------------------------------- # Holdout regimes # --------------------------------------------------------------------------------------------- def split_indices(df: pd.DataFrame, regime: str) -> tuple[np.ndarray, np.ndarray]: """Return (train_idx, test_idx) positional indices for the requested holdout regime.""" n = len(df) positions = np.arange(n) if regime == 'origin': groups = df['origin'].astype(str).to_numpy() splitter = GroupShuffleSplit(n_splits=1, test_size=0.2, random_state=RANDOM_STATE) return next(splitter.split(positions, groups=groups)) if regime == 'time': if 'yyyymm' not in df.columns: raise SystemExit("time holdout needs a 'yyyymm' column") months = np.sort(df['yyyymm'].unique()) if len(months) < 2: raise SystemExit('time holdout needs at least two distinct months in the panel') cutoff = months[int(len(months) * 0.8)] if len(months) > 4 else months[-1] train = positions[df['yyyymm'].to_numpy() < cutoff] test = positions[df['yyyymm'].to_numpy() >= cutoff] return train, test if regime == 'technology': if 'technology_stack' not in df.columns: raise SystemExit("technology holdout needs a 'technology_stack' column") techs = df['technology_stack'].astype(str) counts = techs.value_counts() if len(counts) < 2: raise SystemExit('technology holdout needs at least two technology segments') held_out = counts.index[-1] # smallest segment is held out entirely train = positions[(techs != held_out).to_numpy()] test = positions[(techs == held_out).to_numpy()] return train, test raise SystemExit(f'unknown holdout regime: {regime}') # --------------------------------------------------------------------------------------------- # Evaluation # --------------------------------------------------------------------------------------------- def evaluate_binary(df: pd.DataFrame, numeric: list[str], train_idx, test_idx) -> dict: features = [c for c in numeric if c in df.columns] cats = [c for c in CATEGORICAL_FEATURES if c in df.columns] if not features: return {'skipped': 'no features from this tier present in panel'} X = df[features + cats] y = df[BINARY_TARGET].astype(int).to_numpy() groups = df['origin'].astype(str).to_numpy() pipe = classifier_pipeline(features, cats) # Origin-aware calibration folds. An integer cv here would use StratifiedKFold and leak # sibling rows of the same origin across the estimator/calibration split. n_groups = len(np.unique(groups[train_idx])) if n_groups >= 3: inner = list(GroupKFold(n_splits=min(3, n_groups)).split( X.iloc[train_idx], y[train_idx], groups[train_idx])) clf = CalibratedClassifierCV(pipe, method='isotonic', cv=inner) else: clf = CalibratedClassifierCV(pipe, method='sigmoid', cv=2) clf.fit(X.iloc[train_idx], y[train_idx]) prob = clf.predict_proba(X.iloc[test_idx])[:, 1] y_test = y[test_idx] slope, intercept = calibration_slope_intercept(y_test, prob) both_classes = len(np.unique(y_test)) > 1 return { 'n_train': int(len(train_idx)), 'n_test': int(len(test_idx)), 'n_features': len(features), 'brier': float(brier_score_loss(y_test, prob)), 'auroc': float(roc_auc_score(y_test, prob)) if both_classes else None, 'auprc': float(average_precision_score(y_test, prob)) if both_classes else None, 'ece_15_bins': expected_calibration_error(y_test, prob, bins=15), 'calibration_slope': slope, 'calibration_intercept': intercept, 'base_rate_test': float(y_test.mean()), } def evaluate_p75_regression(df: pd.DataFrame, numeric: list[str], target: str, train_idx, test_idx) -> dict: """Per-metric p75 regression with split-conformal intervals and realised coverage.""" features = [c for c in numeric if c in df.columns and c != target] cats = [c for c in CATEGORICAL_FEATURES if c in df.columns] if target not in df.columns or not features: return {'skipped': f'{target} or its features absent'} mask = df[target].notna().to_numpy() train = np.array([i for i in train_idx if mask[i]]) test = np.array([i for i in test_idx if mask[i]]) if len(train) < 20 or len(test) < 5: return {'skipped': 'insufficient rows for regression'} # Split-conformal: hold out a calibration slice of the TRAINING data, grouped by origin. groups = df['origin'].astype(str).to_numpy() inner = GroupShuffleSplit(n_splits=1, test_size=0.3, random_state=RANDOM_STATE) fit_pos, calib_pos = next(inner.split(train, groups=groups[train])) fit_idx, calib_idx = train[fit_pos], train[calib_pos] X = df[features + cats] y = df[target].to_numpy(dtype=float) model = regressor_pipeline(features, cats) model.fit(X.iloc[fit_idx], y[fit_idx]) calib_residuals = np.abs(y[calib_idx] - model.predict(X.iloc[calib_idx])) q = float(np.quantile(calib_residuals, 1 - CONFORMAL_ALPHA, method='higher')) pred = model.predict(X.iloc[test_idx_valid := test]) lo, hi = pred - q, pred + q covered = float(np.mean((y[test_idx_valid] >= lo) & (y[test_idx_valid] <= hi))) return { 'n_fit': int(len(fit_idx)), 'n_calibration': int(len(calib_idx)), 'n_test': int(len(test_idx_valid)), 'mae': float(mean_absolute_error(y[test_idx_valid], pred)), 'conformal_alpha': CONFORMAL_ALPHA, 'conformal_halfwidth': q, 'nominal_coverage': 1 - CONFORMAL_ALPHA, 'realised_coverage': covered, } def evaluate_divergence_multiclass(df: pd.DataFrame, numeric: list[str], train_idx, test_idx) -> dict: if DIVERGENCE_TARGET not in df.columns: return {'skipped': f'{DIVERGENCE_TARGET} absent — run the divergence classifier first'} features = [c for c in numeric if c in df.columns] cats = [c for c in CATEGORICAL_FEATURES if c in df.columns] labels = df[DIVERGENCE_TARGET].astype(str) if labels.nunique() < 2 or not features: return {'skipped': 'fewer than two divergence classes present'} X = df[features + cats] y = labels.to_numpy() model = classifier_pipeline(features, cats) model.fit(X.iloc[train_idx], y[train_idx]) pred = model.predict(X.iloc[test_idx]) return { 'n_classes': int(labels.nunique()), 'accuracy': float((pred == y[test_idx]).mean()), 'macro_f1': float(f1_score(y[test_idx], pred, average='macro', zero_division=0)), } def main() -> None: ap = argparse.ArgumentParser(description=__doc__) ap.add_argument('--panel', required=True, help='Leakage-reviewed CSV panel') ap.add_argument('--out', required=True, help='Output metrics JSON path') ap.add_argument('--holdout', default='origin', choices=['origin', 'time', 'technology'], help='Holdout regime (default: origin-blocked)') args = ap.parse_args() panel_path = Path(args.panel) df = pd.read_csv(panel_path) required = [BINARY_TARGET, 'origin'] missing = [c for c in required if c not in df.columns] if missing: raise SystemExit(f'Missing required columns: {missing}') df = df.dropna(subset=[BINARY_TARGET, 'origin']).reset_index(drop=True) train_idx, test_idx = split_indices(df, args.holdout) if len(train_idx) == 0 or len(test_idx) == 0: raise SystemExit(f'{args.holdout} holdout produced an empty split') results = { 'manifest': { 'protocol_version': PROTOCOL_VERSION, 'panel_file': panel_path.name, 'panel_sha256': hashlib.sha256(panel_path.read_bytes()).hexdigest(), 'n_rows': int(len(df)), 'n_origins': int(df['origin'].nunique()), 'holdout_regime': args.holdout, 'random_state': RANDOM_STATE, 'months': sorted(int(m) for m in df['yyyymm'].unique()) if 'yyyymm' in df.columns else None, 'leakage_note': 'No CrUX target metric may be used as a feature for itself; ' 'p75 regressions drop their own target from the feature list.', }, 'binary_cwv_pass': {name: evaluate_binary(df, num, train_idx, test_idx) for name, num in MODELS.items()}, 'p75_regression': {t: evaluate_p75_regression(df, MODEL_C_NUM, t, train_idx, test_idx) for t in P75_TARGETS}, 'divergence_multiclass': evaluate_divergence_multiclass(df, MODEL_C_NUM, train_idx, test_idx), } Path(args.out).write_text(json.dumps(results, indent=2), encoding='utf-8') print(f'wrote {args.out}') if __name__ == '__main__': main()