Analyze Stats

Statistical Analysis Skill

Generate reproducible code (Python preferred, R when necessary) for medical research analyses, run it, and produce publication-ready tables and figures from its output.

Data Privacy Check

Before reading any data file, check whether it might contain Protected Health Information (PHI):

1. If *_deidentified.* files exist in the working directory, use those.
2. If only raw CSV/Excel files exist, ask the user whether the data contain patient identifiers (names, national ID / RRN, contact details); if so, have them de-identify it first with /deidentify.
3. Proceed once the user confirms the data are de-identified or contain no PHI.
4. NEVER display raw PHI values (names, phone numbers, RRN) in your output, because console output is copied into transcripts and drafts. If you encounter them, warn the user and suggest /deidentify.

Reference Files

• ${CLAUDE_SKILL_DIR}/references/templates/ — reusable analysis scripts; read the matching one before generating code.
• ${CLAUDE_SKILL_DIR}/references/analysis_guides/ — methodology guides; the Phase 2 table says which to read.
• ${CLAUDE_SKILL_DIR}/references/table-standards/ — table-standards.md (universal and AMA rules, footnote order, gtsummary pipeline), journal-profiles/ (YAML per journal: radiology, jama, nejm, lancet, european_radiology, ajr), table-types/ (one template per table type), tool-comparison.md (R/Python table tools).
• ${CLAUDE_SKILL_DIR}/references/style/figure_style.mplstyle — figure style.

Workflow

Phase 1: Data Assessment

1. Read the data file (CSV, Excel, TSV, or other tabular format).
2. Report to the user: shape (rows x columns); column names and inferred types (continuous, categorical, ordinal, binary, datetime); missing values per column (count and percentage); first 5 rows; unique value counts for categorical columns.
3. Identify the analysis unit: patient, exam, lesion, image, rater, study, etc.
4. If a variable name or coding is uncertain, write [VERIFY: variable_name] and ask the user to confirm it against the data dictionary; never guess a column name or coding.

Phase 2: Analysis Plan

Observational designs (cohort, case-control, cross-sectional, registry, survey): before planning, look for a variable_operationalization.md from /define-variables (or an equivalent codebook-backed definition table). If none exists, warn the user and recommend running /define-variables first, because exposure/outcome/covariate definitions and cutoffs invented ad hoc from the data dictionary are a common reason reviewers reject observational work. This is a warning, not a block: proceed on explicit user confirmation and record that the artifact was not available.

Propose an analysis plan; the user decides the research question and endpoints. Mark any clinical definition, cutoff, diagnostic criterion or guideline claim you could not confirm [VERIFY].

1. Detect the analysis type from the table in Analysis-Specific Guidelines, or accept the user's specification.

2. List the specific tests to be performed.

3. Identify primary and secondary endpoints.

4. State the assumptions that will be checked (normality, homogeneity, independence).

5. Note data cleaning needed (recoding, outlier handling, missing-data strategy).

6. Anchor the estimand to the research question.

   • Interaction, synergy or effect modification: the primary estimand is the interaction parameter itself (a likelihood-ratio test of the interaction term, or the interaction OR/HR on one stated scale) — not a main-effect OR whose CI is then read as "no synergy". A public-health or biological synergy claim is additive-scale: report RERI, AP or S, each with a CI (R interactionR / epiR; Knol & VanderWeele reporting), not only a multiplicative product term, because a non-significant multiplicative interaction is compatible with a large additive one (and vice versa). Joint categories (high/high vs low/low) show joint association, not interaction; "stronger in A than in B" from separate stratum estimates is the difference-in-significance fallacy — report the formal interaction term.
   • Equivalence or non-inferiority: declare the margin up front (a TOST procedure, or the CI compared against a pre-stated MCID); a non-significant difference is not equivalence without a margin.

