Analyze Stats

Statistical Analysis Skill

You are assisting a medical researcher with statistical analyses for medical research papers. Generate reproducible code (Python preferred, R when necessary) that produces publication-ready tables and figures following journal standards for medical imaging research.

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 preferentially. 2. If only raw CSV/Excel files exist (no *_deidentified.* counterpart), warn the user (ask in the user's preferred language):

"Does this data contain patient identifiers (names, national ID / RRN, contact details, etc.)?

If so, please de-identify it first with the /deidentify skill."

3. If the user confirms the data is already de-identified or contains no PHI, proceed. 4. NEVER display raw PHI values (names, phone numbers, RRN) in your output. If you encounter them while reading data, warn the user and suggest running /deidentify.

Reference Files

• Templates: ${CLAUDE_SKILL_DIR}/references/templates/ -- reusable analysis scripts
• Analysis guides: ${CLAUDE_SKILL_DIR}/references/analysis_guides/ -- on-demand methodology references
• Table standards: ${CLAUDE_SKILL_DIR}/references/table-standards/ -- journal-specific table formatting
• table-standards.md -- universal rules, AMA rules, footnote system, mistakes checklist
• journal-profiles/ -- YAML profiles per journal (radiology, jama, nejm, lancet, eur_rad, ajr)
• table-types/ -- templates per table type (Table 1, diagnostic accuracy, regression, survival/Cox, agreement/reliability, meta-analysis, model comparison, incremental value, reader study (MRMC))
• tool-comparison.md -- R/Python tool comparison and recommended pipelines
• Figure style: ${CLAUDE_SKILL_DIR}/references/style/figure_style.mplstyle
• 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

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 preview
• Unique value counts for categorical columns

3. Identify the analysis unit: patient, exam, lesion, image, rater, study, etc.

Phase 2: Analysis Plan

Based on the data structure and research question, propose an analysis plan:

1. Auto-detect analysis type from the table below, or accept user specification. 2. List specific tests to be performed. 3. Identify primary and secondary endpoints. 4. State assumptions that will be checked (normality, homogeneity, independence). 5. Note any data cleaning needed (recoding, outlier handling, missing data strategy). 6. 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.

Present the plan and wait for user approval before executing.

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

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 all outputs to the same directory as the input data, or to a user-specified output directory. 4. Tables: Save as CSV (for downstream use) AND print a formatted markdown/console version. 5. Figures: Save as both PDF (vector) and PNG (300 DPI). 6. 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-review check_cohort_arithmetic.py PARTITION_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

np.random.seed / set.seed / random_state= / default_rng. Non-reproducible.

• HARDCODED_DATA_LITERAL — a hand-typed, table-shaped numeric literal instead of

read_csv()/read.csv() + subset. This is the data-integrity rule "never hand-type CSV data into scripts."

• HARDCODED_ABS_PATH — an absolute path literal (/Users/, /home/, C:\, ~/Documents).

Non-portable and a PII risk.

• INPLACE_SOURCE_OVERWRITE — writing to the same path read as input; this overwrites raw

data. Write derived outputs to a new path ("never modify raw data").

Flags (fix when tidying): DEBUG_LEFTOVER (a breakpoint() / browser() / debug print / TODO marker left in) and UNUSED_IMPORT (a dead Python dependency).

The gate is conservative on the Major checks — it fires HARDCODED_DATA_LITERAL only on genuinely table-shaped literals and MISSING_SEED only on a real randomness call — so it stays quiet on legitimate analysis code. It is the analysis-side mirror of the data-integrity and reproducibility checks /self-review is built to catch downstream.

Phase 4: Report

After execution, generate manuscript-ready text:

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 used, suitable for the Methods section.

Statistical Reporting Rules (Always Enforced)

These rules apply to ALL analyses without exception:

1. Exact p-values: Report exact values (e.g., p = 0.034), not inequalities. Exception: report as p < 0.001 when the value is below 0.001. 2. Confidence intervals: Always report 95% CIs for primary endpoints. 3. Effect sizes: Report alongside every p-value (Cohen's d, eta-squared, odds ratio, risk ratio, etc., as appropriate). 4. Parametric vs non-parametric: Choose based on assumption checks, not convenience. Report the assumption test results. 5. Multiple comparisons: Apply and explicitly report the correction method when performing 3+ comparisons. 6. Sample size reporting: Always state n for each group/analysis. 7. Missing data: Report how many cases were excluded and why. 8. Decimal places: p-values to 3 decimals, proportions to 1 decimal, means/SDs to appropriate precision for the measurement. 9. Design/power statistics are code outputs, never hand-computed. Any minimum detectable effect (MDE), a-priori or post-hoc power, or required sample size that will appear in the manuscript MUST be emitted by this committed script — printed with its method and inputs (n per arm, alpha, power, allocation ratio, one/two-sided) — not computed in a side tool (GPower, an online calculator) and pasted in. Use one method family consistently (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. A value that exists only in the manuscript with no script that reproduces it is the failure mode `/self-review` Phase 2.5a-2 is built to catch. 10. 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 together with its 95% CI. In the output CSV, carry the interval as explicit columns (estimate, ci_lower, ci_upper) or as a single text column in est (lo–hi) form; never emit a point estimate with no interval in an adjacent column. Round ORs/HRs/sHRs to 2 decimals and AUC/C-statistic to 3. This is the output side of the /self-review §C assertion that "all primary metrics have 95% CIs."

Effect-Size Real-World Translation

Whenever a primary result is a correlation, a standardized coefficient, a regression slope, an OR/HR/RR, or a Cohen's d, also report it as a plain-language unit shift a non-statistician can act on. The coefficient answers "is there an association"; the translation answers "how much, in units I use". This complements rule 3 above (report effect sizes) — it does not replace it.

When to apply

• Any continuous-exposure to continuous-outcome association reported as Spearman's rho, Pearson's r,

or a standardized slope.

• Any OR/HR/RR where the audience needs an absolute-risk feel.
• Reader / expert-elicitation studies, clinical-utility framing, abstracts, and figure captions.

Procedure 1. Pick an anchored contrast on the exposure, not a 1-unit step. Default: 25th to 75th percentile (IQR). State both endpoints in native units. 2. Translate to the outcome scale.

• For a rank/standardized association (Spearman's rho or a per-SD slope) under an approximately

monotonic-linear assumption: delta_outcome ~= ((x_p75 - x_p25) / SD_x) * |rho| * SD_outcome. Report as: "going from {x_p25} to {x_p75} {units} is associated with about {delta_outcome} {outcome units} on average."

• For a regression slope b: delta_outcome = b * (x_p75 - x_p25) (cleaner; no monotonicity caveat).
• State the assumption explicitly; the IQR translation is a more defensible verbal guide than an

SD-scaled one. 3. For OR/HR/RR, accompany the relative measure with an absolute one at a stated baseline risk: the absolute risk difference, and NNT = 1 / ARR (or NNH = 1 / ARI). Always state the baseline risk used. 4. Bound the claim: report the contrast, the assumption, and a CI on the coefficient; do not imply causation from a crude or unadjusted estimate.

Worked example (synthetic) rho = 0.39 between a fasting marker (IQR 0.6 to 3.5 units, SD 3.05) and an index (SD 2.13): ((3.5 - 0.6) / 3.05) * 0.39 * 2.13 ~= 0.8 -> "Going from the 25th to the 75th percentile of the marker is associated with about 0.8 index units higher on average (monotonic-linear approximation; crude, unadjusted)."

Output contract: add a "Real-world translation" line beneath each primary effect size in the results table or its footnote. For OR/HR/RR primary outcomes, add an NNT/NNH line with the baseline risk stated.

Error Handling

• If a script fails to execute, report the error in one line, diagnose the likely cause

(missing package, data format mismatch, wrong column name), and present a fix.

• Do NOT retry the same script more than once without modifying it or asking the user.
• If an R package is unavailable, suggest install.packages() and wait for user confirmation.
• For prediction models: always include calibration assessment (Brier score, calibration plot,

or calibration slope/intercept) alongside discrimination metrics. AUC alone is insufficient.

Output Conventions

Tables

Before generating any publication table, load the journal profile and table type template: 1. Load ${CLAUDE_SKILL_DIR}/references/table-standards/journal-profiles/{journal}.yaml if a target journal is known 2. Load ${CLAUDE_SKILL_DIR}/references/table-standards/table-types/{type}.md for the relevant table type 3. If no journal specified, default to AMA style (Radiology profile)

Output formats (always generate all three):

• CSV file (for downstream use and archival)
• Console markdown rendering (for user review)
• R gtsummary code (for publication-quality Word/LaTeX export)

Universal rules (enforced regardless of journal):

• No vertical lines — horizontal rules only (top, below header, bottom)
• Binary variables: show only one level (e.g., Male only, not Male + Female)
• Units in column headers, not repeated in cells
• Consistent decimal places within each column
• All abbreviations defined in footnotes, self-contained per table
• Exact P values always (never "NS" or "significant")
• Name the statistical test in footnote or general note
• Variability measure always stated: mean (SD) or median (IQR)

Journal-specific parameters (from loaded YAML profile):

• Footnote markers: letters (AMA) vs symbols (NEJM/Lancet)
• P value format: case, leading zero, italic
• CI separator: comma (Radiology) vs "to" (JAMA/NEJM/Lancet)
• Title format: period (AMA) vs colon (Lancet)
• Abbreviation order: appearance (Radiology) vs alphabetical (JAMA)

Footnote placement order (universal): 1. General note (no marker) — e.g., "Data are mean (SD) unless noted" 2. Abbreviations — in order per journal convention 3. Specific notes (superscript markers) — per-cell explanations 4. Probability notes — significance thresholds (if applicable)

gtsummary pipeline (recommended for R table generation):

theme_gtsummary_journal("{journal}")  # "jama", "lancet", "nejm"
theme_gtsummary_compact()
# ... build table ...
tbl %>% as_flex_table() %>% flextable::save_as_docx(path = "table.docx")

Validation checklist (run before finalizing any table):

• [ ] Binary variables show only one level
• [ ] Units in headers, not cells
• [ ] Consistent decimal places per column
• [ ] Statistical test named (footnote or general note)
• [ ] Effect sizes per clinically meaningful unit (per 10 years, not per 1 year)
• [ ] Reference category stated for categorical predictors
• [ ] No "NS" — exact P values only
• [ ] Abbreviations defined in footnotes

Figures

• Format: PDF (vector, for journal) + PNG (300 DPI, for review)
• Style: Use figure_style.mplstyle for consistent appearance
• Font: Arial, 8-10pt
• Colors: Colorblind-safe palette
• Size: 3.5 inches (single column) or 7.0 inches (double column) width
• Always include axis labels with units

Console Output

• Formatted for direct copy-paste into the Results section of a manuscript
• Include all numbers that would appear in the text
• Use the reporting format conventions above

Analysis-Specific Guidelines

Table 1 (Demographics)

• Template: references/templates/table1_demographics.py
• 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

• Template: references/templates/diagnostic_accuracy.py
• Always report: sensitivity, specificity, PPV, NPV, accuracy, AUC
• CIs: Wilson score for proportions, DeLong for AUC
• 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-figures references/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-figures references/exemplar_plots/mrmc_roc.md.

Inter-rater Agreement

• 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)
• Template: references/templates/agreement_analysis.py
• 2 raters + categorical: Cohen's kappa
• 2+ raters + categorical: Fleiss' kappa (or Krippendorff's alpha)
• Continuous: ICC (specify model: one-way, two-way random/mixed; type: single/average)
• Always report interpretation labels (Landis & Koch or Cicchetti)
• Bland-Altman plot for continuous paired measurements
• Bootstrap CIs (1000 iterations, seed=42)

Meta-analysis

• Prefer R (meta/metafor packages) for meta-analysis
• Comparative: metabin() for binary outcomes (OR/RR), metagen() for continuous
• Use method = "Inverse", method.tau = "DL", method.random.ci = "HK"
• Avoid deprecated args: comb.fixed → common, hakn → method.random.ci
• 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)
• Sensitivity analysis: leave-one-out (metainf())
• Subgroup: update(res, subgroup = variable)

