Statistical analysis for medical research papers. Generates reproducible Python/R code with publication-ready tables and figures. Supports diagnostic accuracy, inter-rater agreement, meta-analysis, survival analysis, survey data, group comparisons, regression, propensity score, and repeated measures.
123k tokens
context cost
the whole folder, loaded on every use
63
files
ships runnable scripts
0
copies elsewhere
how many repositories repackaged it
230
stars on the repo
on the repository, not the skill itself
Install
one command, takes just this skill from the repository
Project data: See CLAUDE.md for data locations under 2_Data/
Read relevant templates before generating analysis code. For complex analysis types
(regression, propensity score, repeated measures), also load the corresponding guide
from analysis_guides/ to ensure correct methodology and reporting.
Workflow
Phase 1: Data Assessment
Read the data file (CSV, Excel, TSV, or other tabular format).
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 preview
Unique value counts for categorical columns
Identify the analysis unit: patient, exam, lesion, image, rater, study, etc.
Phase 2: Analysis Plan
Precondition (observational studies). Before proposing an analysis plan for an observational design (cohort, case-control, cross-sectional, registry, or survey), confirm that a literature-grounded variable operationalization exists — 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, so exposure / outcome / covariate definitions and cutoffs are citation-backed rather than invented ad hoc from the data dictionary (ad-hoc phenotype/cutoff definitions are a common reviewer-rejection trigger for observational work — see the dictionary-first discipline). This is a WARN, not a hard block: proceed on explicit user confirmation, recording that the operationalization artifact was not available. For stricter projects, treat the missing artifact as a hard stop until /define-variables has run. (This mirrors the same precondition already enforced in /write-protocol before drafting Methods.)
Based on the data structure and research question, propose an analysis plan:
Auto-detect analysis type from the table below, or accept user specification.
List specific tests to be performed.
Identify primary and secondary endpoints.
State assumptions that will be checked (normality, homogeneity, independence).
Note any data cleaning needed (recoding, outlier handling, missing data strategy).
Anchor the estimand to the research question. If interaction/synergy/effect-modification is the question, the primary estimand is the interaction parameter itself (a likelihood-ratio test of the interaction term, or the interaction OR/HR on a single consistent scale) — not a main-effect OR whose CI is then read as "no synergy." If the claim is 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.
Screen every categorical/binary predictor for separation — before fitting anything.
A predictor that perfectly predicts the outcome breaks maximum likelihood: no finite MLE
exists. The failure is silent — glm does not error, it returns an odds ratio near 0 (or
enormous), *p* ≈ 0.99, and an AUC that then gets written into a table. This is routine in
diagnostic imaging, because the good signs are the pathognomonic ones (T2-FLAIR mismatch,
the string sign, a halo sign): 100% specificity means an empty cell by construction.
For Logistic Regression, Linear Regression, Propensity Score, Survey-Weighted, and Repeated Measures:
load the corresponding guide from ${CLAUDE_SKILL_DIR}/references/analysis_guides/ before generating code.
For Survey-Weighted analysis, also load survey_weighted.md. For NHIS claims-based studies, load nhis_icd10_mapping.md.
For test selection guidance, load ${CLAUDE_SKILL_DIR}/references/analysis_guides/test_selection.md.
Phase 3: Execute
Generate and run a Python (preferred) or R script following these rules:
Script Structure
Every script MUST start 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
Random seed: Always np.random.seed(42) or set.seed(42).
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)
Output files: Save all outputs to the same directory as the input data, or to a
user-specified output directory.
Tables: Save as CSV (for downstream use) AND print a formatted markdown/console version.
Figures: Save as both PDF (vector) and PNG (300 DPI).
Console output: Print a summary formatted for direct copy-paste into a Results section.
Assumption Checking
Before running parametric tests, always check and report:
Normality: Shapiro-Wilk test (n < 50) or Kolmogorov-Smirnov (n >= 50), plus visual QQ plot
Homogeneity of variance: Levene's test
If assumptions violated: Use non-parametric alternatives and report why
Multiple Comparisons
If running 3+ tests on the same dataset, apply Bonferroni or Benjamini-Hochberg correction.
Always report both uncorrected and corrected p-values.
State the correction method used.
Stratified & Ordinal-Trend Reporting
Strata disjointness gate (before any ordinal trend test). Before running a Cochran-Armitage trend test (or any analysis that treats tiers as an ordered partition), assert 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. Emit the per-stratum N/event table and the reconciliation in the output (this is the analysis-side mirror of /self-reviewcheck_cohort_arithmetic.pyPARTITION_OVERLAP).
Secondary stratum-HR validation checklist. Every secondary stratum hazard/odds ratio must be reported with (a) its reference contrast (which category is the referent), (b) the event count in each stratum, and (c) a sparse-stratum caveat when any stratum has a low event count (a rule of thumb: < 10 events makes the estimate unstable). A bare "HR 1.55 in lean participants" without the referent and the events is uninterpretable.
Proportion CI lower-bound clamp. Clamp every proportion confidence-interval lower bound to max(0, lower); a zero-event Wilson/score interval can emit a negative or absurd tiny-exponent lower bound (e.g., 3.47e-16) that is a display artifact, not a real bound. Report 0 (or 0.0%) instead, and prefer an exact (Clopper-Pearson) interval for zero/near-zero cells.
Output Manifest
After all analyses complete, save a manifest file _analysis_outputs.md in the output directory:
# Analysis Outputs
Generated: {YYYY-MM-DD}
Study type: {detected or user-specified type}
## Tables
- `table1_demographics.csv` -- Baseline characteristics
- `diagnostic_accuracy_table.csv` -- Performance metrics with 95% CIs
## Figures
- `roc_curve.pdf` / `roc_curve.png` -- ROC curves (vector / 300 DPI)
## Data
- `predictions.csv` -- Per-subject model predictions with ground truth
This manifest enables downstream skills (/make-figures, /write-paper) to auto-discover analysis outputs without user intervention.
Phase 3.5: Generated-Code Quality Gate
Before reporting any script as final, lint every emitted .py/.R file for the
reproducibility-hygiene "slop" that AI-generated analysis code recurrently carries:
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
Major findings (fix before reporting the script):
MISSING_SEED — randomness used (sampling, bootstrap, train/test split, rng) with no
Table type guide: references/table-standards/table-types/table1_demographics.md
Continuous variables: mean +/- SD if normal, median (IQR) if skewed
Categorical variables: n (%)
Binary variables: show only one level (e.g., Male n (%), not both Male and Female)
Compare groups: t-test/Mann-Whitney for continuous, chi-square/Fisher for categorical
Report standardized mean differences (SMD) if requested (preferred over P for PS-matched studies)
RCTs: P values in Table 1 are usually unnecessary per CONSORT
gtsummary tbl_summary() with journal theme for R pipeline
Diagnostic Accuracy
Methodology guide: references/analysis_guides/diagnostic_accuracy.md (load before generating code — every metric with a CI on a stated analysis unit; the confidence-weighted trap [unweighted-baseline AUC + monotonic-encoding check, produce-side of probe D9]; paired DeLong vs MRMC for reader-generalising claims; per-stratum admissibility [D10]; one-scale-per-comparison [D11])
ROC curve: include diagonal reference line, AUC in legend
If comparing models: DeLong test for AUC comparison
Youden's index for optimal threshold when applicable
Include calibration assessment (Brier score, calibration plot) for prediction models
NRI/IDI: When comparing two models (e.g., base model vs model + AI score), report:
Category-based NRI (with clinically defined risk categories)
Continuous NRI (note: tends to be inflated — report alongside category-based)
IDI (Integrated Discrimination Improvement)
Bootstrap 95% CIs (1000+ iterations)
These supplement, not replace, DeLong AUC comparison
Table type guide (added value beyond a baseline): references/table-standards/table-types/incremental_value.md (paired ΔAUC + DeLong CI, continuous NRI with event/non-event split, IDI, net benefit at a prespecified threshold, same-patient/calibrated-first discipline). Pairs the decision-curve exemplar make-figuresreferences/exemplar_plots/decision_curve.md.
Reader study (MRMC): references/table-standards/table-types/reader_study.md (per-reader + reader-averaged AUC with an Obuchowski–Rockette/DBM reader+case CI, per-patient vs per-lesion unit, superiority vs non-inferiority margin). Use an MRMC method (not a fixed-reader DeLong CI) for a claim that generalises to readers. Pairs make-figuresreferences/exemplar_plots/mrmc_roc.md.
Inter-rater Agreement
Methodology guide: references/analysis_guides/agreement_reliability.md (load before generating code — the pseudoreplication trap for clustered/repeated measurements + the pseudoreplication-safe per-subject / mixed-effects code, ICC model/type selection, agreement-vs-reliability distinction; pairs with self-review probe O18)
Table type guide: references/table-standards/table-types/agreement.md (ICC with model/type + CI, weighted κ for ordinal, Bland–Altman bias + LoA, reliability-vs-agreement distinction, common errors)
Single-arm pooled proportion: metaprop() with sm = "PLOGIT", method.ci = "CP"
Small-study test branch: do not use Egger's regression for a single-arm proportion meta-analysis — funnel-asymmetry tests assume an effect-size-vs-SE relationship that does not hold for raw proportions. If a small-study assessment is needed, use Peters' test or an arcsine-based variant, and only when k >= 10 (note underpowered otherwise)
Standard output: report tau-squared on the logit scale and a 95% prediction interval (metaprop(..., prediction = TRUE)) in addition to the pooled estimate; the PI conveys where a future study's proportion is expected to fall under the random-effects model
Nested observation units: if the proportion's unit is nested within study (e.g., per-lesion within study, per-image within patient), do not report a naive Wilson/binomial CI that ignores clustering — use a cluster-bootstrap or a GLMM with a random intercept per study so the CI reflects the design
Heterogeneity: I-squared, Q test, tau-squared, and a 95% prediction interval for the random-effects pooled estimate
Forest plot: individual studies + pooled estimate
Funnel plot + Egger's test for publication bias (comparative effect sizes only; note: underpowered k<10)
Threshold effect: Spearman correlation between logit(Se) and logit(FPR)
If significant: interpret single pooled Se/Sp with caution, emphasize SROC curve
Forest plots: Paired (sensitivity + specificity side by side)
Publication bias: Deeks' funnel plot asymmetry test (NOT standard funnel plot)
Standard funnel plots are inappropriate for DTA studies
Note: underpowered for k < 10
Dual approach (comparative + single-arm):
Primary: metabin() for comparative studies (OR/RR)
Secondary: metaprop() with sm = "PLOGIT" for single-arm pooled proportion
Use method = "Inverse", method.tau = "DL", method.random.ci = "HK"
Small studies (k < 10): bivariate model may not converge; consider narrative synthesis
Alternative: If mada unavailable, use metafor::rma.mv() with bivariate structure
Network Meta-Analysis
Guide: Load analysis_guides/network_meta_analysis.md before generating code
For ≥3 interventions via combined direct + indirect evidence (incl. component NMA); pairwise machinery (search/screening/random-effects model) via the Meta-analysis section above
Assess transitivity before pooling: compare effect-modifier distributions across comparisons (box plots / table) and/or network meta-regression — it is a clinical judgment, not a test
Test consistency globally (design-by-treatment) AND locally (node-split / back-calculation); a star network (no closed loops) cannot be checked — state it; investigate the source of any inconsistency (often one trial)
R netmeta (frequentist: netsplit, decomp.design, netheat, netrank P-scores, comparison-adjusted funnel) or Bayesian gemtc / multinma / BUGSnet (node-split, SUCRA, DIC)
Present a network plot (node ∝ sample size, edge ∝ #trials); report global τ²; ranking (SUCRA/P-score) is not a superiority test — report it with the league table, intervals, and certainty
Certainty per estimate via CINeMA / GRADE-NMA (downgrade indirect-only); component NMA assumes additivity (state/check it). Report against PRISMA-NMA; risk of bias via RoB-NMA. Review-side probes: NM1–NM8 in network_meta_analysis.md
Health Economic Evaluation
Guide: Load analysis_guides/health_economic_evaluation.md before generating code
For cost-effectiveness (CEA), cost-utility (CUA, QALY), cost-benefit (CBA), cost-minimisation, or budget-impact analyses; trial-based or decision-model-based (decision tree, Markov/state-transition, discrete-event simulation)
Compute incremental cost ΔC, incremental effect ΔE, and the ICER = ΔC/ΔE; with ≥3 options remove dominated / extended-dominated strategies before sequential ICERs; prefer net benefit (INMB = λΔE − ΔC) for regression/probabilistic summaries
State and justify the perspective, time horizon (lifetime for chronic disease), discount rate (both costs and outcomes), currency + price year; QALYs from a named preference-based instrument + value set
Uncertainty is the analytic core: one-way / tornado for drivers, probabilistic sensitivity analysis (PSA) with justified parameter distributions (beta for probabilities/utilities, gamma/log-normal for costs) → cost-effectiveness plane + CEAC; scenario analyses for structural choices
R heemod / dampack / hesim / BCEA (state-transition + PSA + CEAC + EVPI), flexsurv for survival extrapolation. Report against CHEERS 2022; make the "cost-effective" conclusion conditional on a stated willingness-to-pay threshold. Review-side probes: HE1–HE8 in health_economic_evaluation.md
Survey/Likert
Descriptive: median, IQR, frequency distribution per item
Internal consistency: Cronbach's alpha with item-total correlations
Reverse-coding guard (run before reliability): a negatively-worded scale item must be recoded (min+max) - x before computing the scale total or Cronbach's alpha. An un-recoded reverse item produces a *negative* item-rest correlation and a negative alpha — which is a coding bug, not evidence of a multidimensional construct (do not defend it as such; you lose a review round). likert_summary.py prints the per-item item-rest correlations, flags negative ones as reverse-code suspects, warns loudly on a negative alpha, and accepts --reverse-items E3 ... to apply the recode before scoring. To screen at cleaning time, run /clean-datascripts/check_reverse_coding.py. See the global rule survey-scale-reliability.md.
If comparing groups: Mann-Whitney or Kruskal-Wallis (ordinal data)
Visualization: diverging stacked bar chart
Survival Analysis
Methodology guide: references/analysis_guides/survival.md (load before generating code — competing risks first [naive 1−KM overestimates → produce the Aalen–Johansen/Fine–Gray CIF; cause-specific vs subdistribution for which question, produce-side of probe S3]; PH check → RMST when violated; reverse-KM follow-up + C-index variant [S6]; estimand provenance [S8])
Table type guide: references/table-standards/table-types/survival_results.md (Cox results table: events/person-time, reverse-KM median follow-up, univariable + adjusted HR with CI, PH-assumption footnote, EPV/sparse-stratum and RMST-when-PH-violated rules)
Kaplan-Meier curves with number-at-risk table
Log-rank test for group comparison
Cox proportional hazards: report HR (95% CI)
Events-per-variable (EPV) gate: check events / n_covariates >= 10 before fitting Cox (mirror of the logistic EPV rule). Warn if violated and fall back to a Firth/penalized Cox or profile-likelihood CIs; do not report Wald CIs from a sparse-event model as if stable
Nested observation units (cluster-robust CI): when a subject contributes more than one analysed unit (multiple lesions, both eyes, repeated episodes), pass a subject id so the HR CIs use a robust cluster-sandwich variance (coxph(..., cluster = id) / robust = TRUE in R, cluster_col= in lifelines, e.g. survival_analysis.py --cluster <id>). Treating correlated rows as independent understates the standard errors and narrows the CI artificially
PH violation → do not report a single time-averaged HR. If the Schoenfeld global test is significant (or a covariate's residual trends with time), a single Cox HR averages a changing effect and is misleading. Report a piecewise / time-stratified HR (split follow-up at a clinically sensible cut, or tt() time-transform), or switch to RMST difference at a fixed horizon, and state the violation explicitly
Horizon vs follow-up. Do not read a KM/CIF estimate at a horizon beyond the data: if a reported time point (e.g., a 15-year cumulative incidence) exceeds the reverse-KM median follow-up, either restrict the horizon to where the risk set is non-trivial or report the number-at-risk at that horizon so the reader can judge the extrapolation
Report median survival with 95% CI
Warranty period / quantile estimands (T25 etc.): Time to a fixed cumulative incidence. Use quantile() from the KM/survfit object and always emit the 95% CI (the lower/upper from quantile(km, conf.int=TRUE), or a log-transformed / bootstrap CI) alongside the events/n that define it. A quantile point estimate reported without its CI is incomplete. If the event rate is below the target quantile, report "not reached" and consider Weibull parametric extrapolation (also with an interval)
Interval-Censored Survival
When exact event times are unknown (e.g., health screening cohorts where status changes are detected at periodic visits), standard KM underestimates time-to-event. Use interval-censored methods:
R packages: icenReg (parametric/semi-parametric IC regression), interval (NPMLE/Turnbull), survival (Surv type "interval2")
Turnbull estimator: Non-parametric MLE for interval-censored data — analogous to KM but accounts for the interval between last negative and first positive observation
Parametric IC models: Weibull or log-logistic via icenReg::ic_par(). Report shape/scale parameters and compare AIC across distributions
Mid-point imputation: Simple approximation — event time = midpoint of (last negative, first positive). Acceptable as sensitivity analysis but NOT as primary method
When to use: Serial measurement cohorts (e.g., health screening databases), cancer screening intervals, repeated biomarker assessments
Auto-trigger: if the event date is defined by a periodic visit / scheduled re-examination (the event is detected *at* a visit, not observed exactly), interval-censoring is not optional — make an IC model the primary analysis, or at minimum a mandatory pre-specified sensitivity analysis, and do not present a right-censored Cox coxph() on visit-dated events as if the times were exact
Multistate / transition models: for repeated transitions (e.g., msm), account for subject-level clustering with a subject random effect or a sandwich (robust) variance, and check the time-homogeneity assumption (constant transition intensities) before trusting a single rate
Reporting: State the interval-censored nature of the data explicitly in Methods. Report both standard KM (for comparability with prior literature) and IC estimates (as primary or sensitivity)
Competing Risks
When death or other events preclude the outcome of interest, standard KM overestimates cumulative incidence (treats competing events as censored). Use competing risk methods:
R packages: cmprsk (Fine-Gray), tidycmprsk (tidy interface), survival (cause-specific Cox)
Cumulative incidence function (CIF): cmprsk::cuminc() — replaces 1-KM for each event type. Gray's test for group comparison
Fine-Gray subdistribution hazard: cmprsk::crr() or tidycmprsk::crr() — reports subdistribution HR (sHR) with 95% CI. Interpretable as effect on CIF directly. Check the subdistribution-PH assumption the same way you check it for Cox (a time-interaction term on the subdistribution scale, or inspection of scaled-residual analogues); a constant sHR is an assumption, not a given. Report the cause-specific HR alongside it so the etiologic and prognostic readings are both visible
Cause-specific Cox: Standard Cox censoring competing events — reports cause-specific HR. Better for etiology; Fine-Gray better for prognosis/prediction
When to use: Mortality studies with multiple causes of death, cardiovascular events when non-CV death is frequent, any outcome where competing events are common (>5% of total events)
Reporting: Present CIF plots (NOT 1-KM) when competing risks exist. Report both cause-specific HR and subdistribution HR when the research question is etiologic. State which competing events were defined. When a CIF is quoted at a horizon beyond the median follow-up, report the number-at-risk at that horizon (or restrict the horizon) — a CIF extrapolated past the data is not a stable estimate
Group Comparison
2 independent groups: t-test or Mann-Whitney U
2 paired groups: paired t-test or Wilcoxon signed-rank
3+ independent groups: ANOVA or Kruskal-Wallis, with post-hoc
3+ paired groups: repeated measures ANOVA or Friedman, with post-hoc
Always report: test statistic, degrees of freedom, p-value, effect size
Correlation
Pearson r (if bivariate normal) or Spearman rho (if not)
Report: coefficient, 95% CI, p-value
Scatter plot with regression line and CI band
For multiple variables: correlation matrix heatmap
Logistic Regression
Guide: Load analysis_guides/regression.md before generating code
Run univariable analysis first, then multivariable with clinically selected variables
Required outputs: OR table (univariable + multivariable), C-statistic (95% CI), and calibration (intercept + slope + flexible plot — not Hosmer–Lemeshow, which is deprecated; see the calibration guide)
Prediction-model calibration guide: references/analysis_guides/calibration.md (load before generating code for any model that outputs a risk used for a decision — the apparent slope of exactly 1.00 is the in-sample tell, so produce the bootstrap optimism-corrected slope/intercept; Van Calster's calibration levels; scaled Brier; why Hosmer–Lemeshow is dropped; produce-side of probe S7)
Check VIF < 5, EPV >= 10 (warn if violated)
Nested observation units: when rows are clustered within subjects (multiple lesions/visits per patient), use cluster-robust standard errors (cov_type="cluster", cov_kwds={"groups": id} in statsmodels) or a mixed-effects logistic model — a naive logit CI assumes independent rows and is too narrow
Box-Tidwell test for continuous predictor linearity
Forest plot of adjusted ORs
NRI/IDI if comparing models (incremental value assessment)
Linear Regression
Guide: Load analysis_guides/regression.md before generating code
For KNHANES/NHANES/KCHS and similar complex survey designs
Always declare survey design (strata, cluster/PSU, weight) before analysis
Use correct weight variable (interview vs exam vs nutrition)
R survey package strongly recommended over Python for publication
Sequential model building: Model 1 (age+sex) → Model 2 (full adjustment)
Report weighted odds ratios (wOR) with 95% CI
Cross-national: analyze each country separately, never pool
Subgroup analysis: exclude the stratification variable from covariates
Mediation Analysis
Guide: Load analysis_guides/mediation.md before generating code
Bootstrapped product-of-coefficients (a×b) indirect effect (R mediation / CMAverse / PROCESS); ≥2000 resamples, bias-corrected percentile CI — not the Sobel test
Binary outcome: counterfactual / natural-effects decomposition (CMAverse, regmedint), not the naive OR product
Report total, direct, indirect effects each with a bootstrap CI; proportion mediated only with uncertainty and only when the total effect is well-estimated (unstable / can exceed 100% when total is near-null)
Identification, not the bootstrap, is the issue: mediation needs no unmeasured mediator–outcome confounding (sequential ignorability) → report an E-value for the indirect effect (or ρ-based sensitivity). A cross-sectional design cannot order X→M→Y — frame as association-level (review probe O13)
Report against AGReMA
Interaction & Effect Modification
Choose and state the scale. A public-health / biological synergy claim is an additive-scale statement → report RERI, AP (attributable proportion), or S (synergy index), each with a CI — not only a multiplicative OR/HR product term. A non-significant multiplicative interaction is compatible with a large additive one (and vice versa)
"Joint association" via a combined multi-level exposure (high/high vs low/low) shows joint *categories*, not interaction — add the product term (multiplicative) and/or RERI (additive) to claim interaction
Stratified-only "stronger in A than B" is the difference-in-significance fallacy — report the formal interaction term, not two separate stratum estimates
R interactionR / epiR for RERI/AP/S with CIs; follow Knol & VanderWeele interaction-reporting recommendations. Review-side probe: O14 in observational_confounding.md
Multiple Testing & High-Dimensional Screening
Guide: Load analysis_guides/multiplicity.md before generating code
For agnostic many-exposure scans (ExWAS / EWAS / MWAS / proteome-/nutrient-wide) and any "screen N predictors, report the significant ones" pass
Match the correction to the claim: FWER (Bonferroni / Holm / permutation-based study-wide threshold) for a confirmatory single hit; FDR (Benjamini–Hochberg q-value) for discovery — then frame as hypothesis-generating
Report the correction method and the number of tests m (the denominator), applied to the whole tested set — never shrink m to the winners
Replication is the real safeguard: split-half / second cohort / cross-cycle with directional concordance and a reported replication rate; a single-cohort FDR-significant scan is exploratory
Correlated exposures → raw Bonferroni is over-conservative (permutation or effective-number-of-tests via poolr::meff()); a univariate hit may be a marker for a correlated cause (consider WQS / quantile g-computation / BKMR before causal reading)
Report full results (all effect sizes + p/q), not only winners; complex surveys combine design-based SEs (survey_weighted.md) WITH the correction. Review-side probe: O17 in observational_confounding.md
Mendelian Randomization
Guide: Load analysis_guides/mendelian_randomization.md before generating code
State and evidence the 3 IV assumptions: relevance (F-statistic / R²), independence (ancestry + confounder scan), exclusion restriction (no horizontal pleiotropy — the untestable one)
Pre-specify the full sensitivity suite, not IVW alone: IVW + MR-Egger (intercept = directional pleiotropy) + weighted median + weighted mode + MR-PRESSO; Cochran's Q + leave-one-out; concordance across methods is the robustness claim (R TwoSampleMR / MendelianRandomization)
Address reverse causation (Steiger / bidirectional), sample overlap (report fraction; overlap + weak instruments inflate type-1 error), winner's curse (select instruments in an independent GWAS), and ancestry matching
Drug-target / cis-MR: GLS-IVW for correlated cis variants + colocalization (coloc) + positive control + adverse-effect phenome scan. Non-linear MR: residual/doubly-ranked shapes can be artefactual → require negative/positive controls + extreme-stratum sensitivity
Interpret as a lifelong genetic-proxy effect direction, not a clinical-intervention magnitude; report against STROBE-MR. Review-side probes: MR1–MR8 in mendelian_randomization.md
Polygenic Risk Score (PRS / PGS)
Guide: Load analysis_guides/polygenic_risk_score.md before generating code
For developing/validating/applying a genome-wide polygenic score as a predictor or risk-stratifier (distinct from MR: PRS is prediction, MR is causal inference)
Base (discovery GWAS) and target/validation samples must be independent; tune (P+T / LDpred2 / PRS-CS shrinkage / quantile cut) on a separate tuning set and evaluate out-of-sample (avoid overfitting / winner's curse). Tools: PRSice-2, LDpred2 (bigsnpr), PRS-CS / PRS-CSx, BridgePRS
Ancestry portability is the central issue: report performance separately per target ancestry (within-group differences can rival between-group); prefer ancestry-matched / multi-ancestry discovery + PCs; do not extend a European-derived score to other ancestries without per-ancestry validation
Report OR/HR per SD (CI) + quantile absolute risk; discrimination (C/AUC) and calibration (plot + slope/intercept) in the target population — discrimination ≠ calibration
Incremental value is the clinical crux: report PRS on top of the guideline clinical model (SCORE2/QRISK3/PCE/Tyrer-Cuzick) — ΔC-statistic (CI), NRI/IDI, net benefit — not PRS-alone AUC. A screening claim needs detection-rate-at-fixed-FPR / likelihood ratio, not AUC
Prefer prospective/incident validation (prevalent case–control overstates utility); report against PGS-RS / TRIPOD+AI. Review-side probes: PG1–PG8 in polygenic_risk_score.md
NHIS Claims-Based Studies
Guide: Load analysis_guides/nhis_icd10_mapping.md for disease definition patterns
Claims-based algorithms: N-claim rule, claim+medication, look-back period
Always specify ICD-10 code ranges, claim count requirement, and time windows
Charlson comorbidity index: cite Quan 2005 adaptation
Anchor covariates to most recent data prior to index date
Sensitivity analysis: test stricter/looser disease definitions
Burden of Disease, Decomposition & Forecasting
Guide: Load analysis_guides/burden_decomposition_forecasting.md before generating code
For a burden-of-disease estimate, attributable-risk (PAF / comparative risk assessment), temporal-trend (joinpoint / AAPC), decomposition (Das Gupta; Arriaga life-expectancy), or forecast (BAPC / age-period-cohort)
The value-add-layer playbook — a descriptive rate is rarely publishable alone; bolt on ONE layer: decomposition (*why* the rate changed — aging vs population growth vs epidemiological change), PAF (*how much* is modifiable), joinpoint/AAPC pre-vs-post (did a datable policy/shock bend the trend), forecast (*where* it is going), Arriaga (which ages/causes drove ΔLE)
Uncertainty intervals, not CIs: report a draw-based 95% UI (2.5th–97.5th percentile of 250–500 draws propagated end-to-end); a UI crossing the null means insufficient evidence for direction, not a non-significant test
Single-center cohort adaptation: three layers port onto existing follow-up without new data — trend-break (reslice around a datable guideline/scanner-era change), forecast (project a serial imaging trajectory), and global framing (place the individual-level effect next to the published GBD burden from GHDx — contextualization, not re-estimation, so no GATHER trigger). The ecological "UI-replaces-confounding-control" shortcut does not port: individual-level data still needs the DAG / E-value / negative-control toolkit
Report the estimate against GATHER (/check-reporting); keep burden/attribution/decomposition/forecast descriptive or associational unless a causal design (natural experiment, MR — analysis_guides/mendelian_randomization.md) is in place
Repeated Measures
Guide: Load analysis_guides/repeated_measures.md before generating code
Applies to any multivariable adjustment (logistic / linear / Cox / propensity-score / survey-weighted). Two coupled failure modes around a dose/duration variable anchored to a categorical exposure (pack-years under smoking status, grams/week under alcohol use, cessation-duration under former-smoker):
Structural-zero guard (do not impute): a never-smoker's pack_years is a *structural zero*, not missing-at-random — the value is known to be 0 by definition of the category. Feeding it to MICE/MNAR imputation as if it were missing fabricates a non-zero dose for unexposed subjects and corrupts the exposure contrast. Before imputing any dose/duration column, set the implied zero explicitly (IF status == 'never' THEN dose = 0) and impute only the genuinely-missing residual among the exposed. /clean-data flags categorical-implied-zero contradictions (a never row with a NULL dose) and ships scripts/check_structural_zero.py.
Complete-case collapse warning (use status, not dose, for adjustment): when a dose/duration variable enters a *complete-case* multivariable model, the unexposed stratum — which carries structural zeros often stored as NULL — is dropped wholesale, collapsing n (commonly 40–60%) and distorting subgroup estimates (a small stratum can shrink to a handful of subjects). For confounder adjustment use the categorical status variable (never/former/current); reserve the continuous dose for an exposed-only (e.g., ever-smoker-restricted) *secondary* analysis. Always report n before and after model fitting and confirm the denominator did not silently collapse.
Covariate Selection: Over-adjustment in a Cross-Sectional Outcome Model
Applies to any cross-sectional / single-visit outcome regression (the exposure and outcome are measured at one time point, so temporal order is not observed). The selection rule is causal, not statistical:
Do not adjust for a consequence or mediator of the outcome. A covariate that the outcome physiologically *drives* sits on or after the causal path; adjusting for it is over-adjustment / collider bias and removes part of the effect under study. The signature case is a renal-function outcome: with eGFR as the outcome, serum uric acid is renally excreted (a lower eGFR mechanically raises urate), so uric acid is an outcome-consequence, not a confounder; blood pressure and HbA1c are often similarly downstream. Classify each candidate covariate against a DAG as confounder / mediator / outcome-consequence / collider, and keep only confounders in the primary model.
"It differs in Table 1" is not a confounder-selection criterion. Baseline imbalance by exposure justifies *considering* a variable, but a mediator or outcome-consequence stays out regardless of how imbalanced it is. A kitchen-sink "adjust for everything that differs" model is over-adjusted by construction.
Report the suspect-covariate sensitivity + VIF. Make a parsimonious, history-/design-based model the primary one; report the fuller model as a sensitivity analysis that drops the suspect covariate (or adds it, if you start parsimonious), and show whether the headline estimate moves. Always print VIF (collinearity between an outcome-consequence and the outcome's other correlates is common) and the n actually fitted. If dropping the covariate materially changes the estimate, propagate to the abstract and conclusion.
Compare adjusted-vs-unadjusted on the SAME frame (extended-adjustment missingness trap). When an extended-adjustment model adds covariates that carry missingness, the analytic n shrinks (e.g. 84 → 49 events). Comparing that adjusted estimate to the full-frame unadjusted/base estimate confounds *adjustment* with *case-concentrated missingness* — it can look as if "adjustment inflated the estimate" when the drift is who-was-dropped. The fair anchor is the unadjusted estimate refit on the reduced complete-case frame (the same rows the adjusted model used): report unadjusted-and-adjusted on the reduced frame alongside the full-frame estimate, and never describe "adjustment changed the estimate" from a comparison across different frames. (Equivalently, use multiple imputation so all models share one frame.)
Language
Code and output: English
Communication with user: Match user's preferred language
Medical terms: English only
What This Skill Does NOT Do
Does not fabricate or simulate data to fill gaps
Does not choose analysis endpoints -- the user decides the research question
Does not interpret clinical significance -- only statistical results
Does not replace biostatistician review for complex designs (e.g., adaptive trials)
Anti-Hallucination
Never fabricate variable names, dataset column names, or variable codings. If a variable mapping is uncertain, output [VERIFY: variable_name] and ask the user to confirm against the data dictionary.
Never fabricate statistical results — no invented p-values, effect sizes, confidence intervals, or sample sizes. All numbers must come from executed code output.
Never generate references from memory. Use /search-lit for all citations.
If a function, package, or API does not exist or you are unsure, say so explicitly rather than guessing.
How to use it
Copy the folder
Take aperivue/analyze-stats from the repository into ~/.claude/skills for personal
use, or into .claude/skills inside a project.
Check the name does not clash
The agent identifies a skill by the name field in its header. Two skills with the
same name cannot sit side by side — one of them will be ignored.