7. Screen every categorical/binary predictor for separation before fitting anything. A predictor that perfectly predicts the outcome has no finite MLE, and the failure is silent: glm returns an odds ratio near 0 (or enormous), p ≈ 0.99, and an AUC that ends up in a table. Pathognomonic imaging signs (T2-FLAIR mismatch, the string sign, a halo sign) cause it routinely: 100% specificity means an empty cell by construction.

   python3 "${CLAUDE_SKILL_DIR}/scripts/check_separation.py" \
     --data cohort.csv --outcome idh_mutant --auto --strict

   COMPLETE_SEPARATION (an empty cell) and QUASI_SEPARATION (a cell below the sparsity floor) both halt the plan. With two or more predictors the script also tests them jointly (a linear combination can separate the outcome when no single predictor does); that test needs scipy, and without it the report says the joint check was not run. The remedy is a design decision, made in the plan: Firth's penalised likelihood keeps one model, while a two-stage rule — classify the sign-positive cases directly, model only the sign-negative remainder — is usually the clinically meaningful choice for a pathognomonic sign, because a sign-positive patient is already diagnosed.

Present the plan and wait for user approval before executing.

Phase 3: Execute

Generate and run a Python (preferred) or R script. Every number you report MUST come from its executed output, never estimated or hand-typed, because a number no script reproduces cannot be verified.

Script Structure

Start every script with a reproducibility header:

"""
Analysis: {description}
Date: {YYYY-MM-DD}
Random seed: 42
Python: {version}
Key packages: {package==version, ...}
"""
import numpy as np
import pandas as pd
np.random.seed(42)

Execution Rules

1. Random seed: always np.random.seed(42) or set.seed(42).
2. Figure style: always load the matplotlib style file:

   import matplotlib.pyplot as plt
   style_path = os.path.join(os.environ.get('CLAUDE_SKILL_DIR', '.'), 'references/style/figure_style.mplstyle')
   if os.path.exists(style_path):
       plt.style.use(style_path)

3. Output files: save next to the input data, or in a user-specified output directory.
4. Tables as CSV plus a printed markdown version (see Output Conventions); figures as PDF (vector) and PNG (300 DPI).
5. Console output: a summary with every number that would appear in the text, formatted for direct copy-paste into the Results section.

Assumption Checking

Choose the summary and test from the design and the distribution, not a preliminary test:

• Shape: QQ plot / histogram and skewness; for Table 1, |skewness| > 1 → median (IQR) with a rank test, else mean (SD) (table-types/table1_demographics.md). Never gate on a Shapiro-Wilk or KS P value: it flags trivial departures at large n, misses real ones at small n, and test-then-test inflates the type I error (Rochon et al. 2012).
• Unequal variances: default to Welch's t / Welch ANOVA, not a Levene gate; Mann-Whitney tests a different hypothesis and does not fix them. State the choice and reason in Methods.

Stratified & Ordinal-Trend Reporting

• Strata disjointness gate (before any ordinal trend test). Before a Cochran-Armitage trend test (or any analysis that treats tiers as an ordered partition), assert that the strata are mutually exclusive and exhaustive: sum(n per stratum) == unique N and sum(events per stratum) == total events. A trend test on overlapping or non-exhaustive strata is invalid. Print the per-stratum N/event table and the reconciliation (the analysis-side mirror of /self-review check_cohort_arithmetic.py PARTITION_OVERLAP).
• Secondary stratum HR/OR: report each with (a) its reference contrast (the referent category), (b) the event count in each stratum, and (c) a sparse-stratum caveat when any stratum has few events (< 10 is unstable). "HR 1.55 in lean participants" without the referent and the events is uninterpretable.
• Proportion CI lower-bound clamp: clamp every lower bound to max(0, lower). A zero-event Wilson/score interval can print a negative or tiny-exponent bound (e.g. 3.47e-16) that is a display artifact; report 0 (or 0.0%), and prefer an exact (Clopper-Pearson) interval for zero/near-zero cells.

Output Manifest

After all analyses complete, save _analysis_outputs.md in the output directory, in the format given in references/analysis_run_workflow.md, so /make-figures and /write-paper can find the outputs without asking.

For prespecified binary predictions on independent units, use the bundled scripts/run_analysis.py run workflow in that file: it embeds data/configuration/code/output hashes, exact counts, metric-specific denominators and the reproduction command in the same manifest (audit checks recorded versions without rewriting them; compare separates declared context and numeric equality from byte drift). It does not select thresholds or establish study validity, privacy clearance or reuse rights. Synthetic example: python3 ${CLAUDE_SKILL_DIR}/scripts/demo_analysis_run.py --out demo-project.

