Preliminary testing of a prespecified liability architecture for autism: theory-guided pathogenetic triad models outperform strength-matched alternatives.

Code ↔ Paper

12 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 12 matches

[1] § Materials and methods › Analytic overview ↔ scripts/_04_pt_triadindex.py, lines 615–724 · score 0.82 · outer training, log loss, model fitting, logistic regression, leakage free, stratified
[2] § Materials and methods › Constructing the triplet multiverse space ↔ scripts/_05_pt_multiverse.py, lines 4–91 · score 0.82 · ridge logistic regression, repeated domains, distinct domains, full predictor, domain restriction, multiverse
[3] § Materials and methods › Analytic overview ↔ scripts/_05_pt_multiverse.py, lines 511–571 · score 0.79 · outer training, log loss, leakage free, model fitting, stratified, tuned
[4] § Results › Does a prespecified PT operationalization show out-of-sample separation? › Primary PT model in triad space ↔ scripts/_03_pt_primary_classifier.py, lines 970–1078 · score 0.69 · bootstrap envelope, ROC curve, Triad space, decision boundary, LR, clustering
[5] § Materials and methods › Constructing the triplet multiverse space › Matched paired comparisons of PT vs non-PT triplets ↔ scripts/_05_pt_multiverse.py, lines 4–91 · score 0.67 · comparator triplets, PT triplet, univariate AUC, nearest, profile, strength
[6] § Materials and methods › Data processing and univariate characterization ↔ scripts/_02_pt_preprocessing.py, lines 4–70 · score 0.65 · right skewed, log transformed, oriented, strength, univariate, variables
[7] § Results › Does a prespecified PT operationalization show out-of-sample separation? › TriadIndex: collapsing triad space onto a single axis ↔ scripts/_03_pt_primary_classifier.py, lines 970–1078 · score 0.60 · calibration curve, ROC curve, decision boundary, sensitivity, triad, space
[8] § Materials and methods › High-dimensional landscape ↔ scripts/_06_pt_kitchen_sink.py, lines 4–77 · score 0.59 · exhaustive enumeration, 2–3, PT models, landscape, predictor
[9] § Materials and methods › The Triadindex as a continuous PT risk score ↔ scripts/_04_pt_triadindex.py, lines 1290–1431 · score 0.59 · TriadIndex, equal weight, raw, risk, TI, axis
[10] § Results › Does a prespecified PT operationalization show out-of-sample separation? › Permutation test and unsupervised clustering ↔ scripts/_03_pt_primary_classifier.py, lines 1294–1345 · score 0.53 · ROC curve, triad space, unsupervised, clustering, BA, model
[11] § Materials and methods › Measures › Non-PT predictors ↔ scripts/_02_pt_preprocessing.py, lines 4–70 · score 0.52 · domain assignments, preprocessing, MEG, pipelines, variables, NB
[12] § Results › Does a prespecified PT operationalization show out-of-sample separation? › TriadIndex: collapsing triad space onto a single axis ↔ scripts/_04_pt_triadindex.py, lines 899–1031 · score 0.52 · calibration curve, ROC curve, sensitivity, BA, TriadIndex, AUC

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 1,435 lines · 49 KB · no license · 3 matches