DTA Meta-Analysis

• Template: references/templates/dta_meta_analysis.R
• Prefer R (mada, meta, metafor packages) for DTA meta-analysis
• Bivariate model (Reitsma): mada::reitsma() — recommended over separate pooling of Se/Sp
• Accounts for correlation between sensitivity and specificity
• Produces SROC curve with confidence + prediction regions
• Key outputs: Pooled Se/Sp (95% CI), positive/negative LR, DOR, SROC AUC
• 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

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-data scripts/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

• 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
• Check proportional hazards assumption (Schoenfeld residuals)
• 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
• Template: references/templates/regression.py (set regression_type = "logistic")
• Run univariable analysis first, then multivariable with clinically selected variables
• Required outputs: OR table (univariable + multivariable), C-statistic (95% CI), Hosmer-Lemeshow
• 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
• Template: references/templates/regression.py (set regression_type = "linear")
• Required outputs: coefficient table (β, 95% CI, P), R²/adjusted R², VIF
• Always generate 4-panel diagnostic plot (residuals vs fitted, Q-Q, scale-location, leverage)
• Check assumptions: normality of residuals, homoscedasticity, multicollinearity
• Report both unstandardized β (primary) and standardized β (for effect size comparison)

Propensity Score

• Guide: Load analysis_guides/propensity_score.md before generating code
• Template: references/templates/propensity_score.py
• Step 1: PS estimation (logistic regression)
• Step 2: Apply method (matching with caliper = 0.2 × SD logit PS, IPTW/SIPTW with stabilized weights, or overlap weighting)
• Step 3: Balance assessment — SMD < 0.10 for all covariates, Love plot mandatory
• Step 4: Weighted/matched outcome analysis with robust SE
• Step 5: Sensitivity analysis (E-value for unmeasured confounding)
• Always state the estimand (ATE/ATT/ATO) explicitly
• Recommend overlap weighting as default (no extreme weight issues)
• SIPTW: Stabilized IPTW variant used in emulated target trial frameworks; report effective sample size

Survey-Weighted Analysis

• Guide: Load analysis_guides/survey_weighted.md before generating code
• Template: references/templates/survey_weighted_analysis.py
• 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

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

Repeated Measures

• Guide: Load analysis_guides/repeated_measures.md before generating code
• Template: references/templates/repeated_measures.py
• Default method: LMM (handles missing data, no sphericity assumption)
• RM ANOVA only if: no missing data AND few time points AND sphericity met
• GEE for: population-averaged effects or non-normal outcomes
• Always convert wide → long format first
• Time × Group interaction is the key result — always report and interpret
• Generate spaghetti plot (individual trajectories) + group mean trajectory plot
• For LMM: report random effects structure, covariance structure (CS/AR1/UN), AIC/BIC
• For RM ANOVA: report Mauchly's test, epsilon, correction method (Greenhouse-Geisser)
• If missing > 5%: load analysis_guides/missing_data.md and apply MICE before analysis

Covariate Pitfalls: Structural Zeros & Dose/Duration Variables

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.

Mediation Analysis Guide

Decomposing an exposure–outcome association into a path through a mediator (indirect effect) and a path not through it (direct effect). The estimate is easy; the identification is hard — that is where mediation analyses fail review.

---

When to Use

• A pre-specified hypothesis that exposure X affects outcome Y partly via mediator M
• M is measured, plausibly on the causal path, and (ideally) measured after X and before Y
• NOT for: a single-timepoint design used to assert a causal X→M→Y chain (see the design caveat

below — that is the most common reason mediation papers are rejected); not a substitute for a longitudinal or experimental test of the mechanism

---

Identification comes first (this, not the bootstrap, is the issue)

A bootstrapped indirect-effect CI quantifies the sampling variability of the a×b product. It says nothing about whether the indirect effect is identified. Causal mediation requires, beyond no exposure–outcome and no exposure–mediator confounding:

• No unmeasured mediator–outcome confounding (sequential ignorability). This must hold *even in

a trial that randomizes X, because M is not randomized. In observational data it essentially never holds unmodelled — so a sensitivity analysis is mandatory*, not optional.

• No mediator–outcome confounder affected by the exposure (an exposure-induced M–Y confounder

breaks the standard product/difference method; use a method that allows it, e.g. an interventional or G-computation estimator).

• Correct temporal order. X precedes M precedes Y. A cross-sectional design measures all

three together and cannot establish this — report the indirect effect as consistent with mediation, state plainly that the design cannot order X/M/Y, and reserve the causal claim for a two-wave / longitudinal design (review probe O13 in observational_confounding.md).

