# Radiomics / classical-ML — how to build a pipeline that passes review

Companion to `radiomics-ml`. *Produce* knowledge for a clinician building a radiomics or tabular
clinical-ML model without an engineer. It wires pyradiomics / scikit-learn / xgboost by name; it does
not reimplement them.

## 1. Feature extraction (radiomics) — reproducibly

Use **pyradiomics** with settings recorded for reproducibility (CLEAR / IBSI):

- **Resample** to a fixed voxel spacing; state interpolator.
- **Intensity discretisation** — a fixed **bin width** (preferred for CT/PET) or fixed bin count;
  state which and the value.
- **Normalisation** for MR (z-score) — fit per image or on train only (never on test).
- **Segmentation source** — who/what drew the ROI; single vs multi-rater (feeds stability, §2).
- Extract the standard classes (first-order, shape, GLCM/GLRLM/GLSZM/GLDM/NGTDM) ± wavelet/LoG; record
  the pyradiomics version and parameter file.

For **tabular clinical** data, the same pipeline applies from §2 onward — the "features" are labs,
demographics, and measurements instead of radiomic descriptors.

## 2. Feature stability (radiomics-specific)

Radiomic features drift with acquisition and segmentation. With test-retest or multi-rater ROIs,
compute **ICC** per feature and keep the stable ones (commonly ICC ≥ 0.75) *before* modelling. Report
how many features survived. No stability step → `NO_FEATURE_STABILITY`.

## 3. The events-per-feature problem (the classic trap)

Radiomics yields hundreds-to-thousands of features on tens of patients. Fitting a flexible model in
that regime overfits. Control it:

- **Dimensionality reduction / regularisation** — LASSO (embedded selection), a stability + redundancy
  filter (drop |r| > 0.9 pairs), or PCA.
- Keep an eye on **events per candidate feature** (the limiting count is the minority class). Report
  the ratio honestly. Features ≥ events with no reduction → `HIGH_DIM_LOW_EVENTS`.

## 4. Nested cross-validation (the non-negotiable)

If you tune hyperparameters and report performance from the **same** CV, the performance is
optimistic. Two acceptable designs:

- **Nested CV** — outer folds estimate performance; an inner CV inside each outer training fold tunes
  hyperparameters **and** does feature selection + scaling. Nothing from the outer test fold touches
  fitting.
- **Held-out test set** — tune with CV on the training split, evaluate once on an untouched test split.

**Everything data-driven goes inside the fold**: imputation, scaling, feature selection, class-imbalance
resampling. Selection on the whole dataset → `SELECTION_OUTSIDE_CV` (and, in prose,
`self-review/check_cv_leakage`). Flat CV or no validation → `NO_NESTED_CV`.

```python
# sklearn nested-CV skeleton (integrate, don't reimplement)
from sklearn.pipeline import Pipeline
from sklearn.feature_selection import SelectKBest
from sklearn.model_selection import GridSearchCV, StratifiedKFold, cross_val_score
pipe = Pipeline([("scale", StandardScaler()),
                 ("select", SelectKBest(k=20)),        # selection INSIDE the fold
                 ("clf", XGBClassifier(eval_metric="logloss"))])
inner = StratifiedKFold(5, shuffle=True, random_state=42)
outer = StratifiedKFold(5, shuffle=True, random_state=42)
grid  = GridSearchCV(pipe, param_grid, cv=inner, scoring="roc_auc")
auc   = cross_val_score(grid, X, y, cv=outer, scoring="roc_auc")  # outer = unbiased estimate
```

## 5. Models — the full classical family (not just RF / XGBoost)

Always report a **simple baseline** (penalised logistic) alongside any complex learner — a model that
does not beat penalised logistic on your data is a finding, not a failure. Fix the seed. Pick by task
and sample size; the pipeline rigor in §2–§4 is identical across all of them (the gate is
learner-agnostic).

