#!/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. # # This is the script that produced the published results, with exactly two # changes, both confined to how the database is reached: # # 1. `DB_CONFIG` read literal connection parameters. Here it reads them from # the environment, so that no credentials are distributed. # 2. The "Run as" line below has been updated to match. # # Nothing else has been altered. The query, the SAP §8 freeze, every exclusion # and every computation are byte-for-byte as they ran. # # The individual responses this script reads are NOT published, and running it # will not reproduce them: it reads Creative Touch's own database, and we # publish aggregated findings only. See the article's data and methods note. # The script is published so that the exclusions, the freeze and the derivation # of every variable can be checked line by line. # --------------------------------------------------------------------------- """ YFAS 2.0 criterion-level study — data extraction (step 1 of 2). Implements the following sections of the committed pre-registration (`../02-statistical-analysis-plan.md`, SAP v1.0, committed 2026-07-26): * SAP §6 — data pull from the `creativetouch` MariaDB `calculators` table, read-only, exclusions applied in code, timestamped CSV snapshot. * SAP §7 — pre-specified exclusions, applied in the order given there, with the running sample size logged at every step. * SAP §8 — data freeze. `date_create <= '2026-07-25 23:59:59'` is hard-coded as a literal WHERE clause below. It is deliberately NOT parameterised, NOT configurable and NOT derived from wall-clock time, so that the freeze is auditable from the source file alone. * SAP §5 — the first descriptive step only: per-criterion endorsement rates, BMI/WHO distribution, sex and age breakdowns, discordant-cell sizes. These are explicitly non-inferential and explicitly permitted before modelling (SAP §5 and §8). THIS SCRIPT COMPUTES NOTHING INFERENTIAL. No correlation, no regression, no model fit, no p-value, no power calculation. Those belong to the analysis script and running them here would breach the pre-registration's no-peeking rule (SAP §8). --------------------------------------------------------------------------- Reproducibility — resolved package versions in the WSL venv `~/yfas-env` (recorded here rather than in a requirements.txt, which would imply a reinstall; these are the versions the committed run used): Python 3.10.12 pandas 2.3.3 numpy 2.2.6 scipy 1.15.3 statsmodels 0.14.6 PyMySQL 1.2.0 girth 0.8.0 patsy 1.0.2 Run as (see the published-copy note above — the credentials were literals in the run that produced the published results): YFAS_DB_USER=... YFAS_DB_PASSWORD=... ~/yfas-env/bin/python extract.py --------------------------------------------------------------------------- Ordering note (not a deviation). SAP §7 lists the `under_18` exclusion as step 1 and the complete-case cohort restriction as step 5, but §7 step 1 states its own count as "99 of 1,898 in the research cohort (BMI + age_range + sex all present)" — i.e. the SAP itself defines the under-18 exclusion as operating *within* the complete-case cohort. `age_range` is one of the three complete-case fields, so the under-18 rule is not evaluable before the cohort is defined. The cohort restriction is therefore applied first and the exclusions then run in the SAP's stated order. This reproduces the SAP's own arithmetic exactly and is a reading of §7, not a change to it. Corrected data facts (established by this script's own structural assertions, 2026-07-26 — they contradict the briefing note this script was written from, and the correction is recorded here so no later reader re-inherits the error): * The item-level data is NOT a full 35-item response matrix. `questions` is non-empty if and only if its criterion was met, and lists only the items that crossed threshold. Verified 29,023 met-with-items / 38,132 unmet-with-none across the frozen set, zero off-diagonal. * `answer_index` is present on only ~53% of recorded item entries; older submissions carry the `answer` text without an index. * Where present, `answer_index` only takes values 2-7, never 0 or 1 — sub-threshold responses are never written. * Item 11 never appears anywhere; item 1 is shared by criteria 1 and 7. None of this obstructs the pre-registration: SAP §3 specifies the IRT/DIF arm on the eleven binary criteria, not on the 35 items. No SAP §9 deviation arises. British English throughout. """ from __future__ import annotations import json import os import sys from collections import Counter from datetime import datetime from pathlib import Path from typing import Any, Iterable import pandas as pd import pymysql # --------------------------------------------------------------------------- # Constants — pre-registered, not configurable # --------------------------------------------------------------------------- #: SAP §8. Hard-coded data freeze. Do not parameterise. Do not move. FREEZE_SQL = ( "SELECT id, code, series_code, date_create, valid, deletion_requested, " " country, campaign, result " "FROM calculators " "WHERE module = 'food_addiction' " " AND date_create <= '2026-07-25 23:59:59' " "ORDER BY id" ) #: Connection parameters. In the run that produced the published results these #: were literals in this file; in this published copy they are read from the #: environment so that no credentials are distributed. A missing user or #: password raises KeyError naming the variable. Nothing downstream of the #: connection differs. DB_CONFIG: dict[str, Any] = { "host": os.environ.get("YFAS_DB_HOST", "127.0.0.1"), "user": os.environ["YFAS_DB_USER"], "password": os.environ["YFAS_DB_PASSWORD"], "database": os.environ.get("YFAS_DB_NAME", "creativetouch"), "charset": "utf8mb4", "cursorclass": pymysql.cursors.DictCursor, } #: The eleven DSM-5-adapted criteria. Key "12" in `dsmCriteriaMet` is the #: clinical-significance criterion and is NOT one of the eleven — it must never #: enter `total_symptom_count` or any rest-score. This is the single most #: dangerous trap in the dataset. CRITERIA = tuple(range(1, 12)) CLINICAL_SIGNIFICANCE_KEY = "12" EXPECTED_CRITERION_KEYS = {str(k) for k in range(1, 13)} N_ITEMS = 35 #: SAP §2.2 criterion sets. CONTENT_CONFOUNDED = {4, 5, 8, 9, 10} PHARMACOLOGICAL_CORE = {6, 7, 11} NEITHER = {1, 2, 3} #: SAP §7 step 1. EXCLUDED_AGE_RANGE = "under_18" VALID_AGE_RANGES = { "under_18", "18-24", "25-34", "35-44", "45-54", "55-64", "65-74", "75+", } VALID_SEX = {"female", "male", "prefer_not_to_say", "other"} #: SAP §7 step 6 — excluded from the primary model only, retained descriptively. NON_PRIMARY_SEX = {"other", "prefer_not_to_say"} #: SAP §7 step 3. BMI_PLAUSIBLE_MIN = 12.0 BMI_PLAUSIBLE_MAX = 80.0 #: SAP §7 step 4 — non-overlapping WHO cut-points. 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), ) BASE_DIR = Path(__file__).resolve().parent SNAPSHOT_DIR = BASE_DIR / "snapshots" # --------------------------------------------------------------------------- # Logging + failure helpers # --------------------------------------------------------------------------- class Logger: """Tee to stdout and to the audit log 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("=" * 78) if title: self(title) self("=" * 78) def flush(self) -> None: self.path.write_text("\n".join(self._lines) + "\n", encoding="utf-8") 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. """ def require(condition: bool, message: str) -> None: """Explicit, non-strippable assertion (survives `python -O`).""" if not condition: raise SapDeparture(message) # --------------------------------------------------------------------------- # Parsing # --------------------------------------------------------------------------- def who_category(bmi: float) -> str: """SAP §7 step 4 — non-overlapping WHO cut-points, lower bound inclusive.""" 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") def parse_row(row: dict[str, Any]) -> dict[str, Any]: """Flatten one `calculators` row into a respondent record. Raises `ValueError` for malformed/partially-missing JSON so the caller can count it into the funnel rather than skipping it silently. """ raw = row["result"] if raw is None or (isinstance(raw, str) and not raw.strip()): raise ValueError("result column is NULL or empty") try: payload = json.loads(raw) except (TypeError, ValueError) as exc: raise ValueError(f"result is not valid JSON: {exc}") from exc if not isinstance(payload, dict): raise ValueError("result JSON is not an object") criteria = payload.get("dsmCriteriaMet") if not isinstance(criteria, dict): raise ValueError("dsmCriteriaMet missing or not an object") # Structural assertion: keys are exactly "1".."12". keys = set(criteria.keys()) if keys != EXPECTED_CRITERION_KEYS: raise ValueError( f"dsmCriteriaMet keys are {sorted(keys)}, expected '1'..'12'" ) record: dict[str, Any] = { "respondent_id": row["id"], "code": row["code"], "date_create": row["date_create"], "country": row["country"], "campaign": row["campaign"], "_valid": row["valid"], "_deletion_requested": row["deletion_requested"], } # --- criteria 1..11 ----------------------------------------------------- # # IMPORTANT — the item-level data is ENDORSEMENT-CONDITIONAL, not a full # 35-item response matrix. Verified across all 6,105 frozen rows: # * `questions` is non-empty if and only if that criterion was met # (29,023 met-with-items vs 38,132 unmet-with-none; zero off-diagonal), # and lists only the items that crossed the criterion's threshold; # * `answer_index` is present on only ~53% of recorded item entries # (older submissions store `answer` text only); # * where present, `answer_index` only ever takes values 2-7, never 0/1, # because sub-threshold responses are never written at all; # * item 11 never appears anywhere in the dataset, and item 1 is shared # by criteria 1 and 7 (values always agree where duplicated). # # An empty q-cell therefore means "not recorded", NEVER "answered 0". See # the caveat block emitted in the log and `cohort_summary.json`. items: dict[int, int] = {} n_items_recorded = 0 for k in CRITERIA: block = criteria[str(k)] if not isinstance(block, dict) or "met" not in block: raise ValueError(f"criterion {k} block malformed") met = block["met"] if met not in (0, 1, "0", "1", True, False): raise ValueError(f"criterion {k} 'met' is {met!r}, expected 0/1") met = int(met) record[f"c{k}"] = met questions = block.get("questions") if not isinstance(questions, list): raise ValueError(f"criterion {k} 'questions' is not a list") # Verified structural invariant across the whole frozen set. if bool(questions) != bool(met): raise ValueError( f"criterion {k}: met={met} but {len(questions)} item(s) recorded " "(invariant 'items recorded iff criterion met' broken)" ) for q in questions: if not isinstance(q, dict): raise ValueError(f"criterion {k} question entry is not an object") if "question" not in q: raise ValueError(f"criterion {k} item entry has no 'question' key") try: qnum = int(q["question"]) except (TypeError, ValueError) as exc: raise ValueError(f"criterion {k} 'question' not an integer: {exc}") from exc if not 1 <= qnum <= N_ITEMS: raise ValueError(f"item number {qnum} outside 1..{N_ITEMS}") n_items_recorded += 1 if "answer_index" not in q: continue # older schema: `answer` text only, no index try: aidx = int(q["answer_index"]) except (TypeError, ValueError) as exc: raise ValueError(f"item {qnum} answer_index not an integer: {exc}") from exc if not 0 <= aidx <= 7: raise ValueError(f"item {qnum} answer_index {aidx} outside 0-7") if qnum in items and items[qnum] != aidx: raise ValueError( f"item {qnum} recorded twice with conflicting answer_index " f"({items[qnum]} vs {aidx})" ) items[qnum] = aidx for qnum, aidx in items.items(): record[f"q{qnum}"] = aidx record["n_items_recorded"] = n_items_recorded record["n_items_with_answer_index"] = len(items) # --- clinical significance (criterion key "12", never a symptom) -------- cs_block = criteria[CLINICAL_SIGNIFICANCE_KEY] if not isinstance(cs_block, dict) or "met" not in cs_block: raise ValueError("clinical-significance block (key '12') malformed") record["clinical_significance"] = int(cs_block["met"]) record["total_symptom_count"] = sum(record[f"c{k}"] for k in CRITERIA) record["_stored_symptom_count"] = payload.get("symptomCount") # --- demographics ------------------------------------------------------- demo = payload.get("demographics") if not isinstance(demo, dict): record["_has_demographics"] = False return record record["_has_demographics"] = True height = demo.get("height") or {} weight = demo.get("weight") or {} bmi_block = demo.get("bmi") or {} record["age_range"] = demo.get("age_range") record["sex"] = demo.get("sex") record["height_cm"] = height.get("cm") record["weight_kg"] = weight.get("kg") record["stored_bmi_value"] = bmi_block.get("value") record["stored_bmi_category"] = bmi_block.get("category") h_unit = height.get("unit") w_unit = weight.get("unit") record["unit_system"] = h_unit # Brief-verified assumption A: units are never mixed within a row. if h_unit is not None and w_unit is not None and h_unit != w_unit: raise ValueError( f"mixed unit systems: height={h_unit!r} weight={w_unit!r} " "(brief asserted rows are metric/metric or imperial/imperial)" ) # Brief-verified assumption B: cm and kg are always populated whenever a # unit is set, including for imperial rows (already converted upstream). # We therefore never convert; we recompute directly from cm/kg. if h_unit is not None and record["height_cm"] in (None, ""): raise ValueError(f"unit_system={h_unit!r} but height.cm is empty") if w_unit is not None and record["weight_kg"] in (None, ""): raise ValueError(f"unit_system={w_unit!r} but weight.kg is empty") return record def fetch_rows() -> list[dict[str, Any]]: """Read-only pull under the hard-coded SAP §8 freeze.""" conn = pymysql.connect(**DB_CONFIG) try: with conn.cursor() as cur: cur.execute(FREEZE_SQL) return list(cur.fetchall()) finally: conn.close() # --------------------------------------------------------------------------- # Descriptives (SAP §5 first step — non-inferential only) # --------------------------------------------------------------------------- def counts_table(log: Logger, series: pd.Series, title: str) -> dict[str, int]: log("") log(f"{title}:") vc = series.value_counts(dropna=False) total = int(vc.sum()) out: dict[str, int] = {} for key, n in vc.sort_index().items(): label = "MISSING" if pd.isna(key) else str(key) pct = 100.0 * n / total if total else 0.0 log(f" {label:<22} {int(n):>6} ({pct:5.1f}%)") out[label] = int(n) return out 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" # --------------------------------------------------------------------------- # main # --------------------------------------------------------------------------- def main() -> int: SNAPSHOT_DIR.mkdir(parents=True, exist_ok=True) stamp = datetime.now().strftime("%Y%m%d-%H%M%S") log = Logger(SNAPSHOT_DIR / f"extraction_log_{stamp}.txt") summary: dict[str, Any] = {"snapshot_timestamp": stamp} log.rule("YFAS 2.0 criterion-level study — extraction (SAP §5, §6, §7, §8)") log("Pre-registration: 02-statistical-analysis-plan.md v1.0 (2026-07-26)") log("This script is NON-INFERENTIAL. No model, no p-value, no correlation.") log("") log("Resolved environment (WSL venv ~/yfas-env):") log(f" Python {sys.version.split()[0]}") for pkg in ("pandas", "numpy", "scipy", "statsmodels", "PyMySQL", "girth", "patsy"): try: from importlib.metadata import version as _v log(f" {pkg:<13} {_v(pkg)}") except Exception: # pragma: no cover - reporting only log(f" {pkg:<13} (not resolvable)") log("") log("SAP §8 data freeze, hard-coded WHERE clause:") log(" module = 'food_addiction' AND date_create <= '2026-07-25 23:59:59'") # -- Step 0: freeze pull ------------------------------------------------- log.rule("EXCLUSION FUNNEL (SAP §7, applied in the order stated there)") rows = fetch_rows() n_frozen = len(rows) log(f"[0] Rows under the freeze n = {n_frozen}") summary["n_frozen"] = n_frozen require(n_frozen > 0, "Freeze returned zero rows — the freeze clause or DB is wrong.") max_dc = max(r["date_create"] for r in rows) log(f" max(date_create) in frozen set {max_dc}") summary["max_date_create"] = str(max_dc) # -- Step 0b: JSON parse ------------------------------------------------- parsed: list[dict[str, Any]] = [] malformed: list[tuple[int, str]] = [] for row in rows: try: parsed.append(parse_row(row)) except ValueError as exc: malformed.append((row["id"], str(exc))) n_parsed = len(parsed) log(f"[1] Parsed `result` JSON successfully n = {n_parsed}" f" (dropped {len(malformed)})") summary["n_malformed_result_json"] = len(malformed) summary["malformed_result_ids"] = [i for i, _ in malformed][:50] if malformed: log(" Malformed/partially-missing rows (id: reason) — reported, not silent:") reasons = Counter(reason for _, reason in malformed) for reason, n in reasons.most_common(): log(f" {n:>5} x {reason}") for rid, reason in malformed[:20]: log(f" id={rid}: {reason}") if len(malformed) > 20: log(f" ... and {len(malformed) - 20} more") summary["malformed_reason_counts"] = dict( Counter(reason for _, reason in malformed) ) require(bool(parsed), "Every row failed to parse — the parser is wrong, not the data.") df = pd.DataFrame(parsed) for i in range(1, N_ITEMS + 1): if f"q{i}" not in df.columns: df[f"q{i}"] = pd.NA # -- Item-level data caveat (a data fact, not a SAP departure) ----------- log("") log(" ITEM-LEVEL DATA CAVEAT — read before using q1..q35:") log(" The `questions` arrays are ENDORSEMENT-CONDITIONAL. Items are") log(" recorded if and only if their criterion was met, and only the") log(" threshold-crossing items are stored. There is NO full 35-item") log(" response matrix in this dataset. An empty q-cell means") log(" 'not recorded', never 'answered 0'.") item_cov = { f"q{i}": int(df[f"q{i}"].notna().sum()) for i in range(1, N_ITEMS + 1) } n_any_item = int(df["n_items_recorded"].gt(0).sum()) n_any_idx = int(df["n_items_with_answer_index"].gt(0).sum()) log(f" rows with >=1 item recorded {n_any_item}") log(f" rows with >=1 usable `answer_index` {n_any_idx}") log(f" total item entries recorded " f"{int(df['n_items_recorded'].sum())}") log(f" of which carry `answer_index` " f"{int(df['n_items_with_answer_index'].sum())}") log(f" items with no usable answer_index anywhere " f"{[k for k, v in item_cov.items() if v == 0]}") log(" (q11 never appears in the dataset at all; q16/q17 belong solely") log(" to criterion key '12', the clinical-significance criterion,") log(" which is deliberately not harvested as a symptom criterion.)") log(" Consequence: SAP §3's IRT/DIF arm is specified on the ELEVEN") log(" BINARY CRITERIA, not the 35 items, so this does not obstruct") log(" the pre-registered analysis. Any 35-item psychometric use would") log(" be structurally missing-not-at-random and is out of scope.") summary["item_level_data_caveat"] = { "endorsement_conditional": True, "full_35_item_matrix_available": False, "rows_with_any_item_recorded": n_any_item, "rows_with_any_answer_index": n_any_idx, "total_item_entries": int(df["n_items_recorded"].sum()), "total_item_entries_with_answer_index": int( df["n_items_with_answer_index"].sum() ), "per_item_answer_index_coverage_frozen_set": item_cov, "items_never_observed": [k for k, v in item_cov.items() if v == 0], } # symptomCount integrity check (data-integrity signal, not an exclusion) stored = pd.to_numeric(df["_stored_symptom_count"], errors="coerce") mismatch_mask = stored.notna() & (stored != df["total_symptom_count"]) n_stored_missing = int(stored.isna().sum()) mismatch_ids = df.loc[mismatch_mask, "respondent_id"].tolist() log("") log(f" Integrity: stored `symptomCount` vs sum(c1..c11)") log(f" mismatches {len(mismatch_ids)}") log(f" stored symptomCount missing/non-numeric {n_stored_missing}") if mismatch_ids: log(f" example ids: {mismatch_ids[:20]}") ex = df.loc[mismatch_mask, ["respondent_id", "total_symptom_count"]].head(10) for _, r in ex.iterrows(): sid = stored.loc[r.name] log(f" id={int(r['respondent_id'])} stored={sid} recomputed=" f"{int(r['total_symptom_count'])}") summary["n_symptomcount_mismatches"] = len(mismatch_ids) summary["symptomcount_mismatch_ids"] = mismatch_ids[:100] summary["n_stored_symptomcount_missing"] = n_stored_missing # -- SAP §7 step 5: complete-case cohort --------------------------------- has_demo = df["_has_demographics"].fillna(False) log("") log(f"[2] Rows carrying a `demographics` object n = {int(has_demo.sum())}") for col in ("age_range", "sex", "height_cm", "weight_kg", "stored_bmi_value", "stored_bmi_category", "unit_system"): if col not in df.columns: df[col] = pd.NA df["height_cm"] = pd.to_numeric(df["height_cm"], errors="coerce") df["weight_kg"] = pd.to_numeric(df["weight_kg"], errors="coerce") df["stored_bmi_value"] = pd.to_numeric(df["stored_bmi_value"], errors="coerce") complete = ( df["stored_bmi_value"].notna() & df["age_range"].notna() & df["sex"].notna() & df["height_cm"].notna() & df["weight_kg"].notna() ) cohort = df.loc[complete].copy() n_cohort = len(cohort) log(f"[3] SAP §7 step 5 — complete-case cohort n = {n_cohort}") log(" (BMI value + age_range + sex all present; no imputation)") summary["n_complete_case_cohort"] = n_cohort # Vocabulary assertions bad_age = sorted(set(cohort["age_range"].unique()) - VALID_AGE_RANGES) bad_sex = sorted(set(cohort["sex"].unique()) - VALID_SEX) require(not bad_age, f"Unexpected age_range values in cohort: {bad_age}") require(not bad_sex, f"Unexpected sex values in cohort: {bad_sex}") # Unit-system homogeneity across the cohort (brief-verified assumption) unit_counts = cohort["unit_system"].value_counts(dropna=False).to_dict() log(f" unit systems in cohort: " f"{ {str(k): int(v) for k, v in unit_counts.items()} }") summary["unit_system_counts_cohort"] = { str(k): int(v) for k, v in unit_counts.items() } # -- Data-integrity guard (NOT a new exclusion rule) --------------------- guard_fail = cohort.loc[ (cohort["_valid"] != 1) | (cohort["_deletion_requested"] != 0) ] log("") log(f"[3a] Data-integrity guard: valid = 1 AND deletion_requested = 0") log(f" rows removed from the complete-case cohort {len(guard_fail)}") if len(guard_fail) > 0: ids = guard_fail["respondent_id"].tolist() raise SapDeparture( "DATA-INTEGRITY GUARD TRIPPED. The filter " "`valid = 1 AND deletion_requested = 0` removed " f"{len(guard_fail)} row(s) from the complete-case cohort " f"(ids: {ids[:50]}). The SAP §7 exclusion list contains no such " "exclusion, so dropping these rows would constitute a DEPARTURE " "FROM THE PRE-REGISTRATION and requires an entry in SAP §9 " "(Deviations Log) — date, section, change, reason — before the " "analysis may proceed. Aborting rather than silently dropping." ) summary["n_removed_by_integrity_guard"] = 0 # -- SAP §7 step 1: under_18 exclusion ----------------------------------- n_under18_cohort = int((cohort["age_range"] == EXCLUDED_AGE_RANGE).sum()) n_under18_frozen = int((df["age_range"] == EXCLUDED_AGE_RANGE).sum()) cohort = cohort.loc[cohort["age_range"] != EXCLUDED_AGE_RANGE].copy() log("") log(f"[4] SAP §7 step 1 — exclude age_range='under_18' n = {len(cohort)}" f" (removed {n_under18_cohort})") log(f" under_18 in the whole frozen parsed set {n_under18_frozen}") log(f" under_18 within the complete-case cohort {n_under18_cohort}") log(" SAP §7 step 1 quotes 99 as the within-cohort figure.") summary["n_under18_in_frozen_set"] = n_under18_frozen summary["n_under18_in_cohort"] = n_under18_cohort summary["n_after_under18_exclusion"] = len(cohort) # -- SAP §7 step 2: BMI recomputation ------------------------------------ require( bool((cohort["height_cm"] > 0).all()), "Non-positive height_cm present in cohort — cannot recompute BMI.", ) require( bool((cohort["weight_kg"] > 0).all()), "Non-positive weight_kg present in cohort — cannot recompute BMI.", ) cohort["bmi"] = cohort["weight_kg"] / (cohort["height_cm"] / 100.0) ** 2 log("") log(f"[5] SAP §7 step 2 — BMI recomputed from height.cm / weight.kg") log(" Formula: kg / (cm/100)^2 for every row, metric and imperial alike") log(" (imperial rows already carry converted cm/kg upstream — asserted,") log(" no conversion branch exists in this script).") log(f" recomputed BMI range min={cohort['bmi'].min():.2f} " f"max={cohort['bmi'].max():.2f}") summary["recomputed_bmi_min_pre_screen"] = round(float(cohort["bmi"].min()), 4) summary["recomputed_bmi_max_pre_screen"] = round(float(cohort["bmi"].max()), 4) # Audit the recomputation against the stored value/category. cohort["bmi_who_category"] = cohort["bmi"].map(who_category) stored_cat = cohort["stored_bmi_category"].astype("string") recomputed_cat = cohort["bmi_who_category"].astype("string") # Stored categories use a coarser vocabulary in places; compare on a # normalised basis so the count is meaningful rather than cosmetic. def normalise(c: object) -> str: s = str(c).strip().lower().replace(" ", "_").replace("-", "_") aliases = { "obese_class_1": "obese_i", "obese_class_i": "obese_i", "obese_1": "obese_i", "obese1": "obese_i", "obese_i": "obese_i", "obese_class_2": "obese_ii", "obese_class_ii": "obese_ii", "obese_2": "obese_ii", "obese2": "obese_ii", "obese_ii": "obese_ii", "obese_class_3": "obese_iii", "obese_class_iii": "obese_iii", "obese_3": "obese_iii", "obese3": "obese_iii", "obese_iii": "obese_iii", "obese": "obese_i", "normal_weight": "normal", "healthy": "normal", "normal": "normal", "under_weight": "underweight", "underweight": "underweight", "over_weight": "overweight", "overweight": "overweight", } return aliases.get(s, s) unmapped = sorted( { str(c) for c in cohort["stored_bmi_category"].unique() if normalise(c) not in {w[0] for w in WHO_CATEGORIES} } ) require( not unmapped, "Unrecognised stored_bmi_category value(s) — the audit comparison would " f"be meaningless: {unmapped}", ) cohort["_stored_cat_norm"] = cohort["stored_bmi_category"].map(normalise) cat_changed = cohort["_stored_cat_norm"] != cohort["bmi_who_category"] n_cat_changed = int(cat_changed.sum()) bmi_delta = (cohort["bmi"] - cohort["stored_bmi_value"]).abs() log("") log(f" Audit vs stored values (stored values are NEVER used for") log(f" classification — SAP §7 step 2; retained in the CSV for audit only):") log(f" rows changing WHO category under recomputation {n_cat_changed}" f" ({100.0 * n_cat_changed / len(cohort):.2f}%)") log(f" |recomputed BMI - stored BMI| max={bmi_delta.max():.3f} " f"mean={bmi_delta.mean():.3f}") # Quantify exposure to the documented stored-category boundary-overlap bug # (00-STATE.md §3 issue 1): how many rows sit exactly on a WHO cut-point, # where the stored classifier's overlapping ranges are ambiguous? boundary_counts: dict[str, int] = {} for cut in (18.5, 25.0, 30.0, 35.0, 40.0): n_on = int((cohort["stored_bmi_value"] == cut).sum()) n_near = int((cohort["bmi"] - cut).abs().le(0.05).sum()) boundary_counts[f"{cut}"] = n_on log(f" rows at WHO cut-point {cut:<5} stored=={cut:<5} {n_on:>4}" f" |recomputed-cut|<=0.05 {n_near:>4}") log(f" total rows sitting exactly on a cut-point " f"{sum(boundary_counts.values())}") summary["rows_on_who_cutpoint_stored"] = boundary_counts if n_cat_changed: pairs = Counter( zip( cohort.loc[cat_changed, "_stored_cat_norm"], cohort.loc[cat_changed, "bmi_who_category"], ) ) log(" stored -> recomputed transitions:") for (a, b), n in pairs.most_common(): log(f" {a:<14} -> {b:<14} {n:>5}") summary["who_category_transitions"] = { f"{a}->{b}": int(n) for (a, b), n in pairs.items() } summary["n_who_category_changed"] = n_cat_changed summary["max_abs_bmi_delta_vs_stored"] = round(float(bmi_delta.max()), 4) # -- SAP §7 step 3: plausibility screen ---------------------------------- implausible = (cohort["bmi"] < BMI_PLAUSIBLE_MIN) | (cohort["bmi"] > BMI_PLAUSIBLE_MAX) n_implausible = int(implausible.sum()) implausible_ids = cohort.loc[implausible, "respondent_id"].tolist() implausible_vals = [round(float(v), 2) for v in cohort.loc[implausible, "bmi"]] cohort = cohort.loc[~implausible].copy() log("") log(f"[6] SAP §7 step 3 — plausibility screen (<12 or >80) n = {len(cohort)}" f" (removed {n_implausible})") if n_implausible: log(f" excluded ids: {implausible_ids}") log(f" excluded recomputed BMI values: {implausible_vals}") summary["n_implausible_bmi"] = n_implausible summary["implausible_bmi_ids"] = implausible_ids summary["implausible_bmi_values"] = implausible_vals # -- SAP §7 step 4: WHO categories (already assigned, re-assert) --------- cohort["bmi_who_category"] = cohort["bmi"].map(who_category) log("") log("[7] SAP §7 step 4 — non-overlapping WHO categories applied to the") log(" recomputed BMI (underweight <18.5; normal 18.5-<25; overweight") log(" 25-<30; obese I 30-<35; obese II 35-<40; obese III >=40).") # -- SAP §7 step 6: sex exclusion (primary model only) ------------------- non_primary = cohort["sex"].isin(NON_PRIMARY_SEX) n_non_primary = int(non_primary.sum()) cohort["primary_model_eligible"] = (~non_primary).astype(int) log("") log(f"[8] SAP §7 step 6 — sex in {sorted(NON_PRIMARY_SEX)} excluded from the") log(f" PRIMARY MODEL ONLY, retained descriptively. count = {n_non_primary}") log(f" n eligible for the primary confirmatory GEE n = " f"{int((~non_primary).sum())}") log(f" n in the descriptive analytic cohort n = {len(cohort)}") summary["n_sex_other_or_prefer_not_to_say"] = n_non_primary summary["n_primary_model_eligible"] = int((~non_primary).sum()) # -- SAP §7 step 7: no deduplication ------------------------------------- dup_codes = int(cohort["code"].duplicated().sum()) log("") log(f"[9] SAP §7 step 7 — NO deduplication is performed (declared") log(f" limitation). Duplicate `code` values in cohort: {dup_codes}") summary["n_duplicate_codes_in_cohort"] = dup_codes n_analytic = len(cohort) log("") log(f"FINAL ANALYTIC n (descriptive cohort) n = {n_analytic}") log(f"FINAL ANALYTIC n (primary confirmatory model) n = " f"{summary['n_primary_model_eligible']}") summary["n_analytic"] = n_analytic # ----------------------------------------------------------------------- # SAP §5 first descriptive step — non-inferential only # ----------------------------------------------------------------------- log.rule("SAP §5 FIRST DESCRIPTIVE STEP (non-inferential; no test is run)") log("") log("Per-criterion endorsement rates in the analytic cohort:") log(f" {'crit':<6}{'set':<22}{'n met':>8}{'rate %':>10}") endorsement: dict[str, dict[str, float]] = {} for k in CRITERIA: n_met = int(cohort[f"c{k}"].sum()) rate = 100.0 * n_met / n_analytic cset = criterion_set_of(k) log(f" C{k:<5}{cset:<22}{n_met:>8}{rate:>10.1f}") endorsement[f"c{k}"] = { "criterion_set": cset, "n_met": n_met, "rate_pct": round(rate, 2), } summary["criterion_endorsement"] = endorsement log("") log(f"Clinical significance (criterion key '12', not a symptom):") n_cs = int(cohort["clinical_significance"].sum()) log(f" met: {n_cs} ({100.0 * n_cs / n_analytic:.1f}%)") summary["n_clinical_significance_met"] = n_cs log("") log("Total symptom count (sum of c1..c11) distribution:") tsc = cohort["total_symptom_count"].value_counts().sort_index() for val, n in tsc.items(): log(f" {int(val):>2} symptoms {int(n):>6} ({100.0 * n / n_analytic:5.1f}%)") log(f" mean = {cohort['total_symptom_count'].mean():.3f} " f"median = {cohort['total_symptom_count'].median():.1f}") summary["total_symptom_count_distribution"] = { str(int(v)): int(n) for v, n in tsc.items() } summary["total_symptom_count_mean"] = round( float(cohort["total_symptom_count"].mean()), 4 ) who_ordered = pd.Series( pd.Categorical( cohort["bmi_who_category"], categories=[c[0] for c in WHO_CATEGORIES], ordered=True, ), index=cohort.index, ) summary["bmi_who_distribution"] = counts_table( log, who_ordered, "BMI distribution across recomputed WHO categories" ) log(f" recomputed BMI: mean = {cohort['bmi'].mean():.2f} " f"median = {cohort['bmi'].median():.2f} " f"min = {cohort['bmi'].min():.2f} max = {cohort['bmi'].max():.2f}") summary["bmi_mean"] = round(float(cohort["bmi"].mean()), 4) summary["bmi_median"] = round(float(cohort["bmi"].median()), 4) summary["bmi_min"] = round(float(cohort["bmi"].min()), 4) summary["bmi_max"] = round(float(cohort["bmi"].max()), 4) summary["sex_distribution"] = counts_table(log, cohort["sex"], "Sex breakdown") summary["age_range_distribution"] = counts_table( log, cohort["age_range"], "Age-range breakdown" ) summary["country_top15"] = { str(k): int(v) for k, v in cohort["country"].value_counts().head(15).items() } # H2 discordant cells (SAP §4) — counts only, no comparison performed. normal = (cohort["bmi"] >= 18.5) & (cohort["bmi"] < 25.0) obese = cohort["bmi"] >= 30.0 high = cohort["total_symptom_count"] >= 6 low = cohort["total_symptom_count"] <= 1 cells = { "normal_bmi_high_symptom": int((normal & high).sum()), "obese_high_symptom": int((obese & high).sum()), "obese_low_symptom": int((obese & low).sum()), } log("") log("H2 discordant cell sizes (SAP §4) — counts only, no test performed:") log(f" normal BMI (18.5-<25) with >=6 symptoms {cells['normal_bmi_high_symptom']}" f" (SAP §4 quotes 182)") log(f" obese (>=30) with >=6 symptoms {cells['obese_high_symptom']}" f" (SAP §4 quotes 407)") log(f" obese (>=30) with <=1 symptoms {cells['obese_low_symptom']}" f" (SAP §4 quotes 99)") summary["h2_discordant_cells"] = cells # ----------------------------------------------------------------------- # Outputs # ----------------------------------------------------------------------- log.rule("SNAPSHOTS") wide_cols = ( ["respondent_id", "code", "date_create", "bmi", "bmi_who_category", "stored_bmi_value", "stored_bmi_category", "height_cm", "weight_kg", "unit_system", "age_range", "sex", "country", "campaign"] + [f"c{k}" for k in CRITERIA] + ["total_symptom_count", "clinical_significance", "primary_model_eligible", "n_items_recorded", "n_items_with_answer_index"] + [f"q{i}" for i in range(1, N_ITEMS + 1)] ) missing_cols = [c for c in wide_cols if c not in cohort.columns] require(not missing_cols, f"Wide output missing columns: {missing_cols}") wide = cohort[wide_cols].copy() require( int(wide["respondent_id"].duplicated().sum()) == 0, "Wide output has duplicate respondent_id values.", ) require( bool((wide[[f"c{k}" for k in CRITERIA]].sum(axis=1) == wide["total_symptom_count"]).all()), "total_symptom_count does not equal sum(c1..c11) in the wide output.", ) long_records: list[dict[str, Any]] = [] for rec in wide.to_dict("records"): total = rec["total_symptom_count"] for k in CRITERIA: met = rec[f"c{k}"] long_records.append({ "respondent_id": rec["respondent_id"], "criterion_id": k, "criterion_met": met, "criterion_set": criterion_set_of(k), "rest_score": total - met, "total_symptom_count": total, "bmi": rec["bmi"], "bmi_who_category": rec["bmi_who_category"], "age_range": rec["age_range"], "sex": rec["sex"], }) long = pd.DataFrame(long_records) require( len(long) == 11 * len(wide), f"Long output has {len(long)} rows, expected {11 * len(wide)}.", ) require( bool(long["rest_score"].between(0, 10).all()), "rest_score outside the SAP §2.1 range 0-10.", ) wide_ts = SNAPSHOT_DIR / f"analytic_wide_{stamp}.csv" long_ts = SNAPSHOT_DIR / f"analytic_long_{stamp}.csv" wide.to_csv(wide_ts, index=False) long.to_csv(long_ts, index=False) wide.to_csv(SNAPSHOT_DIR / "analytic_wide_latest.csv", index=False) long.to_csv(SNAPSHOT_DIR / "analytic_long_latest.csv", index=False) log(f" wide (immutable) {wide_ts.name} {len(wide)} rows x {len(wide.columns)} cols") log(f" long (immutable) {long_ts.name} {len(long)} rows x {len(long.columns)} cols") log(f" stable copies analytic_wide_latest.csv / analytic_long_latest.csv") summary["files"] = { "wide_timestamped": wide_ts.name, "long_timestamped": long_ts.name, "wide_latest": "analytic_wide_latest.csv", "long_latest": "analytic_long_latest.csv", "log": log.path.name, } summary["n_long_rows"] = len(long) summary["freeze_clause"] = ( "module = 'food_addiction' AND date_create <= '2026-07-25 23:59:59'" ) summary["environment"] = {"python": sys.version.split()[0]} for pkg in ("pandas", "numpy", "scipy", "statsmodels", "PyMySQL", "girth", "patsy"): try: from importlib.metadata import version as _v summary["environment"][pkg] = _v(pkg) except Exception: # pragma: no cover summary["environment"][pkg] = None summary_path = SNAPSHOT_DIR / "cohort_summary.json" summary_path.write_text( json.dumps(summary, indent=2, sort_keys=False, default=str), encoding="utf-8" ) (SNAPSHOT_DIR / f"cohort_summary_{stamp}.json").write_text( json.dumps(summary, indent=2, sort_keys=False, default=str), encoding="utf-8" ) log(f" summary cohort_summary.json (+ timestamped copy)") log.rule("EXTRACTION COMPLETE — nothing inferential was computed.") log.flush() print(f"\nAudit log written to: {log.path}") return 0 if __name__ == "__main__": sys.exit(main())