Method

• Estimator: bootstrapped product-of-coefficients (a×b) is standard for continuous M and Y

(R mediation, CMAverse, or Hayes PROCESS; Python pingouin/statsmodels). For a binary outcome use the counterfactual / natural-effects decomposition (e.g. CMAverse, regmedint), not the naive product on the odds-ratio scale (non-collapsibility distorts it).

• Bootstrap the indirect effect (≥ 2000, ideally 5000 resamples) and report the **bias-corrected

percentile CI** — not the Sobel test (Sobel assumes normality of a×b and is underpowered).

• Exposure–mediator interaction: with the counterfactual approach, report the natural direct and

indirect effects allowing X×M interaction rather than assuming it away.

Reporting (AGReMA)

Report against AGReMA (A Guideline for Reporting Mediation Analyses). Minimum:

• The full mediation model (every path), the estimator, and the confounder set for each of the

X→Y, X→M, and M→Y relationships (they differ — name them separately).

• Total, direct, and indirect effects, each with a bootstrap CI.
• Proportion mediated only with uncertainty, and only when the total effect is well-estimated.

Proportion mediated = indirect / total is unstable when the total effect is small or near-null — it can exceed 100% or flip sign. Do not headline a proportion-mediated when the total-effect CI is wide or crosses the null.

• The sensitivity analysis for unmeasured mediator–outcome confounding: an **E-value for the

indirect effect** (or a ρ-/correlation-based mediation sensitivity, e.g. mediation::medsens), stating how strong unmeasured M–Y confounding would have to be to null the indirect effect.

• The temporal structure of X, M, Y (and, if cross-sectional, the explicit caveat).

Common failures (flag at review)

• Causal "X affects Y through M" from cross-sectional data with no temporal caveat → reframe to

association-level + add the sensitivity analysis (O13).

• A significant bootstrap CI treated as evidence of identification.
• No sensitivity analysis for unmeasured M–Y confounding.
• Proportion-mediated quoted as a stable figure when the total effect is small.
• Binary-outcome mediation via the naive OR product instead of a counterfactual decomposition.

Anti-Hallucination

• Never fabricate path coefficients, indirect effects, or bootstrap CIs — all from executed code.
• Do not assert a causal direction the design cannot support; state the identification assumptions.
• Generate references via /search-lit (AGReMA, VanderWeele).

Missing Data Guide

Strategies for handling missing data in medical research analyses. This is a preprocessing step, not an independent analysis type.

---

Missing Data Mechanisms

MNAR cannot be verified statistically — always consider as a possibility.

---

Decision by Missing Rate

---

Multiple Imputation via Chained Equations (MICE)

Procedure

1. Specify imputation model for each variable with missing data 2. Generate m imputed datasets (m = 5 minimum; m = 20 if missing > 20%) 3. Analyze each dataset independently (same analysis model) 4. Pool results using Rubin's rules: estimate = mean of m estimates; variance = within + between variance

Imputation Model Rules

• Include ALL variables from the analysis model in the imputation model
• Include the outcome variable (in imputation model only)
• Match variable type to method:
• Continuous → predictive mean matching (PMM) or regression
• Binary → logistic regression
• Multinomial → polytomous regression
• Ordinal → proportional odds

Methods NOT Recommended

• Mean imputation: underestimates variance — never use
• LOCF (Last Observation Carried Forward): biased in most settings — avoid
• Single imputation: does not account for imputation uncertainty

---

Reporting Template

"Missing data ranged from X% (variable A) to Y% (variable B). Little's MCAR test suggested data were not MCAR (P = 0.003). Multiple imputation using chained equations (MICE) was performed with m = 20 imputed datasets. Results were pooled using Rubin's rules. Sensitivity analyses using complete case analysis yielded consistent results (Supplementary Table X)."

---

When to Trigger in analyze-stats

During Phase 1 (Data Assessment), if any analysis variable has > 5% missing: 1. Report missing counts and percentages 2. Suggest MICE before proceeding 3. Generate imputation code as a preprocessing step 4. Run the primary analysis on imputed data 5. Run complete case analysis as sensitivity analysis

---

Python Implementation

from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
import numpy as np
import pandas as pd

m = 20  # number of imputations
results = []
for i in range(m):
    imputer = IterativeImputer(max_iter=10, random_state=i, sample_posterior=True)
    df_imputed = pd.DataFrame(imputer.fit_transform(df), columns=df.columns)
    # Run analysis on each imputed dataset
    # ... append results
# Pool using Rubin's rules
pooled_estimate = np.mean(results, axis=0)

R Implementation

library(mice)
imp <- mice(df, m = 20, method = 'pmm', seed = 42)
fit <- with(imp, lm(y ~ x1 + x2))
pooled <- pool(fit)
summary(pooled)

---

Common Reviewer Flags

1. Not reporting missing data counts per variable 2. Using listwise deletion without justification 3. Mean imputation or LOCF without acknowledging limitations 4. Not performing sensitivity analysis (complete case vs MI comparison) 5. Not stating the assumed missing mechanism (MCAR/MAR/MNAR)

NHIS Claims-Based Disease Definition Guide

Methodological patterns for defining exposures, outcomes, and comorbidities using ICD-10 codes in Korean NHIS claims data. This guide covers the validated algorithm patterns -- NOT specific code lists for individual diseases.

---

When to Use

• Defining study populations, exposures, or outcomes in NHIS cohort studies
• Building claims-based algorithms for disease ascertainment
• Validating disease definitions using multiple data sources within NHIS
• NOT for: survey data (KNHANES/NHANES), clinical registry data

---

Core Principle: Claims-Based Algorithms

A single ICD-10 claim is insufficient for disease ascertainment due to rule-out diagnoses, coding errors, and provisional codes. Validated algorithms combine multiple data elements to improve positive predictive value.

Algorithm Components

---

Validated Algorithm Patterns

Pattern 1: N-Claim Rule (Most Common)

Require N or more claims with the target ICD-10 code within a defined time window to establish a diagnosis.

Disease = (ICD-10 code X appears >= N times within T months)

Reporting template: "[DISEASE] was identified by at least [N] claims with ICD-10 codes [CODE RANGE] within [T] months during the observation period."

Pattern 2: Claim + Medication (Recommended for Chronic Diseases)

Combine diagnostic codes with disease-specific medication prescriptions.

Disease = (ICD-10 code X >= 1 time) AND (related medication dispensed >= 1 time)

Reporting template: "[DISEASE] was defined as at least one claim with ICD-10 codes [CODE RANGE] combined with a prescription for [MEDICATION CLASS] (ATC codes [RANGE])."

Common medication validation pairs:

Pattern 3: Look-Back Period for Incident Cases

To identify incident (new-onset) cases, require a disease-free washout period before the index date.

Incident case = (ICD-10 code X after index date)
              AND (NO claims with code X during look-back period)

Reporting template: "To ensure the identification of incident cases, we applied a minimum [N]-year look-back period and excluded individuals with any diagnostic claims for [DISEASE] during this period."

Pattern 4: Outcome Validation with Time Windows

For complications or sequelae, require temporal proximity to the index event.

Complication = (ICD-10 code Y within T1-T2 days after event X)

Examples of time-window patterns:

Pattern 5: Hierarchical Disease Classification

When exposure or outcome categories overlap, use a hierarchical scheme.

Level 1 (broadest): Any code in ICD-10 block (e.g., E10-E14 for all diabetes)
Level 2: Specific subcategory (e.g., E11 for type 2 diabetes only)
Level 3 (strictest): Subcategory + medication + lab confirmation

Sensitivity analyses should test multiple levels to assess robustness.

---

ICD-10-KM vs ICD-10 Differences

Korean NHIS uses ICD-10-KM (Korean Modification), which is based on ICD-10 but includes Korea-specific extensions:

For most major disease categories (cardiovascular, metabolic, respiratory, neurological), the ICD-10-KM codes are identical to international ICD-10. Country-specific extensions primarily affect rare or culture-bound conditions.

---

Charlson Comorbidity Index in NHIS

The CCI is computed from claims using the Quan adaptation for ICD-10:

Always specify which CCI adaptation was used (Quan 2005 is standard for ICD-10 claims data).

---

Covariate Definitions from NHIS Data Sources

NHIS integrates three data sources, each providing different covariate types:

Source 1: Claims Database