| Family | Learner (scikit-learn / library) | When |
|---|---|---|
| Penalised regression | `LogisticRegression(penalty=l1/l2/elasticnet)` — LASSO / ridge / elastic-net | baseline; small n; interpretable coefficients |
| Margin / kernel | `SVC` (linear / RBF) | moderate n, clear margin; scale features first |
| Instance-based | `KNeighborsClassifier` | small feature set, local structure |
| Probabilistic / discriminant | `GaussianNB`, `LinearDiscriminantAnalysis`, `QuadraticDiscriminantAnalysis` | fast baselines; LDA when classes ~Gaussian |
| Single tree | `DecisionTreeClassifier` | interpretability demo (rarely final) |
| Bagging | `RandomForestClassifier`, `ExtraTreesClassifier` | robust default, low tuning |
| Boosting | `XGBClassifier`, `LGBMClassifier`, `CatBoostClassifier`, `HistGradientBoostingClassifier`, `AdaBoost` | usually top tabular performance; tune with the inner CV |
| Shallow neural | `MLPClassifier` | non-linear, enough samples; scale + early-stop |
| Meta / ensemble | `StackingClassifier`, `VotingClassifier` | squeeze marginal gain; guard overfitting |
| Survival ML | random survival forest, Cox-net, DeepSurv | time-to-event outcome (+ `/analyze-stats` survival) |

**Unsupervised (upstream, not the endpoint):** PCA / UMAP for dimensionality reduction (fit inside the
fold, §4), and k-means / hierarchical / Gaussian-mixture for phenotype discovery — report cluster
stability, never as a supervised performance claim.

For imbalanced outcomes prefer probability-calibrated models + threshold analysis over resampling that
distorts prevalence; calibrate with Platt / isotonic (`/analyze-stats` calibration).

## 6. Report discrimination AND calibration AND utility

- **Discrimination** — AUROC (+ AUPRC at low prevalence) with bootstrap CIs.
- **Calibration** — slope + intercept and a flexible calibration curve (not decile bins); use the
  `/analyze-stats` calibration guide. Discrimination-only → `NO_CALIBRATION`.
- **Clinical utility** — decision-curve net benefit / NNT at a stated threshold (`/analyze-stats`).
- **Interpretation** — SHAP for feature contributions (global + a few local), framed as association.

## 7. Validation and reporting

- **External / temporal validation** for a clinical claim; single-cohort development → report as
  development-only (`NO_EXTERNAL_VALIDATION`).
- **Reporting** — CLEAR (radiomics), TRIPOD+AI (prediction model), PROBAST-AI (risk of bias) via
  `/check-reporting`.

## 8. Common reviewer objections this pre-empts

1. "Was performance from nested CV or the same folds you tuned on?" → §4.
2. "Features ≥ patients — how did you avoid overfitting?" → §3.
3. "Were features selected inside the CV?" → §4.
4. "Are the features reproducible (ICC)?" → §2.
5. "Calibration, not just AUC?" → §6.
6. "External validation?" → §7.

#!/usr/bin/env python3
"""Radiomics / classical-ML pipeline-rigor gate (radiomics-ml).

Radiomics + tree-ensemble studies (features -> random forest / XGBoost -> clinical
outcome) are the most common solo-doable clinical-ML workflow, and the most commonly
over-optimistic: hundreds-to-thousands of features on tens of patients, hyperparameters
tuned on the same folds performance is reported from, features selected on the whole
dataset, unstable features never filtered, and discrimination reported without
calibration (Lambin 2017; Park & Kim radiomics-quality; CLEAR; TRIPOD+AI; PROBAST-AI).

This gate reads a declarative **pipeline manifest** (JSON — the artifact this skill
emits, or one the researcher writes) and decides each rigor requirement by rule. It
complements `self-review/check_cv_leakage` (which greps a finished manuscript's prose):
this one audits the pipeline spec at build time.

CHECKS (verdicts):
  1. NO_NESTED_CV          (Major)  hyperparameters are tuned and performance reported on
                                    the same CV (flat k-fold) or with no validation at all;
                                    nested CV or a held-out test set is required.
  2. HIGH_DIM_LOW_EVENTS   (Major)  at least as many features as events (p >= events) with no
                                    dimensionality reduction / regularisation — the classic
                                    radiomics overfitting trap.
  3. SELECTION_OUTSIDE_CV  (Major)  feature selection is fit outside the CV fold (on the whole
                                    dataset), leaking the held-out folds into selection.
  4. NO_FEATURE_STABILITY  (Minor)  no test-retest / ICC feature-stability filtering; radiomics
                                    features are notoriously unstable across acquisition.
  5. NO_CALIBRATION        (Minor)  a clinical prediction model reported by discrimination only
                                    (no calibration slope/intercept or flexible curve).
  6. NO_EXTERNAL_VALIDATION(Minor)  single cohort, no external / temporal validation for a
                                    clinical claim.

MANIFEST (JSON)
  {
    "task": "classification",
    "n_features": 1200,
    "n_samples": 140,
    "n_events": 40,                       // minority-class count (for events-per-feature)
    "cv_scheme": "nested",                // nested / single_split / held_out_test / flat / loocv / none
    "feature_selection_stage": "inside_cv", // inside_cv / outside_cv / none
    "dimensionality_reduction": true,     // LASSO / PCA / regularisation applied
    "feature_stability": "icc",           // icc / test_retest / none
    "calibration_reported": true,
    "external_validation": "temporal",    // external / temporal / none
    "model": "xgboost"
  }

INPUTS
  --manifest  radiomics/classical-ML pipeline manifest JSON (required).

OUTPUT
  A reconciliation table (stdout) and, with --out, a JSON artifact:
    {manifest, model, n_features, n_samples, n_events, cv_scheme, claims[...], summary}
  NO_NESTED_CV / HIGH_DIM_LOW_EVENTS / SELECTION_OUTSIDE_CV are Major.

Stdlib-only (json / argparse / pathlib). Exit codes: 0 clean (or report-only),
1 Major claim(s) found (with --strict), 2 input/usage error.
"""