# Copyright (C) 2025 Darko Sarovic
# SPDX-License-Identifier: AGPL-3.0-or-later only
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
04_pt_triadindex.py
===================
Purpose
----------------------------------------------------------------------------------------
TriadIndex (TI) construction and leakage-free 1-SE nested-CV evaluation (ridge logistic)
AUTOMATIC MULTI-MODE RUN:
- specified weights (default 2:1:1)
- equal weights (1:1:1)
- Exploratory (not implemented here because of small sample): 3D-derived weights (|beta| from 3-predictor logistic), computed leakage-free within CV
Conceptual separation (important for Methods):
- Panel A (TI distribution) uses full-sample standardization of TI_raw -> TI_z to define a stable TI axis.
This is descriptive only.
- Predictive performance uses leakage-free nested CV:
outer LOOCV + inner 5-fold CV to tune C by minimizing log loss.
All standardization and TI construction are performed within training partitions only.
Outputs: <OUTPUT_ROOT>/04_TriadIndex/
- ti_models_summary.csv
- ti_weight_info.csv
- For each model:
ti_predictions_<model>.csv
ti_calibration_<model>.csv
ti_C_diagnostics_<model>.csv
ti_outer_fold_summary_<model>.csv
- Figures:
TriadIndex_TI_Calibration_ROC_combined_<model>.png/.pdf/.eps
"""
from __future__ import annotations
import argparse
import math
from dataclasses import dataclass
from pathlib import Path
from typing import Dict, List, Optional, Tuple
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, roc_curve, brier_score_loss, log_loss
from sklearn.model_selection import StratifiedKFold
from _01_pt_utils import set_global_random_seed, RANDOM_SEED
# -----------------------------------------------------------------------------
# Defaults
# -----------------------------------------------------------------------------
DEFAULT_OUTPUT_ROOT = Path("/User/Desktop/PT_project_folder/Output_folder")
DEFAULT_INPUT_XLSX = DEFAULT_OUTPUT_ROOT / "derived_data" / "triad_data_logfixed.xlsx"
DIAG_COL = "Diagnosis" # 0 = non-autism, 1 = autism
AUTISM_CODE = 1
CONTROL_CODE = 0
AP_CANDS = ["AQ", "AQ_total", "AP_AQ"]
CC_CANDS = ["WMIQ", "CC_WMIQ"]
NB_CANDS = ["SD1", "NB_SD1"]
ID_CANDS = ["ID", "Id", "Participant", "ParticipantID", "Subject", "SubjectID", "Case", "CaseID"]
# Match your primary classifier grid (10^-4 ... 10^4, 17 values)
C_GRID = np.logspace(-4, 4, 17)
# -----------------------------------------------------------------------------
# Styling
# -----------------------------------------------------------------------------
def set_plot_style() -> None:
plt.rcParams["font.family"] = "serif"
plt.rcParams["font.serif"] = ["Times New Roman"]
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.color"] = "0.82"
plt.rcParams["grid.linewidth"] = 1.0
plt.rcParams["grid.linestyle"] = "-"
plt.rcParams["axes.edgecolor"] = "0.7"
plt.rcParams["axes.linewidth"] = 1.0
plt.rcParams["savefig.dpi"] = 300
plt.rcParams["figure.dpi"] = 150
# -----------------------------------------------------------------------------
# Helpers
# -----------------------------------------------------------------------------
# -------------------------
# Positive-class convention
# -------------------------
POS_LABEL = AUTISM_CODE # autism/case must be coded as 1
def proba_of_label(model: LogisticRegression, X: np.ndarray, label: int = POS_LABEL) -> np.ndarray:
"""Return P(y==label) robustly even if class ordering differs."""
probs = model.predict_proba(X)
classes = list(model.classes_)
if label not in classes:
raise ValueError(f"Label {label} not in model.classes_={classes}")
return probs[:, classes.index(label)]
def centered_decision_score(model: LogisticRegression, X: np.ndarray) -> np.ndarray:
"""Centered decision score: decision_function - intercept = X @ coef (binary)."""
s = np.asarray(model.decision_function(X), dtype=float).reshape(-1)
try:
intercept = float(np.asarray(model.intercept_).reshape(-1)[0])
except Exception:
intercept = 0.0
return s - intercept
def sigmoid(z: float) -> float:
return float(1.0 / (1.0 + math.exp(-float(z))))
def best_ba_point(y_true: np.ndarray, scores: np.ndarray) -> Tuple[float, float, float, float]:
"""
Maximize BA over thresholds on `scores`.
Returns: (best_thr, best_BA, FPR_at_best, TPR_at_best)
"""
y_true = y_true.astype(int)
scores = scores.astype(float)
thr_unique = np.unique(scores)
best_thr = float(thr_unique[0])
best_ba = -np.inf
best_fpr = float("nan")
best_tpr = float("nan")
for thr in thr_unique:
pred = (scores >= thr).astype(int)
tp = np.sum((pred == 1) & (y_true == 1))
tn = np.sum((pred == 0) & (y_true == 0))
fp = np.sum((pred == 1) & (y_true == 0))
fn = np.sum((pred == 0) & (y_true == 1))
if (tp + fn) == 0 or (tn + fp) == 0:
continue
sens = tp / (tp + fn)
spec = tn / (tn + fp)
ba = 0.5 * (sens + spec)
if ba > best_ba:
best_ba = float(ba)
best_thr = float(thr)
best_fpr = float(1.0 - spec)
best_tpr = float(sens)
return best_thr, float(best_ba), float(best_fpr), float(best_tpr)
def score_to_prob_threshold(thr_score: float, intercept: float) -> float:
"""Probability at the decision threshold (intercept + centered_score)."""
return sigmoid(float(intercept) + float(thr_score))
def score_to_ti_threshold(thr_score: float, beta: float) -> float:
"""Map centered-score threshold to TI threshold: thr_score = beta * TI."""
if not np.isfinite(beta) or np.isclose(beta, 0.0):
return float("nan")
return float(thr_score / beta)
def ensure_dir(p: Path) -> None:
p.mkdir(parents=True, exist_ok=True)
def find_first_present(df: pd.DataFrame, candidates: List[str]) -> str:
for c in candidates:
if c in df.columns:
return c
raise KeyError(f"None of these columns were found: {candidates}")
def find_optional_id_col(df: pd.DataFrame) -> Optional[str]:
for c in ID_CANDS:
if c in df.columns:
return c
return None
def _mean_sd(x: np.ndarray) -> Tuple[float, float]:
mu = float(np.nanmean(x))
sd = float(np.nanstd(x, ddof=0))
if (not np.isfinite(sd)) or sd <= 0:
sd = 1.0
return mu, sd
def z_from_train(x: np.ndarray, mu: float, sd: float) -> np.ndarray:
return (x - mu) / sd
def best_ba_threshold_from_probs(y_true: np.ndarray, y_prob: np.ndarray) -> Tuple[float, float]:
fpr, tpr, thr = roc_curve(y_true, y_prob)
ba = 0.5 * (tpr + (1 - fpr))
k = int(np.argmax(ba))
return float(thr[k]), float(ba[k])
def prob_to_ti_threshold(thr_prob: float, intercept: float, beta: float) -> float:
if thr_prob <= 0.0 or thr_prob >= 1.0:
raise ValueError("thr_prob must be in (0,1).")
if np.isclose(beta, 0.0):
return float("nan")
logit = math.log(thr_prob / (1.0 - thr_prob))
return float((logit - intercept) / beta)
def calibration_table_quantile(y_true: np.ndarray, y_prob: np.ndarray, n_bins: int = 5) -> pd.DataFrame:
df = pd.DataFrame({"y": y_true.astype(int), "p": y_prob.astype(float)}).sort_values("p").reset_index(drop=True)
qs = np.linspace(0, 1, n_bins + 1)
edges = df["p"].quantile(qs).to_numpy()
for i in range(1, len(edges)):
if edges[i] <= edges[i - 1]:
edges[i] = edges[i - 1] + 1e-8
rows = []
for b in range(n_bins):
lo, hi = float(edges[b]), float(edges[b + 1])
if b == n_bins - 1:
m = (df["p"] >= lo) & (df["p"] <= hi)
else:
m = (df["p"] >= lo) & (df["p"] < hi)
sub = df.loc[m]
rows.append(
{
"bin": b + 1,
"lower": lo,
"upper": hi,
"mean_prob": float(sub["p"].mean()) if len(sub) else float("nan"),
"obs_rate": float(sub["y"].mean()) if len(sub) else float("nan"),
"n": int(len(sub)),
}
)
return pd.DataFrame(rows)
def calibration_table_fixed_edges(y_true: np.ndarray, y_prob: np.ndarray, bin_edges: List[float]) -> pd.DataFrame:
df = pd.DataFrame({"y": y_true.astype(int), "p": y_prob.astype(float)})
edges = np.asarray(bin_edges, dtype=float)
if edges.ndim != 1 or edges.size < 3:
raise ValueError("bin_edges must be a 1D list/array with at least 3 values (>=2 bins).")
# Ensure [0,1] coverage and strict monotonicity
edges[0] = min(edges[0], 0.0)
edges[-1] = max(edges[-1], 1.0)
for i in range(1, len(edges)):
if edges[i] <= edges[i - 1]:
edges[i] = edges[i - 1] + 1e-8
df["bin"] = pd.cut(df["p"], bins=edges, include_lowest=True, right=True)
rows = []
for interval, g in df.groupby("bin", observed=True):
if g.shape[0] == 0:
continue
rows.append(
{
"bin": str(interval),
"lower": float(interval.left),
"upper": float(interval.right),
"mean_prob": float(g["p"].mean()),
"obs_rate": float(g["y"].mean()),
"n": int(g.shape[0]),
}
)
return pd.DataFrame(rows).sort_values("mean_prob").reset_index(drop=True)
def bootstrap_auc_ci(
y_true: np.ndarray,
y_score: np.ndarray,
n_boot: int = 2000,
seed: int = 12345,
alpha: float = 0.05,
) -> Tuple[float, float, int]:
"""
Percentile bootstrap CI for ROC AUC on paired (y_true, y_score).
Skips samples with only one class.
"""
y_true = np.asarray(y_true).astype(int)
y_score = np.asarray(y_score).astype(float)
n = y_true.shape[0]
if n_boot <= 0:
return float("nan"), float("nan"), 0
rng = np.random.default_rng(int(seed))
aucs: List[float] = []
for _ in range(int(n_boot)):
idx = rng.integers(0, n, size=n)
ys = y_true[idx]
if np.unique(ys).size < 2:
continue
aucs.append(float(roc_auc_score(ys, y_score[idx])))
if len(aucs) == 0:
return float("nan"), float("nan"), 0
lo, hi = np.quantile(np.asarray(aucs, float), [alpha / 2.0, 1.0 - alpha / 2.0])
return float(lo), float(hi), int(len(aucs))
def bootstrap_brier_ci(
y_true: np.ndarray,
y_prob: np.ndarray,
n_boot: int = 2000,
seed: int = 12345,
alpha: float = 0.05,
) -> Tuple[float, float, int]:
"""
Percentile bootstrap CI for Brier score on paired (y_true, y_prob).
Returns (ci_low, ci_high, n_valid_boot).
"""
y_true = np.asarray(y_true).astype(int)
y_prob = np.asarray(y_prob).astype(float)
n = y_true.shape[0]
if n_boot <= 0:
return float("nan"), float("nan"), 0
rng = np.random.default_rng(int(seed))
vals = []
for _ in range(int(n_boot)):
idx = rng.integers(0, n, size=n)
vals.append(float(brier_score_loss(y_true[idx], y_prob[idx])))
lo, hi = np.quantile(np.asarray(vals, float), [alpha / 2.0, 1.0 - alpha / 2.0])
return float(lo), float(hi), int(len(vals))
def bootstrap_ba_ci(
y_true: np.ndarray,
y_score: np.ndarray,
n_boot: int = 2000,
seed: int = 12345,
alpha: float = 0.05,
) -> Tuple[float, float, int]:
"""
Percentile bootstrap CI for best BA from thresholds on paired (y_true, y_score).
Skips samples with only one class.
"""
y_true = np.asarray(y_true).astype(int)
y_score = np.asarray(y_score).astype(float)
n = y_true.shape[0]
if n_boot <= 0:
return float("nan"), float("nan"), 0
rng = np.random.default_rng(int(seed))
vals: List[float] = []
for _ in range(int(n_boot)):
idx = rng.integers(0, n, size=n)
ys = y_true[idx]
if np.unique(ys).size < 2:
continue
_, ba, _, _ = best_ba_point(ys, y_score[idx])
vals.append(float(ba))
if len(vals) == 0:
return float("nan"), float("nan"), 0
lo, hi = np.quantile(np.asarray(vals, float), [alpha / 2.0, 1.0 - alpha / 2.0])
return float(lo), float(hi), int(len(vals))
def save_figure_triplet(fig: plt.Figure, out_base: Path) -> None:
out_base.parent.mkdir(parents=True, exist_ok=True)
fig.savefig(out_base.with_suffix(".png"), bbox_inches="tight")
fig.savefig(out_base.with_suffix(".pdf"), bbox_inches="tight")
fig.savefig(out_base.with_suffix(".eps"), bbox_inches="tight")
def normalize_weights(w_ap: float, w_cc: float, w_nb: float, method: str = "l1") -> Tuple[float, float, float]:
w = np.array([w_ap, w_cc, w_nb], dtype=float)
if method == "none":
return float(w[0]), float(w[1]), float(w[2])
if method == "sum":
s = float(np.sum(w))
if np.isclose(s, 0.0):
raise ValueError("Sum of weights is zero; cannot normalize by sum.")
return float(w[0] / s), float(w[1] / s), float(w[2] / s)
s = float(np.sum(np.abs(w)))
if np.isclose(s, 0.0):
raise ValueError("All weights are zero; provide nonzero weights.")
return float(w[0] / s), float(w[1] / s), float(w[2] / s)
def format_model_tag(tag: str) -> str:
return tag.replace(" ", "_").replace(":", "").replace("/", "_").replace("|", "_")
def safe_kfold_k(y: np.ndarray, desired_k: int) -> int:
"""Choose a safe K for label-independent KFold in small, imbalanced samples.
Logistic regression requires both classes in each *training* fold. With KFold, rare edge cases
can produce an inner split whose training side has only one class if the minority class is tiny.
We cap K by the per-class counts to reduce that risk, while still using label-independent splits.
"""
y = y.astype(int)
n0 = int(np.sum(y == 0))
n1 = int(np.sum(y == 1))
if n0 == 0 or n1 == 0:
raise ValueError("Cannot do CV: only one class present.")
k = min(int(desired_k), n0, n1, len(y))
return max(2, k)
def choose_C_one_se_smallest(
Cs: np.ndarray,
mean_losses: np.ndarray,
se_losses: np.ndarray,
) -> float:
"""1-SE rule for a minimization metric (log loss): pick smallest C within 1 SE of the best mean."""
Cs = np.asarray(Cs, dtype=float)
mean_losses = np.asarray(mean_losses, dtype=float)
se_losses = np.asarray(se_losses, dtype=float)
finite = np.isfinite(mean_losses) & np.isfinite(se_losses) & np.isfinite(Cs)
if not np.any(finite):
return float(Cs[0])
best_idx = int(np.nanargmin(np.where(finite, mean_losses, np.nan)))
best_mean = float(mean_losses[best_idx])
best_se = float(se_losses[best_idx])
threshold = best_mean + best_se
eligible = finite & (mean_losses <= threshold + 1e-12)
if not np.any(eligible):
return float(Cs[best_idx])
return float(np.min(Cs[eligible]))
def ridge_logistic(C: float) -> LogisticRegression:
return LogisticRegression(
penalty="l2",
C=float(C),
solver="lbfgs",
max_iter=6000,
)
# -----------------------------------------------------------------------------
# Data construction
# -----------------------------------------------------------------------------
@dataclass(frozen=True)
class TriadData:
df: pd.DataFrame
ap_var: str
cc_var: str
nb_var: str
def build_triad_fullsample(df: pd.DataFrame) -> TriadData:
ap_var = find_first_present(df, AP_CANDS)
cc_var = find_first_present(df, CC_CANDS)
nb_var = find_first_present(df, NB_CANDS)
id_col = find_optional_id_col(df)
cols = [DIAG_COL, ap_var, cc_var, nb_var]
if id_col is not None:
cols = [id_col] + cols
tri = df[cols].copy().dropna(subset=[DIAG_COL, ap_var, cc_var, nb_var]).reset_index(drop=True)
if id_col is not None:
tri = tri.rename(columns={id_col: "participant_id"})
else:
tri["participant_id"] = np.arange(tri.shape[0], dtype=int)
tri = tri.rename(columns={ap_var: "AP_raw", cc_var: "CC_raw", nb_var: "NB_raw"})
# Full-sample z's for descriptive TI axis (Panel A only)
ap = tri["AP_raw"].to_numpy(float)
cc = tri["CC_raw"].to_numpy(float)
nb = tri["NB_raw"].to_numpy(float)
ap_mu, ap_sd = _mean_sd(ap)
cc_mu, cc_sd = _mean_sd(cc)
nb_mu, nb_sd = _mean_sd(nb)
tri["z_AP"] = z_from_train(ap, ap_mu, ap_sd)
tri["z_CC_risk"] = -z_from_train(cc, cc_mu, cc_sd)
tri["z_NB_risk"] = -z_from_train(nb, nb_mu, nb_sd)
return TriadData(df=tri, ap_var=ap_var, cc_var=cc_var, nb_var=nb_var)
# -----------------------------------------------------------------------------
# TI construction (leakage-free) utilities
# -----------------------------------------------------------------------------
def compute_ti_from_train(
ap: np.ndarray,
cc: np.ndarray,
nb: np.ndarray,
train_idx: np.ndarray,
apply_idx: np.ndarray,
w_ap: float,
w_cc: float,
w_nb: float,
) -> Tuple[np.ndarray, float, float]:
"""
Compute TI_z for apply_idx using parameters learned from train_idx only:
- z-score AP/CC/NB using train means/SDs
- risk-orient CC and NB by sign reversal
- compute TI_raw = w_ap*z_AP + w_cc*z_CC_risk + w_nb*z_NB_risk
- standardize TI_raw using train TI_raw mean/SD to yield TI_z
Returns
-------
ti_apply : array of TI_z for apply_idx
ti_mu, ti_sd : mean and SD of TI_raw in train set (used for standardization)
"""
ap_tr = ap[train_idx]
cc_tr = cc[train_idx]
nb_tr = nb[train_idx]
ap_mu, ap_sd = _mean_sd(ap_tr)
cc_mu, cc_sd = _mean_sd(cc_tr)
nb_mu, nb_sd = _mean_sd(nb_tr)
# z for train (needed to compute TI_raw train stats)
z_ap_tr = z_from_train(ap_tr, ap_mu, ap_sd)
z_cc_tr = z_from_train(cc_tr, cc_mu, cc_sd)
z_nb_tr = z_from_train(nb_tr, nb_mu, nb_sd)
z_cc_risk_tr = -z_cc_tr
z_nb_risk_tr = -z_nb_tr
ti_raw_tr = w_ap * z_ap_tr + w_cc * z_cc_risk_tr + w_nb * z_nb_risk_tr
ti_mu, ti_sd = _mean_sd(ti_raw_tr)
# apply
ap_ap = ap[apply_idx]
cc_ap = cc[apply_idx]
nb_ap = nb[apply_idx]
z_ap = z_from_train(ap_ap, ap_mu, ap_sd)
z_cc = z_from_train(cc_ap, cc_mu, cc_sd)
z_nb = z_from_train(nb_ap, nb_mu, nb_sd)
ti_raw = w_ap * z_ap + w_cc * (-z_cc) + w_nb * (-z_nb)
ti_z = z_from_train(ti_raw, ti_mu, ti_sd)
return ti_z.astype(float), float(ti_mu), float(ti_sd)
def derive_weights_from_3d_logistic_trainonly(
ap: np.ndarray,
cc: np.ndarray,
nb: np.ndarray,
y: np.ndarray,
train_idx: np.ndarray,
) -> Tuple[float, float, float, float, float, float]:
"""
Derive positive L1-normalized weights from a 3-predictor logistic model fit
on standardized risk-oriented components within train_idx only.
Returns: (w_ap, w_cc, w_nb, beta_ap, beta_cc, beta_nb)
"""
y_tr = y[train_idx].astype(int)
ap_tr = ap[train_idx]
cc_tr = cc[train_idx]
nb_tr = nb[train_idx]
ap_mu, ap_sd = _mean_sd(ap_tr)
cc_mu, cc_sd = _mean_sd(cc_tr)
nb_mu, nb_sd = _mean_sd(nb_tr)
z_ap = z_from_train(ap_tr, ap_mu, ap_sd)
z_cc_risk = -z_from_train(cc_tr, cc_mu, cc_sd)
z_nb_risk = -z_from_train(nb_tr, nb_mu, nb_sd)
X = np.column_stack([z_ap, z_cc_risk, z_nb_risk]).astype(float)
# Use (near-)unpenalized as a descriptive 3D fit; fall back if it fails.
try:
m = LogisticRegression(penalty=None, solver="lbfgs", max_iter=10000)
m.fit(X, y_tr)
except Exception:
m = LogisticRegression(penalty="l2", C=1e6, solver="lbfgs", max_iter=10000)
m.fit(X, y_tr)
beta = m.coef_[0].astype(float)
s = float(np.sum(np.abs(beta)))
if np.isclose(s, 0.0):
w = np.array([1.0, 1.0, 1.0], float) / 3.0
else:
w = np.abs(beta) / s
return float(w[0]), float(w[1]), float(w[2]), float(beta[0]), float(beta[1]), float(beta[2])
# -----------------------------------------------------------------------------
# Nested tuning: inner CV chooses C by minimizing log loss
# -----------------------------------------------------------------------------
def tune_C_inner_cv_for_ti(
ap: np.ndarray,
cc: np.ndarray,
nb: np.ndarray,
y: np.ndarray,
outer_train_idx: np.ndarray,
w_mode: str,
fixed_w: Tuple[float, float, float],
inner_k_desired: int,
rng_seed: int,
) -> Tuple[float, pd.DataFrame]:
"""Tune C on outer_train_idx using inner CV + 1-SE rule (smallest C within 1 SE of the best mean log loss).
Inner splitting is stratified (KFold with shuffle). For each inner split:
- TI is constructed using inner-train only (leakage-free)
- ridge logistic is fit on TI(inner-train)
- log loss is evaluated on TI(inner-val)
w_mode:
- "fixed": use fixed_w for all inner splits
- "derived": derive weights from a 3D logistic model fit on inner-train only (leakage-free)
Returns:
best_C (float): selected by the 1-SE rule
diagnostics (DataFrame): per-C mean/se log loss and folds used
"""
outer_idx = np.asarray(outer_train_idx, dtype=int)
y_outer = y[outer_idx].astype(int)
k_inner = safe_kfold_k(y_outer, inner_k_desired)
kf = StratifiedKFold(n_splits=k_inner, shuffle=True, random_state=int(rng_seed))
rows = []
mean_losses = []
se_losses = []
for C in C_GRID:
fold_losses = []
folds_used = 0
for fold_id, (tr_pos, va_pos) in enumerate(kf.split(np.zeros(len(outer_idx)), y_outer), start=1):
inner_tr_idx = outer_idx[tr_pos]
inner_va_idx = outer_idx[va_pos]
y_tr = y[inner_tr_idx].astype(int)
y_va = y[inner_va_idx].astype(int)
# Guard: training must contain both classes for logistic regression
if np.unique(y_tr).size < 2:
continue
# Weights for TI (fixed or derived on inner-train only)
if w_mode == "derived":
w_ap, w_cc, w_nb, *_ = derive_weights_from_3d_logistic_trainonly(
ap=ap, cc=cc, nb=nb, y=y, train_idx=inner_tr_idx
)
else:
w_ap, w_cc, w_nb = fixed_w
# Leakage-free TI construction: stats learned on inner-train only
ti_tr, _, _ = compute_ti_from_train(
ap=ap, cc=cc, nb=nb,
train_idx=inner_tr_idx, apply_idx=inner_tr_idx,
w_ap=w_ap, w_cc=w_cc, w_nb=w_nb
)
ti_va, _, _ = compute_ti_from_train(
ap=ap, cc=cc, nb=nb,
train_idx=inner_tr_idx, apply_idx=inner_va_idx,
w_ap=w_ap, w_cc=w_cc, w_nb=w_nb
)
X_tr = ti_tr.reshape(-1, 1)
X_va = ti_va.reshape(-1, 1)
m = ridge_logistic(float(C))
m.fit(X_tr, y_tr)
p = proba_of_label(m, X_va, POS_LABEL)
fold_losses.append(float(log_loss(y_va, p, labels=[0, 1])))
folds_used += 1
if folds_used == 0:
mean_ll = float("nan")
se_ll = float("nan")
else:
mean_ll = float(np.mean(fold_losses))
se_ll = float(np.std(fold_losses, ddof=1) / np.sqrt(folds_used)) if folds_used >= 2 else 0.0
rows.append({
"C": float(C),
"k_inner": int(k_inner),
"folds_used": int(folds_used),
"mean_logloss": mean_ll,
"se_logloss": se_ll,
"w_mode": str(w_mode),
})
mean_losses.append(mean_ll)
se_losses.append(se_ll)
df = pd.DataFrame(rows)
if np.all(np.isnan(np.asarray(mean_losses, float))):
raise RuntimeError("Inner CV failed: all candidate C values produced zero usable folds (single-class training folds).")
best_C = choose_C_one_se_smallest(
np.asarray(C_GRID, float),
np.asarray(mean_losses, float),
np.asarray(se_losses, float),
)
return float(best_C), df
def loocv_ti_probs_nested_ridge(
triad: pd.DataFrame,
w_mode: str,
fixed_w: Tuple[float, float, float],
inner_k: int,
random_seed: int,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, pd.DataFrame]:
"""
Outer LOOCV:
- within each outer fold: tune C via inner CV (log loss)
- construct TI using outer-train only
- fit ridge logistic with tuned C on TI_train
- predict held-out
Returns:
ti_cv (held-out foldwise TI_z),
prob_cv (held-out probabilities),
fold_df (per-fold diagnostics: C, and (if derived) weights, betas)
"""
y = triad[DIAG_COL].to_numpy(int)
ap = triad["AP_raw"].to_numpy(float)
cc = triad["CC_raw"].to_numpy(float)
nb = triad["NB_raw"].to_numpy(float)
n = len(y)
ti_cv = np.full(n, np.nan, float)
prob_cv = np.full(n, np.nan, float)
score_cv = np.full(n, np.nan, float)
fold_rows = []
for i in range(n):
outer_train_idx = np.array([j for j in range(n) if j != i], dtype=int)
# Tune C (inner CV) using leakage-free TI construction within inner splits
best_C, diag = tune_C_inner_cv_for_ti(
ap=ap, cc=cc, nb=nb, y=y,
outer_train_idx=outer_train_idx,
w_mode=w_mode,
fixed_w=fixed_w,
inner_k_desired=inner_k,
rng_seed=random_seed + 1000 + i,
)
# For derived mode: derive final weights on full outer training only (leakage-free w.r.t. held-out)
if w_mode == "derived":
w_ap, w_cc, w_nb, b_ap, b_cc, b_nb = derive_weights_from_3d_logistic_trainonly(
ap=ap, cc=cc, nb=nb, y=y, train_idx=outer_train_idx
)
else:
w_ap, w_cc, w_nb = fixed_w
b_ap = b_cc = b_nb = float("nan")
# Compute outer-train TI and held-out TI using outer-train stats only
ti_tr, _, _ = compute_ti_from_train(
ap=ap, cc=cc, nb=nb,
train_idx=outer_train_idx, apply_idx=outer_train_idx,
w_ap=w_ap, w_cc=w_cc, w_nb=w_nb
)
ti_te, _, _ = compute_ti_from_train(
ap=ap, cc=cc, nb=nb,
train_idx=outer_train_idx, apply_idx=np.array([i], dtype=int),
w_ap=w_ap, w_cc=w_cc, w_nb=w_nb
)
Xtr = ti_tr.reshape(-1, 1)
ytr = y[outer_train_idx].astype(int)
m = ridge_logistic(best_C)
m.fit(Xtr, ytr)
p_te = float(proba_of_label(m, ti_te.reshape(1, 1), POS_LABEL)[0])
s_te = float(centered_decision_score(m, ti_te.reshape(1, 1))[0])
ti_cv[i] = float(ti_te[0])
prob_cv[i] = p_te
score_cv[i] = s_te
fold_rows.append({
"heldout_index": int(i),
"heldout_id": str(triad.loc[i, "participant_id"]),
"heldout_y": int(y[i]),
"C_best": float(best_C),
"inner_k": int(safe_kfold_k(y[outer_train_idx], inner_k)),
"w_mode": str(w_mode),
"w_ap": float(w_ap),
"w_cc": float(w_cc),
"w_nb": float(w_nb),
"beta_ap_3d": float(b_ap),
"beta_cc_3d": float(b_cc),
"beta_nb_3d": float(b_nb),
})
# Store per-fold inner-CV diagnostics (one file per model; append with fold index)
diag2 = diag.copy()
# Safe column insertion: avoid crashing if the diagnostics already include these fields.
for _loc, _col, _val in [
(0, "heldout_index", int(i)),
(1, "C_best_fold", float(best_C)),
(2, "w_mode", str(w_mode)),
]:
if _col in diag2.columns:
diag2[_col] = _val
else:
diag2.insert(_loc, _col, _val)
# attach later by caller
fold_df = pd.DataFrame(fold_rows)
if np.any(~np.isfinite(prob_cv)) or np.any(~np.isfinite(ti_cv)) or np.any(~np.isfinite(score_cv)):
raise RuntimeError("Non-finite LOOCV results encountered; check data and model stability.")
return ti_cv, prob_cv, score_cv, fold_df
# -----------------------------------------------------------------------------
# In-sample ridge (with 5-fold tuning of C) on full-sample TI axis (descriptive)
# -----------------------------------------------------------------------------
def tune_C_on_fixed_feature(
X: np.ndarray,
y: np.ndarray,
desired_k: int,
seed: int,
) -> Tuple[float, pd.DataFrame]:
"""Tune C on a fixed feature matrix using stratified KFold and a 1-SE rule (smallest C).
This is used for descriptive/full-sample models (not the leakage-free nested LOOCV outputs).
"""
y = y.astype(int)
k = safe_kfold_k(y, desired_k)
kf = StratifiedKFold(n_splits=k, shuffle=True, random_state=int(seed))
rows = []
mean_losses = []
se_losses = []
for C in C_GRID:
losses = []
folds_used = 0
for tr, va in kf.split(X, y):
y_tr = y[tr]
if np.unique(y_tr).size < 2:
continue
m = ridge_logistic(C)
m.fit(X[tr], y_tr)
p = m.predict_proba(X[va])[:, 1]
losses.append(float(log_loss(y[va], p, labels=[0, 1])))
folds_used += 1
if folds_used == 0:
mean_ll = np.nan
se_ll = np.nan
else:
mean_ll = float(np.mean(losses))
if folds_used >= 2:
se_ll = float(np.std(losses, ddof=1) / np.sqrt(folds_used))
else:
se_ll = 0.0
rows.append({
"C": float(C),
"k": int(k),
"folds_used": int(folds_used),
"mean_logloss": float(mean_ll) if np.isfinite(mean_ll) else np.nan,
"se_logloss": float(se_ll) if np.isfinite(se_ll) else np.nan,
})
mean_losses.append(mean_ll)
se_losses.append(se_ll)
df = pd.DataFrame(rows)
best_C = choose_C_one_se_smallest(np.asarray(C_GRID, float), np.asarray(mean_losses, float), np.asarray(se_losses, float))
return float(best_C), df
def make_joint_figure(
ti_full: np.ndarray,
y: np.ndarray,
score_in: np.ndarray,
score_cv: np.ndarray,
cal_df: pd.DataFrame,
auc_in: float,
auc_cv: float,
thr_ti_for_plot: float,
title_tag: str,
use_boxplot: bool,
jitter_seed: int = 12345,
) -> plt.Figure:
set_plot_style()
fig = plt.figure(figsize=(8, 9))
gs = GridSpec(2, 2, width_ratios=[0.6, 1.0], height_ratios=[1.0, 1.0], figure=fig)
ax_ti = fig.add_subplot(gs[:, 0])
ax_roc = fig.add_subplot(gs[0, 1])
ax_cal = fig.add_subplot(gs[1, 1])
diag = y.astype(int)
TI = np.asarray(ti_full, float)
mask_ctrl = diag == CONTROL_CODE
mask_aut = diag == AUTISM_CODE
x_ctrl_pos, x_aut_pos = 0.4, 0.6
if use_boxplot:
bp = ax_ti.boxplot(
[TI[mask_ctrl], TI[mask_aut]],
positions=[x_ctrl_pos, x_aut_pos],
widths=0.075,
patch_artist=True,
showfliers=False,
whis=1.5,
manage_ticks=False
)
for k in bp:
for art in bp[k]:
try:
art.set_zorder(1)
except Exception:
pass
for box in bp["boxes"]:
box.set_facecolor("white")
box.set_edgecolor("0.35")
box.set_linewidth(1.0)
box.set_alpha(0.35)
for whisk in bp["whiskers"]:
whisk.set_color("0.35")
whisk.set_linewidth(1.0)
whisk.set_alpha(0.35)
for cap in bp["caps"]:
cap.set_color("0.35")
cap.set_linewidth(1.0)
cap.set_alpha(0.35)
for med in bp["medians"]:
med.set_color("0.20")
med.set_linewidth(1.2)
med.set_alpha(0.55)
med.set_zorder(1)
rng = np.random.default_rng(jitter_seed)
x_ctrl = x_ctrl_pos + rng.normal(scale=0.015, size=int(mask_ctrl.sum()))
x_aut = x_aut_pos + rng.normal(scale=0.015, size=int(mask_aut.sum()))
ax_ti.scatter(x_ctrl, TI[mask_ctrl], marker="o", s=90, alpha=0.7, edgecolors="none", zorder=3)
ax_ti.scatter(x_aut, TI[mask_aut], marker="x", s=90, alpha=0.7, zorder=3)
ax_ti.axhline(thr_ti_for_plot, color="black", linestyle="--", linewidth=1.5,
label="BA-optimal cutoff (descriptive)")
ax_ti.set_xlim(0.3, 0.7)
ax_ti.set_xticks([x_ctrl_pos, x_aut_pos])
ax_ti.set_xticklabels(["Non-autism", "Autism"], fontsize=15)
ax_ti.set_xlabel("Diagnosis", fontsize=15)
ax_ti.xaxis.grid(False)
ax_ti.yaxis.grid(True, linestyle="-", color="0.9", linewidth=1, alpha=0.6, zorder=1)
ax_ti.set_ylabel("TriadIndex (standardized)", fontsize=15)
ax_ti.set_title("A. TriadIndex by diagnosis", fontweight="bold", fontsize=18)
ax_ti.set_ylim(-2, 3.5)
ax_ti.legend(loc="lower right", fontsize=9.4, framealpha=0.9)
fpr_in, tpr_in, _ = roc_curve(y, score_in)
fpr_cv, tpr_cv, _ = roc_curve(y, score_cv)
ax_roc.plot(fpr_cv, tpr_cv, color="k", linewidth=1.6, label=f"LOOCV (AUC = {auc_cv:.3f})")
ax_roc.plot(fpr_in, tpr_in, color="k", linewidth=1.6, linestyle="--", label=f"In-sample (AUC = {auc_in:.3f})")
ax_roc.plot([0, 1], [0, 1], "k:", linewidth=1)
ax_roc.set_xlim(-0.001, 1.001)
ax_roc.set_ylim(-0.001, 1.001)
ax_roc.set_xticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax_roc.set_xlabel("False positive rate (1 - specificity)", fontsize=15)
ax_roc.set_ylabel("True positive rate (sensitivity)", fontsize=15)
ax_roc.set_title("B. ROC curve", fontweight="bold", fontsize=18)
ax_roc.grid(True, linestyle="-", color="0.9", linewidth=1, alpha=0.6)
ax_roc.set_aspect("equal", adjustable="box")
ax_roc.legend(loc="lower right", fontsize=12, framealpha=0.9)
mean_pred = cal_df["mean_prob"].to_numpy(float)
obs_rate = cal_df["obs_rate"].to_numpy(float)
n_bin = cal_df["n"].to_numpy(int)
ok = np.isfinite(mean_pred) & np.isfinite(obs_rate)
mean_pred, obs_rate, n_bin = mean_pred[ok], obs_rate[ok], n_bin[ok]
ax_cal.plot([0, 1], [0, 1], linestyle="--", linewidth=1, color="0.7")
ax_cal.plot(mean_pred, obs_rate, marker="o", linestyle="-", linewidth=1.2, markersize=8)
for x, yv, n in zip(mean_pred, obs_rate, n_bin):
ax_cal.text(x, yv + 0.03, f"{int(n)}", ha="center", va="bottom", fontsize=8)
ax_cal.set_xlim(0.0, 1.0)
ax_cal.set_ylim(0.0, 1.0)
ax_cal.set_xticks([0, 0.2, 0.4, 0.6, 0.8, 1])
ax_cal.set_title("C. Calibration curve", fontweight="bold", fontsize=18)
ax_cal.set_xlabel("Mean predicted probability of autism", fontsize=15)
ax_cal.set_ylabel("Observed proportion of autism", fontsize=15)
ax_cal.grid(True, linestyle="-", color="0.9", linewidth=1, alpha=0.6)
ax_cal.set_aspect("equal", adjustable="box")
fig.suptitle(title_tag, fontsize=1) # keep layout stable; no visible suptitle
fig.tight_layout()
return fig
# -----------------------------------------------------------------------------
# Model runner
# -----------------------------------------------------------------------------
@dataclass
class TIModelResult:
model_tag: str
w_mode: str
w_ap: float
w_cc: float
w_nb: float
C_in: float
C_median_outer: float
auc_in: float
auc_cv: float
auc_cv_ci_low: float
auc_cv_ci_high: float
auc_cv_ci_nboot: int
brier_in: float
brier_cv: float
brier_cv_ci_low: float
brier_cv_ci_high: float
brier_cv_ci_nboot: int
ba_in: float
ba_cv: float
ba_cv_ci_low: float
ba_cv_ci_high: float
ba_cv_ci_nboot: int
thr_prob_in: float
thr_prob_cv: float
thr_ti_in: float
thr_ti_cv_for_plot: float
def run_one_ti_model(
triad: pd.DataFrame,
model_tag: str,
w_mode: str,
w_ap_in: float,
w_cc_in: float,
w_nb_in: float,
normalization: str,
n_cal_bins: int,
inner_k: int,
n_boot_auc: int,
out_dir: Path,
fig_dir: Path,
boxplot: bool,
random_seed: int,
) -> TIModelResult:
y = triad[DIAG_COL].to_numpy(int)
ap = triad["AP_raw"].to_numpy(float)
cc = triad["CC_raw"].to_numpy(float)
nb = triad["NB_raw"].to_numpy(float)
# Fixed weights (for specified/equal). For derived mode, these are just for the descriptive axis.
w_ap, w_cc, w_nb = normalize_weights(w_ap_in, w_cc_in, w_nb_in, method=normalization)
# -------------------------
# Full-sample TI axis (descriptive Panel A)
# -------------------------
z_ap = triad["z_AP"].to_numpy(float)
z_cc = triad["z_CC_risk"].to_numpy(float)
z_nb = triad["z_NB_risk"].to_numpy(float)
ti_raw_full = w_ap * z_ap + w_cc * z_cc + w_nb * z_nb
mu_full, sd_full = _mean_sd(ti_raw_full)
ti_full = z_from_train(ti_raw_full, mu_full, sd_full).astype(float)
# -------------------------
# In-sample ridge on full-sample TI axis (C tuned by 5-fold CV)
# -------------------------
X_in = ti_full.reshape(-1, 1)
C_in, C_diag = tune_C_on_fixed_feature(X_in, y, desired_k=5, seed=random_seed + 42)
m_in = ridge_logistic(C_in)
m_in.fit(X_in, y)
prob_in = proba_of_label(m_in, X_in, POS_LABEL)
score_in = centered_decision_score(m_in, X_in)
auc_in = float(roc_auc_score(y, score_in)) # Route 1 AUC
auc_in_prob = float(roc_auc_score(y, prob_in)) # optional diagnostic
brier_in = float(brier_score_loss(y, prob_in))
thr_score_in, ba_in, _, _ = best_ba_point(y, score_in) # Route 1 BA
intercept_in = float(m_in.intercept_[0])
beta_in = float(m_in.coef_[0, 0])
thr_prob_in = score_to_prob_threshold(thr_score_in, intercept_in) # descriptive
thr_ti_in = score_to_ti_threshold(thr_score_in, beta_in)
intercept_in = float(m_in.intercept_[0])
beta_in = float(m_in.coef_[0, 0])
thr_ti_in = prob_to_ti_threshold(thr_prob_in, intercept_in, beta_in)
# Save in-sample C tuning diagnostics (optional but helpful)
C_diag.to_csv(out_dir / f"ti_C_in_sample_{format_model_tag(model_tag)}.csv", index=False)
# -------------------------
# LOOCV nested ridge (leakage-free)
# -------------------------
ti_cv, prob_cv, score_cv, fold_df = loocv_ti_probs_nested_ridge(
triad=triad,
w_mode=w_mode,
fixed_w=(w_ap, w_cc, w_nb),
inner_k=inner_k,
random_seed=random_seed,
)
auc_cv = float(roc_auc_score(y, score_cv)) # Route 1 AUC
auc_cv_prob = float(roc_auc_score(y, prob_cv)) # optional diagnostic
auc_cv_ci_low, auc_cv_ci_high, auc_cv_ci_nboot = bootstrap_auc_ci(
y_true=y,
y_score=score_cv,
n_boot=int(n_boot_auc),
seed=int(random_seed) + 777,
alpha=0.05,
)
brier_cv = float(brier_score_loss(y, prob_cv))
thr_score_cv, ba_cv, _, _ = best_ba_point(y, score_cv) # Route 1 BA
thr_prob_cv = score_to_prob_threshold(thr_score_cv, intercept_in) # descriptive
brier_cv_ci_low, brier_cv_ci_high, brier_cv_ci_nboot = bootstrap_brier_ci(
y_true=y,
y_prob=prob_cv,
n_boot=int(n_boot_auc),
seed=int(random_seed) + 888,
alpha=0.05,
)
ba_cv_ci_low, ba_cv_ci_high, ba_cv_ci_nboot = bootstrap_ba_ci(
y_true=y,
y_score=score_cv,
n_boot=int(n_boot_auc),
seed=int(random_seed) + 999,
alpha=0.05,
)
# Map LOOCV best-BA probability threshold onto the *full-sample TI axis* for Panel A display
thr_ti_cv_for_plot = score_to_ti_threshold(thr_score_cv, beta_in)
# Calibration (LOOCV)
cal_df = calibration_table_fixed_edges(y, prob_cv, bin_edges=[0.0, 0.25, 0.5, 0.75, 1.0])
# Save per-model predictions
pred_df = triad.copy()
pred_df["TI_fullsample_axis"] = ti_full
pred_df["TI_LOOCV_foldwise"] = ti_cv
pred_df["prob_in_sample"] = prob_in
pred_df["prob_LOOCV"] = prob_cv
pred_df["score_in_sample_centered"] = score_in
pred_df["score_LOOCV_centered"] = score_cv
# Attach per-subject C (the C chosen for the fold where that subject was held out)
# (one row per subject; fold_df is indexed by heldout_index)
pred_df = pred_df.merge(
fold_df[["heldout_index", "C_best"]].rename(columns={"heldout_index": "row_index"}),
left_index=True, right_on="row_index", how="left"
).drop(columns=["row_index"])
pred_df.to_csv(out_dir / f"ti_predictions_{format_model_tag(model_tag)}.csv", index=False)
# Save calibration
cal_df.to_csv(out_dir / f"ti_calibration_{format_model_tag(model_tag)}.csv", index=False)
# Save outer-fold summary (C, and for derived also fold weights/betas)
fold_df.to_csv(out_dir / f"ti_outer_fold_summary_{format_model_tag(model_tag)}.csv", index=False)
# Summary of per-fold C
C_median_outer = float(np.median(fold_df["C_best"].to_numpy(float)))
cdiag_out = fold_df[["heldout_index", "heldout_id", "heldout_y", "C_best", "inner_k", "w_mode", "w_ap", "w_cc", "w_nb"]].copy()
cdiag_out.to_csv(out_dir / f"ti_C_diagnostics_{format_model_tag(model_tag)}.csv", index=False)
# Figure
fig = make_joint_figure(
ti_full=ti_full,
y=y,
score_in=score_in,
score_cv=score_cv,
cal_df=cal_df,
auc_in=auc_in,
auc_cv=auc_cv,
thr_ti_for_plot=thr_ti_cv_for_plot,
title_tag=model_tag,
use_boxplot=boxplot,
jitter_seed=12345,
)
fig_base = fig_dir / f"TriadIndex_TI_Calibration_ROC_combined_{format_model_tag(model_tag)}"
save_figure_triplet(fig, fig_base)
plt.close(fig)
return TIModelResult(
model_tag=model_tag,
w_mode=w_mode,
w_ap=float(w_ap), w_cc=float(w_cc), w_nb=float(w_nb),
C_in=float(C_in),
C_median_outer=float(C_median_outer),
auc_in=float(auc_in), auc_cv=float(auc_cv),
auc_cv_ci_low=float(auc_cv_ci_low),
auc_cv_ci_high=float(auc_cv_ci_high),
auc_cv_ci_nboot=int(auc_cv_ci_nboot),
brier_in=float(brier_in), brier_cv=float(brier_cv),
brier_cv_ci_low=float(brier_cv_ci_low),
brier_cv_ci_high=float(brier_cv_ci_high),
brier_cv_ci_nboot=int(brier_cv_ci_nboot),
ba_in=float(ba_in), ba_cv=float(ba_cv),
ba_cv_ci_low=float(ba_cv_ci_low),
ba_cv_ci_high=float(ba_cv_ci_high),
ba_cv_ci_nboot=int(ba_cv_ci_nboot),
thr_prob_in=float(thr_prob_in), thr_prob_cv=float(thr_prob_cv),
thr_ti_in=float(thr_ti_in),
thr_ti_cv_for_plot=float(thr_ti_cv_for_plot),
)
# -----------------------------------------------------------------------------
# CLI
# -----------------------------------------------------------------------------
def parse_args() -> argparse.Namespace:
p = argparse.ArgumentParser(description="Script 04 (1-SE nested-CV) — TriadIndex evaluation with nested LOOCV + inner C tuning (ridge).")
p.add_argument("--input-xlsx", type=str, default=str(DEFAULT_INPUT_XLSX))
p.add_argument("--output-root", type=str, default=str(DEFAULT_OUTPUT_ROOT))
p.add_argument("--random-seed", type=int, default=RANDOM_SEED)
p.add_argument("--n-boot-auc", type=int, default=2000,
help="Bootstrap replicates for LOOCV metric 95% CIs (AUC, Brier, BA); 0 disables.")
# User-specified weights (main model)
p.add_argument("--w-ap", type=float, default=2.0)
p.add_argument("--w-cc", type=float, default=1.0)
p.add_argument("--w-nb", type=float, default=1.0)
p.add_argument("--weight-normalization", type=str, choices=["l1", "sum", "none"], default="l1")
# Which models to run
p.add_argument("--run-specified", action="store_false", help="Run specified-weight TI (default ON if none selected).")
p.add_argument("--run-equal", action="store_false", help="Run equal-weight TI (1:1:1).")
p.add_argument("--run-derived", action="store_false", help="Run 3D-derived-weight TI (exploratory).")
# CV / calibration / plots
p.add_argument("--inner-k", type=int, default=5, help="Inner folds for C tuning (safe-adjusted if class counts are small).")
p.add_argument("--n-cal-bins", type=int, default=4)
p.add_argument("--boxplot", dest="boxplot", action="store_true", help="Show boxplots behind points (default: ON).")
p.add_argument("--no-boxplot", dest="boxplot", action="store_false", help="Disable boxplots.")
p.set_defaults(boxplot=True)
return p.parse_args()
# -----------------------------------------------------------------------------
# Main
# -----------------------------------------------------------------------------
def main() -> None:
args = parse_args()
# If user didn’t specify any model flags, run all three by default.
if not (args.run_specified or args.run_equal or args.run_derived):
args.run_specified = True
args.run_equal = True
args.run_derived = True
output_root = Path(args.output_root)
out_dir = output_root / "04_TriadIndex"
fig_dir = out_dir / "figures"
ensure_dir(fig_dir)
input_path = Path(args.input_xlsx)
if not input_path.exists():
raise FileNotFoundError(f"Input file not found: {input_path}")
df = pd.read_excel(input_path, sheet_name=0)
triad_data = build_triad_fullsample(df)
triad = triad_data.df.copy()
ensure_dir(out_dir)
# Full-sample derived weights (for reporting / descriptive axis only)
ap = triad["AP_raw"].to_numpy(float)
cc = triad["CC_raw"].to_numpy(float)
nb = triad["NB_raw"].to_numpy(float)
y = triad[DIAG_COL].to_numpy(int)
full_idx = np.arange(len(y), dtype=int)
w_ap_der, w_cc_der, w_nb_der, b_ap, b_cc, b_nb = derive_weights_from_3d_logistic_trainonly(
ap=ap, cc=cc, nb=nb, y=y, train_idx=full_idx
)
# Save weight info (always)
weight_info = pd.DataFrame([
{
"type": "specified_input",
"w_ap_input": float(args.w_ap), "w_cc_input": float(args.w_cc), "w_nb_input": float(args.w_nb),
"normalization": str(args.weight_normalization),
},
{
"type": "equal_weights",
"w_ap_input": 1.0, "w_cc_input": 1.0, "w_nb_input": 1.0,
"normalization": "l1",
},
{
"type": "derived_from_fullsample_3d_logistic_trainonly",
"beta_AP": float(b_ap),
"beta_CC_risk": float(b_cc),
"beta_NB_risk": float(b_nb),
"w_AP_l1": float(w_ap_der),
"w_CC_l1": float(w_cc_der),
"w_NB_l1": float(w_nb_der),
"note": "These full-sample derived weights are for descriptive reporting; LOOCV uses train-only derived weights per fold.",
}
])
weight_info.to_csv(out_dir / "ti_weight_info.csv", index=False)
results: List[TIModelResult] = []
if args.run_specified:
results.append(
run_one_ti_model(
triad=triad,
model_tag=f"TI specified ({args.w_ap:g}:{args.w_cc:g}:{args.w_nb:g})",
w_mode="fixed",
w_ap_in=float(args.w_ap),
w_cc_in=float(args.w_cc),
w_nb_in=float(args.w_nb),
normalization=str(args.weight_normalization),
n_cal_bins=int(args.n_cal_bins),
inner_k=int(args.inner_k),
out_dir=out_dir,
fig_dir=fig_dir,
boxplot=bool(args.boxplot),
n_boot_auc=int(args.n_boot_auc),
random_seed=int(args.random_seed),
)
)
if args.run_equal:
results.append(
run_one_ti_model(
triad=triad,
model_tag="TI equal (1:1:1)",
w_mode="fixed",
w_ap_in=1.0,
w_cc_in=1.0,
w_nb_in=1.0,
normalization="l1",
n_cal_bins=int(args.n_cal_bins),
inner_k=int(args.inner_k),
out_dir=out_dir,
fig_dir=fig_dir,
boxplot=bool(args.boxplot),
n_boot_auc=int(args.n_boot_auc),
random_seed=int(args.random_seed),
)
)
if args.run_derived:
# Use full-sample derived weights for the descriptive axis/in-sample part,
# but LOOCV will derive weights leakage-free per fold and inner split (w_mode="derived").
results.append(
run_one_ti_model(
triad=triad,
model_tag="TI derived (3D logistic |beta|, L1-normalized)",
w_mode="derived",
w_ap_in=float(w_ap_der),
w_cc_in=float(w_cc_der),
w_nb_in=float(w_nb_der),
normalization="none",
n_cal_bins=int(args.n_cal_bins),
inner_k=int(args.inner_k),
out_dir=out_dir,
fig_dir=fig_dir,
boxplot=bool(args.boxplot),
n_boot_auc=int(args.n_boot_auc),
random_seed=int(args.random_seed),
)
)
# Summary table
summary_df = pd.DataFrame([r.__dict__ for r in results])
summary_df.insert(0, "ap_var", triad_data.ap_var)
summary_df.insert(1, "cc_var", triad_data.cc_var)
summary_df.insert(2, "nb_var", triad_data.nb_var)
summary_df.insert(3, "n", int(triad.shape[0]))
summary_df.to_csv(out_dir / "ti_models_summary.csv", index=False)
print("------------------------------------------------------------")
print("TriadIndex multi-model analysis complete")
print(f"Output directory: {out_dir}")
print("Models run:")
for r in results:
print(
f" - {r.model_tag}: LOOCV AUC={r.auc_cv:.3f}, LOOCV Brier={r.brier_cv:.3f}, "
f"LOOCV BA(best)={r.ba_cv:.3f}, C_median_outer={r.C_median_outer:.5g}"
)
print("------------------------------------------------------------")
if __name__ == "__main__":
main()

_04_pt_triadindex.py, no license · at the source

Overview

Authors: Darko Sarovic^1,2,3,4,5

Institutions: Uppsala University (Sweden); Harvard University (United States); Sahlgrenska University Hospital (Sweden); Massachusetts General Hospital (United States); Athinoula A. Martinos Center for Biomedical Imaging (United States); University of Gothenburg (Sweden)

Journal: Frontiers in psychiatry, volume 17, article 1837909

Dates: received 24 March 2026; accepted 19 May 2026; published online 6 July 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.3389/fpsyt.2026.1837909 · PMID 42519192 · PMCID PMC13381483 · OpenAlex W7167459730

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: structural MRI / diffusion (modality), autism (population)

Methods: Statistics, Machine learning, Preprocessing, Physiology & signal measures

Keywords: autism, autistic traits, heart rate variability (HRV), liability architecture, multilevel framework, multiverse analyses, predictive modeling, pathogenetic triad

Topic: Autism Spectrum Disorder Research (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: Hjärnfonden (PS2023-0041); Göteborgs Läkaresällskap; Märta Lundqvists Stiftelse; Fredrik och Ingrid Thurings Stiftelse (2018-00419, 2020-00581, 2022-00704)

Citations: not cited yet (Europe PMC); 26 references in the paper

Abstract

Background: Neuropsychiatric conditions are heterogeneous and mechanistically diverse, yet predictive modeling studies commonly aggregate features without evaluating prespecified liability architectures. The Pathogenetic Triad (PT) is a multilevel framework proposing that diagnostic outcomes reflect the joint configuration of a trait-related domain, cognitive capacity (CC), and neuropathological burden (NB). In autism, the trait-related domain corresponds to autistic personality (AP). We evaluated this framework using out-of-sample prediction as a structured empirical test of architectural coherence rather than as a purely data-driven exercise in classification.

Methods: We analyzed a case–comparison cohort (n = 42, 21 autistic) with dense multimodal characterization, including behavioral phenotyping, psychometric measures, autonomic physiology, structural MRI morphometry, and magnetoencephalographic indices. AP was indexed by the Autism-Spectrum Quotient, CC by Wechsler’s scales of intelligence, and NB by heart-rate variability as a proxy indicator. This multimodal structure enabled direct comparison of theory-constrained PT models with strength-matched, domain-restricted atheoretical combinations. Predictive discrimination was evaluated using leakage-free nested cross-validation with permutation inference.

Results: Across multiverse specifications, low-dimensional PT models ranked among the strongest models of comparable size and achieved discrimination broadly comparable to higher-dimensional models within this dataset. Matched comparisons indicated that including all three PT domains was associated with systematic advantages relative to alternatives with similar univariate input strength, consistent with complementary configurational information relating to case status.

Conclusions: These findings provide preliminary evidence supporting the Pathogenetic Triad as a multilevel architecture of autism liability. Although based on a small and demographically restricted cohort, this densely characterized dataset permitted explicit comparison of theory-guided and atheoretical model spaces. The results illustrate how prespecified multilevel frameworks can be operationalized and empirically evaluated in neuropsychiatric samples.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 12 matches between paragraphs and lines of code.

OSF sg8ad

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Languages: Python (7)

Size: 12 files, 7 scripts

Software Heritage: not checked

Found in: “Data availability statement”

Holds: environment (docs/requirements.txt), documentation

Not found: README, license file, CITATION.cff, tests, continuous integration

Tools: NumPy (7 files), pandas (6 files), scikit-learn (6 files), Matplotlib (5 files), SciPy (3 files), seaborn (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)

27 September 2026: the link answers (HTTP 200)

7 files

scripts/_00_pt_create_synthetic_data.py, Python, 1,143 lines
scripts/_01_pt_utils.py, Python, 240 lines
scripts/_02_pt_preprocessing.py, Python, 739 lines, 2 matches
scripts/_03_pt_primary_classifier.py, Python, 1,603 lines, 3 matches
scripts/_04_pt_triadindex.py, Python, 1,435 lines, 3 matches
scripts/_05_pt_multiverse.py, Python, 1,583 lines, 3 matches
scripts/_06_pt_kitchen_sink.py, Python, 2,543 lines, 1 match

At the source: doi.org/10.17605/OSF.IO/SG8AD

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
7 scripts, each with its path and the digest of its content;
12 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The datasets presented in this article are not readily available because of ethical/consent constraints but are available from the corresponding author on reasonable request. To enable computational reproducibility, analysis scripts and a fully synthetic dataset are available at https://doi.org/10.17605/OSF.IO/SG8AD. Requests to access the datasets should be directed to DS, .

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 1 author, 8 keywords, 4 funders, 20 references.

Cite

This paper

Sarovic, D. (2026). Preliminary testing of a prespecified liability architecture for autism: theory-guided pathogenetic triad models outperform strength-matched alternatives. Frontiers in psychiatry, 17, 1837909. https://doi.org/10.3389/fpsyt.2026.1837909

BibTeX

@article{sarovic2026preliminary, author = {Sarovic, Darko}, title = {{Preliminary testing of a prespecified liability architecture for autism: theory-guided pathogenetic triad models outperform strength-matched alternatives}}, journal = {Frontiers in psychiatry}, year = {2026}, month = jul, volume = {17}, pages = {1837909}, publisher = {Frontiers Media SA}, issn = {1664-0640}, doi = {10.3389/fpsyt.2026.1837909}, url = {https://doi.org/10.3389/fpsyt.2026.1837909}, pmid = {42519192}, pmcid = {PMC13381483} }

RIS

TY - JOUR AU - Sarovic, Darko TI - Preliminary testing of a prespecified liability architecture for autism: theory-guided pathogenetic triad models outperform strength-matched alternatives T2 - Frontiers in psychiatry J2 - Front Psychiatry PY - 2026 DA - 2026/07/06 VL - 17 SP - 1837909 SN - 1664-0640 PB - Frontiers Media SA DO - 10.3389/fpsyt.2026.1837909 UR - https://doi.org/10.3389/fpsyt.2026.1837909 LA - en ER -

CSL-JSON

{ "id": "10.3389/fpsyt.2026.1837909", "type": "article-journal", "title": "Preliminary testing of a prespecified liability architecture for autism: theory-guided pathogenetic triad models outperform strength-matched alternatives", "container-title": "Frontiers in psychiatry", "author": [ { "family": "Sarovic", "given": "Darko" } ], "container-title-short": "Front Psychiatry", "volume": "17", "page": "1837909", "DOI": "10.3389/fpsyt.2026.1837909", "PMID": "42519192", "PMCID": "PMC13381483", "ISSN": "1664-0640", "publisher": "Frontiers Media SA", "URL": "https://doi.org/10.3389/fpsyt.2026.1837909", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 6 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41598-026-55163-y [code]: Autism spectrum disorder identification using machine learning models on MRI data.
Journal: Scientific reports
In common: seaborn, scikit-learn, pandas, 3 other tools, autism, structural MRI / diffusion
[2] doi:10.1007/s44192-026-00521-5 [code]: Associations between anterior hypothalamic subunits and ADHD and autistic traits revealed by deep learning MRI segmentation.
Journal: Discover mental health
In common: scikit-learn, SciPy, Matplotlib, 1 other tool, autism, structural MRI / diffusion, 1 reference
[3] doi:10.1038/s41467-026-76675-1 [code]: Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, autism
[4] doi:10.3389/fgene.2026.1799530 [code]: Sex-dependent prediction of autism.
Journal: Frontiers in genetics
In common: seaborn, scikit-learn, pandas, 3 other tools, autism
[5] doi:10.1016/j.celrep.2026.117590 [code]: Impaired behavioral inhibition in Fmr1 KO mice is linked to disrupted visual cortex theta oscillations.
Journal: Cell reports
In common: seaborn, scikit-learn, pandas, 3 other tools, autism
[6] doi:10.1038/s41586-026-10679-1 [code]: Cortical development dynamics across autism spectrum disorder mouse models.
Journal: Nature
In common: seaborn, scikit-learn, pandas, 3 other tools, autism
[7] doi:10.1162/imag.a.1220 [code]: Brain functional network connectivity interpolation characterizes the neuropsychiatric continuum and heterogeneity.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: seaborn, scikit-learn, pandas, 3 other tools, autism
[8] doi:10.1038/s42003-026-10094-2 [code]: E/I imbalance and internal noise cause weak neural representations and face recognition challenges in ASD.
Journal: Communications biology
In common: seaborn, scikit-learn, pandas, 3 other tools, autism
[9] doi:10.1002/advs.202519479 [code]: Diminished Signal-to-Noise Ratio Disrupts Somatosensory Population Encoding and Drives Tactile Hyposensitivity in the Fmr1<sup>-/y</sup> Autism Model.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: seaborn, scikit-learn, pandas, 3 other tools, autism
[10] doi:10.1002/hbm.70496 [code]: Transdiagnostic Profiles of BOLD Signal Variability in Autism and Schizophrenia Spectrum Disorders: Associations With Cognition and Functioning.
Journal: Human brain mapping
In common: seaborn, scikit-learn, pandas, 3 other tools, autism

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.