• Comorbidities (via ICD-10 codes)
• Medication history (via dispensing records)
• Healthcare utilization (visit counts, hospitalization days)

Source 2: National Health Examination

• BMI, blood pressure, fasting glucose, GFR
• Liver enzymes (AST, ALT, gamma-GTP)
• Hemoglobin, total cholesterol
• Urinalysis

Source 3: Health Interview / Questionnaire

• Smoking status (never / former / current)
• Alcohol consumption (frequency categories)
• Physical activity (MET-based sufficiency)
• Region of residence, household income level

Baseline alignment: Covariates must be anchored to the most recent data available prior to the index date. Do not use post-index measurements.

---

Sensitivity Analysis Patterns

Stricter Disease Definition

Run the primary analysis with the standard algorithm, then repeat with a more stringent definition (e.g., 2-claim → 3-claim, or adding medication confirmation).

Negative Control Outcomes

Select outcomes with no plausible biological link to the exposure (e.g., neuromuscular junction disorders, self-harm) to detect residual confounding or coding artifacts.

Negative Control Exposures

Use an exposure with no expected effect on the outcome as a falsification test.

Multiple Look-Back Periods

Test 1-year, 3-year, and 5-year look-back periods and compare incidence rates to assess the impact of prevalent case misclassification.

---

Reporting Checklist

When reporting NHIS claims-based definitions in a manuscript:

• [ ] ICD-10 code ranges specified (main text or supplementary table)
• [ ] Number of claims required (1 vs 2+ vs with medication)
• [ ] Time window for claims (within 12 months, any time, etc.)
• [ ] Look-back period for incident case definition
• [ ] Medication confirmation specified with ATC codes (if used)
• [ ] Data source specified (claims vs exam vs questionnaire)
• [ ] CCI adaptation cited (Quan 2005 or other)
• [ ] Sensitivity analysis with alternative definition included

---

Common Reviewer Flags

1. Single-claim definition without justification (low PPV concern) 2. No look-back period for "incident" cases 3. Mixing ICD-10 and ICD-9 codes without noting transition dates 4. Not specifying whether inpatient, outpatient, or both claims were used 5. CCI version not cited 6. Covariate timing not anchored to index date 7. Missing sensitivity analysis with stricter/looser definitions 8. Rule-out diagnoses not addressed (diagnostic codes without treatment)

---

References (Methodological)

• Quan H et al. Coding algorithms for defining comorbidities in ICD-9-CM

and ICD-10 administrative data. Med Care. 2005;43(11):1130-9.

• Choi JY et al. Validation of administrative data for the identification

of chronic diseases in South Korea. J Korean Med Sci. (methodology reference)

• NHIS data user guide: nhiss.nhis.or.kr (Korean National Health Insurance

Sharing Service)

Propensity Score Analysis Guide

Methods for estimating causal treatment effects in observational studies by balancing confounders between treatment groups.

---

When to Use

• Observational study comparing treatment vs control (or two interventions)
• Multiple confounders to adjust for
• Goal: estimate causal effect analogous to an RCT
• NOT for: randomized trials (already balanced), single-arm studies

---

Estimands — Choose Before Analysis

Comparing PSM and IPTW results directly is inappropriate — they estimate different estimands.

---

Step-by-Step Workflow

Step 1: PS Estimation

• Model: logistic regression (standard)
• Dependent variable: treatment assignment (binary)
• Covariates: all variables that affect the outcome (confounders)
• Do NOT include: instrumental variables (affect treatment but not outcome)
• Individual PS model coefficients have no clinical meaning — only PS distribution matters

Step 2: Apply PS Method

Option A — PS Matching (PSM)

• Nearest-neighbor matching with caliper = 0.2 x SD(logit PS)
• 1:1 matching is standard; 1:N or full matching available
• Estimand: ATT (typically)
• Drawback: unmatched subjects are excluded → sample size reduction

Option B — IPTW (Inverse Probability of Treatment Weighting)

• Weights: treated = 1/PS, control = 1/(1-PS) for ATE
• Always use stabilized weights to prevent extreme values
• Stabilized: treated = P(T=1)/PS, control = P(T=0)/(1-PS)
• All subjects included (no exclusion)
• Flag extreme weights > 10

Option C — SIPTW (Stabilized Inverse Probability of Treatment Weighting)

• Weights: treated = P(T=1)/PS, control = P(T=0)/(1-PS)
• Mathematically equivalent to stabilized IPTW (Option B with stabilized=True)
• Named explicitly as SIPTW in emulated target trial literature
• Maintains entire cohort sample size (no exclusion)
• Allows appropriate variance estimation of main effect
• Estimand: ATE
• Increasingly used in large-scale NHIS cohort studies
• Report effective sample size (ESS) alongside raw N

Option D — Overlap Weighting (Recommended for most cases)

• Weights: treated = (1-PS), control = PS
• Naturally down-weights subjects at PS extremes
• No extreme weight problem (advantage over IPTW)
• Estimand: ATO
• Increasingly recommended in recent guidelines (JAMA 2020, AJE 2024)

Step 3: Assess Balance

• SMD < 0.10 for all covariates (Austin, 2009)
• SMD < 0.25 is acceptable but suboptimal
• Use SMD, NOT p-values (SMD is sample-size independent)
• Love plot: pre/post-matching SMD comparison (mandatory figure)
• Variance ratios: should be within 0.5-2.0
• PS distribution overlap: histogram comparing treated vs control

Step 4: Outcome Analysis

• After PSM: paired analysis (paired t-test, conditional logistic, stratified Cox)
• After IPTW/OW: weighted regression using survey methods
• Always use robust/sandwich standard errors for weighted analyses

Step 5: Sensitivity Analysis

• E-value: quantifies how strong unmeasured confounding would need to be to explain away the result
• Report E-value for the point estimate and lower CI bound

---

Balance Table (Required Output)

---

Reporting Templates

PSM: "Propensity scores were estimated using logistic regression with the following covariates: [list]. PS matching was performed using 1:1 nearest-neighbor matching with a caliper of 0.2 SD of the logit PS. After matching, all SMDs were below 0.10 (Figure X). In the matched cohort (n = X pairs), ..."

IPTW/OW: "Inverse probability of treatment weighting (or overlap weighting) was applied using stabilized weights. Covariate balance was assessed using SMDs (all < 0.10; Figure X). The weighted analysis showed ..."

SIPTW: "Stabilized inverse probability of treatment weighting was used to balance covariate distributions between the [exposed] and [unexposed] groups. This approach maintains the sample size of the entire cohort and allows for appropriate estimation of the variance of the main effect. Covariate balance was assessed using SMDs (all < 0.10; Figure X)."

---

Common Reviewer Flags

1. SMD not reported (using p-values instead) 2. Caliper width not specified 3. Estimand (ATE/ATT/ATO) not stated 4. No sensitivity analysis for unmeasured confounding 5. Number of unmatched subjects not reported (for PSM) 6. Extreme weights not assessed (for IPTW) 7. Individual PS model coefficients interpreted clinically

---

Python Packages

• sklearn.linear_model.LogisticRegression — PS estimation
• causalinference — matching (limited)
• Manual implementation for IPTW/OW (see template)
• statsmodels — weighted regression

R Packages

• MatchIt — PS matching
• WeightIt — IPTW, overlap weighting
• cobalt — balance assessment, Love plot
• survey — weighted outcome analysis
• tableone — baseline table with SMD
• EValue — sensitivity analysis

Regression Analysis Guide

Covers both logistic regression (binary outcome) and multiple linear regression (continuous outcome).

---

Logistic Regression

When to Use

Binary outcome variable (0/1, event/no event) with one or more predictors. Used for risk factor identification and prediction model building.

Assumptions

1. Binary outcome 2. Linear relationship between continuous predictors and log-odds (Box-Tidwell test) 3. Independent observations (if repeated → GEE or mixed logistic) 4. No multicollinearity: VIF < 5 5. Sufficient sample: EPV >= 10 (minimum), >= 20 (recommended)

Variable Selection

• Clinical rationale first — avoid purely data-driven stepwise selection
• Include: variables significant at P < 0.10 in univariable analysis + known confounders
• STROBE/TRIPOD guideline compliance required

Model Assessment

Calibration:

• Hosmer-Lemeshow test: P > 0.05 = adequate fit (overpowered with large N)
• Calibration plot: predicted probability vs observed frequency