from __future__ import annotations

import argparse
import json
import sys
from pathlib import Path

VALID_CV = {"nested", "nested_cv", "single_split", "held_out_test", "holdout_test",
            "held-out", "train_test_holdout"}
NONE_VALUES = {"", "none", "no", "na", "n/a", "false", "0"}


def _norm(s) -> str:
    return str(s).strip().lower() if s is not None else ""


def check(m: dict) -> list[dict]:
    claims: list[dict] = []
    n_features = m.get("n_features")
    n_events = m.get("n_events")
    n_samples = m.get("n_samples")
    cv = _norm(m.get("cv_scheme"))
    sel = _norm(m.get("feature_selection_stage"))
    dimred = m.get("dimensionality_reduction")
    stability = _norm(m.get("feature_stability"))
    calib = m.get("calibration_reported")
    extval = _norm(m.get("external_validation"))

    # 1. No nested CV / no held-out validation.
    if cv not in VALID_CV:
        detail = ("no validation scheme (`none`)" if cv in NONE_VALUES
                  else f"flat CV ('{cv or 'missing'}') tunes and reports on the same folds")
        claims.append({
            "verdict": "NO_NESTED_CV", "severity": "Major",
            "detail": (f"{detail}; use nested cross-validation or a held-out test set so tuning "
                       f"does not inflate the reported performance"),
            "where": "cv_scheme",
        })

    # 2. High dimensionality vs events, no reduction.
    denom = n_events if isinstance(n_events, (int, float)) else n_samples
    if isinstance(n_features, (int, float)) and isinstance(denom, (int, float)) and denom > 0:
        if n_features >= denom and dimred is not True:
            unit = "events" if isinstance(n_events, (int, float)) else "samples"
            claims.append({
                "verdict": "HIGH_DIM_LOW_EVENTS", "severity": "Major",
                "detail": (f"{n_features} features vs {int(denom)} {unit} (p >= {unit}) with no "
                           f"dimensionality reduction / regularisation; radiomics overfits badly "
                           f"in this regime — apply LASSO / PCA / a stability+redundancy filter"),
                "where": "n_features",
            })

    # 3. Feature selection outside the CV fold.
    if sel in {"outside_cv", "outside", "whole_dataset", "before_cv", "pre_cv", "global"}:
        claims.append({
            "verdict": "SELECTION_OUTSIDE_CV", "severity": "Major",
            "detail": ("feature selection is fit outside the CV fold (on the whole dataset), so the "
                       "held-out folds leak into selection; nest selection inside each training fold"),
            "where": "feature_selection_stage",
        })

    # 4. No feature-stability filtering.
    if stability in NONE_VALUES:
        claims.append({
            "verdict": "NO_FEATURE_STABILITY", "severity": "Minor",
            "detail": ("no test-retest / ICC feature-stability filtering; radiomics features are "
                       "unstable across acquisition and segmentation — filter to reproducible features"),
            "where": "feature_stability",
        })

    # 5. No calibration.
    if calib is not True:
        claims.append({
            "verdict": "NO_CALIBRATION", "severity": "Minor",
            "detail": ("a clinical prediction model is reported without calibration (slope/intercept "
                       "or a flexible calibration curve), only discrimination"),
            "where": "calibration_reported",
        })

    # 6. No external / temporal validation.
    if extval in NONE_VALUES:
        claims.append({
            "verdict": "NO_EXTERNAL_VALIDATION", "severity": "Minor",
            "detail": ("single-cohort development with no external / temporal validation; a clinical "
                       "claim needs validation beyond the development sample"),
            "where": "external_validation",
        })

    return claims


