#!/usr/bin/env python3 # --------------------------------------------------------------------------- # PUBLISHED COPY — supplementary material for: # "Does BMI Predict All 11 Food Addiction Criteria?" # https://creativetouchrotherham.co.uk/lab/insights/food-addiction-criteria-bmi # # Creative Touch, Rotherham. Released under CC BY 4.0 # (https://creativecommons.org/licenses/by/4.0/), the same licence as the # article and its figures. # # Published verbatim. This is the file as it was committed and run; no line has # been changed for publication. # # It reads the analytic dataset produced by `extract.py`, which is derived from # individual quiz responses. Those responses are not published, so this script # cannot be re-run against our data by a third party. It is published so that # every model specification, every pre-specified decision rule and every logged # deviation can be checked against the analysis plan and against the numbers # printed in the article. # --------------------------------------------------------------------------- """ YFAS 2.0 criterion-level study — pre-registered modelling (step 2 of 2). Implements the committed pre-registration `../02-statistical-analysis-plan.md` (SAP v1.0, committed 2026-07-26). Sections implemented, in the order run: * SAP §5 — a priori power table, per criterion, minimum detectable odds ratio per 5 BMI units at 80% power, alpha 0.05 two-sided, using the analytic cohort's own endorsement rates. [A] * SAP §2.2 — THE primary confirmatory GEE. One contrast, one p-value, one confidence interval. H1 lives or dies on it. [B] * SAP §2.3 — secondary per-criterion GEE, eleven BMI slopes, BH-FDR q=0.05. Secondary and descriptive. Not a test of H1. [C] * SAP §2.4 — unconditioned per-criterion logistic regressions. Descriptive only, retained as the didactic contrast against [C]. [D] * Locked decision "BMI continuous primary, linearity tested not assumed" — linearity of the BMI logit tested, never used to switch the primary specification. [E] * SAP §3 — IRT/DIF convergent-methods check. 2PL on the eleven binary criteria, then Swaminathan-Rogers logistic-regression DIF. [F] * SAP §4 — H2 discordant-phenotype confirmatory secondary, with the mandatory total_symptom_count adjustment and an exact power recomputation against the real cell sizes. Fitted AS AMENDED by SAP §9 'Deviation 1' — see below. [G] * SAP §1 — H3 (age, exploratory, no direction) and H4 (sex, C9-specific replication attempt of Saffari et al. 2022). [H] * Per-criterion heterogeneity within the H1 criterion sets — each criterion's conditioned slope against its set's fitted direction. Structural description, no interpretation. [I] --------------------------------------------------------------------------- SAP §9 DEVIATION 1 (2026-07-26) — the one amendment this script implements SAP §4's H2 model, as originally pre-registered, does not converge on the real analytic sample: the outcome is compositional within respondent, within-cluster dependence is negative by construction, and the exchangeable dependence parameter diverges. The project owner authorised changing the working correlation to `Independence()` AND NOTHING ELSE — same GEE, same logit/binomial family, same formula including the mandatory `total_symptom_count` adjustment, same cluster-robust sandwich SEs on `respondent_id`, same population-averaged estimand. The full record, including the disclosure that the ruling was made with the Independence result already visible, is at `02-statistical-analysis- plan.md` §9. [G] re-derives the original failure on every run and prints it as the audit trail, then reports the amended fit under an inline amendment notice. No other model in this script is affected. THE SAP IS BINDING. Every model below is specified in advance. This script is a faithful implementation, not a place for statistical judgement. Where the SAP proved ambiguous, the reading taken is stated in-line and reported in the output under "SPECIFICATION NOTES" rather than resolved silently. --------------------------------------------------------------------------- DATA SOURCE — no database connection The frozen snapshots written by `extract.py` are the analytic sample by design (SAP §6 and §8). This script reads them and never contacts MariaDB, so the freeze cannot be re-opened by re-running the analysis: snapshots/analytic_wide_latest.csv 1,796 rows, one per respondent snapshots/analytic_long_latest.csv 19,756 rows, respondent x criterion `rest_score` is taken as given from the long snapshot (verified there as total_symptom_count - criterion_met, range 0-10). It is not recomputed. The q1..q35 columns in the wide snapshot are UNUSABLE and are never read. Item responses were recorded if and only if the criterion was met, so their missingness is deterministic on the outcome. Every model here runs on the eleven binary criteria, exactly as SAP §2 and §3 specify. --------------------------------------------------------------------------- OPERATIONALISATION DECISIONS (coordinator rulings, recorded here as SAP §3 and §6 require implementation choices to be recorded in the script header) 1. AGE. `age_range` is banded (18-24 ... 75+), not continuous. The SAP writes "age" without specifying a functional form, so it is modelled as a categorical `C(age_range)` throughout. This is assumption-light and faithful: it is an implementation choice, not a §9 deviation. Band midpoints are NOT invented. 2. BMI SCALING. BMI enters every model on its natural (per-1-unit) scale. SAP §5 expresses effects "per 5-unit BMI increase", so every odds ratio printed by this script is exponentiated at 5 units and labelled "per 5 BMI units" on the line where it appears. Per-1-unit log-coefficients are printed alongside so the two scales can never be confused. 3. IRT/DIF METHOD (SAP §3 requires this choice be recorded here). `girth` 0.8.0 is available in `~/yfas-env`, so the arm gets a genuine unidimensional 2PL (`girth.twopl_mml`) with EAP ability estimates (`girth.ability_eap`). DIF is then Swaminathan-Rogers logistic regression conditioned on that ability estimate: uniform DIF = the group main effect, non-uniform DIF = the group x ability interaction, both by likelihood-ratio test. Grouping is BMI split at the analytic cohort's median, with a continuous-BMI sensitivity check as SAP §3 permits. The Mantel-Haenszel fallback in SAP §3 is therefore NOT used and no §9 deviation arises from this arm. --------------------------------------------------------------------------- RESULT-BLIND DEBUGGING (SAP §8 no-peeking) `--permute SEED` shuffles the `bmi` column across respondents before anything is fitted, keeping every other column intact. That fits the real model structure against a null association, so the script could be debugged to convergence without any real coefficient being seen. Permuted runs are stamped PERMUTED in the banner, the footer, every output filename and the JSON payload, so a permuted run can never be mistaken for the real one. There is no other source of randomness in this script. LIMIT OF THE BLIND, stated honestly: permuting BMI blinds only the BMI-dependent output — [A], [B], [C], [D], [E], the BMI arm of [F], and [G]. H3 (age) and H4 (sex DIF) do not involve BMI, so a permuted run reproduces their real values exactly. Those two sections were therefore not blinded during development. They are exploratory (H3) and a pre-registered directional replication attempt whose verdict is mechanical (H4), so no analytic choice in either could be steered by having seen them; the fact is recorded here rather than glossed. --------------------------------------------------------------------------- Reproducibility — resolved package versions in the WSL venv `~/yfas-env`: Python 3.10.12 pandas 2.3.3 numpy 2.2.6 scipy 1.15.3 statsmodels 0.14.6 girth 0.8.0 patsy 1.0.2 Run as: wsl bash -lc "~/yfas-env/bin/python \\ /mnt/c/xampp/htdocs/creative-touch/modules/calculators/_docs/phase-3/analysis/analyse.py" British English throughout. """ from __future__ import annotations import argparse import json import sys import warnings from datetime import datetime from pathlib import Path from typing import Any import numpy as np import pandas as pd import patsy import statsmodels.api as sm import statsmodels.formula.api as smf from scipy import stats as sp_stats from statsmodels.stats.multitest import multipletests from statsmodels.stats.power import NormalIndPower, normal_power from statsmodels.stats.proportion import proportion_effectsize # --------------------------------------------------------------------------- # Constants — pre-registered, not configurable # --------------------------------------------------------------------------- CRITERIA = tuple(range(1, 12)) CONTENT_CONFOUNDED = (4, 5, 8, 9, 10) PHARMACOLOGICAL_CORE = (6, 7, 11) NEITHER = (1, 2, 3) #: SAP §2.2 — the confirmatory model is fitted on these two sets only. CONFIRMATORY_SETS = ("content_confounded", "pharmacological_core") #: SAP §2.2 — reference level, so the contrast reads "content-confounded minus core". SET_REFERENCE = "pharmacological_core" CRITERION_LABELS = { 1: "C1 larger amounts / longer than intended", 2: "C2 persistent desire / failed cut-down", 3: "C3 great deal of time spent", 4: "C4 activities given up", 5: "C5 use despite knowledge of harm", 6: "C6 tolerance", 7: "C7 withdrawal", 8: "C8 interpersonal problems", 9: "C9 role obligation failure", 10: "C10 hazardous use", 11: "C11 craving", } #: Verified sample facts (00-STATE.md "Step 3a results"). Asserted, not re-derived. EXPECTED_N_DESCRIPTIVE = 1796 EXPECTED_N_PRIMARY = 1773 EXPECTED_N_LONG = 19756 EXPECTED_H2_NORMAL_HIGH = 155 EXPECTED_H2_OBESE_HIGH = 394 EXPECTED_H2_OBESE_LOW = 98 ALPHA = 0.05 POWER_TARGET = 0.80 FDR_Q = 0.05 BMI_UNIT = 5.0 # SAP §5 reporting scale BASE_DIR = Path(__file__).resolve().parent SNAPSHOT_DIR = BASE_DIR / "snapshots" RESULTS_DIR = BASE_DIR / "results" # --------------------------------------------------------------------------- # Infrastructure # --------------------------------------------------------------------------- class SapDeparture(RuntimeError): """Raised when an assumption underpinning the SAP fails. Any occurrence requires an entry in SAP §9 (Deviations Log) before the analysis may proceed. We abort rather than continue with numbers that would look publishable but not be. """ def require(condition: bool, message: str) -> None: """Explicit, non-strippable assertion (survives `python -O`).""" if not condition: raise SapDeparture(message) class Logger: """Tee to stdout and to the timestamped results file.""" def __init__(self, path: Path) -> None: self.path = path self._lines: list[str] = [] def __call__(self, message: str = "") -> None: print(message) self._lines.append(message) def rule(self, title: str = "") -> None: self("") self("=" * 100) if title: self(title) self("=" * 100) def sub(self, title: str) -> None: self("") self("-" * 100) self(title) self("-" * 100) def flush(self) -> None: self.path.write_text("\n".join(self._lines) + "\n", encoding="utf-8") def jsonable(obj: Any) -> Any: if isinstance(obj, (np.integer,)): return int(obj) if isinstance(obj, (np.floating,)): return float(obj) if isinstance(obj, (np.bool_,)): return bool(obj) if isinstance(obj, np.ndarray): return obj.tolist() return str(obj) def criterion_set_of(k: int) -> str: if k in CONTENT_CONFOUNDED: return "content_confounded" if k in PHARMACOLOGICAL_CORE: return "pharmacological_core" return "neither" def or_per_5(beta: float, se: float) -> dict[str, float]: """Odds ratio per 5 BMI units from a per-1-unit log-odds coefficient. Every OR this script prints goes through here, so the per-1 / per-5 mix-up the coordinator warned about cannot happen silently. """ z = sp_stats.norm.ppf(1.0 - ALPHA / 2.0) b5 = beta * BMI_UNIT se5 = se * BMI_UNIT return { "beta_per_1": float(beta), "se_per_1": float(se), "beta_per_5": float(b5), "se_per_5": float(se5), "or_per_5": float(np.exp(b5)), "ci_lo_per_5": float(np.exp(b5 - z * se5)), "ci_hi_per_5": float(np.exp(b5 + z * se5)), } def check_converged(res: Any, label: str, log: Logger) -> None: """Assert convergence and say so. A silently non-converged GEE is the worst failure mode available here, so this aborts rather than warns.""" converged = getattr(res, "converged", None) if converged is None: mle = getattr(res, "mle_retvals", None) converged = None if mle is None else mle.get("converged") if converged is not None: converged = bool(converged) log(f" convergence [{label}]: {converged}") require( converged is True, f"MODEL DID NOT CONVERGE: {label} (converged={converged!r}). Aborting — " "a non-converged fit must never be reported as a result.", ) def assert_full_rank(formula: str, data: pd.DataFrame, label: str, log: Logger) -> dict[str, Any]: """SAP §2.2 requires the design-matrix rank be asserted in the script rather than trusting the formula parser to drop the right term.""" y, X = patsy.dmatrices(formula, data, return_type="dataframe") rank = int(np.linalg.matrix_rank(np.asarray(X, dtype=float))) ncols = int(X.shape[1]) log(f" design matrix [{label}]: {X.shape[0]} rows x {ncols} cols, rank = {rank}") require( rank == ncols, f"DESIGN MATRIX IS RANK-DEFICIENT for {label}: rank {rank} < {ncols} columns. " "SAP §2.2 requires an explicit rank assertion here; a rank-deficient fit " "would silently alias the interaction contrast. Aborting.", ) return {"n_rows": int(X.shape[0]), "n_cols": ncols, "rank": rank, "columns": list(X.columns)} # --------------------------------------------------------------------------- # Data loading # --------------------------------------------------------------------------- def load_data(permute_seed: int | None, log: Logger) -> tuple[pd.DataFrame, pd.DataFrame]: wide_path = SNAPSHOT_DIR / "analytic_wide_latest.csv" long_path = SNAPSHOT_DIR / "analytic_long_latest.csv" require(wide_path.exists(), f"Missing snapshot: {wide_path}") require(long_path.exists(), f"Missing snapshot: {long_path}") # q1..q35 are structurally unusable (missingness deterministic on the # outcome) and are deliberately not loaded at all. wide = pd.read_csv(wide_path, usecols=lambda c: not c.startswith("q")) long = pd.read_csv(long_path) log(f" wide snapshot: {wide_path.name} {len(wide)} rows x {len(wide.columns)} cols " "(q1..q35 deliberately not loaded)") log(f" long snapshot: {long_path.name} {len(long)} rows x {len(long.columns)} cols") require(len(wide) == EXPECTED_N_DESCRIPTIVE, f"Descriptive n is {len(wide)}, expected {EXPECTED_N_DESCRIPTIVE}.") require(len(long) == EXPECTED_N_LONG, f"Long n is {len(long)}, expected {EXPECTED_N_LONG}.") require(int(wide["primary_model_eligible"].sum()) == EXPECTED_N_PRIMARY, f"Primary-eligible n is {int(wide['primary_model_eligible'].sum())}, " f"expected {EXPECTED_N_PRIMARY}.") require(int(wide["respondent_id"].duplicated().sum()) == 0, "Duplicate respondent_id in the wide snapshot.") # rest_score is taken as given (verified in extraction), only range-checked. require(bool(long["rest_score"].between(0, 10).all()), "rest_score outside the SAP §2.1 range 0-10 in the long snapshot.") require(bool((long["rest_score"] == long["total_symptom_count"] - long["criterion_met"]).all()), "rest_score != total_symptom_count - criterion_met in the long snapshot.") # criterion_set vocabulary check against SAP §2.2. for k in CRITERIA: sets = set(long.loc[long["criterion_id"] == k, "criterion_set"].unique()) require(sets == {criterion_set_of(k)}, f"criterion {k} carries criterion_set {sets}, expected " f"{criterion_set_of(k)!r} per SAP §2.2.") if permute_seed is not None: rng = np.random.default_rng(permute_seed) order = rng.permutation(len(wide)) wide = wide.copy() wide["bmi"] = wide["bmi"].to_numpy()[order] wide["bmi_who_category"] = [who_category(b) for b in wide["bmi"]] bmi_map = dict(zip(wide["respondent_id"], wide["bmi"])) long = long.copy() long["bmi"] = long["respondent_id"].map(bmi_map) long["bmi_who_category"] = long["bmi"].map(who_category) log("") log(f" *** BMI PERMUTED across respondents, seed = {permute_seed}. " "THESE ARE NOT REAL RESULTS. ***") # Formulas below are written to read exactly as the SAP writes them, so the # model frames carry a `BMI` column alias alongside the snapshot's `bmi`. wide["BMI"] = wide["bmi"] long["BMI"] = long["bmi"] long = long.sort_values(["respondent_id", "criterion_id"]).reset_index(drop=True) elig = dict(zip(wide["respondent_id"], wide["primary_model_eligible"])) long["primary_model_eligible"] = long["respondent_id"].map(elig) require(long["primary_model_eligible"].notna().all(), "Some long rows have no matching respondent in the wide snapshot.") return wide, long WHO_CATEGORIES = ( ("underweight", None, 18.5), ("normal", 18.5, 25.0), ("overweight", 25.0, 30.0), ("obese_i", 30.0, 35.0), ("obese_ii", 35.0, 40.0), ("obese_iii", 40.0, None), ) def who_category(bmi: float) -> str: for label, low, high in WHO_CATEGORIES: if (low is None or bmi >= low) and (high is None or bmi < high): return label raise SapDeparture(f"BMI {bmi!r} fell through the WHO category ladder") # --------------------------------------------------------------------------- # [A] SAP §5 — a priori power table # --------------------------------------------------------------------------- def section_a_power(wide: pd.DataFrame, log: Logger) -> dict[str, Any]: log.rule("[A] SAP §5 — A PRIORI POWER, PER CRITERION") log("Minimum detectable odds ratio per 5 BMI units at 80% power, alpha = 0.05") log("two-sided, using the analytic cohort's OWN endorsement rates (SAP §5 step 1),") log("adjusted for the covariates already in the model (rest_score, age, sex) via") log("Hsieh's variance-inflation factor 1/(1-R^2) — the 'logistic-regression-") log("appropriate power approximation' SAP §5 step 2 permits and prefers over a") log("bare two-proportion approximation. Every returned OR is verified back") log("through statsmodels.stats.power.normal_power to reproduce 80% power.") log("") log(f"n used for the power table = {len(wide)} (SAP §5 specifies the analytic cohort).") log("") n = len(wide) z_a = sp_stats.norm.ppf(1.0 - ALPHA / 2.0) z_b = sp_stats.norm.ppf(POWER_TARGET) sd_bmi = float(wide["BMI"].std(ddof=1)) log(f"BMI SD in the analytic cohort = {sd_bmi:.4f}") log("") log(f" {'crit':<5}{'set':<21}{'rate %':>8}{'R^2 adj':>9}{'SE(logOR/1)':>13}" f"{'min det OR/5':>14}{'power check':>13}") rows: list[dict[str, Any]] = [] for k in CRITERIA: p = float(wide[f"c{k}"].mean()) rest = wide["total_symptom_count"] - wide[f"c{k}"] frame = pd.DataFrame({ "BMI": wide["BMI"].to_numpy(), "rest_score": rest.to_numpy(), "age_range": wide["age_range"].to_numpy(), "sex": wide["sex"].to_numpy(), }) # R^2 of BMI on the other covariates — a descriptive quantity, not a test. r2 = float(smf.ols("BMI ~ rest_score + C(age_range) + sex", data=frame).fit().rsquared) se_per_1 = 1.0 / np.sqrt(n * p * (1.0 - p) * (sd_bmi ** 2) * (1.0 - r2)) mdes_log_per_1 = (z_a + z_b) * se_per_1 # Verification through statsmodels: power at the returned effect must be 0.80. achieved = float(normal_power(mdes_log_per_1 / se_per_1, 1, ALPHA, alternative="two-sided")) require(abs(achieved - POWER_TARGET) < 1e-6, f"Power verification failed for C{k}: {achieved} != {POWER_TARGET}") est = or_per_5(mdes_log_per_1, se_per_1) cset = criterion_set_of(k) log(f" C{k:<4}{cset:<21}{100 * p:>8.1f}{r2:>9.4f}{se_per_1:>13.5f}" f"{est['or_per_5']:>14.3f}{achieved:>13.3f}") rows.append({ "criterion_id": k, "label": CRITERION_LABELS[k], "criterion_set": cset, "endorsement_rate_pct": round(100 * p, 2), "r2_bmi_on_covariates": round(r2, 5), "se_log_or_per_1": se_per_1, "min_detectable_or_per_5": est["or_per_5"], "verified_power": achieved, }) lowest = sorted(rows, key=lambda r: r["endorsement_rate_pct"])[:2] log("") log("SAP §5 step 3 — REQUIRED LIMITATION STATEMENT, reported regardless of outcome:") for r in lowest: log(f" C{r['criterion_id']} ({r['endorsement_rate_pct']}% endorsement, " f"{r['criterion_set']}) is among the two lowest-powered criteria; " f"min detectable OR per 5 BMI units = {r['min_detectable_or_per_5']:.3f}") flagged = {r["criterion_id"] for r in lowest} log(f" Both lowest-powered criteria are in the content-confounded set: " f"{all(criterion_set_of(k) == 'content_confounded' for k in flagged)}") log(" SAP §5 pre-flagged C8 and C9 specifically; this is confirmed above and must") log(" be stated as a design limitation in the write-up whatever the result.") return {"n_used": n, "bmi_sd": sd_bmi, "criteria": rows, "two_lowest_powered": sorted(flagged)} # --------------------------------------------------------------------------- # [B] SAP §2.2 — the primary confirmatory model # --------------------------------------------------------------------------- PRIMARY_FORMULA = ( "criterion_met ~ C(criterion_id) " f"+ BMI:C(criterion_set, Treatment(reference='{SET_REFERENCE}')) " "+ rest_score + C(age_range) + sex" ) #: Algebraically identical reparameterisation (same column space, same fit) in #: which the interaction coefficient IS directly "content-confounded minus core". #: Used only to verify the contrast computed from PRIMARY_FORMULA. PRIMARY_FORMULA_REPARAM = ( "criterion_met ~ C(criterion_id) + BMI " f"+ BMI:C(criterion_set, Treatment(reference='{SET_REFERENCE}')) " "+ rest_score + C(age_range) + sex" ) def fit_gee(formula: str, data: pd.DataFrame, label: str, log: Logger, cov_struct: Any | None = None, allow_failure: bool = False) -> Any: """SAP §2.2 / §4 GEE: logit link, Binomial family, Exchangeable working correlation, clustered on respondent_id. `cov_struct` is a parameter ONLY so that a non-pre-registered structure can be fitted under an explicit quarantine label in [G]; it defaults to the pre-registered Exchangeable everywhere else and is never varied silently. """ model = sm.GEE.from_formula( formula, groups="respondent_id", data=data, family=sm.families.Binomial(), cov_struct=cov_struct if cov_struct is not None else sm.cov_struct.Exchangeable(), ) with warnings.catch_warnings(): warnings.simplefilter("ignore") res = model.fit(maxiter=200) if allow_failure: converged = bool(getattr(res, "converged", False)) log(f" convergence [{label}]: {converged}") return res if converged else None check_converged(res, label, log) return res def section_b_primary(long: pd.DataFrame, log: Logger) -> dict[str, Any]: log.rule("[B] SAP §2.2 — PRIMARY CONFIRMATORY MODEL (the single test of H1)") df = long[ long["criterion_set"].isin(CONFIRMATORY_SETS) & (long["primary_model_eligible"] == 1) ].copy() n_resp = int(df["respondent_id"].nunique()) log(f"Rows entering the fit: {len(df)} respondents: {n_resp}") log("C1-C3 rows are excluded from the FIT ENTIRELY (SAP §2.2), not merely from") log("interpretation. Sex 'other'/'prefer_not_to_say' excluded per SAP §7.6 —") log("primary model only.") require(n_resp == EXPECTED_N_PRIMARY, f"Primary model respondents = {n_resp}, expected {EXPECTED_N_PRIMARY}.") require(len(df) == 8 * EXPECTED_N_PRIMARY, f"Primary model rows = {len(df)}, expected {8 * EXPECTED_N_PRIMARY} " "(8 criteria x 1,773 respondents).") log(f"criteria in the fit: {sorted(int(k) for k in df['criterion_id'].unique())}") log(f"sex levels in the fit: {sorted(df['sex'].unique())}") log("") log("Model, exactly as SAP §2.2 writes it (criterion_set MAIN EFFECT OMITTED —") log("it is a deterministic function of criterion_id and therefore collinear with") log("the criterion intercepts):") log(f" {PRIMARY_FORMULA}") log(" GEE, logit link, Binomial family, Exchangeable working correlation,") log(" clustered on respondent_id.") log("") rank_info = assert_full_rank(PRIMARY_FORMULA, df, "primary", log) log(" fitted parameter names (so a reader can verify the right term was dropped):") for c in rank_info["columns"]: log(f" {c}") log("") cs_main = [c for c in rank_info["columns"] if "criterion_set" in c and "BMI" not in c] log(f" criterion_set MAIN-EFFECT columns present: {cs_main if cs_main else 'none'}" " <- must be none") require(not cs_main, "A criterion_set main effect entered the design matrix. SAP §2.2 requires " "it be omitted. Aborting.") bmi_set_cols = [c for c in rank_info["columns"] if c.startswith("BMI:")] log(f" BMI x criterion_set columns present: {bmi_set_cols}") require(len(bmi_set_cols) == 2, f"Expected 2 BMI x criterion_set columns, got {bmi_set_cols}.") res = fit_gee(PRIMARY_FORMULA, df, "primary GEE", log) log("") log(str(res.summary())) # --- the confirmatory contrast ----------------------------------------- # # SPECIFICATION NOTE. SAP §2.2 asks that pharmacological_core be the # reference level "so the interaction term is read directly as # content-confounded minus core". With the BMI main effect absent from the # formula (as §2.2 writes it), patsy codes C(criterion_set) inside the # interaction with FULL dummies, so the fit yields two set-specific slopes # rather than one difference. The contrast is therefore formed explicitly # below, and cross-checked against an algebraically identical # reparameterisation in which it IS a single coefficient. Same fit, same # number; the reference level is honoured in the sign convention. # NB the reference level's name appears inside the term string itself # (`Treatment(reference='pharmacological_core')`), so match on the trailing # level bracket rather than on substring presence. names = list(res.params.index) cc_col = [c for c in names if c.startswith("BMI:") and c.endswith("[content_confounded]")] core_col = [c for c in names if c.startswith("BMI:") and c.endswith(f"[{SET_REFERENCE}]")] require(len(cc_col) == 1 and len(core_col) == 1, f"Could not identify the two BMI x set slopes: cc={cc_col} core={core_col}") cc_col, core_col = cc_col[0], core_col[0] contrast = np.zeros(len(names)) contrast[names.index(cc_col)] = 1.0 contrast[names.index(core_col)] = -1.0 ct = res.t_test(contrast) beta = float(np.ravel(ct.effect)[0]) se = float(np.ravel(ct.sd)[0]) pval = float(np.ravel(ct.pvalue)[0]) zval = float(np.ravel(ct.statistic)[0]) est = or_per_5(beta, se) cc_slope = or_per_5(float(res.params[cc_col]), float(res.bse[cc_col])) core_slope = or_per_5(float(res.params[core_col]), float(res.bse[core_col])) # Verification via the reparameterisation. assert_full_rank(PRIMARY_FORMULA_REPARAM, df, "primary (reparameterised)", log) res2 = fit_gee(PRIMARY_FORMULA_REPARAM, df, "primary GEE (reparameterised)", log) rep_col = [c for c in res2.params.index if c.startswith("BMI:") and c.endswith("[T.content_confounded]")] require(len(rep_col) == 1, f"Reparameterised interaction term not found: {rep_col}") rep_beta = float(res2.params[rep_col[0]]) rep_se = float(res2.bse[rep_col[0]]) log("") log(" Contrast verification against the algebraically identical") log(" reparameterisation `BMI + BMI:C(criterion_set)` (same column space):") log(f" contrast from SAP formula beta = {beta:+.6f} se = {se:.6f}") log(f" single coefficient, reparam beta = {rep_beta:+.6f} se = {rep_se:.6f}") require(abs(beta - rep_beta) < 1e-6 and abs(se - rep_se) < 1e-6, "The explicit contrast and the reparameterised coefficient disagree — " "one of the two parameterisations is not what the SAP specifies.") log(" agreement to 1e-6: yes") log.sub("THE CONFIRMATORY TEST — SAP §2.2. One contrast. One p-value. One CI.") log("H1 lives or dies on this line and nothing else in this script.") log("Contrast = (BMI slope, content-confounded) - (BMI slope, pharmacological-core),") log("conditioned on rest_score, with criterion-specific intercepts.") log("NOT subject to multiplicity correction (SAP summary point 7).") log("") log(f" set-specific slope, content_confounded : OR per 5 BMI units = " f"{cc_slope['or_per_5']:.4f} (95% CI {cc_slope['ci_lo_per_5']:.4f} to " f"{cc_slope['ci_hi_per_5']:.4f}) [log-odds per 1 unit = " f"{cc_slope['beta_per_1']:+.6f}]") log(f" set-specific slope, pharmacological_core: OR per 5 BMI units = " f"{core_slope['or_per_5']:.4f} (95% CI {core_slope['ci_lo_per_5']:.4f} to " f"{core_slope['ci_hi_per_5']:.4f}) [log-odds per 1 unit = " f"{core_slope['beta_per_1']:+.6f}]") log("") log(f" CONFIRMATORY CONTRAST ratio of odds ratios, PER 5 BMI UNITS = " f"{est['or_per_5']:.4f}") log(f" 95% CI {est['ci_lo_per_5']:.4f} to " f"{est['ci_hi_per_5']:.4f} (per 5 BMI units)") log(f" z = {zval:+.4f} p = {pval:.6g}") log(f" log-contrast per 1 BMI unit = {beta:+.6f} " f"(se {se:.6f})") log("") direction = "content-confounded > core" if beta > 0 else "content-confounded < core" sig = pval < ALPHA log(f" Direction of the point estimate: {direction}") log(f" Significant at alpha = {ALPHA}: {sig}") log(" H1 predicted a contrast > 1 on the OR-ratio scale (content-confounded") log(" slope steeper than the pharmacological-core slope).") return { "n_rows": len(df), "n_respondents": n_resp, "formula": PRIMARY_FORMULA, "rank_assertion": rank_info, "slope_content_confounded": cc_slope, "slope_pharmacological_core": core_slope, "confirmatory_contrast": { **est, "z": zval, "p_value": pval, "significant_at_0_05": bool(sig), "direction": direction, "scale": "ratio of odds ratios per 5 BMI units", "multiplicity_corrected": False, }, "reparameterisation_check": {"beta": rep_beta, "se": rep_se, "agrees": True}, } # --------------------------------------------------------------------------- # [C] SAP §2.3 — secondary per-criterion model # --------------------------------------------------------------------------- SECONDARY_FORMULA = ( "criterion_met ~ C(criterion_id) + BMI:C(criterion_id) " "+ rest_score + C(age_range) + sex" ) def section_c_secondary(long: pd.DataFrame, log: Logger) -> dict[str, Any]: log.rule("[C] SAP §2.3 — SECONDARY PER-CRITERION MODEL (descriptive, NOT a test of H1)") log("Eleven BMI slopes, one per criterion, rest-score conditioned. Benjamini-") log("Hochberg FDR at q = 0.05 across the eleven tests (SAP §2.3). These do NOT") log("constitute the test of H1 and must not be reported as though eleven small") log("p-values confirm H1 in the absence of the §2.2 contrast.") log("") log("SPECIFICATION NOTE: SAP §7.6 excludes sex 'other'/'prefer_not_to_say' from") log("the PRIMARY model only, and §2.3 does not restate the exclusion, so this") log("secondary model is fitted on the full descriptive cohort (all sex levels).") df = long.copy() n_resp = int(df["respondent_id"].nunique()) log("") log(f"Rows entering the fit: {len(df)} respondents: {n_resp} " f"criteria: {sorted(int(k) for k in df['criterion_id'].unique())}") require(n_resp == EXPECTED_N_DESCRIPTIVE, f"Secondary model respondents = {n_resp}, expected {EXPECTED_N_DESCRIPTIVE}.") require(len(df) == EXPECTED_N_LONG, f"Secondary model rows = {len(df)}, expected {EXPECTED_N_LONG}.") log(f" {SECONDARY_FORMULA}") rank_info = assert_full_rank(SECONDARY_FORMULA, df, "secondary", log) res = fit_gee(SECONDARY_FORMULA, df, "secondary GEE", log) rows: list[dict[str, Any]] = [] for k in CRITERIA: col = [c for c in res.params.index if c.startswith("BMI:") and c.endswith(f"[{k}]")] require(len(col) == 1, f"Could not uniquely identify the BMI x criterion_id slope for C{k}: {col}") name = col[0] est = or_per_5(float(res.params[name]), float(res.bse[name])) rows.append({ "criterion_id": k, "label": CRITERION_LABELS[k], "criterion_set": criterion_set_of(k), "term": name, "p_value": float(res.pvalues[name]), **est, }) pvals = [r["p_value"] for r in rows] reject, qvals, _, _ = multipletests(pvals, alpha=FDR_Q, method="fdr_bh") for r, q, rej in zip(rows, qvals, reject): r["q_value_bh"] = float(q) r["significant_bh_q05"] = bool(rej) log("") log("Eleven secondary BMI slopes (all odds ratios PER 5 BMI UNITS):") log(f" {'crit':<5}{'set':<21}{'OR/5':>9}{'95% CI (per 5 BMI units)':>28}" f"{'p':>12}{'q (BH)':>12}{'sig q<.05':>11}") for r in rows: ci = f"{r['ci_lo_per_5']:.4f} to {r['ci_hi_per_5']:.4f}" log(f" C{r['criterion_id']:<4}{r['criterion_set']:<21}{r['or_per_5']:>9.4f}" f"{ci:>28}{r['p_value']:>12.5g}{r['q_value_bh']:>12.5g}" f"{str(r['significant_bh_q05']):>11}") return {"n_rows": len(df), "n_respondents": n_resp, "formula": SECONDARY_FORMULA, "rank_assertion": rank_info, "criteria": rows, "fdr_method": "fdr_bh", "fdr_q": FDR_Q} # --------------------------------------------------------------------------- # [D] SAP §2.4 — unconditioned descriptive models # --------------------------------------------------------------------------- def section_d_unconditioned(wide: pd.DataFrame, secondary: dict[str, Any], log: Logger) -> dict[str, Any]: log.rule("[D] SAP §2.4 — UNCONDITIONED MODELS (DESCRIPTIVE, NOT A TEST OF H1)") log("Eleven separate logistic regressions, `criterion_met ~ BMI + C(age_range) + sex`,") log("with NO rest-score conditioning. Retained for transparency and as the") log("didactic contrast in SAP §2.4: the gap between these and the conditioned") log("slopes in [C] is itself the reportable demonstration of the circularity") log("problem SAP §2 point 1 asserts.") log("") rows: list[dict[str, Any]] = [] for k in CRITERIA: frame = pd.DataFrame({ "y": wide[f"c{k}"].to_numpy(), "BMI": wide["BMI"].to_numpy(), "age_range": wide["age_range"].to_numpy(), "sex": wide["sex"].to_numpy(), }) with warnings.catch_warnings(): warnings.simplefilter("ignore") res = smf.logit("y ~ BMI + C(age_range) + sex", data=frame).fit(disp=0) check_converged(res, f"unconditioned Logit C{k}", log) est = or_per_5(float(res.params["BMI"]), float(res.bse["BMI"])) rows.append({ "criterion_id": k, "criterion_set": criterion_set_of(k), "n": int(res.nobs), "p_value": float(res.pvalues["BMI"]), **est, }) by_k = {r["criterion_id"]: r for r in secondary["criteria"]} log("") log("SIDE BY SIDE — unconditioned (descriptive) vs rest-score-conditioned " "(secondary [C]).") log("All odds ratios PER 5 BMI UNITS.") log(f" {'crit':<5}{'set':<21}{'uncond OR/5':>13}{'uncond p':>12}" f"{'cond OR/5':>12}{'cond p':>12}{'shift':>10}") for r in rows: c = by_k[r["criterion_id"]] shift = r["or_per_5"] - c["or_per_5"] r["conditioned_or_per_5"] = c["or_per_5"] r["or_shift_uncond_minus_cond"] = float(shift) log(f" C{r['criterion_id']:<4}{r['criterion_set']:<21}{r['or_per_5']:>13.4f}" f"{r['p_value']:>12.5g}{c['or_per_5']:>12.4f}{c['p_value']:>12.5g}" f"{shift:>+10.4f}") n_pos_uncond = sum(1 for r in rows if r["or_per_5"] > 1.0) n_sig_uncond = sum(1 for r in rows if r["p_value"] < ALPHA) log("") log(f" Unconditioned BMI odds ratios above 1: {n_pos_uncond} of 11") log(f" Unconditioned BMI slopes with p < 0.05: {n_sig_uncond} of 11") log(" SAP §2 point 1 predicted the unconditioned picture would be uniformly") log(" positive by construction. The counts above are the empirical check on") log(" that assertion, and are descriptive only.") return {"criteria": rows, "n_unconditioned_or_above_1": n_pos_uncond, "n_unconditioned_p_below_0_05": n_sig_uncond} # --------------------------------------------------------------------------- # [E] Linearity of the BMI logit # --------------------------------------------------------------------------- def rcs_basis(x: np.ndarray, knots: np.ndarray) -> np.ndarray: """Restricted cubic spline basis (Harrell), k knots -> k-2 basis columns beyond the linear term. Returns only the non-linear columns.""" k = len(knots) t = knots out = np.empty((len(x), k - 2)) denom = (t[-1] - t[0]) ** 2 for j in range(k - 2): term = ( np.maximum(x - t[j], 0) ** 3 - np.maximum(x - t[-2], 0) ** 3 * (t[-1] - t[j]) / (t[-1] - t[-2]) + np.maximum(x - t[-1], 0) ** 3 * (t[-2] - t[j]) / (t[-1] - t[-2]) ) out[:, j] = term / denom return out def section_e_linearity(long: pd.DataFrame, log: Logger) -> dict[str, Any]: log.rule("[E] LINEARITY OF THE BMI LOGIT — tested, not assumed") log("Locked decision (00-STATE.md §4): 'BMI continuous primary, linearity tested") log("not assumed.' The primary specification in [B] is NEVER switched on the basis") log("of this test. If BMI is materially non-linear, that is REPORTED as a") log("limitation and the pre-registered primary model stands.") log("") df = long[ long["criterion_set"].isin(CONFIRMATORY_SETS) & (long["primary_model_eligible"] == 1) ].copy() bmi_mean = float(df["BMI"].mean()) df["bmi_c"] = df["BMI"] - bmi_mean df["bmi_c2"] = df["bmi_c"] ** 2 knots = np.percentile(df["BMI"].to_numpy(), [5, 35, 65, 95]) basis = rcs_basis(df["BMI"].to_numpy(), knots) df["rcs1"] = basis[:, 0] df["rcs2"] = basis[:, 1] log(f"RCS knots (5/35/65/95th centiles of BMI in the primary sample): " f"{np.round(knots, 3).tolist()}") log("") out: dict[str, Any] = {"bmi_mean": bmi_mean, "rcs_knots": knots.tolist(), "tests": []} set_term = f"C(criterion_set, Treatment(reference='{SET_REFERENCE}'))" variants = [ ("quadratic, common across sets", PRIMARY_FORMULA + " + bmi_c2", ["bmi_c2"]), ("quadratic, set-specific", PRIMARY_FORMULA + f" + bmi_c2:{set_term}", ["bmi_c2:"]), ("restricted cubic spline (4 knots), common across sets", PRIMARY_FORMULA + " + rcs1 + rcs2", ["rcs1", "rcs2"]), ] for label, formula, prefixes in variants: log(f"Variant: {label}") assert_full_rank(formula, df, label, log) res = fit_gee(formula, df, label, log) names = list(res.params.index) added = [c for c in names if any(c == p or c.startswith(p) for p in prefixes)] require(bool(added), f"No added non-linear terms found for variant {label}.") contrast = np.zeros((len(added), len(names))) for i, a in enumerate(added): contrast[i, names.index(a)] = 1.0 wt = res.wald_test(contrast, scalar=True) stat = float(np.ravel(wt.statistic)[0]) pv = float(np.ravel(wt.pvalue)[0]) log(f" added terms: {added}") log(f" joint Wald test of non-linearity: chi2 = {stat:.4f}, " f"df = {len(added)}, p = {pv:.6g}") log(f" materially non-linear at alpha=0.05: {pv < ALPHA}") out["tests"].append({ "variant": label, "added_terms": added, "chi2": stat, "df": len(added), "p_value": pv, "non_linear_at_0_05": bool(pv < ALPHA), }) log("") any_nl = any(t["non_linear_at_0_05"] for t in out["tests"]) out["any_non_linearity_flagged"] = bool(any_nl) log(f"Any variant flagging non-linearity: {any_nl}") log("The primary model in [B] is UNCHANGED regardless of the above. This is a") log("reported diagnostic, not a model-selection step.") return out # --------------------------------------------------------------------------- # [F] SAP §3 — IRT / DIF convergent check # --------------------------------------------------------------------------- def sr_dif(y: np.ndarray, theta: np.ndarray, group: np.ndarray, label: str) -> dict[str, Any]: """Swaminathan-Rogers logistic-regression DIF. M0: y ~ theta (matched on ability only) M1: y ~ theta + group (uniform DIF = M1 vs M0) M2: y ~ theta + group*theta (non-uniform DIF = M2 vs M1) """ frame = pd.DataFrame({"y": y, "theta": theta, "grp": group}) with warnings.catch_warnings(): warnings.simplefilter("ignore") m0 = smf.logit("y ~ theta", data=frame).fit(disp=0) m1 = smf.logit("y ~ theta + grp", data=frame).fit(disp=0) m2 = smf.logit("y ~ theta * grp", data=frame).fit(disp=0) for m, nm in ((m0, "M0"), (m1, "M1"), (m2, "M2")): conv = m.mle_retvals.get("converged") if conv is not True: raise SapDeparture(f"DIF model {nm} did not converge for {label}.") lr_u = 2.0 * (m1.llf - m0.llf) p_u = float(sp_stats.chi2.sf(lr_u, 1)) lr_nu = 2.0 * (m2.llf - m1.llf) p_nu = float(sp_stats.chi2.sf(lr_nu, 1)) int_name = [c for c in m2.params.index if ":" in c][0] return { "uniform_beta": float(m1.params["grp"]), "uniform_se": float(m1.bse["grp"]), "uniform_lr_chi2": float(lr_u), "uniform_p": p_u, "nonuniform_beta": float(m2.params[int_name]), "nonuniform_se": float(m2.bse[int_name]), "nonuniform_lr_chi2": float(lr_nu), "nonuniform_p": p_nu, } def section_f_irt_dif(wide: pd.DataFrame, secondary: dict[str, Any], log: Logger) -> dict[str, Any]: import girth log.rule("[F] SAP §3 — IRT / DIF CONVERGENT-METHODS CHECK") log("Unidimensional 2PL on the ELEVEN BINARY CRITERIA (not the 35 items — those") log("are structurally unusable). girth 0.8.0 `twopl_mml`, marginal maximum") log("likelihood; abilities by EAP (`girth.ability_eap`).") log("") log("DIF method, recorded as SAP §3 requires: Swaminathan-Rogers logistic") log("regression conditioned on the 2PL ability estimate. Uniform DIF = the group") log("main effect; non-uniform DIF = the group x ability interaction; both by") log("likelihood-ratio test on 1 df. Grouping = BMI split at the analytic cohort's") log("median, plus a continuous-BMI sensitivity check. girth is available, so SAP") log("§3's Mantel-Haenszel fallback is NOT used and no §9 deviation arises.") log("") resp = wide[[f"c{k}" for k in CRITERIA]].to_numpy().astype(int) # respondents x items dataset = resp.T # girth wants items x participants log(f"2PL fitted on {dataset.shape[1]} respondents x {dataset.shape[0]} criteria.") est = girth.twopl_mml(dataset) disc = np.asarray(est["Discrimination"], dtype=float) diff = np.asarray(est["Difficulty"], dtype=float) require(np.isfinite(disc).all() and np.isfinite(diff).all(), "2PL returned non-finite item parameters — the IRT arm cannot be reported.") theta = np.asarray(girth.ability_eap(dataset, diff, disc), dtype=float) require(len(theta) == len(wide), "Ability vector length does not match the cohort.") log(f" 2PL converged to finite parameters for all 11 criteria: True") log(f" ability (EAP): mean = {theta.mean():.4f} sd = {theta.std(ddof=1):.4f} " f"range {theta.min():.3f} to {theta.max():.3f}") log(f" corr(ability, total_symptom_count) = " f"{np.corrcoef(theta, wide['total_symptom_count'])[0, 1]:.4f}") log("") log("2PL item parameters:") log(f" {'crit':<5}{'set':<21}{'difficulty':>12}{'discrimination':>16}") items: list[dict[str, Any]] = [] for i, k in enumerate(CRITERIA): log(f" C{k:<4}{criterion_set_of(k):<21}{diff[i]:>12.4f}{disc[i]:>16.4f}") items.append({"criterion_id": k, "criterion_set": criterion_set_of(k), "difficulty": float(diff[i]), "discrimination": float(disc[i])}) bmi = wide["BMI"].to_numpy(dtype=float) median_bmi = float(np.median(bmi)) group = (bmi >= median_bmi).astype(int) log("") log(f"DIF grouping: BMI >= median ({median_bmi:.4f}) = 1, below = 0. " f"n high = {int(group.sum())}, n low = {int((1 - group).sum())}") dif_rows: list[dict[str, Any]] = [] for i, k in enumerate(CRITERIA): y = resp[:, i] med = sr_dif(y, theta, group.astype(float), f"C{k} median-split") cont = sr_dif(y, theta, bmi, f"C{k} continuous BMI") # Continuous-BMI DIF effects are reported per 5 BMI units for consistency. cont_u = or_per_5(cont["uniform_beta"], cont["uniform_se"]) dif_rows.append({ "criterion_id": k, "criterion_set": criterion_set_of(k), "median_split": med, "continuous_bmi": {**cont, "uniform_or_per_5": cont_u["or_per_5"], "uniform_ci_lo_per_5": cont_u["ci_lo_per_5"], "uniform_ci_hi_per_5": cont_u["ci_hi_per_5"]}, }) for scope in ("median_split", "continuous_bmi"): for kind in ("uniform_p", "nonuniform_p"): ps = [r[scope][kind] for r in dif_rows] _, qs, _, _ = multipletests(ps, alpha=FDR_Q, method="fdr_bh") for r, q in zip(dif_rows, qs): r[scope][kind.replace("_p", "_q_bh")] = float(q) by_k = {r["criterion_id"]: r for r in secondary["criteria"]} log("") log("DIF TABLE, set alongside the [C] per-criterion GEE coefficients so") log("convergence or divergence between the two arms is directly visible.") log("(BH q-values on the DIF columns are supplied for transparency; SAP §3 does") log("not pre-specify multiplicity control for this arm, so raw p is the") log("pre-registered quantity and q is supplementary.)") log("") log(f" {'crit':<5}{'set':<21}{'unif p':>10}{'unif q':>10}{'nonunif p':>11}" f"{'nonunif q':>11}{'GEE OR/5':>10}{'GEE q':>10}{'agree':>8}") for r in dif_rows: m = r["median_split"] g = by_k[r["criterion_id"]] dif_flag = m["uniform_p"] < ALPHA or m["nonuniform_p"] < ALPHA agree = dif_flag == g["significant_bh_q05"] r["dif_flagged_median_split"] = bool(dif_flag) r["agrees_with_gee_secondary"] = bool(agree) log(f" C{r['criterion_id']:<4}{r['criterion_set']:<21}{m['uniform_p']:>10.4g}" f"{m['uniform_q_bh']:>10.4g}{m['nonuniform_p']:>11.4g}" f"{m['nonuniform_q_bh']:>11.4g}{g['or_per_5']:>10.4f}" f"{g['q_value_bh']:>10.4g}{str(agree):>8}") log("") log("Continuous-BMI sensitivity check (uniform DIF, per 5 BMI units):") log(f" {'crit':<5}{'set':<21}{'OR/5':>10}{'95% CI (per 5 BMI units)':>28}{'p':>12}") for r in dif_rows: c = r["continuous_bmi"] ci = f"{c['uniform_ci_lo_per_5']:.4f} to {c['uniform_ci_hi_per_5']:.4f}" log(f" C{r['criterion_id']:<4}{r['criterion_set']:<21}" f"{c['uniform_or_per_5']:>10.4f}{ci:>28}{c['uniform_p']:>12.4g}") n_agree = sum(1 for r in dif_rows if r["agrees_with_gee_secondary"]) log("") log(f"Criteria where the DIF flag and the [C] BH-significant flag agree: " f"{n_agree} of 11") log("SAP §10: a convergence failure between the GEE and IRT-DIF arms is itself") log("reported as a finding, not resolved by preferring whichever arm is tidier.") return {"model": "2PL (girth.twopl_mml), abilities by EAP", "dif_method": "Swaminathan-Rogers logistic regression on 2PL EAP ability", "median_bmi": median_bmi, "items": items, "dif": dif_rows, "n_arms_agreeing": n_agree} # --------------------------------------------------------------------------- # [G] SAP §4 — H2 discordant phenotypes # --------------------------------------------------------------------------- H2_FORMULA = ("endorsed_criterion_is_content_confounded ~ group " "+ total_symptom_count + C(age_range) + sex") H2_FORMULA_UNADJ = "endorsed_criterion_is_content_confounded ~ group" def h2_exchangeable_failure(df: pd.DataFrame, log: Logger) -> dict[str, Any]: """SAP §4's pre-registered Exchangeable H2 fit failed. Diagnose it in full. This diagnosis stays in the output permanently. SAP §9 Deviation 1 authorises the amended working correlation, but the failure that forced the amendment is part of the study's audit trail and is never reduced to a footnote: the `converged = False`, the empirical rho, the exchangeable lower bound and the maximum cluster size are all re-derived and re-printed on every run. """ log.sub("H2 STEP 1 — THE PRE-REGISTERED EXCHANGEABLE FIT FAILED (audit trail)") log("SAP §4 mandates a GEE with an EXCHANGEABLE working correlation, clustered") log("on respondent, 'mirroring the primary model's approach for methodological") log("consistency'. On the real analytic sample that fit does not converge: the") log("estimated exchangeable dependence parameter goes NEGATIVE and the") log("iteration then diverges to NaN at higher iteration limits.") log("") # Empirical within-cluster dependence, from a pooled logit's Pearson residuals. with warnings.catch_warnings(): warnings.simplefilter("ignore") pooled = smf.logit(H2_FORMULA.replace( "endorsed_criterion_is_content_confounded", "y"), data=df.assign(y=df["endorsed_criterion_is_content_confounded"]) ).fit(disp=0) r = pd.Series(pooled.resid_pearson, index=df.index) a: list[float] = [] b: list[float] = [] for _, gg in r.groupby(df["respondent_id"]): v = gg.to_numpy() for i in range(len(v)): for j in range(i + 1, len(v)): a.append(v[i]) b.append(v[j]) emp_rho = float(np.corrcoef(a, b)[0, 1]) sizes = df.groupby("respondent_id").size() log("WHY. The H2 outcome is COMPOSITIONAL within a respondent. Each cluster is") log("that respondent's endorsed criteria, partitioned into content-confounded") log("(at most 5 available) and pharmacological-core (at most 3 available). One") log("row being content-confounded makes another row in the same cluster less") log("likely to be. The within-cluster dependence is therefore NEGATIVE BY") log("CONSTRUCTION, and bounded below by -1/(m-1) for cluster size m — a region") log("in which the exchangeable working correlation is not a valid correlation") log("matrix and the GEE iteration is not guaranteed to converge.") log("") log(f" empirical within-cluster pairwise correlation of Pearson residuals = " f"{emp_rho:+.4f} ({len(a)} pairs)") log(f" cluster sizes: min {int(sizes.min())}, median {int(sizes.median())}, " f"max {int(sizes.max())}, mean {sizes.mean():.2f}") log(f" exchangeable lower bound at the largest cluster (-1/(m-1)) = " f"{-1.0 / (int(sizes.max()) - 1):+.4f}") log("") log("This is NOT the same dependence structure as the primary model's. There,") log("the outcome is criterion endorsement and within-respondent dependence is") log("strongly POSITIVE (all eleven criteria load on one severity trait), which") log("is exactly what exchangeable is for. SAP §4 carried the structure across") log("by analogy without the analogy holding.") log("") log("############################################################################") log("# SAP §4 SPECIFICATION DEFECT — AMENDED UNDER SAP §9 DEVIATION 1") log("# The originally pre-registered H2 fit above CANNOT BE FITTED as written.") log("# The amendment authorised on 2026-07-26 by the project owner changes the") log("# working correlation to Independence() AND NOTHING ELSE. This failure") log("# block is retained permanently as the audit trail for that amendment.") log("############################################################################") return { "prespecified_exchangeable_converged": False, "failure": "SAP §4's originally pre-registered Exchangeable GEE did not converge", "diagnosis": ( "The H2 outcome is compositional within respondent, so within-cluster " "dependence is negative by construction; the exchangeable working " "correlation's dependence parameter goes negative and the iteration " "diverges." ), "empirical_within_cluster_rho": emp_rho, "n_residual_pairs": len(a), "cluster_size_min": int(sizes.min()), "cluster_size_max": int(sizes.max()), "cluster_size_mean": float(sizes.mean()), "exchangeable_lower_bound_at_max_cluster": -1.0 / (int(sizes.max()) - 1), "resolved_by": "SAP §9 Deviation 1 (2026-07-26)", } def section_g_h2(wide: pd.DataFrame, long: pd.DataFrame, log: Logger, permuted: bool) -> dict[str, Any]: log.rule("[G] SAP §4 — H2 DISCORDANT PHENOTYPES (confirmatory secondary)") normal_high = (wide["BMI"] >= 18.5) & (wide["BMI"] < 25.0) & \ (wide["total_symptom_count"] >= 6) obese_high = (wide["BMI"] >= 30.0) & (wide["total_symptom_count"] >= 6) obese_low = (wide["BMI"] >= 30.0) & (wide["total_symptom_count"] <= 1) n_nh, n_oh, n_ol = int(normal_high.sum()), int(obese_high.sum()), int(obese_low.sum()) log(f"Cell sizes in the frozen analytic cohort:") log(f" normal BMI (18.5-<25), >=6 symptoms : n = {n_nh} (SAP §4 quoted 182)") log(f" obese (>=30), >=6 symptoms : n = {n_oh} (SAP §4 quoted 407)") log(f" obese (>=30), <=1 symptoms : n = {n_ol} (SAP §4 quoted 99)") log("The SAP's figures were pre-exclusion counts; the cells SHRANK after the §7") log("exclusions were applied. This shrinkage is reported, not absorbed.") if permuted: log(" (PERMUTED RUN — cells are re-cut on shuffled BMI, so the verified") log(" 155/394/98 assertion is deliberately not enforced here.)") else: require( (n_nh, n_oh, n_ol) == (EXPECTED_H2_NORMAL_HIGH, EXPECTED_H2_OBESE_HIGH, EXPECTED_H2_OBESE_LOW), f"H2 cells are {(n_nh, n_oh, n_ol)}, expected " f"{(EXPECTED_H2_NORMAL_HIGH, EXPECTED_H2_OBESE_HIGH, EXPECTED_H2_OBESE_LOW)} " "(00-STATE.md step 3a). Aborting rather than reporting a cell the " "coordinator has not verified.", ) groups = pd.Series(pd.NA, index=wide.index, dtype="object") groups[obese_high] = "obese_high_symptom" groups[normal_high] = "normal_bmi_high_symptom" wide = wide.assign(_h2_group=groups) mean_nh = float(wide.loc[normal_high, "total_symptom_count"].mean()) mean_oh = float(wide.loc[obese_high, "total_symptom_count"].mean()) log("") log("Mean total_symptom_count by group (SAP §4 requires these be reported") log("alongside, so the reader sees how little severity separation the mandatory") log("adjustment absorbs):") log(f" normal-BMI high-symptom : {mean_nh:.4f}") log(f" obese high-symptom : {mean_oh:.4f}") log(f" difference : {mean_oh - mean_nh:+.4f} symptoms") ids = wide.loc[normal_high | obese_high, ["respondent_id", "_h2_group"]] df = long.merge(ids, on="respondent_id", how="inner") df = df[(df["criterion_met"] == 1) & df["criterion_set"].isin(CONFIRMATORY_SETS)].copy() df["endorsed_criterion_is_content_confounded"] = ( df["criterion_set"] == "content_confounded").astype(int) # Reference = obese_high_symptom, so the group coefficient reads # "normal-BMI high-symptom minus obese high-symptom", the direction H2 states. df["group"] = pd.Categorical( df["_h2_group"], categories=["obese_high_symptom", "normal_bmi_high_symptom"]) n_resp = int(df["respondent_id"].nunique()) log("") log(f"Endorsed-criterion rows entering the H2 fit: {len(df)} respondents: {n_resp}") log("Restricted to ENDORSED criteria in the content-confounded and") log("pharmacological-core sets only (SAP §4). C1-C3 endorsements are not rows here.") log(f" respondents per group: " f"{df.groupby('_h2_group')['respondent_id'].nunique().to_dict()}") require(n_resp == n_nh + n_oh, f"H2 respondents {n_resp} != {n_nh} + {n_oh}.") raw_share = df.groupby("_h2_group")[ "endorsed_criterion_is_content_confounded"].mean().to_dict() log("") log("Raw content-confounded share of endorsed criteria (descriptive):") for g, v in raw_share.items(): log(f" {g:<26} {v:.4f}") log("") log(f"Adjusted (CONFIRMATORY) model — the total_symptom_count adjustment is") log("MANDATORY per SAP §4, not optional:") log(f" {H2_FORMULA}") log(" GEE, logit link, Binomial family, Exchangeable, clustered on respondent_id.") rank_adj = assert_full_rank(H2_FORMULA, df, "H2 adjusted", log) z = sp_stats.norm.ppf(1 - ALPHA / 2) def group_effect(res: Any) -> dict[str, Any]: gn = [c for c in res.params.index if c.startswith("group")] require(len(gn) == 1, f"Could not identify the H2 group term: {gn}") gn = gn[0] b, se = float(res.params[gn]), float(res.bse[gn]) return {"term": gn, "beta": b, "se": se, "odds_ratio": float(np.exp(b)), "ci_lo": float(np.exp(b - z * se)), "ci_hi": float(np.exp(b + z * se)), "z": float(res.tvalues[gn]), "p_value": float(res.pvalues[gn])} # --- Step 1: attempt the ORIGINALLY pre-registered Exchangeable fit ----- # # This attempt is retained even though SAP §9 Deviation 1 already authorises # the amendment. Re-deriving the failure on every run is what keeps the # amendment auditable rather than asserted. res_x = fit_gee(H2_FORMULA, df, "H2 adjusted GEE (Exchangeable — ORIGINAL pre-registration)", log, allow_failure=True) out_h2: dict[str, Any] = {} if res_x is None: out_h2.update(h2_exchangeable_failure(df, log)) else: out_h2["prespecified_exchangeable_converged"] = True out_h2["exchangeable_fit"] = group_effect(res_x) # --- Step 2: the ADOPTED model, per SAP §9 Deviation 1 ------------------ log.sub("H2 STEP 2 — THE ADOPTED CONFIRMATORY MODEL (SAP §9 DEVIATION 1)") log("AMENDMENT NOTICE — this result is NOT the fit originally pre-registered in") log("SAP §4. It is the amended fit authorised by the project owner on 2026-07-26") log("and recorded in full at SAP §9, 'Deviation 1'. Any reader, table or article") log("carrying the numbers below must carry this notice with them.") log("") log(" Amendment (SAP §9 Deviation 1): working correlation changed from") log(" Exchangeable() to Independence(). NOTHING ELSE CHANGED — same GEE, same") log(" logit link and Binomial family, same formula including the mandatory") log(" total_symptom_count adjustment, same cluster-robust sandwich standard") log(" errors on respondent_id, same population-averaged estimand.") log(" Reason: forced by the non-convergence documented in step 1, not chosen.") log(" GEE point estimates are consistent under ANY working correlation; the") log(" working correlation affects efficiency only, and inference rests on the") log(" sandwich estimator, which stays valid under negative within-cluster") log(" dependence.") log(" SAP §9 also discloses that the Independence fit was computed, and held") log(" under quarantine, BEFORE the amendment was authorised. That ordering is") log(" disclosed there rather than concealed.") log("") res = fit_gee(H2_FORMULA, df, "H2 adjusted GEE (Independence — ADOPTED)", log, cov_struct=sm.cov_struct.Independence()) log("") log(str(res.summary())) adj = group_effect(res) adj["working_correlation"] = "Independence" adj["amendment"] = "SAP §9 Deviation 1 (2026-07-26)" log("") log("Unadjusted comparison (DESCRIPTIVE, labelled as such per SAP §4 — this is") log("NOT the confirmatory result and must never be substituted for it):") log(f" {H2_FORMULA_UNADJ}") assert_full_rank(H2_FORMULA_UNADJ, df, "H2 unadjusted", log) res_u = fit_gee(H2_FORMULA_UNADJ, df, "H2 unadjusted GEE (Independence — descriptive)", log, cov_struct=sm.cov_struct.Independence()) unadj = group_effect(res_u) unadj["working_correlation"] = "Independence" unadj["status"] = "DESCRIPTIVE ONLY — not the confirmatory result" out_h2["adjusted_confirmatory"] = adj out_h2["unadjusted_descriptive"] = unadj log.sub("H2 RESULT — SAP §4 confirmatory secondary, as amended by SAP §9 Deviation 1") log("Outcome: is an endorsed criterion content-confounded rather than") log("pharmacological-core? Group coefficient = normal-BMI-high-symptom relative") log("to obese-high-symptom (reference). H2 predicted OR < 1.") log("") log(f" ADJUSTED (CONFIRMATORY, SAP §9 Deviation 1)") log(f" OR = {adj['odds_ratio']:.4f} 95% CI {adj['ci_lo']:.4f} to " f"{adj['ci_hi']:.4f} z = {adj['z']:+.4f} p = {adj['p_value']:.6g}") log(f" UNADJUSTED (DESCRIPTIVE, NOT the confirmatory result)") log(f" OR = {unadj['odds_ratio']:.4f} 95% CI {unadj['ci_lo']:.4f} to " f"{unadj['ci_hi']:.4f} z = {unadj['z']:+.4f} p = {unadj['p_value']:.6g}") log("") supported = bool(adj["ci_hi"] < 1.0) log(" **** H2 IS NOT SUPPORTED. ****" if not supported else " **** H2 IS SUPPORTED. ****") log(f" The pre-specified direction IS the direction observed " f"(OR = {adj['odds_ratio']:.4f} < 1: the normal-BMI high-symptom group shows") log(" a lower content-confounded share), but the 95% confidence interval") log(f" ({adj['ci_lo']:.4f} to {adj['ci_hi']:.4f}) INCLUDES 1 and p = " f"{adj['p_value']:.4f} > {ALPHA}.") log(" Per SAP §10 this null is reported with the same prominence a positive") log(" result would have received.") log("") log(" The unadjusted comparison IS significant (OR " f"{unadj['odds_ratio']:.4f}, p = {unadj['p_value']:.5g}). The gap between") log(" adjusted and unadjusted is the finding SAP §4's adjustment argument") log(" predicted, and is reported as such — it is NOT the confirmatory result.") out_h2["h2_supported"] = supported out_h2["h2_conclusion"] = ( "H2 is NOT supported. Direction as predicted (OR < 1) but the 95% CI " "includes 1." ) log("") # --- exact power recomputation against the REAL cells ------------------- log.sub("H2 power, recomputed EXACTLY against the real cells (SAP §4 commits to this)") log(f"SAP §4's a priori statement was written against n = 182 / 407. The real") log(f"cells after the §7 exclusions are n = {n_nh} / {n_oh}. Both are computed") log("below so the shrinkage is visible and quantified, not just noted.") power_calc = NormalIndPower() base = float(df["endorsed_criterion_is_content_confounded"].mean()) log("") log(f"Pooled content-confounded share used as the null baseline: {base:.4f}") log("Two independent proportions, alpha = 0.05 two-sided, Cohen's h effect size.") log("NOTE: this is a respondent-level normal approximation. It ignores the") log("within-respondent clustering the fitted GEE accounts for, so it is an") log("approximation to the power of the model actually fitted, exactly as SAP §4") log("frames it.") log("") log(f" {'cells':<14}{'ratio':>8}{'min det. diff at 80%':>24}" f"{'power @12pp':>13}{'power @15pp':>13}") power_rows: list[dict[str, Any]] = [] for label, n1, n2 in ((f"real {n_nh}/{n_oh}", n_nh, n_oh), ("SAP 182/407", 182, 407)): ratio = n2 / n1 h = float(power_calc.solve_power(effect_size=None, nobs1=n1, alpha=ALPHA, power=POWER_TARGET, ratio=ratio, alternative="two-sided")) # Convert the detectable Cohen's h back to an absolute difference at `base`. phi1 = 2 * np.arcsin(np.sqrt(base)) p2 = float(np.sin((phi1 + h) / 2.0) ** 2) mdd = p2 - base pw12 = float(power_calc.solve_power( effect_size=proportion_effectsize(min(base + 0.12, 0.999), base), nobs1=n1, alpha=ALPHA, ratio=ratio, alternative="two-sided")) pw15 = float(power_calc.solve_power( effect_size=proportion_effectsize(min(base + 0.15, 0.999), base), nobs1=n1, alpha=ALPHA, ratio=ratio, alternative="two-sided")) log(f" {label:<14}{ratio:>8.3f}{mdd * 100:>21.2f} pp{pw12:>13.4f}{pw15:>13.4f}") power_rows.append({"cells": label, "n1": n1, "n2": n2, "detectable_cohens_h_at_80pc": h, "min_detectable_abs_diff_pp": mdd * 100, "power_at_12pp": pw12, "power_at_15pp": pw15}) log("") log("Achieved power against the observed adjusted effect is not a pre-registered") log("quantity (post-hoc observed power is uninformative), so it is not reported.") log("The pre-registered quantity is the detectable difference at 80% power above.") real_row = power_rows[0] sap_row = power_rows[1] under = real_row["power_at_12pp"] < POWER_TARGET log("") log("**** UNDER-POWERING, ATTACHED TO THE H2 NULL ****") log("This must travel with the null result above; the null must NOT be read as a") log("confident absence of effect.") log(f" The H2 cells SHRANK after the SAP §7 exclusions were applied: the") log(f" pre-registration was designed against 182/407, the realised cells are") log(f" {n_nh}/{n_oh}.") log(f" Minimum detectable difference at 80% power rose from " f"{sap_row['min_detectable_abs_diff_pp']:.2f} pp (pre-registered cells) to " f"{real_row['min_detectable_abs_diff_pp']:.2f} pp (realised cells).") log(f" Power at a 12 pp difference fell from {sap_row['power_at_12pp']:.4f} to " f"{real_row['power_at_12pp']:.4f}, against the {POWER_TARGET:.2f} the design") log(f" assumed.") log(f" H2 was therefore {'UNDER-POWERED' if under else 'adequately powered'} " "relative to its pre-registered design. A non-significant result at this") log(" power is consistent with a real effect of the hypothesised size going") log(" undetected, and is reported as an inconclusive null rather than as") log(" evidence of no difference.") out_h2["under_powered_vs_design"] = bool(under) out_h2["power_caveat_on_null"] = ( f"Cells shrank from the pre-registered 182/407 to {n_nh}/{n_oh}. Minimum " f"detectable difference at 80% power rose from " f"{sap_row['min_detectable_abs_diff_pp']:.2f} pp to " f"{real_row['min_detectable_abs_diff_pp']:.2f} pp, and power at a 12 pp " f"difference fell from {sap_row['power_at_12pp']:.4f} to " f"{real_row['power_at_12pp']:.4f} against the {POWER_TARGET:.2f} assumed. " "The H2 null must not be read as a confident absence of effect." ) # --- obese, <=1 symptom: descriptive only ------------------------------- log.sub(f"Obese with <=1 symptom (n = {n_ol}) — DESCRIPTIVE ONLY, no test (SAP §4)") log("SAP §4 pre-registers no confirmatory hypothesis for this cell.") sub = wide.loc[obese_low] log(f" {'crit':<5}{'set':<21}{'n met':>8}{'rate %':>10}") ol_rows = [] for k in CRITERIA: nm = int(sub[f"c{k}"].sum()) log(f" C{k:<4}{criterion_set_of(k):<21}{nm:>8}{100 * nm / len(sub):>10.1f}") ol_rows.append({"criterion_id": k, "criterion_set": criterion_set_of(k), "n_met": nm, "rate_pct": round(100 * nm / len(sub), 2)}) log(f" mean total_symptom_count = {sub['total_symptom_count'].mean():.4f}") return { "cells": {"normal_bmi_high_symptom": n_nh, "obese_high_symptom": n_oh, "obese_low_symptom": n_ol}, "sap_quoted_cells": {"normal_bmi_high_symptom": 182, "obese_high_symptom": 407, "obese_low_symptom": 99}, "mean_total_symptom_count": {"normal_bmi_high_symptom": mean_nh, "obese_high_symptom": mean_oh}, "raw_content_confounded_share": {k: float(v) for k, v in raw_share.items()}, "n_endorsed_rows": len(df), "n_respondents": n_resp, "rank_assertion": rank_adj, **out_h2, "power": power_rows, "obese_low_symptom_profile": ol_rows, } # --------------------------------------------------------------------------- # [H] H3 (age, exploratory) and H4 (sex / C9 replication) # --------------------------------------------------------------------------- def section_i_heterogeneity(primary: dict[str, Any], secondary: dict[str, Any], log: Logger) -> dict[str, Any]: """Per-criterion heterogeneity within the H1 criterion sets. Structural description only. This section states WHAT the per-criterion slopes do relative to their set's fitted direction; it does not interpret what that means for H1. That is the write-up's job, and doing it here would be this script exceeding its remit. 'Set direction' is taken from the primary model's own fitted set-level slope in [B], not from the hypothesis — so the flag is data-defined, not a scored judgement against what H1 wanted to see. """ log.rule("[I] PER-CRITERION HETEROGENEITY WITHIN THE H1 CRITERION SETS") log("The confirmatory contrast in [B] is a difference between two SET-LEVEL BMI") log("slopes. This section sets each criterion's own conditioned slope from [C]") log("against the fitted direction of the set it was assigned to, so that") log("within-set heterogeneity is visible rather than averaged away.") log("") log("Reported neutrally as a structural feature of the result. No interpretation") log("of its consequences for H1 is offered here.") log("") set_slopes = { "content_confounded": primary["slope_content_confounded"]["or_per_5"], "pharmacological_core": primary["slope_pharmacological_core"]["or_per_5"], } log("Set-level BMI slopes from the primary model [B] (OR per 5 BMI units):") for s, v in set_slopes.items(): log(f" {s:<22} {v:.4f} (direction: " f"{'above 1' if v > 1 else 'below 1'})") log(f" confirmatory contrast (ratio of the two) = " f"{primary['confirmatory_contrast']['or_per_5']:.4f}, " f"p = {primary['confirmatory_contrast']['p_value']:.6g}") log("") rows: list[dict[str, Any]] = [] for r in secondary["criteria"]: cset = r["criterion_set"] set_or = set_slopes.get(cset) if set_or is None: # C1-C3, 'neither' — no set slope exists in [B] counter = None else: counter = (r["or_per_5"] > 1.0) != (set_or > 1.0) rows.append({ "criterion_id": r["criterion_id"], "label": r["label"], "criterion_set": cset, "set_level_or_per_5": set_or, "criterion_or_per_5": r["or_per_5"], "ci_lo_per_5": r["ci_lo_per_5"], "ci_hi_per_5": r["ci_hi_per_5"], "q_value_bh": r["q_value_bh"], "significant_bh_q05": r["significant_bh_q05"], "runs_counter_to_set_direction": counter, "counter_and_significant": bool(counter and r["significant_bh_q05"]), }) log("All odds ratios PER 5 BMI UNITS. 'set OR/5' is that criterion's SET-level") log("slope from [B]; 'counter' means the criterion's own slope sits on the") log("opposite side of 1 from its set's slope. C1-C3 are in the 'neither' set,") log("which the confirmatory model does not fit, so no set slope exists for them.") log("") log(f" {'crit':<5}{'set':<21}{'set OR/5':>10}{'crit OR/5':>11}" f"{'95% CI (per 5 BMI units)':>28}{'q (BH)':>11}{'counter':>9}{'both':>7}") for r in rows: so = " n/a" if r["set_level_or_per_5"] is None \ else f"{r['set_level_or_per_5']:>8.4f}" cf = "n/a" if r["runs_counter_to_set_direction"] is None \ else str(r["runs_counter_to_set_direction"]) ci = f"{r['ci_lo_per_5']:.4f} to {r['ci_hi_per_5']:.4f}" log(f" C{r['criterion_id']:<4}{r['criterion_set']:<21}{so:>10}" f"{r['criterion_or_per_5']:>11.4f}{ci:>28}{r['q_value_bh']:>11.4g}" f"{cf:>9}{str(r['counter_and_significant']):>7}") counter_sig = [r for r in rows if r["counter_and_significant"]] counter_any = [r for r in rows if r["runs_counter_to_set_direction"]] log("") log(f" Criteria running counter to their set's fitted direction: " f"{[('C' + str(r['criterion_id'])) for r in counter_any] or 'none'}") log(f" Of those, BH-significant at q<0.05: " f"{[('C' + str(r['criterion_id'])) for r in counter_sig] or 'none'}") log("") log(" Within-set spread (OR per 5 BMI units), confirmatory sets only:") summary_sets: dict[str, Any] = {} for s in CONFIRMATORY_SETS: vals = [r["criterion_or_per_5"] for r in rows if r["criterion_set"] == s] ids = [r["criterion_id"] for r in rows if r["criterion_set"] == s] lo_i = ids[int(np.argmin(vals))] hi_i = ids[int(np.argmax(vals))] n_above = sum(1 for v in vals if v > 1.0) log(f" {s:<22} n={len(vals)} min {min(vals):.4f} (C{lo_i}) " f"max {max(vals):.4f} (C{hi_i}) above 1: {n_above}/{len(vals)}") summary_sets[s] = { "n": len(vals), "min_or_per_5": float(min(vals)), "min_criterion": lo_i, "max_or_per_5": float(max(vals)), "max_criterion": hi_i, "n_above_1": n_above, "set_level_or_per_5": set_slopes[s], } log("") log(" Neither set is homogeneous in the direction of its own set-level slope.") log(" This is stated as a structural feature of the result, for the write-up to") log(" confront. This script offers no reading of it.") return { "note": ("Structural description only. Set direction is taken from the " "primary model's fitted set-level slopes in [B], not from the " "hypothesis. No interpretation is offered here."), "set_level_slopes_or_per_5": set_slopes, "confirmatory_contrast_or_per_5": primary["confirmatory_contrast"]["or_per_5"], "criteria": rows, "criteria_counter_to_set_direction": [r["criterion_id"] for r in counter_any], "criteria_counter_and_bh_significant": [r["criterion_id"] for r in counter_sig], "within_set_spread": summary_sets, } H3_FORMULA = ("criterion_met ~ C(criterion_id) * C(age_range) + rest_score + sex") def section_h_h3_h4(wide: pd.DataFrame, long: pd.DataFrame, log: Logger) -> dict[str, Any]: out: dict[str, Any] = {} log.rule("[H1 of 2] H3 — AGE x CRITERION PROFILE (EXPLORATORY, NO DIRECTION)") log("SAP §1: no directional hypothesis is pre-registered for H3 — the literature") log("audit found no prior age-by-criterion finding to anchor one. Reported as") log("exploratory. No confirmatory claim attaches to anything below.") log("") log("Descriptive endorsement rate (%) by age band and criterion:") bands = [b for b in ["18-24", "25-34", "35-44", "45-54", "55-64", "65-74", "75+"] if b in set(wide["age_range"].unique())] header = f" {'band':<8}{'n':>6}" + "".join(f"{'C' + str(k):>7}" for k in CRITERIA) log(header) prof: list[dict[str, Any]] = [] for b in bands: sub = wide[wide["age_range"] == b] rates = [100 * float(sub[f"c{k}"].mean()) for k in CRITERIA] log(f" {b:<8}{len(sub):>6}" + "".join(f"{r:>7.1f}" for r in rates)) prof.append({"age_range": b, "n": len(sub), "rates_pct": {f"c{k}": round(r, 2) for k, r in zip(CRITERIA, rates)}}) out["age_criterion_profile"] = prof log("") log("Joint test of the criterion x age interaction block (exploratory):") log(f" {H3_FORMULA}") df = long.copy() rank_h3 = assert_full_rank(H3_FORMULA, df, "H3", log) res = fit_gee(H3_FORMULA, df, "H3 GEE", log) names = list(res.params.index) inter = [c for c in names if "criterion_id" in c and "age_range" in c] contrast = np.zeros((len(inter), len(names))) for i, a in enumerate(inter): contrast[i, names.index(a)] = 1.0 wt = res.wald_test(contrast, scalar=True) stat, pv = float(np.ravel(wt.statistic)[0]), float(np.ravel(wt.pvalue)[0]) log(f" interaction terms: {len(inter)}") log(f" joint Wald chi2 = {stat:.4f}, df = {len(inter)}, p = {pv:.6g}") log(f" criterion profile varies by age band at alpha=0.05: {pv < ALPHA}") out["h3_joint_test"] = {"formula": H3_FORMULA, "n_interaction_terms": len(inter), "chi2": stat, "df": len(inter), "p_value": pv, "significant_at_0_05": bool(pv < ALPHA), "rank_assertion": rank_h3, "status": "exploratory, no directional hypothesis"} # ---------------------------------------------------------------- H4 --- log.rule("[H2 of 2] H4 — SEX DIF ON C9 (replication attempt, Saffari et al. 2022)") log("SAP §1 H4: C9 (role obligation failure) is pre-registered to show a") log("sex-differential pattern, replicating Saffari 2022's gender-DIF finding on") log("this specific criterion. A sex effect on a DIFFERENT criterion is a FAILURE") log("TO REPLICATE, not a new discovery to chase.") log("") log("Method: the same Swaminathan-Rogers DIF machinery as [F], with sex as the") log("grouping variable, so the test is trait-conditioned exactly as Saffari's") log("Rasch-DIF analysis was. Restricted to female/male (the two levels Saffari") log("compared); sex 'other'/'prefer_not_to_say' cannot support a DIF contrast.") log("") import girth sub = wide[wide["sex"].isin(["female", "male"])].copy() log(f"n for the sex-DIF analysis: {len(sub)} " f"({(sub['sex'] == 'female').sum()} female, {(sub['sex'] == 'male').sum()} male)") resp = sub[[f"c{k}" for k in CRITERIA]].to_numpy().astype(int) est = girth.twopl_mml(resp.T) disc = np.asarray(est["Discrimination"], dtype=float) diff = np.asarray(est["Difficulty"], dtype=float) theta = np.asarray(girth.ability_eap(resp.T, diff, disc), dtype=float) grp = (sub["sex"] == "female").astype(float).to_numpy() # 1 = female log("2PL re-fitted on this subsample; group coded 1 = female, 0 = male.") rows: list[dict[str, Any]] = [] for i, k in enumerate(CRITERIA): r = sr_dif(resp[:, i], theta, grp, f"C{k} sex-DIF") r["criterion_id"] = k r["criterion_set"] = criterion_set_of(k) rows.append(r) ps_u = [r["uniform_p"] for r in rows] _, qs_u, _, _ = multipletests(ps_u, alpha=FDR_Q, method="fdr_bh") ps_n = [r["nonuniform_p"] for r in rows] _, qs_n, _, _ = multipletests(ps_n, alpha=FDR_Q, method="fdr_bh") for r, qu, qn in zip(rows, qs_u, qs_n): r["uniform_q_bh"] = float(qu) r["nonuniform_q_bh"] = float(qn) r["flagged"] = bool(r["uniform_p"] < ALPHA or r["nonuniform_p"] < ALPHA) log("") log("Sex DIF across all eleven criteria (uniform = female vs male main effect on") log("the log-odds of endorsement at equal ability; positive = female higher):") log(f" {'crit':<5}{'set':<21}{'unif beta':>11}{'unif p':>11}{'unif q':>11}" f"{'nonunif p':>11}{'flag':>7}") for r in rows: log(f" C{r['criterion_id']:<4}{r['criterion_set']:<21}" f"{r['uniform_beta']:>+11.4f}{r['uniform_p']:>11.4g}" f"{r['uniform_q_bh']:>11.4g}{r['nonuniform_p']:>11.4g}" f"{str(r['flagged']):>7}") c9 = next(r for r in rows if r["criterion_id"] == 9) others = [r for r in rows if r["criterion_id"] != 9 and r["flagged"]] replicates = bool(c9["flagged"]) log("") log("REPLICATION VERDICT — stated plainly as SAP §10 requires:") if replicates: log(f" C9 SHOWS sex DIF in this sample (uniform p = {c9['uniform_p']:.6g}, " f"non-uniform p = {c9['nonuniform_p']:.6g}).") log(" Saffari et al. 2022's C9 gender-DIF finding REPLICATES here.") else: log(f" C9 does NOT show sex DIF in this sample (uniform p = " f"{c9['uniform_p']:.6g}, non-uniform p = {c9['nonuniform_p']:.6g}).") log(" Saffari et al. 2022's C9 gender-DIF finding DOES NOT REPLICATE here.") if others: log(f" Sex DIF appears instead on: " f"{', '.join('C' + str(r['criterion_id']) for r in others)}.") log(" Per SAP §1, a sex effect on a criterion OTHER than C9 is a FAILURE TO") log(" REPLICATE, not a new finding. It is recorded, not chased.") else: log(" No other criterion shows a sex DIF flag either.") out["h4_sex_dif"] = {"n": len(sub), "criteria": rows, "c9_replicates_saffari_2022": replicates, "other_flagged_criteria": [r["criterion_id"] for r in others], "status": "confirmatory replication attempt, C9-specific"} return out # --------------------------------------------------------------------------- # [J] POST-HOC SENSITIVITY ANALYSIS — SAP §9 DEVIATION 2 # # *** NOT PRE-REGISTERED. NOT CONFIRMATORY. NOT A TEST OF H1. *** # # Added on 2026-07-26, AFTER the pre-registered analysis had been run and # reported, in response to a measurement defect found in the scoring code of the # instrument that generated this dataset: criterion 7 (withdrawal) was scored # from item 1 (a criterion-1 item) instead of item 11 (the affective withdrawal # item, which is scored against a threshold and then discarded). C7 as measured # is therefore a contaminated indicator — missing its affective withdrawal item # and carrying criterion-1 content. Item 11's responses were never persisted, so # the frozen cohort cannot be rescored; the only available response is to # quantify how much the pre-registered conclusion depends on C7. # # The full record, including the fact that the direction of the expected effect # was declared BEFORE this section was run, is at `02-statistical-analysis- # plan.md` §9 'Deviation 2'. # # Two variants, both specified in the SAP before either was run, both reported # whatever they return: # S1 set-membership only (PRIMARY SENSITIVITY). C7's rows are dropped from # the stacked model, so the pharmacological core becomes {C6, C11}. # rest_score and total_symptom_count are left EXACTLY as in [B], still # computed over all eleven criteria. One thing changes, so S1 is directly # comparable to the pre-registered contrast. # S2 full removal (SECONDARY). C7 is dropped entirely and # total_symptom_count and every rest_score are recomputed over the # remaining ten criteria. # # S1 is the primary sensitivity because it is the more conservative of the two, # a designation made in the SAP before either was fitted. Nothing in this # section is promoted to primary. The pre-registered contrast in [B] remains the # headline figure and is NOT re-run, re-fitted or altered here. # --------------------------------------------------------------------------- DEVIATION_2_REF = ("SAP §9 'Deviation 2' (2026-07-26) — criterion 7 measurement " "defect; POST-HOC, NOT PRE-REGISTERED") #: The pharmacological core as it stands once the contaminated indicator is #: removed. Stated as a constant so the sensitivity can never silently drift. PHARMACOLOGICAL_CORE_EX_C7 = (6, 11) C7 = 7 def fit_h1_contrast(df: pd.DataFrame, label: str, log: Logger) -> dict[str, Any]: """Fit PRIMARY_FORMULA on `df` and form the H1 contrast exactly as [B] does. Used ONLY by [J]. [B] is untouched and computes its own result inline; this function deliberately duplicates [B]'s logic rather than refactoring it, so that no edit made for the sensitivity analysis can reach the primary model. Every assertion [B] makes is made here too: design-matrix rank, absence of a criterion_set main effect, exactly two BMI x set slopes, GEE convergence, and agreement with the algebraically identical reparameterisation to 1e-6. """ rank_info = assert_full_rank(PRIMARY_FORMULA, df, label, log) cs_main = [c for c in rank_info["columns"] if "criterion_set" in c and "BMI" not in c] log(f" criterion_set MAIN-EFFECT columns present [{label}]: " f"{cs_main if cs_main else 'none'} <- must be none") require(not cs_main, f"A criterion_set main effect entered the design matrix for {label}. " "The sensitivity must use the primary model's specification exactly.") bmi_set_cols = [c for c in rank_info["columns"] if c.startswith("BMI:")] log(f" BMI x criterion_set columns present [{label}]: {bmi_set_cols}") require(len(bmi_set_cols) == 2, f"Expected 2 BMI x criterion_set columns for {label}, got {bmi_set_cols}.") res = fit_gee(PRIMARY_FORMULA, df, label, log) names = list(res.params.index) cc_col = [c for c in names if c.startswith("BMI:") and c.endswith("[content_confounded]")] core_col = [c for c in names if c.startswith("BMI:") and c.endswith(f"[{SET_REFERENCE}]")] require(len(cc_col) == 1 and len(core_col) == 1, f"Could not identify the two BMI x set slopes for {label}: " f"cc={cc_col} core={core_col}") cc_col, core_col = cc_col[0], core_col[0] contrast = np.zeros(len(names)) contrast[names.index(cc_col)] = 1.0 contrast[names.index(core_col)] = -1.0 ct = res.t_test(contrast) beta = float(np.ravel(ct.effect)[0]) se = float(np.ravel(ct.sd)[0]) pval = float(np.ravel(ct.pvalue)[0]) zval = float(np.ravel(ct.statistic)[0]) est = or_per_5(beta, se) cc_slope = or_per_5(float(res.params[cc_col]), float(res.bse[cc_col])) core_slope = or_per_5(float(res.params[core_col]), float(res.bse[core_col])) assert_full_rank(PRIMARY_FORMULA_REPARAM, df, f"{label} (reparameterised)", log) res2 = fit_gee(PRIMARY_FORMULA_REPARAM, df, f"{label} (reparameterised)", log) rep_col = [c for c in res2.params.index if c.startswith("BMI:") and c.endswith("[T.content_confounded]")] require(len(rep_col) == 1, f"Reparameterised interaction term not found for {label}: {rep_col}") rep_beta = float(res2.params[rep_col[0]]) rep_se = float(res2.bse[rep_col[0]]) log(f" contrast from SAP formula beta = {beta:+.6f} se = {se:.6f}") log(f" single coefficient, reparam beta = {rep_beta:+.6f} se = {rep_se:.6f}") require(abs(beta - rep_beta) < 1e-6 and abs(se - rep_se) < 1e-6, f"The explicit contrast and the reparameterised coefficient disagree " f"for {label}.") log(" agreement to 1e-6: yes") return { "n_rows": int(len(df)), "n_respondents": int(df["respondent_id"].nunique()), "criteria_in_fit": sorted(int(k) for k in df["criterion_id"].unique()), "formula": PRIMARY_FORMULA, "rank_assertion": rank_info, "converged": True, "slope_content_confounded": cc_slope, "slope_pharmacological_core": core_slope, "contrast": { **est, "z": zval, "p_value": pval, "significant_at_0_05": bool(pval < ALPHA), "direction": ("content-confounded > core" if beta > 0 else "content-confounded < core"), "scale": "ratio of odds ratios per 5 BMI units", "multiplicity_corrected": False, }, "reparameterisation_check": {"beta": rep_beta, "se": rep_se, "agrees": True}, } def _report_variant(name: str, title: str, r: dict[str, Any], log: Logger) -> None: log("") log(f" {name} — {title}") log(f" n respondents = {r['n_respondents']} stacked rows = {r['n_rows']} " f"criteria in fit = {r['criteria_in_fit']}") cc, core, c = r["slope_content_confounded"], r["slope_pharmacological_core"], r["contrast"] log(f" set-level slope, content_confounded : OR/5 = {cc['or_per_5']:.4f} " f"(95% CI {cc['ci_lo_per_5']:.4f} to {cc['ci_hi_per_5']:.4f}) " f"[log-odds per 1 unit = {cc['beta_per_1']:+.6f}]") log(f" set-level slope, pharmacological_core: OR/5 = {core['or_per_5']:.4f} " f"(95% CI {core['ci_lo_per_5']:.4f} to {core['ci_hi_per_5']:.4f}) " f"[log-odds per 1 unit = {core['beta_per_1']:+.6f}]") log(f" RATIO OF ODDS RATIOS, PER 5 BMI UNITS = {c['or_per_5']:.4f} " f"95% CI {c['ci_lo_per_5']:.4f} to {c['ci_hi_per_5']:.4f}") log(f" z = {c['z']:+.4f} p = {c['p_value']:.6g} " f"log-contrast per 1 BMI unit = {c['beta_per_1']:+.6f} (se {c['se_per_1']:.6f})") log(f" direction: {c['direction']} significant at alpha = {ALPHA}: " f"{c['significant_at_0_05']}") def section_j_c7_sensitivity(wide: pd.DataFrame, long: pd.DataFrame, primary: dict[str, Any], log: Logger) -> dict[str, Any]: log.rule("[J] POST-HOC SENSITIVITY — C7 MEASUREMENT DEFECT (NOT PRE-REGISTERED)") log("*** THIS SECTION IS NOT PRE-REGISTERED AND IS NOT A TEST OF H1. ***") log(f"Reference: {DEVIATION_2_REF}") log("") log("Added after the pre-registered analysis had been run and reported. The") log("instrument scored criterion 7 (withdrawal) from item 1 — 'I ate much more") log("than planned', a criterion-1 item — instead of item 11, the affective") log("withdrawal item, which was scored against a threshold and then discarded.") log("C7 as measured is therefore a contaminated indicator. Item 11's responses") log("were never persisted, so the frozen cohort cannot be rescored; this section") log("quantifies how much the pre-registered conclusion depends on C7 instead.") log("") log("DECLARED IN THE SAP BEFORE THIS WAS RUN: C6 and C11 are both strongly") log("negative, so dropping C7 from the pharmacological core is EXPECTED to move") log("the contrast further above 1. Any such movement is the mechanical") log("consequence of removing the one core criterion that ran above 1 — it is") log("NOT evidence for H1 and must not be written up as though it were.") log("") log("The pre-registered primary model in [B] is retained unchanged as THE result.") log("Nothing below is promoted to primary. The headline figure remains " f"{primary['confirmatory_contrast']['or_per_5']:.4f}.") # --- descriptive: C7 in the analytic cohort ---------------------------- log.sub("[J.0] DESCRIPTIVE — C7 in the analytic cohort") p_desc = float(wide["c7"].mean()) n_desc_pos = int(wide["c7"].sum()) elig = wide[wide["primary_model_eligible"] == 1] p_prim = float(elig["c7"].mean()) n_prim_pos = int(elig["c7"].sum()) log(f" C7 endorsement, descriptive cohort : {n_desc_pos} / {len(wide)} = " f"{100 * p_desc:.1f}% (the figure reported in [A] and pre-registered " "reporting)") log(f" C7 endorsement, primary-model-eligible : {n_prim_pos} / {len(elig)} = " f"{100 * p_prim:.1f}%") log("") log(" How many C7-positive respondents were positive SOLELY on item 1 cannot be") log(" computed here: item entries are persisted only where they crossed") log(" threshold AND only ~53% carry an `answer_index`, so an empty q-cell in the") log(" wide snapshot means 'not recorded', never 'sub-threshold'. That figure is") log(" computed directly from the stored JSON against the database, read-only,") log(" and reported alongside this run rather than inferred from the snapshot.") # --- S1 ---------------------------------------------------------------- log.sub("[J.1] S1 — SET-MEMBERSHIP ONLY (the primary sensitivity)") log("C7's rows are dropped from the stacked model; the pharmacological core") log(f"becomes {{C{PHARMACOLOGICAL_CORE_EX_C7[0]}, C{PHARMACOLOGICAL_CORE_EX_C7[1]}}}. " "rest_score and total_symptom_count are LEFT EXACTLY") log("AS IN [B], still computed over all eleven criteria. Everything else — the") log("formula, the eligibility filter, the working correlation, the clustering,") log("the contrast — is identical to the primary model.") log("") base = long[ long["criterion_set"].isin(CONFIRMATORY_SETS) & (long["primary_model_eligible"] == 1) ].copy() require(len(base) == 8 * EXPECTED_N_PRIMARY, f"[J] base frame is {len(base)} rows, expected {8 * EXPECTED_N_PRIMARY} — " "the sensitivity must start from exactly the primary model's rows.") require(len(base) == primary["n_rows"], "[J] base frame does not match the row count [B] fitted on.") s1 = base[base["criterion_id"] != C7].copy() require(sorted(int(k) for k in s1["criterion_id"].unique()) == sorted(CONTENT_CONFOUNDED + PHARMACOLOGICAL_CORE_EX_C7), "S1 criteria are not {content-confounded} u {C6, C11}.") require(len(s1) == 7 * EXPECTED_N_PRIMARY, f"S1 rows = {len(s1)}, expected {7 * EXPECTED_N_PRIMARY}.") require(int(s1["respondent_id"].nunique()) == EXPECTED_N_PRIMARY, "S1 respondent count changed — dropping C7's rows must not drop a " "respondent.") # The whole point of S1: the conditioning variable is untouched. merged = s1.merge(base[["respondent_id", "criterion_id", "rest_score", "total_symptom_count"]], on=["respondent_id", "criterion_id"], suffixes=("", "_base")) require(bool((merged["rest_score"] == merged["rest_score_base"]).all()) and bool((merged["total_symptom_count"] == merged["total_symptom_count_base"]).all()), "S1 altered rest_score or total_symptom_count — it must not.") log(" rest_score / total_symptom_count identical to [B] on every retained row: " "asserted") log(f" rest_score range in S1: {int(s1['rest_score'].min())} to " f"{int(s1['rest_score'].max())} (still over all eleven criteria)") log("") r_s1 = fit_h1_contrast(s1, "S1 sensitivity GEE", log) _report_variant("S1", "set-membership only (primary sensitivity)", r_s1, log) # --- S2 ---------------------------------------------------------------- log.sub("[J.2] S2 — FULL REMOVAL (secondary sensitivity)") log("C7 is dropped entirely: total_symptom_count and every rest_score are") log("recomputed over the remaining TEN criteria, so the contaminated variance is") log("removed from the conditioning variable as well. This is what the measurement") log("concern strictly implies, at the cost of no longer being the same model") log("conditioned differently — S2 and [B] are not the same model.") log("") c7_by_resp = (long[long["criterion_id"] == C7] .set_index("respondent_id")["criterion_met"]) require(int(c7_by_resp.index.nunique()) == len(c7_by_resp), "Duplicate C7 row per respondent in the long snapshot.") s2 = base.copy() c7_vec = s2["respondent_id"].map(c7_by_resp) require(bool(c7_vec.notna().all()), "A respondent in the primary frame has no C7 row to remove.") s2["total_symptom_count"] = (s2["total_symptom_count"] - c7_vec).astype(int) s2["rest_score"] = s2["total_symptom_count"] - s2["criterion_met"] s2 = s2[s2["criterion_id"] != C7].copy() require(len(s2) == 7 * EXPECTED_N_PRIMARY, f"S2 rows = {len(s2)}, expected {7 * EXPECTED_N_PRIMARY}.") require(bool(s2["total_symptom_count"].between(0, 10).all()), "S2 total_symptom_count outside 0-10 after removing C7.") require(bool(s2["rest_score"].between(0, 9).all()), "S2 rest_score outside 0-9 — with C7 removed the maximum rest-score over " "ten criteria is 9.") require(bool((s2["rest_score"] == s2["total_symptom_count"] - s2["criterion_met"]).all()), "S2 rest_score != recomputed total_symptom_count - criterion_met.") n_changed = int((s2["rest_score"].to_numpy() != s1["rest_score"].to_numpy()).sum()) log(f" total_symptom_count now over ten criteria; range " f"{int(s2['total_symptom_count'].min())} to " f"{int(s2['total_symptom_count'].max())}") log(f" rest_score range in S2: {int(s2['rest_score'].min())} to " f"{int(s2['rest_score'].max())}") log(f" rows whose rest_score differs from S1: {n_changed} of {len(s2)} " f"({100 * n_changed / len(s2):.1f}%) — S1 and S2 differ in the conditioning") log(" variable and in nothing else.") log("") r_s2 = fit_h1_contrast(s2, "S2 sensitivity GEE", log) _report_variant("S2", "full removal (secondary sensitivity)", r_s2, log) # --- side by side ------------------------------------------------------ log.sub("[J.3] SIDE BY SIDE — pre-registered primary against both sensitivities") pc = primary["confirmatory_contrast"] rows = [ ("[B] PRE-REGISTERED", "pre-registered confirmatory", pc, primary["slope_content_confounded"]["or_per_5"], primary["slope_pharmacological_core"]["or_per_5"], primary["n_rows"], primary["n_respondents"]), ("S1 sensitivity", "post-hoc, NOT pre-registered", r_s1["contrast"], r_s1["slope_content_confounded"]["or_per_5"], r_s1["slope_pharmacological_core"]["or_per_5"], r_s1["n_rows"], r_s1["n_respondents"]), ("S2 sensitivity", "post-hoc, NOT pre-registered", r_s2["contrast"], r_s2["slope_content_confounded"]["or_per_5"], r_s2["slope_pharmacological_core"]["or_per_5"], r_s2["n_rows"], r_s2["n_respondents"]), ] log(f" {'model':<20}{'cc OR/5':>9}{'core OR/5':>11}{'ROR/5':>9}" f"{'95% CI (per 5 BMI units)':>28}{'z':>9}{'p':>12}{'rows':>8}{'resp':>7}") for name, _, c, cc, core, nrows, nresp in rows: ci = f"{c['ci_lo_per_5']:.4f} to {c['ci_hi_per_5']:.4f}" log(f" {name:<20}{cc:>9.4f}{core:>11.4f}{c['or_per_5']:>9.4f}{ci:>28}" f"{c['z']:>+9.4f}{c['p_value']:>12.6g}{nrows:>8}{nresp:>7}") log("") log(" Only the first row is a pre-registered result. The other two are post-hoc") log(f" sensitivity fits under {DEVIATION_2_REF}.") log("") d1 = r_s1["contrast"]["or_per_5"] - pc["or_per_5"] d2 = r_s2["contrast"]["or_per_5"] - pc["or_per_5"] log(f" Movement in the ratio of odds ratios relative to the pre-registered " f"{pc['or_per_5']:.4f}:") log(f" S1 {r_s1['contrast']['or_per_5']:.4f} ({d1:+.4f})") log(f" S2 {r_s2['contrast']['or_per_5']:.4f} ({d2:+.4f})") log(" The SAP declared this direction of movement in advance. It is the") log(" mechanical consequence of removing the one pharmacological-core criterion") log(" whose own slope sat above 1, and it is NOT independent support for H1.") log("") agree_sign = ((r_s1["contrast"]["beta_per_1"] > 0) == (r_s2["contrast"]["beta_per_1"] > 0)) agree_sig = (r_s1["contrast"]["significant_at_0_05"] == r_s2["contrast"]["significant_at_0_05"]) log(f" S1 and S2 agree in sign: {agree_sign}") log(f" S1 and S2 agree in significance at alpha = {ALPHA}: {agree_sig}") log(f" |ROR(S1) - ROR(S2)| = {abs(r_s1['contrast']['or_per_5'] - r_s2['contrast']['or_per_5']):.4f}") log(" Whether that constitutes material disagreement is a reporting judgement") log(" for the write-up; the SAP requires any disagreement be reported rather") log(" than resolved by preferring whichever variant reads better. Neither") log(" variant is selected on the basis of what it returned.") return { "status": "POST-HOC SENSITIVITY ANALYSIS — NOT PRE-REGISTERED, NOT CONFIRMATORY", "deviation_reference": DEVIATION_2_REF, "added": "2026-07-26, after the pre-registered analysis had been run and reported", "reason": ("Criterion 7 was scored from item 1 (a criterion-1 item) instead of " "item 11 (the affective withdrawal item, scored then discarded). C7 " "as measured is a contaminated indicator and item 11 was never " "persisted, so the frozen cohort cannot be rescored."), "expectation_declared_in_advance": ( "C6 and C11 are both strongly negative, so excluding C7 was expected to " "move the contrast further above 1. Any such movement is mechanical and " "is not evidence for H1."), "primary_model_unchanged": True, "pharmacological_core_after_exclusion": list(PHARMACOLOGICAL_CORE_EX_C7), "c7_descriptive": { "n_positive_descriptive_cohort": n_desc_pos, "n_descriptive_cohort": int(len(wide)), "endorsement_rate_pct_descriptive": round(100 * p_desc, 2), "n_positive_primary_eligible": n_prim_pos, "n_primary_eligible": int(len(elig)), "endorsement_rate_pct_primary_eligible": round(100 * p_prim, 2), "positive_solely_via_item_1": ( "not computable from the frozen snapshot — item entries are stored " "only where they crossed threshold and only ~53% carry an " "answer_index, so an empty q-cell means 'not recorded'. Computed " "read-only from the stored JSON and reported alongside this run."), }, "s1_set_membership_only": { "designation": "primary sensitivity (more conservative; designated in the " "SAP before either variant was fitted)", "rest_score_basis": "unchanged — all eleven criteria, exactly as in [B]", **r_s1, }, "s2_full_removal": { "designation": "secondary sensitivity", "rest_score_basis": "recomputed over the remaining ten criteria", "n_rows_with_rest_score_differing_from_s1": n_changed, **r_s2, }, "comparison_with_preregistered": { "preregistered_ror_per_5": pc["or_per_5"], "s1_ror_per_5": r_s1["contrast"]["or_per_5"], "s2_ror_per_5": r_s2["contrast"]["or_per_5"], "s1_minus_preregistered": float(d1), "s2_minus_preregistered": float(d2), "s1_s2_agree_in_sign": bool(agree_sign), "s1_s2_agree_in_significance": bool(agree_sig), "abs_ror_difference_s1_s2": float( abs(r_s1["contrast"]["or_per_5"] - r_s2["contrast"]["or_per_5"])), }, } # --------------------------------------------------------------------------- # main # --------------------------------------------------------------------------- def main(argv: list[str] | None = None) -> int: parser = argparse.ArgumentParser( description="YFAS 2.0 criterion-level pre-registered analysis (SAP v1.0).") parser.add_argument( "--permute", type=int, default=None, metavar="SEED", help="Shuffle the `bmi` column across respondents with this seed before " "fitting anything. For RESULT-BLIND debugging only: the model " "structure is real, the association is null. Output is stamped " "PERMUTED throughout and must never be reported.") args = parser.parse_args(argv) permuted = args.permute is not None RESULTS_DIR.mkdir(parents=True, exist_ok=True) stamp = datetime.now().strftime("%Y%m%d-%H%M%S") tag = f"PERMUTED-seed{args.permute}_" if permuted else "" log = Logger(RESULTS_DIR / f"analysis_{tag}{stamp}.txt") banner = ("*** PERMUTED NULL RUN — NOT REAL RESULTS ***" if permuted else "REAL DATA RUN") log.rule("YFAS 2.0 CRITERION-LEVEL STUDY — PRE-REGISTERED ANALYSIS") log(f"Pre-registration : 02-statistical-analysis-plan.md v1.0 (2026-07-26), BINDING") log(f"Run mode : {banner}") if permuted: log("") log(" ####################################################################") log(f" # BMI HAS BEEN PERMUTED ACROSS RESPONDENTS (seed {args.permute}).") log(" # Every coefficient below is fitted against a NULL association.") log(" # This run exists to debug model structure result-blind (SAP §8).") log(" # DO NOT REPORT, CITE OR INTERPRET ANY NUMBER FROM THIS RUN.") log(" # Limit of the blind: H3 (age) and H4 (sex DIF) do not involve BMI,") log(" # so those two sections reproduce their real values even here.") log(" ####################################################################") log(f"Run timestamp : {stamp}") log("") log("Resolved environment (WSL venv ~/yfas-env):") log(f" Python {sys.version.split()[0]}") versions: dict[str, str | None] = {"python": sys.version.split()[0]} for pkg in ("pandas", "numpy", "scipy", "statsmodels", "girth", "patsy"): try: from importlib.metadata import version as _v versions[pkg] = _v(pkg) except Exception: # pragma: no cover versions[pkg] = None log(f" {pkg:<13} {versions[pkg]}") log("") log("Operationalisation decisions (coordinator rulings, see the module docstring):") log(" 1. age modelled as categorical C(age_range) — banded data, no invented") log(" midpoints. An implementation choice, not a §9 deviation.") log(" 2. BMI on its natural scale in every model; EVERY odds ratio printed is") log(" exponentiated at 5 BMI units and labelled 'per 5 BMI units'.") log(" 3. SAP §3 IRT/DIF: girth 2PL + Swaminathan-Rogers logistic DIF on the EAP") log(" ability estimate. Mantel-Haenszel fallback not needed; no §9 deviation.") log("") log("q1..q35 are never read: item responses were recorded only when the criterion") log("was met, so their missingness is deterministic on the outcome. Every model") log("below runs on the eleven binary criteria, exactly as SAP §2 and §3 specify.") log.rule("DATA — frozen snapshots only, no database connection (SAP §6, §8)") wide, long = load_data(args.permute, log) payload: dict[str, Any] = { "run_timestamp": stamp, "permuted": permuted, "permutation_seed": args.permute, "WARNING": ("PERMUTED NULL RUN — NOT REAL RESULTS" if permuted else "real analytic sample"), "preregistration": "02-statistical-analysis-plan.md v1.0", "environment": versions, "n_descriptive": int(len(wide)), "n_primary_model": int(wide["primary_model_eligible"].sum()), "n_long_rows": int(len(long)), } payload["a_power"] = section_a_power(wide, log) payload["b_primary_confirmatory"] = section_b_primary(long, log) payload["c_secondary_per_criterion"] = section_c_secondary(long, log) payload["d_unconditioned_descriptive"] = section_d_unconditioned( wide, payload["c_secondary_per_criterion"], log) payload["e_linearity"] = section_e_linearity(long, log) payload["f_irt_dif"] = section_f_irt_dif( wide, payload["c_secondary_per_criterion"], log) payload["g_h2"] = section_g_h2(wide, long, log, permuted) payload.update(section_h_h3_h4(wide, long, log)) payload["i_criterion_heterogeneity"] = section_i_heterogeneity( payload["b_primary_confirmatory"], payload["c_secondary_per_criterion"], log) # [J] is POST-HOC and NOT pre-registered — SAP §9 'Deviation 2'. It is run # last, additively, and consumes [B]'s already-computed result read-only so # that it cannot touch the primary model. payload["j_c7_sensitivity_post_hoc"] = section_j_c7_sensitivity( wide, long, payload["b_primary_confirmatory"], log) # ---------------------------------------------------------------- notes log.rule("SPECIFICATION NOTES — where the SAP proved ambiguous") notes = [ "SAP §2.2 asks that pharmacological_core be the reference level 'so the " "interaction term is read directly as content-confounded minus core'. With " "the BMI main effect absent (as §2.2 writes the formula), patsy codes " "C(criterion_set) inside the interaction with FULL dummies, so the fit " "yields two set-specific slopes and no single difference coefficient. The " "contrast is therefore formed explicitly and cross-checked against the " "algebraically identical reparameterisation `BMI + BMI:C(criterion_set)`, " "in which it IS one coefficient. Identical fit, agreement asserted to 1e-6. " "The model is unchanged; only the reporting route is.", "SAP §5 specifies NormalIndPower 'or a logistic-regression-appropriate power " "approximation ... the simulation-based approach is preferred where " "feasible'. NormalIndPower solves two-proportion problems and does not " "accept a continuous predictor, so the per-criterion table uses Hsieh's " "logistic-regression formula with the covariate variance-inflation factor " "1/(1-R^2), and every returned OR is verified back through " "statsmodels.stats.power.normal_power to reproduce exactly 80% power. " "NormalIndPower IS used where it fits the design: the two-proportion H2 " "power recomputation in [G].", "SAP §2.3 and §2.4 do not restate the §7.6 sex exclusion, which §7.6 itself " "scopes to 'the primary confirmatory GEE model'. The secondary and " "unconditioned models are therefore fitted on the full descriptive cohort " "(n=1,796, all sex levels); only the primary model uses n=1,773.", "SAP §4 does not name a reference group for the H2 comparison. " "obese_high_symptom is set as the reference so the coefficient reads " "'normal-BMI-high-symptom relative to obese-high-symptom', matching the " "direction in which §1 states H2 (OR < 1 supports H2).", "SAP §3 does not pre-specify multiplicity control for the DIF arm. Raw " "p-values are the pre-registered quantity; BH q-values are printed " "alongside as supplementary transparency, clearly labelled as such.", "SAP §4's power statement asks for the design's detectable difference at 80% " "power. Post-hoc observed power against the realised effect is not computed, " "because it is a deterministic function of the p-value and carries no " "information; the detectable-difference figure IS reported, for both the " "real cells and the SAP's pre-exclusion cells.", "SAP §1's H3 asks for an age x criterion profile with no direction. The " "joint Wald test on the interaction block plus the descriptive rate table " "is the implementation; no post-hoc contrast is pursued, since a " "significant block test with no pre-registered direction cannot support a " "directional claim.", "SAP §4 WAS NOT IMPLEMENTABLE AS ORIGINALLY WRITTEN — the most serious of " "these notes, now RESOLVED by amendment. §4 mandated an Exchangeable " "working correlation 'mirroring the primary model's approach for " "methodological consistency', but the H2 outcome is compositional within " "respondent and its within-cluster dependence is negative by construction, " "so the exchangeable dependence parameter goes negative and the fit " "diverges. The analogy to the primary model does not hold: there the " "within-respondent dependence is strongly positive. This was a defect in " "the specification, not in the data or the code. It is amended at SAP §9 " "'Deviation 1' (Independence working correlation, nothing else changed), " "authorised by the project owner on 2026-07-26. The failure that forced " "the amendment is retained in [G] step 1 as the audit trail, and the " "amended result in [G] step 2 carries the amendment notice inline.", "LIMIT OF THE RESULT-BLIND DEBUGGING PROTOCOL, recorded honestly: the " "--permute flag shuffles BMI only, so H3 (age) and H4 (sex DIF) reproduce " "their real values in a permuted run and were therefore not blinded during " "development. Neither involves an analytic choice that could have been " "steered by seeing them. Separately, permuting BMI re-cuts the H2 cells, " "which is why the SAP §4 convergence failure above did not surface until " "the real run.", ] for i, note in enumerate(notes, 1): log("") log(f"{i}. {note}") payload["specification_notes"] = notes log.rule("SAP §9 DEVIATIONS LOG STATUS") log("No departure from the pre-registration was MADE by this script. The") log("operationalisation decisions in the header are implementation choices the") log("SAP explicitly leaves open (§3 package choice, §5 power method, age's") log("functional form), not changes to a specified model.") deviations: list[dict[str, str]] = [] if not payload["g_h2"].get("prespecified_exchangeable_converged", True): deviations.append({ "date": "2026-07-26", "section": "SAP §4 (H2 confirmatory secondary)", "status": "CLOSED — authorised by the project owner, recorded at " "SAP §9 'Deviation 1'", "deviation": ( "H2 working correlation changed from Exchangeable() to " "Independence(). Nothing else changed — same GEE, same " "logit/binomial family, same formula including the mandatory " "total_symptom_count adjustment, same cluster-robust sandwich SEs " "on respondent_id, same population-averaged estimand." ), "reason": ( "Forced by non-convergence, not chosen. The H2 outcome is " "compositional within respondent, so within-cluster dependence is " "negative by construction (rho = " f"{payload['g_h2'].get('empirical_within_cluster_rho', float('nan')):.4f} " "against the exchangeable bound " f"{payload['g_h2'].get('exchangeable_lower_bound_at_max_cluster', float('nan')):.4f} " "at the observed maximum cluster size), and the fit diverges." ), }) log("") log("ONE DEVIATION IS ON RECORD, AND IT IS CLOSED:") log(" SAP §9 'Deviation 1' (2026-07-26) — H2 working correlation changed") log(" from Exchangeable() to Independence(), nothing else. Forced by the") log(" non-convergence documented in [G] step 1, authorised by the project") log(" owner, and disclosed there in full including the fact that the ruling") log(" was made with the Independence result already visible.") log(" The amended H2 result is reported in [G] step 2 with that amendment") log(" annotated inline, never as though it were the original fit.") else: log("Nothing here requires a §9 entry.") payload["deviations_required"] = deviations footer = ("*** END OF PERMUTED NULL RUN — NOT REAL RESULTS ***" if permuted else "*** END OF REAL ANALYSIS RUN ***") log.rule(footer) log.flush() json_path = RESULTS_DIR / f"analysis_{tag}{stamp}.json" json_path.write_text(json.dumps(payload, indent=2, default=jsonable), encoding="utf-8") print(f"\nResults written to: {log.path}") print(f"Machine-readable : {json_path}") return 0 if __name__ == "__main__": sys.exit(main())