Discrimination:

• C-statistic (= AUC): 0.7-0.8 acceptable, 0.8-0.9 excellent, > 0.9 outstanding

Multicollinearity:

• VIF > 5 → remove or combine variables

Required Outputs

1. OR table: univariable AND multivariable OR (95% CI), P-value per variable 2. C-statistic with 95% CI 3. Hosmer-Lemeshow test result 4. VIF table (supplementary) 5. Box-Tidwell results for continuous predictors (supplementary)

OR Table Format

Reporting Template

"Multivariable logistic regression was performed to identify independent predictors of [outcome]. Variables with P < 0.10 in univariable analysis and clinically relevant confounders were included. The model showed good discrimination (C-statistic = 0.82, 95% CI 0.78-0.86) and calibration (Hosmer-Lemeshow P = 0.45). [Variable] was independently associated with [outcome] (adjusted OR = 2.15, 95% CI 1.43-3.24; P < 0.001)."

Pitfalls

• OR != RR: when event rate > 10%, OR overestimates RR
• EPV < 10 → overfitting risk. Consider penalized regression (LASSO/Ridge)
• Specify reference category for categorical variables
• Specify units for continuous variables (per 1 year vs per 10 years)
• Complete separation → use Firth's penalized likelihood

---

Multiple Linear Regression

When to Use

Continuous outcome variable with one or more predictors. Used for identifying determinants and estimating adjusted effects.

Assumptions (LINE + No Multicollinearity)

1. Linearity: residuals vs fitted plot 2. Independence: no repeated measures (if repeated → LMM/GEE) 3. Normality of residuals: Q-Q plot, Shapiro-Wilk on residuals 4. Equal variance (homoscedasticity): residuals vs fitted, Scale-Location plot 5. No multicollinearity: VIF < 5 6. No influential outliers: Cook's distance < 4/n

Assumption Violations → Alternatives

Model Evaluation

• R² (coefficient of determination): proportion of variance explained
• Adjusted R²: penalized for number of predictors — use for model comparison
• In medical research, R² = 0.2-0.4 can be meaningful (high biological variability)

Diagnostic Plots (4-panel, always generate)