Phase 3.5: Generated-Code Quality Gate

Before reporting any script as final, lint every emitted .py/.R file:

python3 ${CLAUDE_SKILL_DIR}/scripts/check_generated_code.py {script.py} --strict
# or scan a whole output directory:
python3 ${CLAUDE_SKILL_DIR}/scripts/check_generated_code.py --code-dir {analysis_dir} --strict

Fix every Major before reporting the script: MISSING_SEED (randomness with no seed), HARDCODED_DATA_LITERAL (hand-typed, table-shaped data instead of read_csv()/read.csv() + subset), HARDCODED_ABS_PATH (non-portable and a PII risk), and INPLACE_SOURCE_OVERWRITE (writing to the path read as input — never modify raw data; write derived outputs to a new path). Fix the flags DEBUG_LEFTOVER and UNUSED_IMPORT when tidying. Minor, Python only: API_DEFAULT_STUDENT_T (scipy's ttest_ind without equal_var is Student's t; state Welch or the reason for Student) and API_DEFAULT_PENALIZED_OR (LogisticRegression with none of penalty, C, l1_ratio, in a file that exponentiates a coef_). R's t.test() is Welch by default.

Phase 4: Report

After execution, generate manuscript-ready text from the script's output:

1. Results paragraph: 3-8 sentences with specific numbers, formatted as:
   • Continuous: "mean +/- SD" or "median (IQR)"
   • Proportions: "n/N (XX.X%)"
   • Test results: "statistic = X.XX, p = 0.XXX"
   • Effect sizes: "Cohen's d = X.XX (95% CI: X.XX-X.XX)"
   • AUC: "AUC = 0.XXX (95% CI: 0.XXX-0.XXX)"
2. Table/figure captions: draft captions referencing table/figure numbers.
3. Methods snippet: 2-3 sentences describing the statistical methods, for the Methods section. Cite only references whose DOI/PMID /search-lit confirmed; mark any other [UNVERIFIED - NEEDS MANUAL CHECK].

Report statistical results; leave judgments of clinical significance to the user. For designs this skill does not cover (adaptive or Bayesian trials, complex multilevel or causal-mediation models), say so and recommend biostatistician review before the results are reported.

Statistical Reporting Rules (Always Enforced)

1. Exact p-values: report exact values (e.g., p = 0.034), not inequalities; below 0.001, report p < 0.001.
2. Confidence intervals: always report 95% CIs for primary endpoints.
3. Effect sizes: report one alongside every p-value (Cohen's d, eta-squared, odds ratio, risk ratio, etc., as appropriate).
4. Multiple comparisons: for 3+ tests on the same dataset, apply Bonferroni or Benjamini-Hochberg correction, state the method, and report both uncorrected and corrected p-values.
5. Sample size: state n for each group/analysis.
6. Missing data: report how many cases were excluded and why.
7. Decimal places: p-values to 3 decimals, proportions to 1 decimal, means/SDs to the precision of the measurement.
8. Design/power statistics are code outputs, never hand-computed. Any minimum detectable effect (MDE), a-priori power, or required sample size that will appear in the manuscript MUST be printed by the committed script with its method and inputs (n per arm, alpha, power, allocation ratio, one/two-sided), not computed in a side tool (G*Power, an online calculator) and pasted in, because a manuscript value no script reproduces cannot be checked. Use one method family throughout (e.g. the exact noncentral-t via statsmodels TTestIndPower or scipy's nct); do not mix a normal approximation for some values with exact-t for others. Do not report post-hoc (observed) power: it is a function of the P value and adds nothing (Hoenig & Heisey 2001); report the CI of the effect instead.
9. Estimand & CI output contract. Every primary point estimate — including quantile estimands (T25, median time-to-event), pooled proportions and subdistribution HRs, not just ORs/HRs/AUCs — MUST be emitted with its 95% CI, because /self-review treats a primary metric without one as a defect. In the output CSV carry the interval as explicit columns (estimate, ci_lower, ci_upper) or as one text column in est (lo–hi) form; never a point estimate with no interval beside it. Round ORs/HRs/sHRs to 2 decimals and AUC/C-statistic to 3.