Python Biostats Cheatsheet

Lab reference for statistical computing in clinical research. Python + LLM-assisted vibe-coding is the preferred stack over R for new analyses in the Holubar Lab — see the Vibe Coding primer for the rationale and workflow.

Core stack — install once, use forever

# With uv (recommended) or pip:
uv add pandas numpy scipy statsmodels scikit-learn matplotlib seaborn \
       lifelines pingouin tableone forestplot duckdb polars pyarrow \
       python-docx python-pptx Pillow pydicom

Reporting conventions — match SDH's manuscripts

Every published number in the lab follows these rules. Hardcode them in your analysis scripts so the output is paste-ready:

• P-values: 2 significant non-zero decimals. p = 0.04, p = 0.0003. Never p < 0.001 — show the actual value.
• Effect sizes: aRR 0.63 (95% CI 0.51–0.76, p = 0.001) — en-dash for CI range.
• IQR: 3 (1 - 7) with spaces around dash.
• Continuous: mean ± SD or median (IQR), max 1 decimal.
• Categorical: n (%) with the percent sign — 423 (48.2%).
• Tables: AMA three-line border style — top rule, header-row underline, bottom rule. No interior verticals.
• Figures: TIFF 600 DPI, Arial, color-blind-safe palette. Horizontal text only — never rotate axis labels.

Full conventions live in .claude/rules/statistical-output.md + manuscript-conventions.md. Auto-loaded for every Claude session.

Idiomatic patterns SDH uses

Pattern 1 — Compose, don't overwrite

Never overwrite an original data column. Derived variables get a new name with a suffix:

# WRONG — destroys canonical pexy column
df["pexy"] = df["pexy"] | df["pouch_inlet_transposition"]

# RIGHT — new column with composite name
df["pexy_or_transposition"] = (df["pexy"] | df["pouch_inlet_transposition"]).astype(int)

# Conventions:
df["bmi_n"]       = pd.to_numeric(df["bmi"], errors="coerce")     # _n = numeric
df["dod_dt"]      = pd.to_datetime(df["dod"], errors="coerce")    # _dt = datetime
df["asa_high"]    = df["asa"].str.startswith(("3","4")).astype(int)  # _high = binary high
df["twist_deg_cat"] = pd.cut(df["degree_twisted"], …)             # _cat = categorical

Pattern 2 — GEE with site clustering (default for NSQIP analyses)

import statsmodels.api as sm

cov = ["age", "asa_high", "bmi_n", "mis", "ibd_dx_uc", "year"]
X = sm.add_constant(df[cov])
model = sm.GEE(
    df["organ_space_ssi"], X,
    groups=df["site_id"],
    family=sm.families.Poisson(),
    cov_struct=sm.cov_struct.Exchangeable(),
).fit()
print(model.summary())
arr = np.exp(model.params)
ci  = np.exp(model.conf_int())

Pattern 3 — Kaplan-Meier with at-risk table

from lifelines import KaplanMeierFitter
from lifelines.plotting import add_at_risk_counts

kmf_uc = KaplanMeierFitter().fit(uc.t_years, uc.died, label="UC")
kmf_cd = KaplanMeierFitter().fit(cd.t_years, cd.died, label="CD")

ax = kmf_uc.plot_survival_function(ci_show=True)
kmf_cd.plot_survival_function(ax=ax, ci_show=True)
add_at_risk_counts(kmf_uc, kmf_cd, ax=ax)
plt.savefig("Figures/Figure_2_clean.tiff", dpi=600, format="tiff",
            pil_kwargs={"compression": "tiff_lzw"})

Pattern 4 — DuckDB for big CSVs

import duckdb
# Query the CSV directly — no pandas load needed
out = duckdb.sql("""
    SELECT site, COUNT(*) n, AVG(los_days) mean_los
    FROM 'Data/nsqip_ibd_2025.csv'
    WHERE ibd_diagnosis = 'UC'
    GROUP BY site
    ORDER BY n DESC
""").df()

Pattern 5 — python-docx surgical edits (NEVER rewrite manuscripts from scratch)

from docx import Document
TEMPLATE = "SDH Admin/SDH Letterhead Headshots Templates/Canonical SDH_Letterhead 2026.docx"
doc = Document(TEMPLATE)   # ✓ inherits letterhead + signature
# Walk runs, surgically replace {PLACEHOLDER}, preserve bold/italic
for p in doc.paragraphs:
    for run in p.runs:
        if "{TITLE}" in run.text:
            run.text = run.text.replace("{TITLE}", "Apixaban for VTE in IBD")
doc.save("Letter_v1.docx")

The vibe-coding loop (with Claude / Cursor)

1. Describe the analysis in plain English — "Run a Cox model on this cohort with site clustering, output the HR table as AMA-style python-docx."
2. Let the LLM write the script — paste the data shape and columns. Claude/Cursor produce working Python far more reliably than R for this kind of task.
3. Read it before you run it — every line. If you don't understand it, ask the model to explain. Vibe coding ≠ blind execution.
4. Run, eyeball results, iterate — "the HRs look right but the CI is reversed; can you check the conf_int() output?" — the model fixes it.
5. Cascade through your deliverables — see project-structure.md for the Code → JSON → Notebook → Memo → Manuscript → Figures cascade rule.

See the Vibe Coding primer for the longer treatment, including pitfalls and how to evaluate LLM-generated code for clinical-research correctness.

Exceptions where R remains canonical

• Network meta-analysis — netmeta (R). No equivalent Python package. Drive from Python via rpy2.
• Cost-effectiveness modeling at full Drummond/CHEERS rigor — BCEA, hesim, dampack (R). The Holubar Lab CEA pipeline (/cea-pipeline/) is Python-native but reaches for these via rpy2 for specific PSA / EVPPI work.
• Legacy CCF biostatistician hand-offs — if Dr. Tang or Li Tang is the analyst of record, they may prefer R. Don't fight it.

R Cheatsheets — still useful

The R Cheatsheets card holds 12 RStudio one-page references (tidyverse, lubridate, purrr, shiny, plumber). Useful for reading legacy CCF analyses, debugging rpy2 calls, and reviewing fellow scripts in R. Lab default for new work is still Python.