def analyze(manifest_path: str) -> dict:
    p = Path(manifest_path)
    if not p.is_file():
        sys.stderr.write(f"ERROR: manifest not found: {manifest_path}\n")
        sys.exit(2)
    try:
        m = json.loads(p.read_text(encoding="utf-8"))
    except (json.JSONDecodeError, ValueError) as e:
        sys.stderr.write(f"ERROR: manifest is not valid JSON: {e}\n")
        sys.exit(2)
    if not isinstance(m, dict):
        sys.stderr.write("ERROR: manifest JSON must be an object\n")
        sys.exit(2)

    claims = check(m)
    n_major = sum(1 for c in claims if c["severity"] == "Major")
    return {
        "manifest": str(p),
        "model": m.get("model"),
        "n_features": m.get("n_features"),
        "n_samples": m.get("n_samples"),
        "n_events": m.get("n_events"),
        "cv_scheme": m.get("cv_scheme"),
        "claims": claims,
        "summary": {
            "n_claims": len(claims),
            "n_major": n_major,
            "n_flag": len(claims) - n_major,
            "verdict": "MAJOR_CANDIDATE" if n_major else "OK",
        },
    }


def render(result: dict) -> str:
    lines = ["| Check | Severity | Detail |", "|---|---|---|"]
    for c in result["claims"]:
        lines.append(f"| {c['verdict']} | {c['severity']} | {c['detail']} |")
    if len(lines) == 2:
        lines.append("| (none) | — | radiomics/ML pipeline meets the rigor bar |")
    return "\n".join(lines)


def main() -> int:
    ap = argparse.ArgumentParser(description="Radiomics / classical-ML pipeline-rigor gate.")
    ap.add_argument("--manifest", required=True, help="radiomics/classical-ML pipeline manifest JSON")
    ap.add_argument("--out", help="write JSON artifact to this path")
    ap.add_argument("--strict", action="store_true", help="exit 1 if any Major claim exists")
    ap.add_argument("--quiet", action="store_true", help="suppress stdout table")
    args = ap.parse_args()

    result = analyze(args.manifest)

    if not args.quiet:
        print("=" * 41)
        print(" Radiomics / Classical-ML Gate (radiomics-ml)")
        print("=" * 41)
        print(f"  model={result['model']}  n_features={result['n_features']}  "
              f"n_samples={result['n_samples']}  n_events={result['n_events']}  "
              f"cv_scheme={result['cv_scheme']}")
        print(render(result))
        print()
        s = result["summary"]
        if s["n_major"]:
            print(f"MAJOR candidate: {s['n_major']} radiomics/ML rigor issue(s).")
        elif s["n_flag"]:
            print(f"MINOR flag: {s['n_flag']} radiomics/ML rigor issue(s) (see table).")
        else:
            print("OK: radiomics/ML pipeline meets the rigor bar.")

    if args.out:
        Path(args.out).parent.mkdir(parents=True, exist_ok=True)
        Path(args.out).write_text(json.dumps({"detector": "check_radiomics_ml", **result}, indent=2), encoding="utf-8")
        if not args.quiet:
            print(f"\nwrote {args.out}")

    return 1 if (args.strict and result["summary"]["n_major"]) else 0


if __name__ == "__main__":
    sys.exit(main())

#!/usr/bin/env bash
# Regression test for the radiomics/classical-ML pipeline-rigor gate (radiomics-ml).
# Synthetic, PII-free JSON manifests reproduce each verdict class + the suppressions.
# Stdlib-only (python3).
set -u

HERE="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
SCRIPT="$HERE/../scripts/check_radiomics_ml.py"
CH="$HERE/../scripts/check_radiomics_ml_challenge"
TMP="$(mktemp -d -t radml_XXXX)"
OUT="$TMP/out.json"
trap 'rm -rf "$TMP"' EXIT