1. Residuals vs Fitted → linearity + homoscedasticity 2. Q-Q plot → residual normality 3. Scale-Location → homoscedasticity 4. Residuals vs Leverage → influential outliers (Cook's distance)

Required Outputs

1. Coefficient table: β (95% CI), P-value per variable 2. R² and adjusted R² 3. VIF table 4. 4-panel diagnostic plot 5. Standardized coefficients (optional, for effect size comparison)

Coefficient Table Format

Reporting Template

"Multiple linear regression was performed with [outcome] as the dependent variable. The model explained X% of the variance (adjusted R² = 0.XX). After adjusting for [covariates], [variable] was significantly associated with [outcome] (β = X.XX, 95% CI X.XX to X.XX; P = exact). Model assumptions were verified using diagnostic plots."

Pitfalls

• Always report β units (per 1 year, per 10 kg/m², etc.)
• Standardized β useful for comparing effect sizes but unstandardized is standard in papers
• EPV for continuous outcome: N >= 10-20 per predictor
• Always present diagnostic plots (at minimum in supplementary)

Repeated Measures / Mixed Models Guide

Analysis methods for longitudinal data where the same subjects are measured at multiple time points.

---

When to Use

• Same subjects measured at 2+ time points
• Key research question: does change over time differ between groups? (Time x Group interaction)
• Examples: treatment response over weeks, serial imaging measurements, before/after/follow-up

---

Method Selection

Default recommendation: LMM — handles missing data, does not require sphericity, allows unequal time spacing.

---

Data Format

All methods (LMM, GEE) require long format:

Convert from wide format before analysis.

---

1. RM ANOVA

Sphericity Check (Mandatory)

• Mauchly's test: H0 = sphericity holds
• P >= 0.05 → use standard F-test
• P < 0.05 → apply correction

Corrections for Sphericity Violation

1. Greenhouse-Geisser (G-G): conservative, recommended 2. Huynh-Feldt: slightly liberal 3. Multivariate test (Pillai's trace): no sphericity assumption needed

Limitations

• Complete cases only — any missing time point drops the entire subject
• Difficult to add covariates
• Equal time spacing required

Key Results to Report

• Mauchly's test result (W statistic, P-value, epsilon)
• Correction method used (if sphericity violated)
• Within-subject effect: Time (F, df, P)
• Between-subject effect: Group (F, df, P)
• Interaction: Time x Group (most important)

---

2. Linear Mixed Model (LMM)

Structure

Y = Xβ (fixed effects) + Zb (random effects) + ε

Random Effects Selection

• Random intercept only: subjects differ in baseline (default, start here)
• Random intercept + slope: subjects differ in both baseline and rate of change
• Decision: visualize individual trajectories (spaghetti plot). If slopes vary → add random slope.
• If model does not converge with random slope → simplify to random intercept only

Covariance Structure Selection

Compare structures using AIC/BIC (lower = better).

Advantages over RM ANOVA

• Handles missing data under MAR assumption
• No sphericity requirement
• Allows unequal time intervals
• Easy covariate adjustment
• Estimates individual trajectories

---

3. GEE (Generalized Estimating Equations)

Key Properties

• Estimates population-averaged (marginal) effects — not individual-level
• Specify working correlation: exchangeable, AR(1), unstructured
• Works for non-normal outcomes (binary, count)
• Requires MCAR assumption for missing data (or use weighted GEE for MAR)
• Needs sufficient clusters: subjects >= 30-40

LMM vs GEE Decision

---

Required Outputs

1. Spaghetti plot: individual trajectories by group 2. Model summary table: fixed effects (β, 95% CI, P), random effects variance 3. Time x Group interaction result prominently reported 4. Model fit: AIC/BIC for covariance structure comparison 5. Missing data description (n per time point, mechanism assumed)

---

Reporting Templates

RM ANOVA: "Repeated-measures ANOVA was performed with Greenhouse-Geisser correction for violation of sphericity (Mauchly's test P < 0.001, ε = 0.42). There was a significant time × group interaction (F(X, Y) = Z, P = exact)."

LMM: "A linear mixed-effects model with random intercepts for subjects and [CS/AR(1)] correlation structure was fitted. The time × group interaction was significant (β = -2.34, 95% CI -3.87 to -0.81; P = 0.003), indicating that the rate of change in [outcome] differed between groups."

GEE: "GEE with exchangeable correlation structure was used to estimate population-averaged effects. The time × group interaction was ..."

---

Common Reviewer Flags

1. Sphericity test not reported (for RM ANOVA) 2. Covariance structure selection rationale not stated (for LMM) 3. Time × Group interaction not interpreted 4. Missing data handling not described 5. Using RM ANOVA with substantial missing data (should use LMM) 6. Not reporting individual trajectories (spaghetti plot)

---

Python Packages

• pingouin — RM ANOVA with G-G correction
• statsmodels.formula.api.mixedlm — LMM
• statsmodels.genmod.generalized_estimating_equations.GEE — GEE

R Packages

• lme4 + lmerTest — LMM (standard)
• nlme — LMM with correlation structures (CS, AR1, UN)
• geepack — GEE
• Base R aov() with Error() — RM ANOVA

Survey-Weighted Analysis Guide

Methods for analyzing complex survey data (KNHANES, NHANES, KCHS, and similar nationally representative health surveys) that use stratified, multistage probability sampling designs.

---

When to Use

• Data from a national health survey with sampling weights (e.g., KNHANES, NHANES, KCHS)
• Goal: produce nationally representative prevalence, means, or associations
• Cross-national comparisons using parallel survey datasets
• NOT for: simple random samples, census data, or claims-based cohorts (NHIS, JMDC)

Claims-based databases (NHIS, JMDC) are NOT surveys -- they do not have sampling weights or complex sampling design. Use standard regression for these.

---

Key Concepts

Complex Survey Design Elements

Ignoring these elements produces biased standard errors, incorrect p-values, and non-representative point estimates.

Dataset-Specific Design Variables

Weight Selection Rules

• Interview-only variables: use interview weight
• Exam/lab variables: use exam weight (smaller denominator)
• Nutrition variables (KNHANES): use nutrition weight
• Multi-cycle NHANES: divide weight by number of cycles combined (e.g., 4 cycles: weight/4)
• Single-cycle analysis: use the cycle-specific weight as-is

Subpopulation (domain) analysis — never row-delete

A restricted analysis (adults only, one sex, a disease subgroup) must keep the full design and select the domain, not filter the data frame and refit. Row-deletion discards the strata/PSU structure of the dropped units and gives wrong standard errors and design degrees of freedom.

• R survey: subset(design, age >= 18) on the design object (or svyby), not svydesign(data = df[df$age>=18, ]).
• Stata: svy, subpop(if age>=18): — never keep if age>=18 before svy:.
• Python samplics / R is preferred; statsmodels has no native domain estimator.

Reporting & common errors (these invalidate the inference, flag at review)

• Model-based SEs on weighted points. Applying the weight but computing SEs without strata + PSU (or replicate weights) understates uncertainty. Always declare strata + id/cluster.
• Weighted total ≠ sample size. Report the unweighted n as the analytic sample; the weighted figure is a population estimate, not "n".
• Design effect / effective n. Report DEFF or the effective sample size where precision is load-bearing; a large DEFF means far fewer independent observations than rows.
• Unweighted vs weighted divergence. If they differ materially, that signals weight-dependent selection — discuss it, do not hide it.
• Data-driven thresholds. A restricted-cubic-spline "non-linear/saturation" claim needs a pre-specified non-linearity test (LRT/Wald vs the linear model), and any "inflection point / threshold" from a recursive breakpoint search must carry a confidence interval and a stability check — a searched cutoff is exploratory, not a validated target (review-side probe O12 in observational_confounding.md).

---

Step-by-Step Workflow

Step 1: Declare Survey Design

Always declare the design before any analysis. This ensures correct variance estimation.

Python (statsmodels):

# statsmodels does not have a native survey design object.
# Use linearmodels or manual weight application.
# For publication-quality survey analysis, R is strongly recommended.

R (survey package):

library(survey)

# KNHANES
design_kr <- svydesign(
  id = ~psu,
  strata = ~kstrata,
  weights = ~wt_itvex,
  data = df_kr,
  nest = TRUE
)

# NHANES (2-year cycle)
design_us <- svydesign(
  id = ~SDMVPSU,
  strata = ~SDMVSTRA,
  weights = ~WTMECXYR,
  data = df_us,
  nest = TRUE
)

SAS:

/* KNHANES */
PROC SURVEYLOGISTIC DATA=kr;
  STRATA kstrata;
  CLUSTER psu;
  WEIGHT wt_itvex;
  MODEL outcome(event='1') = exposure covariates;
RUN;

/* NHANES */
PROC SURVEYLOGISTIC DATA=us;
  STRATA SDMVSTRA;
  CLUSTER SDMVPSU;
  WEIGHT WTMECXYR;
  MODEL outcome(event='1') = exposure covariates;
RUN;

Step 2: Weighted Descriptive Statistics

R:

# Weighted means
svymean(~continuous_var, design, na.rm = TRUE)

# Weighted proportions
svymean(~factor(categorical_var), design, na.rm = TRUE)

# Weighted Table 1 by group
library(tableone)
svyCreateTableOne(
  vars = c("age", "sex", "bmi", "income"),
  strata = "exposure_group",
  data = design,
  test = TRUE
)

SAS:

PROC SURVEYMEANS DATA=dataset;
  STRATA strata_var;
  CLUSTER cluster_var;
  WEIGHT weight_var;
  VAR continuous_var1 continuous_var2;
RUN;

PROC SURVEYFREQ DATA=dataset;
  STRATA strata_var;
  CLUSTER cluster_var;
  WEIGHT weight_var;
  TABLES group * categorical_var / CHISQ;
RUN;

Step 3: Weighted Regression — Sequential Model Building

The standard cross-national analysis pattern uses sequential model building:

R (survey-weighted logistic regression):

# Model 1: age + sex
model1 <- svyglm(
  outcome ~ exposure + age + sex,
  design = design,
  family = quasibinomial()
)

# Model 2: full adjustment
model2 <- svyglm(
  outcome ~ exposure + age + sex + income + education +
            smoking + alcohol + bmi + cvd_history,
  design = design,
  family = quasibinomial()
)

# Extract weighted OR (wOR) with 95% CI
extract_wor <- function(model, var) {
  coef_val <- coef(model)[var]
  se_val <- summary(model)$coefficients[var, "Std. Error"]
  or <- exp(coef_val)
  ci_lo <- exp(coef_val - 1.96 * se_val)
  ci_hi <- exp(coef_val + 1.96 * se_val)
  p_val <- summary(model)$coefficients[var, "Pr(>|t|)"]
  data.frame(wOR = or, CI_lower = ci_lo, CI_upper = ci_hi, P = p_val)
}

SAS (PROC SURVEYLOGISTIC):

/* Model 2: full adjustment */
PROC SURVEYLOGISTIC DATA=dataset;
  STRATA strata_var;
  CLUSTER cluster_var;
  WEIGHT weight_var;
  CLASS sex(ref='Male') income(ref='High') smoking(ref='Never') / PARAM=REF;
  MODEL outcome(event='1') = exposure age sex income education smoking alcohol bmi;
  ODDSRATIO exposure / CL=WALD;
RUN;

Step 4: Subgroup / Stratified Analyses

# Stratified by sex
svyglm(
  outcome ~ exposure + age + income + education + smoking + alcohol + bmi,
  design = subset(design, sex == "Male"),
  family = quasibinomial()
)
# Repeat for Female
# Note: exclude the stratification variable from covariates

Reporting pattern: "Weighted odds ratios are adjusted for all covariates except for the stratification variable."

Step 5: Advanced Analyses (Optional)

Restricted cubic spline (dose-response):

library(rms)
design_rms <- svydesign(id = ~psu, strata = ~kstrata,
                         weights = ~wt_itvex, data = df)
model_rcs <- svyglm(
  outcome ~ rcs(continuous_exposure, 3) + age + sex + covariates,
  design = design_rms,
  family = quasibinomial()
)

Weighted quantile sum (WQS) regression:

library(gWQS)
# WQS for composite exposure (e.g., LE8 components)
result_wqs <- gwqs(
  outcome ~ wqs + age + sex + covariates,
  mix_name = c("comp1", "comp2", "comp3", "comp4"),
  data = df,
  q = 4,           # quartiles
  b = 500,         # bootstrap iterations
  seed = 42,
  family = "binomial",
  weights = df$weight_var
)

---

Weighted SMD (Standardized Mean Difference)

For balance assessment in weighted analyses:

library(survey)
library(tableone)

# SMD in survey design
tab <- svyCreateTableOne(
  vars = covariates,
  strata = "treatment",
  data = design,
  smd = TRUE
)
print(tab, smd = TRUE)

Manual calculation (Python):

def weighted_smd(x, treatment, weights, is_binary=False):
    """Calculate weighted standardized mean difference."""
    t_mask = treatment == 1
    c_mask = treatment == 0
    
    w1, w0 = weights[t_mask], weights[c_mask]
    x1, x0 = x[t_mask], x[c_mask]
    
    wm1 = np.average(x1, weights=w1)
    wm0 = np.average(x0, weights=w0)
    
    if is_binary:
        denom = np.sqrt((wm1 * (1 - wm1) + wm0 * (1 - wm0)) / 2)
    else:
        wv1 = np.average((x1 - wm1) ** 2, weights=w1)
        wv0 = np.average((x0 - wm0) ** 2, weights=w0)
        denom = np.sqrt((wv1 + wv0) / 2)
    
    return (wm1 - wm0) / denom if denom > 0 else 0.0

---

Cross-National Analysis Pattern

When comparing two countries using parallel surveys:

1. Never pool raw data across countries into a single regression 2. Analyze each country independently with country-specific survey design 3. Present results side-by-side in the same table 4. Compare effect magnitudes narratively (not via interaction terms)

Standard Output Table Format

---

Reporting Templates

Methods: "To account for the complex survey designs and ensure nationally representative estimates, appropriate stratification, clustering, and sampling weights were applied. Weighted logistic regression models were used to estimate weighted odds ratios (wORs) and 95% confidence intervals (CIs). Model 1 adjusted for age and sex, while Model 2 further included [covariates]. All analyses were performed using SAS version 9.4 (SAS Institute Inc.) [and R version X.X.X (R Foundation for Statistical Computing)]. A two-sided P value of less than 0.05 was considered statistically significant."

Results: "In the weighted analysis of [N] participants from [DATASET], [EXPOSURE] was significantly associated with [OUTCOME] (wOR [X.XX]; 95% CI [X.XX-X.XX]) after adjusting for [covariates] (Model 2)."

---

Common Reviewer Flags

1. Survey weights not applied (unweighted analysis of survey data) 2. Strata/cluster variables not specified (incorrect SE estimation) 3. Wrong weight variable used (interview weight for lab variables) 4. Multi-cycle NHANES weights not adjusted (divided by number of cycles) 5. Data pooled across countries instead of analyzed separately 6. Weighted proportions not reported (using raw counts instead) 7. Subgroup analysis includes the stratification variable as covariate 8. Missing data handling not stated (complete case vs imputation)

---

Python vs R Recommendation

For publication-quality survey analysis, R is strongly recommended over Python. Python's statsmodels supports basic weighted regression but lacks full survey design support (no strata/cluster specification for variance estimation).

---

R Packages

• survey -- core survey design and analysis
• tableone -- survey-weighted baseline tables with SMD
• rms -- restricted cubic splines with survey design
• gWQS -- weighted quantile sum regression
• srvyr -- tidyverse-compatible survey analysis wrapper

SAS Procedures

• PROC SURVEYMEANS -- weighted means and proportions
• PROC SURVEYFREQ -- weighted frequency tables and chi-square
• PROC SURVEYLOGISTIC -- weighted logistic regression (wOR)
• PROC SURVEYREG -- weighted linear regression
• PROC SURVEYPHREG -- weighted Cox proportional hazards

Statistical Test Selection Guide

Decision tree for selecting the appropriate statistical test based on data structure. Reference: Petrie & Sabin flowchart, Kirkwood Summary Guide.

---

Step 1: Outcome Variable Type

---

Step 2A: Continuous Outcome

---

Step 2B: Categorical Outcome

---

Step 3: Confounder Control

---

Normality Assessment

Practical rule: n > 30 → t-test is generally robust (CLT), unless extreme skew or outliers.

---

Common Reviewer Flags

• Using independent test on paired data
• Using chi-squared when expected cell count < 5 (should use Fisher's exact)
• Not reporting assumption check results
• Missing multiple comparison correction for 3+ groups
• Not specifying the test selection rationale in Methods

Publication Table Standards — Knowledge Base

Reference document for medical journal table formatting.

Source: YouTube tutorials, journal author guidelines, AMA Manual of Style, tool documentation.

Last updated: 2026-04-11

---

1. Universal Rules (All Medical Journals)

1. No vertical lines — horizontal rules only (top, below header, bottom) 2. Editable Word tables — use Word Insert > Table. Never submit images, Excel, or PDF 3. Sequential numbering — Table 1, 2, 3... in order of first citation in text 4. Define all abbreviations — in footnotes, independently for each table 5. Explanatory content in footnotes, not headings 6. Display item limit — most journals: tables + figures combined 4-5 7. Exact P values — never just "significant/not significant" 8. Variability measures required — always state mean (SD) or median (IQR) 9. Self-contained titles — table title alone must convey content without reading text 10. No duplication — tables and figures must not repeat the same data

---

2. Journal-Specific Differences

Footnote Markers

P Value Formatting

P Value Decimal Places

Horizontal Lines

95% CI Format

Table Title Format

---

3. AMA Manual of Style (11th Ed) — Table Rules

Structure

• Column headers: Bold
• Row headers (stub): Regular weight, sentence case
• Decking (nested column headers): Maximum 3 levels
• Alignment: Text = left-aligned, Numbers = center or decimal-aligned
• Units: In column header parentheses, not repeated in cells
• Maximum size: ~40 rows x 6-8 columns, fit within 1 page

Footnote System (AMA 11th)

• Table footnotes: Superscript lowercase letters (a, b, c...)
• NOT symbols — symbols (*, †, ‡) are for bottom-of-page text footnotes only
• Assignment order: Top-to-bottom, left-to-right through the table
• Whole-table footnote: Place superscript "a" at end of table title
• Reference numbers + footnotes: Reference first, comma, then footnote letter (e.g., 5,b)

Footnote Placement Order

Below the table, in this sequence:

1. General note (no marker, applies to entire table)
   "Data are presented as median (IQR) unless otherwise noted."

2. Abbreviation note
   "ASA = American Society of Anesthesiologists, BMI = body mass index,
    CI = confidence interval."
   → Listed in order of appearance in table (left→right, top→bottom)
   → Some journals (JAMA) require alphabetical order — check guide

3. Specific notes (superscript letter markers)
   "a Excludes patients with missing follow-up data."
   "b Adjusted for age and sex."

4. Probability notes (asterisk markers, if used)
   "* P < .05; ** P < .01."

P Value Rules

• Uppercase italic: P
• No leading zero: .05, not 0.05
• Thin space around operators: P = .03
• Only P, α, β omit leading zero; all other statistics keep it

---

4. Footnote System Specification

When to Use Symbols vs Letters vs Numbers

Caption vs Footnote — What Goes Where

Statistical Test Footnote Patterns

Pattern A — Individual markers (preferred when tests vary by row):

Age, y                    65.3 (12.1)    62.1 (11.8)    .04ᵃ
Sex, male                 53 (52)        48 (47)        .41ᵇ

ᵃ Wilcoxon rank-sum test.
ᵇ Fisher exact test.

Pattern B — General note (preferred when same test for all rows of a type):

P values were calculated using the Wilcoxon rank-sum test for continuous
variables and the Fisher exact test for categorical variables.

---

5. Common Mistakes Checklist

Structure

• [ ] Binary variables: show only one level (Male 53%, not Male 53% / Female 47%)
• [ ] Remove derivable columns (Total column when groups are shown)
• [ ] Keep tables under 50 rows (excess → supplementary)
• [ ] No sub-part numbering (Table 3A/3B → Table 3, Table 4)
• [ ] Do not repeat all table data in Results text

Formatting

• [ ] Consistent decimal places within each column
• [ ] Units in column headers, not in cells
• [ ] No vertical lines, minimal horizontal lines
• [ ] Numbers right-aligned or center-aligned (not left)
• [ ] Specify variability: "Mean (SD)" or "Median (IQR)" in header or footnote

Statistics

• [ ] RCTs: P values in Table 1 are usually unnecessary (randomization)
• [ ] Always name the statistical test (in footnote or general note)
• [ ] Report effect sizes per clinically meaningful unit (OR per 10-year, not per 1-year)
• [ ] CI notation consistent throughout (parentheses vs "to" — match journal style)
• [ ] Never use "NS" — report exact P values

Submission

• [ ] Editable Word table (not image/screenshot)
• [ ] No color as sole information carrier (use bold/italic/symbols instead)
• [ ] All abbreviations defined in footnotes
• [ ] Each table's footnotes are self-contained (no "see Table 1 footnote")

---

6. Tool Recommendations

Primary Pipeline: R {gtsummary}

Why gtsummary:

• Built-in journal themes: JAMA, Lancet, NEJM
• Auto-selects statistical tests (Wilcoxon, Fisher, Chi-square)
• Auto-generates footnotes (test names, summary statistics)
• Outputs to Word (via flextable), LaTeX (via huxtable), HTML (via gt)

Core API:

library(gtsummary)

# Set journal theme FIRST
theme_gtsummary_journal("jama")  # or "lancet", "nejm"
theme_gtsummary_compact()

# Table 1
tbl <- df %>%
  tbl_summary(
    by = group,
    type = list(age ~ "continuous2"),
    statistic = list(
      all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})"),
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = list(all_continuous() ~ 1),
    missing = "ifany"
  ) %>%
  add_p() %>%
  add_overall() %>%
  add_stat_label() %>%
  bold_labels() %>%
  modify_footnote(all_stat_cols() ~ "Mean (SD); Median (Q1, Q3); n (%)")

# Export
tbl %>% as_flex_table() %>% flextable::save_as_docx(path = "table1.docx")
tbl %>% as_hux_table() %>% huxtable::to_latex() %>% writeLines("table1.tex")

Regression table:

model <- glm(outcome ~ age + sex + bmi, data = df, family = binomial)

tbl_regression(model, exponentiate = TRUE) %>%
  add_global_p() %>%
  bold_p(t = 0.05) %>%
  bold_labels()

Supporting Tools

Python Limitation

Python has no equivalent to gtsummary's journal themes or automatic statistical test selection. For publication tables, R is strongly recommended. Python-only pipelines should use tableone for generation + python-docx for Word formatting, accepting manual style adjustments.

---

7. Format-Specific Notes

Word (DOCX) Submission

• Use Word's Insert > Table function (never tab-separated)
• Font: Times New Roman 10-12pt (table body can be 9-10pt)
• Single spacing within cells
• Remove all vertical borders; keep 3 horizontal lines
• Bold headers, regular body text
• Each table on a separate page, after References

LaTeX Submission

• Use booktabs package (\toprule, \midrule, \bottomrule)
• Never use \hline or | vertical separators
• \siunitx for decimal alignment
• \multirow / \multicolumn for merged cells
• Generate via gtsummary → huxtable → to_latex() for consistency

HTML (Review/Proofing)

• gt produces the highest quality HTML tables
• Useful for co-author review before final Word/LaTeX export
• Can include interactive features (sorting, filtering) for supplementary materials

Table: Reliability / Agreement Results

Reporting Guidelines

• GRRAS: Reliability and agreement studies — report the coefficient, its model/definition,

the 95% CI, and the number of raters/subjects/replicates.

• Distinguish relative reliability (how well raters/methods separate subjects relative to

total variance — a consistency ICC) from absolute agreement (how close the actual values are — an absolute-agreement ICC, or Bland–Altman LoA for two continuous methods). κ measures categorical agreement corrected for chance. Report the metric that matches the question, and both when relevant.

Standard Structure

Table 3. Inter-rater Reliability and Agreement for [measurement]
(n = [N] subjects, [R] raters, [k] replicates)

Measure / metric          Estimate (95% CI)        Notes

ICC (continuous)          0.88 (0.82-0.92)         two-way random, absolute agreement,
                                                   single measures [ICC(2,1)]
Weighted κ (ordinal)      0.74 (0.65-0.82)         quadratic weights
Bland-Altman bias          1.2 (0.6 to 1.8)        mean difference (rater A - rater B)
  Limits of agreement     -4.8 to 7.2              ±1.96 SD; CI on each limit in footnote
Percent agreement         86%                      within ±[clinically acceptable Δ]

ICC = intraclass correlation coefficient. Heteroscedasticity assessed (LoA constant across
the measurement range). Interpretation: ICC >0.90 excellent, 0.75-0.90 good, 0.50-0.75
moderate, <0.50 poor.

Rules

• ICC is uninterpretable without its model and definition — always state: one-way vs

two-way; consistency vs absolute agreement; single vs average measures (the Shrout–Fleiss form, e.g., ICC(2,1)). Report the point estimate with its 95% CI.

• Ordinal categories → weighted κ (state linear or quadratic weights); unweighted κ ignores

the magnitude of disagreement. Nominal → unweighted κ (or Fleiss' κ for >2 raters). Report the CI; note that κ is prevalence-sensitive (report observed agreement alongside it).

• Agreement of two continuous methods → Bland–Altman: mean difference (bias) with its CI,

and the limits of agreement (bias ± 1.96 SD) with a CI on each limit; assess heteroscedasticity (do the LoA widen with magnitude? — if so, log-transform or model the SD). State whether the LoA fall within a pre-defined clinically acceptable difference.

• Never report Pearson/Spearman correlation as agreement — high correlation is compatible

with large systematic bias; correlation measures association, not agreement.

• Sample: report n subjects, n raters/methods, replicates per subject, and a sample-size /

precision justification for the CI width.

• Match the metric to the data and the question; do not report an ICC for nominal categories or

a κ for continuous measurements.

Python / R Code

library(irr); library(psych); library(blandr)

# ICC — specify model & type explicitly (psych::ICC reports all six forms)
psych::ICC(ratings_wide)                 # pick the row matching your design, e.g. ICC2 (2,1)

# Weighted kappa for ordinal (two raters)
irr::kappa2(cbind(rater_a, rater_b), weight = "squared")   # quadratic weights

# Bland-Altman (two continuous methods)
blandr::blandr.statistics(method_a, method_b)   # bias, LoA, and CIs; plot via blandr.draw()

import pingouin as pg
pg.intraclass_corr(data=df, targets="subject", raters="rater", ratings="value")  # all ICC forms + CI
# pg.cohen_kappa(...) for categorical; Bland-Altman bias = mean(diff), LoA = mean ± 1.96*sd

Table: Diagnostic Accuracy Results

Reporting Guidelines

• STARD: Sensitivity, Specificity, PPV, NPV with 95% CIs
• TRIPOD-AI / CLAIM: For AI model performance

Standard Structure

Table 2. Diagnostic Performance of [Model/Test] for [Condition]

                     Sensitivity       Specificity       PPV              NPV              AUC
                     (95% CI)          (95% CI)          (95% CI)         (95% CI)         (95% CI)

Model A              0.92 (0.87-0.96)  0.85 (0.79-0.90)  0.78 (0.71-0.84) 0.95 (0.91-0.97) 0.94 (0.91-0.97)
Model B              0.88 (0.82-0.93)  0.90 (0.85-0.94)  0.83 (0.76-0.88) 0.93 (0.89-0.96) 0.93 (0.89-0.96)
P value              .12               .08               .21              .34              .45

95% CIs were calculated using the Wilson score method. AUC comparison
by DeLong test.

Rules

• Always include 95% CIs for all metrics (STARD requirement)
• CI method: Specify in footnote (Wilson score, Clopper-Pearson, DeLong)
• Threshold: State the decision threshold used (e.g., "at Youden optimal threshold")
• Per-class results: For multi-class, show per-class + macro/micro average
• Comparison P values: DeLong test for AUC, McNemar for sensitivity/specificity
• Decimal places: 2-3 for proportions (0.92), 2-3 for AUC (0.94)
• Reader studies: Include per-reader AND pooled results

Key Footnote Content

• CI calculation method
• Comparison test used (DeLong, McNemar, bootstrap)
• Threshold selection method
• Whether results are per-patient or per-lesion

Tables: Meta-Analysis

Reporting Guidelines

• PRISMA 2020: Study characteristics + pooled results
• PRISMA-DTA: For diagnostic test accuracy meta-analyses

---

Table A: Characteristics of Included Studies

Table 1. Characteristics of Included Studies

Author, Year   Country   Design   N      Population       Index Test        Reference Standard   Quality
Kim 2023       Korea     Retro    450    Suspected PE     CTPA AI (v2.1)    Expert consensus     Low risk
Smith 2024     USA       Prosp    1200   ED patients      CTPA AI (v3.0)    Pulmonary DSA        Some concerns
...

Retro = retrospective, Prosp = prospective, PE = pulmonary embolism,
CTPA = CT pulmonary angiography, ED = emergency department,
DSA = digital subtraction angiography.

Rules

• Column order: Author/Year, Country, Design, N, Population, Index test, Reference standard, Quality/RoB
• Author format: First author surname + year
• Study design: Use standard abbreviations (Retro, Prosp, RCT)
• Quality assessment tool: QUADAS-2 (DTA), RoB 2 (RCT), NOS (observational)
• Quality rating: "Low risk" / "Some concerns" / "High risk" (QUADAS-2 terms)

---

Table B: Pooled Results / Summary Estimates

Table 2. Pooled Diagnostic Accuracy Estimates

                    k    N       Pooled Estimate (95% CI)    I²     P_het

Sensitivity         12   3400    0.91 (0.87-0.94)           78%    <.001
Specificity         12   3400    0.88 (0.83-0.92)           65%    .002
Positive LR         12   3400    7.58 (5.12-11.2)           72%    <.001
Negative LR         12   3400    0.10 (0.07-0.15)           69%    .001
DOR                 12   3400    75.8 (42.1-136.5)          58%    .008

k = number of studies, N = total participants, CI = confidence interval,
LR = likelihood ratio, DOR = diagnostic odds ratio.
Pooled estimates from bivariate random-effects model.
I² = Higgins inconsistency statistic; P_het from Cochran Q test.

Rules

• k and N: Always report number of studies and total participants
• Heterogeneity: I² + P from Cochran Q test (mandatory)
• Model: State random-effects vs fixed-effects in footnote
• Subgroup analyses: Separate rows or separate table
• Prediction interval: Include for random-effects if k >= 3
• Forest plot complement: Table complements but does not replace forest plot