fail=0
check() { local label="$1"; shift
    if "$@" >/dev/null 2>&1; then printf '  PASS  %s\n' "$label"
    else printf '  FAIL  %s\n' "$label"; fail=$((fail+1)); fi
}
has_verdict() { python3 -c "
import json,sys
d=json.load(open('$OUT'))
assert any(c['verdict']=='$1' for c in d['claims']), '$1 not found'
"; }
no_verdict() { python3 -c "
import json,sys
d=json.load(open('$OUT'))
assert not any(c['verdict']=='$1' for c in d['claims']), '$1 unexpectedly present'
"; }

[[ -f "$SCRIPT" ]] || { echo "ENV-ERR: script missing" >&2; exit 2; }

# (1) weak fixture -> 3 Major + exit 1
python3 "$SCRIPT" --manifest "$CH/fixture/pipeline_weak.json" --out "$OUT" --strict --quiet >/dev/null 2>&1
check "exit 1 (weak pipeline)" test "$?" -eq 1
check "NO_NESTED_CV detected" has_verdict NO_NESTED_CV
check "HIGH_DIM_LOW_EVENTS detected" has_verdict HIGH_DIM_LOW_EVENTS
check "SELECTION_OUTSIDE_CV detected" has_verdict SELECTION_OUTSIDE_CV
check "NO_FEATURE_STABILITY detected" has_verdict NO_FEATURE_STABILITY
check "NO_CALIBRATION detected" has_verdict NO_CALIBRATION
check "NO_EXTERNAL_VALIDATION detected" has_verdict NO_EXTERNAL_VALIDATION

# (2) strong fixture -> exit 0, no claims
python3 "$SCRIPT" --manifest "$CH/fixture/pipeline_strong.json" --strict --quiet >/dev/null 2>&1
check "exit 0 (strong pipeline)" test "$?" -eq 0

# (3) dimensionality_reduction=true suppresses HIGH_DIM_LOW_EVENTS even when p >= events
cat > "$TMP/dimred.json" <<'EOF'
{"n_features": 1000, "n_events": 30, "cv_scheme": "nested", "feature_selection_stage": "inside_cv",
 "dimensionality_reduction": true, "feature_stability": "icc", "calibration_reported": true,
 "external_validation": "external", "model": "lasso_logistic"}
EOF
python3 "$SCRIPT" --manifest "$TMP/dimred.json" --out "$OUT" --quiet >/dev/null 2>&1
check "dim-reduction suppresses HIGH_DIM_LOW_EVENTS" no_verdict HIGH_DIM_LOW_EVENTS

# (4) single_split is an acceptable validation scheme (no NO_NESTED_CV)
cat > "$TMP/single.json" <<'EOF'
{"n_features": 50, "n_events": 120, "cv_scheme": "single_split", "feature_selection_stage": "inside_cv",
 "dimensionality_reduction": true, "feature_stability": "icc", "calibration_reported": true,
 "external_validation": "temporal", "model": "random_forest"}
EOF
python3 "$SCRIPT" --manifest "$TMP/single.json" --out "$OUT" --quiet >/dev/null 2>&1
check "single_split does NOT fire NO_NESTED_CV" no_verdict NO_NESTED_CV
python3 "$SCRIPT" --manifest "$TMP/single.json" --strict --quiet >/dev/null 2>&1
check "exit 0 on rigorous single-split pipeline" test "$?" -eq 0

# (5) cv_scheme none -> NO_NESTED_CV (no validation at all)
cat > "$TMP/nocv.json" <<'EOF'
{"n_features": 20, "n_events": 200, "cv_scheme": "none", "feature_selection_stage": "inside_cv",
 "dimensionality_reduction": true, "feature_stability": "icc", "calibration_reported": true,
 "external_validation": "external", "model": "xgboost"}
EOF
python3 "$SCRIPT" --manifest "$TMP/nocv.json" --out "$OUT" --quiet >/dev/null 2>&1
check "cv_scheme=none fires NO_NESTED_CV" has_verdict NO_NESTED_CV

# (6) n_events absent -> falls back to n_samples for the dimensionality check
cat > "$TMP/samplesonly.json" <<'EOF'
{"n_features": 500, "n_samples": 60, "cv_scheme": "nested", "feature_selection_stage": "inside_cv",
 "dimensionality_reduction": false, "feature_stability": "icc", "calibration_reported": true,
 "external_validation": "external", "model": "random_forest"}
EOF
python3 "$SCRIPT" --manifest "$TMP/samplesonly.json" --out "$OUT" --quiet >/dev/null 2>&1
check "HIGH_DIM_LOW_EVENTS uses n_samples fallback when n_events missing" has_verdict HIGH_DIM_LOW_EVENTS

# (7) the shipped challenge card passes
check "challenge verify.sh passes" bash "$CH/verify.sh"

echo "fail=$fail"; [[ "$fail" -eq 0 ]] && echo "ALL PASS" || echo "FAILURES: $fail"
exit "$fail"