Tooluniverse Statistical Modeling

Statistical Modeling for Biomedical Data Analysis

RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:

• *_executed.ipynb → read with tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' and cite its cell outputs as the authoritative answer
• Pre-computed result files (CSV/TSV with names like *results*, *deseq*, *enrich*, *stats*, *_simplified.csv) → read directly and report the requested value
• Canonical analysis scripts (analysis.R, run_*.py, find_*.R, *.Rmd) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if none of the above exist. Re-running from raw data produces different numbers than the published answer and is much slower (often 5-10× turn count).

---

PRIMARY SCRIPTS — use these FIRST

These scripts encode the question-specific gotchas in scripts/ and emit labelled, parseable output. Prefer them over ad-hoc statsmodels / scipy code.

Concrete invocations

Natural-spline regression (R^2, overall F-test p, peak Y + 95% CI):

python skills/tooluniverse-statistical-modeling/scripts/r_natural_spline_regression.py \
  --csv data.csv --y-col Area \
  --ratio-col Ratio --new-x-col Frequency_strain \
  --filter "StrainNumber not in ['1', '98']" \
  --df 4 --workdir /tmp/spline_run

Quadratic vs cubic vs natural-spline comparison + best-model peak:

python skills/tooluniverse-statistical-modeling/scripts/spline_model_compare.py \
  --csv data.csv --y-col Area \
  --ratio-col Ratio --new-x-col Frequency_strain \
  --filter "StrainNumber not in ['1', '98']" \
  --ns-df 4 --workdir /tmp/spline_cmp

*Report the peak location (`x) in the units of the fitted x-variable, not a derived label.** When the model is fit on a frequency/proportion column (e.g. Frequency_strain, a 0–1 value), the answer to "at what ratio/frequency is the maximum" is that fraction (e.g. 0.909), NOT the colon-ratio it was derived from (e.g. 10:1). Convert a colon ratio a:b to the fraction a/(a+b)` when the question expects a 0–1 value or the fitted x-column is a fraction.

Ordinal logistic regression with interaction term (e.g. trial AE severity):

python skills/tooluniverse-statistical-modeling/scripts/logistic_regression_or.py \
  --csv merged.csv --outcome AESEV --outcome-type ordinal --outcome-order "1,2,3,4" \
  --predictors TRTGRP,expect_interact,patients_seen,MHONGO \
  --encode TRTGRP,expect_interact,patients_seen \
  --encode-map "TRTGRP:Placebo=0,BCG=1" \
  --interaction MHONGO:TRTGRP_cat \
  --coef-name TRTGRP_cat

Two-sample power analysis from a pilot CSV:

python skills/tooluniverse-statistical-modeling/scripts/power_analysis.py \
  --csv pilot.csv --value-col MeasuredValue --group-col Group \
  --group-a Treatment --group-b Control \
  --power 0.8 --alpha 0.05

---

Workspace isolation (CRITICAL)

The input data folder for any analysis must remain untouched so re-runs are reproducible. Scripts that write intermediate files (R drivers, prepared CSVs, comparison tables) must write to /tmp/ or to a --workdir you pass in. Both R-based scripts in this skill refuse to run if --workdir resolves to the input CSV's parent directory (or any ancestor of it).

# OK
--workdir /tmp/spline_run

# Refused:
--workdir <path-equal-to-or-containing-the-input-csv>/...

---

CRITICAL — Read before writing any code

1. Clinical trial AE analysis (regression, chi-square, ANY severity test): Use the bundled script (or the clinical_trial_ae_severity_test ToolUniverse tool which wraps it):

   tu run clinical_trial_ae_severity_test '{"dm_file":"DM.csv","ae_file":"AE.csv","test":"chi-square","group_col":"TRTGRP"}'

   # Or directly:
   python skills/tooluniverse-statistical-modeling/scripts/prepare_ae_cohort.py \
     --dm DM.csv --ae AE.csv --test chi-square --group TRTGRP \
     --subgroup "expect_interact=Yes"   # optional

The script/tool handles: encoding='latin1' for SDTM CSVs, max(AESEV) per subject across ALL AEs (no AEPT filtering), inner join with DM, optional subgroup filter, optional ordinal-logistic with covariates.

Why no AEPT filter — AESEV is a protocol-defined severity scale on the AE table. Filtering AE by AEPT (e.g. keeping only AEPT == "COVID-19") drops subjects whose worst severity was recorded under a different AEPT label, drastically changes the contingency table, and can flip the test result. The phrase "COVID-19 severity" describes the OUTCOME, NOT a filter criterion.

• ❌ WRONG: ae[ae['AEPT'].str.contains('COVID-19')].groupby('USUBJID')['AESEV'].max() — filters to COVID-19 events
• ✅ RIGHT: ae.groupby('USUBJID')['AESEV'].max() — uses ALL AE records

2. Expression ANOVA / fold change with multi-feature data (gene × sample matrix): For "the F-statistic" or "a fold change" as a single value, run per-gene then summarize — NEVER pool expr.values.ravel() across all genes.

• For F-statistic: derive a per-sample quantity (like DESeq2 LFC of each gene between two cell types, then ANOVA on those LFCs across groups) OR run on a single target gene.
• For median/mean log2 fold change between two groups: run DESeq2 with design=~<group>, extract per-gene log2FoldChange (with shrinkage if the pipeline uses it), then take median/mean across genes.

❌ WRONG (aggregate): log2(sum_counts_groupA / sum_counts_groupB) per sample then summarize — gives ratio of totals, dominated by high-expression genes. ✅ RIGHT (per-gene): DESeq2 → results_df['log2FoldChange'].median().

Sanity heuristics: F > 50 for biological ANOVA across a few groups means you aggregated (typical biological F is 0.5–10). |median LFC| > 2 between similar groups means you aggregated (typical |median| < 1).

Use the bundled script: python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py (or the expression_anova_per_gene ToolUniverse tool).

3. Spline models — R splines::ns(x, df=K) ≠ Python patsy.dmatrix("cr(x, df=K)"). They produce different design matrices because of internal-knot placement, boundary-knot placement, and basis orthogonalization. For ANY question that references R syntax like lm(y ~ ns(x, df = 4)), run R via Rscript. Use the bundled wrapper:

   python skills/tooluniverse-statistical-modeling/scripts/r_natural_spline_regression.py \
     --csv data.csv --y-col Y --x-col X --df 4 --workdir /tmp/spline_run

For "frequency of strain X" co-culture models, include pure focal strain (freq=1) but exclude non-focal pure strain (freq=0).

4. CSV encoding: Clinical trial CSVs often need encoding='latin1'.

5. Pearson correlation between count-like and length-like variables: when one variable spans orders of magnitude (raw read counts, TPM, gene length, transcript abundance), raw Pearson r is often near 0 even when log-transformed r is moderate. ALWAYS compute and explicitly report ALL FOUR variants in your final answer body: r(x, y), r(log10(x+1), y), r(x, log10(y+1)), r(log10(x+1), log10(y+1)). List as a table; mark one as your primary pick. The published answer can be ANY of the four, and the question text rarely disambiguates which transform combination was used.

   ## Primary answer: r = X.XXX (transform: <name>)

   ## Sensitivity (all 4 transform combinations):
   - r(x, y)                   = ...
   - r(log10(x+1), y)          = ...
   - r(x, log10(y+1))          = ...
   - r(log10(x+1), log10(y+1)) = ...

Background — for any single transform variant:

   import numpy as np
   from scipy.stats import pearsonr
   r_raw, _ = pearsonr(x, y)
   r_log, _ = pearsonr(np.log10(x + 1), y)
   print(f"r_raw={r_raw:.4f}  r_log10={r_log:.4f}")

Defaults:

• Question says "log-transformed" / "log expression" → report r_log10
• Question doesn't specify but the variable is gene expression / RNA count → also report r_log10 as the canonical answer (most published correlations between gene length and expression are log-scale)
• When |r_raw| < 0.1 AND |r_log10| > 0.2, prefer r_log10

❌ WRONG: report only r_raw ≈ 0.05 when log is 0.35 ✅ RIGHT: "r_raw = 0.05; r_log10 = 0.35 (canonical for log-distributed expression)"

---

COMPUTE, DON'T DESCRIBE

Write and run Python code (via Bash) for every statistical analysis. Never describe what you "would do" — do it. Use pandas for data wrangling, statsmodels for regression, scipy for tests, and matplotlib for plots. Execute the code and report actual numbers (β, p-value, CI, N).

LOOK UP, DON'T GUESS

When uncertain about any scientific fact, SEARCH databases first rather than reasoning from memory.

Features

• Linear Regression - OLS for continuous outcomes with diagnostic tests
• Logistic Regression - Binary, ordinal, and multinomial models with odds ratios
• Survival Analysis - Cox proportional hazards and Kaplan-Meier curves
• Mixed-Effects Models - LMM/GLMM for hierarchical/repeated measures data
• ANOVA - One-way/two-way ANOVA, per-feature ANOVA for omics data
• Model Diagnostics - Assumption checking, fit statistics, residual analysis
• Statistical Tests - t-tests, chi-square, Mann-Whitney, Kruskal-Wallis, etc.

When to Use

Apply this skill when user asks:

• "What is the odds ratio of X associated with Y?"
• "What is the hazard ratio for treatment?"
• "Fit a linear regression of Y on X1, X2, X3"
• "Perform ordinal logistic regression for severity outcome"
• "What is the Kaplan-Meier survival estimate at time T?"
• "What is the percentage reduction in odds ratio after adjusting for confounders?"
• "Run a mixed-effects model with random intercepts"
• "Compute the interaction term between A and B"
• "What is the F-statistic from ANOVA comparing groups?"
• "Test if gene/miRNA expression differs across cell types"

Model Selection Decision Tree

START: What type of outcome variable?
|
+-- CONTINUOUS (height, blood pressure, score)
|   +-- Independent observations -> Linear Regression (OLS)
|   +-- Repeated measures -> Mixed-Effects Model (LMM)
|   +-- Count data -> Poisson/Negative Binomial
|
+-- BINARY (yes/no, disease/healthy)
|   +-- Independent observations -> Logistic Regression
|   +-- Repeated measures -> Logistic Mixed-Effects (GLMM/GEE)
|   +-- Rare events -> Firth logistic regression
|
+-- ORDINAL (mild/moderate/severe, stages I/II/III/IV)
|   +-- Ordinal Logistic Regression (Proportional Odds)
|
+-- MULTINOMIAL (>2 unordered categories)
|   +-- Multinomial Logistic Regression
|
+-- TIME-TO-EVENT (survival time + censoring)
    +-- Regression -> Cox Proportional Hazards
    +-- Survival curves -> Kaplan-Meier

Workflow

Phase 0: Data Validation

Goal: Load data, identify variable types, check for missing values.

CRITICAL: Identify the Outcome Variable First

Before any analysis, verify what you're actually predicting:

1. Read the full question - Look for "predict [outcome]", "model [outcome]", or "dependent variable" 2. Examine available columns - List all columns in the dataset 3. Match question to data - Find the column that matches the described outcome 4. Verify outcome exists - Don't create outcome variables from predictors

Common mistake: Question mentions "obesity" -> Assumed outcome = BMI >= 30 (circular logic with BMI predictor). Always check data columns first: print(df.columns.tolist())

import pandas as pd
import numpy as np

df = pd.read_csv('data.csv')
print(f"Observations: {len(df)}, Variables: {len(df.columns)}, Missing: {df.isnull().sum().sum()}")

for col in df.columns:
    n_unique = df[col].nunique()
    if n_unique == 2:
        print(f"{col}: binary")
    elif n_unique <= 10 and df[col].dtype == 'object':
        print(f"{col}: categorical ({n_unique} levels)")
    elif df[col].dtype in ['float64', 'int64']:
        print(f"{col}: continuous (mean={df[col].mean():.2f})")

Phase 1: Model Fitting

Goal: Fit appropriate model based on outcome type.

Use the decision tree above to select model type, then refer to the appropriate reference file for detailed code:

• Linear Regression: references/linear_models.md
• Logistic Regression (binary): references/logistic_regression.md
• Ordinal Logistic: references/ordinal_logistic.md
• Cox Proportional Hazards: references/cox_regression.md
• ANOVA / Statistical Tests: anova_and_tests.md

Quick reference for key models:

import statsmodels.formula.api as smf
import numpy as np

# Linear regression
model = smf.ols('outcome ~ predictor1 + predictor2', data=df).fit()

# Logistic regression (odds ratios)
model = smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0)
ors = np.exp(model.params)
ci = np.exp(model.conf_int())

# Cox proportional hazards
from lifelines import CoxPHFitter
cph = CoxPHFitter()
cph.fit(df[['time', 'event', 'treatment', 'age']], duration_col='time', event_col='event')
hr = cph.hazard_ratios_['treatment']

Phase 1b: ANOVA for Multi-Feature Data

When data has multiple features (genes, miRNAs, metabolites), use per-feature ANOVA (not aggregate). This is the most common pattern in genomics.

See anova_and_tests.md for the full decision tree, both methods, and worked examples.

Default for gene expression data: Per-feature ANOVA (Method B).

Phase 2: Model Diagnostics

Goal: Check model assumptions and fit quality.

Key diagnostics by model type:

• OLS: Shapiro-Wilk (normality), Breusch-Pagan (heteroscedasticity), VIF (multicollinearity)
• Cox: Proportional hazards test via cph.check_assumptions()
• Logistic: Hosmer-Lemeshow, ROC/AUC

See references/troubleshooting.md for diagnostic code and common issues.

Phase 3: Interpretation

Goal: Generate publication-quality summary.

For every result, report: effect size (OR/HR/coefficient), 95% CI, p-value, and model fit statistic. See common_patterns_summary.md for common question-answer patterns.

Common Patterns

See common_patterns_summary.md for solution code for each pattern. See references/common_patterns.md for 15+ detailed question patterns.

Statsmodels vs Scikit-learn

General rule: Use statsmodels for statistical inference questions (p-values, ORs, HRs).

Python Package Requirements

statsmodels>=0.14.0
scikit-learn>=1.3.0
lifelines>=0.27.0
pandas>=2.0.0
numpy>=1.24.0
scipy>=1.10.0

Key Principles

1. Data-first approach - Always inspect and validate data before modeling 2. Model selection by outcome type - Use decision tree above 3. Assumption checking - Verify model assumptions (linearity, proportional hazards, etc.) 4. Complete reporting - Always report effect sizes, CIs, p-values, and model fit statistics 5. Confounder awareness - Adjust for confounders when specified or clinically relevant 6. Reproducible analysis - All code must be deterministic and reproducible 7. Robust error handling - Graceful handling of convergence failures, separation, collinearity 8. Round correctly - Match the precision requested (typically 2-4 decimal places)

Reasoning Framework for Result Interpretation

Evidence Grading

Interpretation Guidance

• Model diagnostics (R-squared): For OLS, R-squared > 0.7 indicates good fit in biomedical data; 0.3-0.7 is moderate. Adjusted R-squared penalizes added predictors. For logistic models, use pseudo-R-squared (McFadden > 0.2 is acceptable) and AUC (> 0.7 = acceptable, > 0.8 = good discrimination).
• AIC/BIC for model comparison: Lower is better. AIC difference > 10 between models is strong evidence for the lower-AIC model. BIC penalizes complexity more heavily than AIC, preferring simpler models. Use AIC for prediction-focused selection, BIC for inference.
• Coefficient significance thresholds: Report exact p-values rather than just significance stars. For multiple predictors, apply Bonferroni or FDR correction. A coefficient with p = 0.049 in a model with 20 predictors is likely a false positive without correction.
• Survival analysis HR interpretation: HR > 1 means increased hazard (worse outcome) for the exposed group. HR = 2.0 means twice the instantaneous risk of the event. Always verify the proportional hazards assumption -- if violated, the HR is an average over time and may be misleading. Report median survival times alongside HRs for clinical interpretability.
• Odds ratio interpretation: OR = 1.0 means no association. OR > 1 indicates increased odds. The 95% CI excluding 1.0 confirms significance. For rare outcomes, OR approximates relative risk; for common outcomes (> 10% prevalence), OR overstates the relative risk.
• Confounding assessment: Compare crude vs adjusted ORs/HRs. A change > 10% in the effect estimate after adjusting for a covariate suggests confounding by that variable.

Synthesis Questions

1. Do the model diagnostics (residual plots, Hosmer-Lemeshow, PH test) support the validity of the chosen model, or do assumption violations require alternative approaches (e.g., robust standard errors, stratified models)? 2. For adjusted models, does the inclusion of confounders change the primary effect estimate by more than 10%, indicating meaningful confounding? 3. Are the reported effect sizes (OR, HR, coefficients) clinically meaningful in addition to being statistically significant, considering the scale of the predictor and outcome? 4. When comparing nested models via AIC/BIC, does the more complex model provide substantially better fit, or is the simpler model preferred by parsimony? 5. For survival analysis, is the proportional hazards assumption met throughout the follow-up period, or do Schoenfeld residuals suggest time-varying effects?

---

Completeness Checklist

Before finalizing any statistical analysis:

• [ ] Outcome variable identified: Verified which column is the actual outcome
• [ ] Data validated: N, missing values, variable types confirmed
• [ ] Multi-feature data identified: If multiple features, use per-feature approach
• [ ] Model appropriate: Outcome type matches model family
• [ ] Assumptions checked: Relevant diagnostics performed
• [ ] Effect sizes reported: OR/HR/Cohen's d with CIs
• [ ] P-values reported: With appropriate correction if needed
• [ ] Model fit assessed: R-squared, AIC/BIC, concordance
• [ ] Results interpreted: Plain-language interpretation
• [ ] Precision correct: Numbers rounded appropriately

Bundled Scripts

These ready-to-run scripts live in skills/tooluniverse-statistical-modeling/scripts/. Use them via the Bash tool — they are the deterministic answer for the recurring question patterns documented above.

r_natural_spline_regression.py — Natural spline regression in R

Shells out to Rscript to fit lm(y ~ ns(x, df=K)) with splines::ns(). Emits R², adj R², F-stat with df1/df2, overall F-test p-value, residual SE, coefficient table (estimate, SE, t, p), and the prediction-grid peak with 95% CI from predict.lm(..., interval='confidence'). Supports a --ratio-col shortcut to convert "a:b" string ratios into a frequency fraction a/(a+b). Refuses to write into the input CSV's parent directory.

spline_model_compare.py — Quadratic vs cubic vs natural spline

Fits all three models on the same x,y in R, ranks by adjusted R², AIC, and BIC, and emits the best model's peak (x, y) with 95% CI. Use for "best-fitting model" questions and "maximum predicted y at optimal x".

logistic_regression_or.py — Binary or ordinal logistic regression with ORs

Fits sm.Logit (binary) or OrderedModel (ordinal proportional-odds) and emits ORs (exp(coef)) plus 95% CIs and p-values for every coefficient. Handles label encoding (--encode A,B,C), explicit value maps (--encode-map TRTGRP:Placebo=0,BCG=1), and interaction columns (--interaction A:B -> creates A_B = A*B). With --coef-name <NAME> also prints a SCALARS block tagged for the requested coefficient.

power_analysis.py — Two-sample required-N for a t-test

Computes Cohen's d (pooled SD) from a CSV given --value-col, --group-col, --group-a, --group-b, then TTestIndPower.solve_power with --alpha, --power, --alternative. Use for "minimum sample size per group" power questions. Returns both the raw and the ceil-ed N.

stat_tests.py — Basic statistical tests (pure stdlib, no scipy)

Implements chi-square goodness-of-fit, Fisher's exact test, and simple linear regression without any external dependencies. All p-values are computed from first principles using the gamma function (chi-square) or hypergeometric enumeration (Fisher's).

# Chi-square goodness-of-fit
python stat_tests.py --type chi_square --observed 100,50,25 --expected 87.5,50,37.5

# Fisher's exact test (2×2 table)
python stat_tests.py --type fisher_exact --a 10 --b 5 --c 3 --d 20
python stat_tests.py --type fisher_exact --a 10 --b 5 --c 3 --d 20 --alternative greater

# Simple linear regression (OLS)
python stat_tests.py --type regression --x "1,2,3,4,5" --y "2.1,4.0,5.9,8.1,10.0"

Key formulas:

• chi_square: χ² = Σ (O−E)²/E; p-value via upper regularized incomplete gamma Q(df/2, χ²/2)
• fisher_exact: hypergeometric PMF; p-value = sum of probabilities ≤ P(observed)
• regression: b1 = Sxy/Sxx; b0 = ȳ − b1x̄; R² = 1 − SSR/SST; SE and t-statistics included

Output includes: full contingency/data table, step-by-step arithmetic, significance statement, and a round-trip verification for each test.

When to use stat_tests.py vs statsmodels:

• Use stat_tests.py when you need a quick sanity check with no imports, or when the

environment lacks scipy/statsmodels.

• Use statsmodels when you need multivariate regression, logistic models, or survival analysis.

format_statistical_output.py — Format results for reporting

Utility functions to format fitted statsmodels results as publication-ready tables. Import and call from analysis scripts; not a standalone CLI tool.

model_diagnostics.py — Automated model diagnostics

Runs assumption checks (normality, heteroscedasticity, multicollinearity) on fitted models. Import and call from analysis scripts; not a standalone CLI tool.

---

File Structure

tooluniverse-statistical-modeling/
+-- SKILL.md                          # This file (workflow guide)
+-- QUICK_START.md                    # 8 quick examples
+-- EXAMPLES.md                       # Legacy examples
+-- TOOLS_REFERENCE.md                # ToolUniverse tool catalog
+-- anova_and_tests.md                # ANOVA decision tree and code
+-- common_patterns_summary.md         # Common solution patterns
+-- test_skill.py                     # Test suite
+-- references/
|   +-- logistic_regression.md        # Detailed logistic examples
|   +-- ordinal_logistic.md           # Ordinal logit guide
|   +-- cox_regression.md             # Survival analysis guide
|   +-- linear_models.md              # OLS and mixed-effects
|   +-- common_patterns.md            # 15+ question patterns
|   +-- troubleshooting.md            # Diagnostic issues
+-- scripts/
    +-- r_natural_spline_regression.py # lm(y ~ ns(x, df=K)) via Rscript
    +-- spline_model_compare.py        # quadratic vs cubic vs natural-spline (Rscript)
    +-- logistic_regression_or.py      # binary / ordinal logistic + ORs + interactions
    +-- power_analysis.py              # TTestIndPower required-N from CSV
    +-- expression_anova.py            # per-gene ANOVA / log2FC summary
    +-- prepare_ae_cohort.py           # SDTM AE/DM cohort prep
    +-- stat_tests.py                  # Chi-square, Fisher's exact, OLS (stdlib)
    +-- format_statistical_output.py   # Format results for reporting
    +-- model_diagnostics.py           # Automated diagnostics

ToolUniverse Integration

While this skill is primarily computational, ToolUniverse tools can provide data:

See TOOLS_REFERENCE.md for complete tool catalog.

References

• statsmodels: https://www.statsmodels.org/
• lifelines: https://lifelines.readthedocs.io/
• scikit-learn: https://scikit-learn.org/
• Ordinal models: statsmodels.miscmodels.ordinal_model.OrderedModel

Analysis conventions

These conventions are validated best practices. Apply when the dataset/question matches.

MANDATORY: Use bundled script for expression ANOVA / fold change

For per-gene ANOVA or median log2FC questions, use the bundled script:

python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py \
  --counts counts.csv --meta meta.csv --group-col cell_type \
  --exclude-groups PBMC --mode anova

python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py \
  --counts counts.csv --meta meta.csv --group-col cell_type \
  --group-a CD14 --group-b CD19 --mode fold_change

Do NOT write your own pandas ANOVA — the aggregation level (per-gene, not per-sample) is critical and easy to get wrong.

CSV encoding

Clinical-trial exports (SDTM) are often latin1. If pd.read_csv() fails with UnicodeDecodeError, retry with encoding='latin1'.

Clinical-trial AE analysis — applies to regression AND chi-square AND any severity test

For any statistical test of a clinical-trial AE severity outcome (chi-square, ordinal/logistic regression, Mann-Whitney, etc.) on trial covariates, fit/test on the cohort that reported any AE (inner-joined to demographics). Do NOT pre-filter AE records by AEPT — NOT even when the question says "COVID-19 severity" or names any specific condition. Use max(AESEV) across all AE records per subject, regardless of AEPT.

Why: AESEV on the AE table reflects the study's protocol-defined severity scale. Pre-filtering to specific AEPT values (e.g., keeping only certain condition labels) drops subjects whose worst severity was recorded under a different AEPT label, which drastically changes the contingency table and can flip results from significant to non-significant.

Do NOT pad subjects with no AEs as AESEV=0 — that dilutes the signal.

dm = pd.read_csv("DM.csv", encoding='latin1')
ae = pd.read_csv("AE.csv", encoding='latin1')
sev = ae.groupby('USUBJID')['AESEV'].max().reset_index()
df = dm.merge(sev, on='USUBJID', how='inner').dropna(subset=['AESEV'])
df['AESEV'] = df['AESEV'].astype(int)
# Now ordinal regression / chi-square on df

Odds-ratio deviation

"Percentage reduction in odds ratio" means (1 − OR) × 100% — deviation from OR=1 (no effect). OR=0.75 → 25% reduction vs reference. Do NOT interpret as (unadjusted_OR − adjusted_OR) / unadjusted_OR; that's almost always ≈0% because adjustment barely moves a well-specified OR.

F-statistic vs p-value

scipy.stats.f_oneway(g1, g2, ...) returns (F, p). If GT looks like (0.76, 0.78) and you computed F=91.6, you returned the F-statistic when the question asked for the p-value (or vice-versa). Always re-read the question for which the answer expects.

ANOVA on expression levels across groups — aggregation matters

When asked for "F-statistic comparing expression levels across cell types/groups":

• Each gene/miRNA is one observation. For N genes across K groups, you have N values per group (mean or median expression of that gene across samples in that group)
• Run f_oneway(group1_values, group2_values, ...) where each group has N gene-level values
• Do NOT sum all genes per sample — that gives total RNA content, a completely different quantity
• The F-stat is typically LOW (0.5–2.0) when most genes don't differ across groups, not HIGH (50+)
• If your F-statistic is >10 but the biological context suggests "no significant difference", you probably aggregated to sample level instead of gene level

Log2 fold change between groups — per-gene, then summarize

When asked for "median log2 fold change between group A and group B":

• Compute log2FC per gene: for each gene, log2(mean_expr_A / mean_expr_B)
• Then take the median across all genes
• Do NOT sum all genes per sample first — that gives a single total-expression ratio, not the median per-gene fold change
• A median log2FC near 0 means most genes have similar expression between groups (expected when groups are similar cell types)

Bundled script for both ANOVA and fold change:

# Per-gene ANOVA across cell types
python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py \
  --counts counts.csv --meta meta.csv --group-col cell_type \
  --exclude-groups PBMC --mode anova

# Per-gene log2FC between two groups
python skills/tooluniverse-statistical-modeling/scripts/expression_anova.py \
  --counts counts.csv --meta meta.csv --group-col cell_type \
  --group-a CD14 --group-b CD19 --mode fold_change

Natural spline regression on strain co-culture data

When fitting models on strain co-culture frequency data:

1. Ratio → frequency conversion: If the Ratio column contains strings like "10:1", convert to a frequency fraction: first / (first + second) = 0.909. Report as a fraction in [0, 1], not as ratio notation.

2. Pure-strain endpoints: Whether to include pure-strain data (freq=0 and freq=1) depends on the model type:

• Cubic/polynomial models: Fit on co-culture rows ONLY (exclude pure strains). The cubic model captures the mixed-population response curve; pure strains are fundamentally different biological regimes and including them typically lowers R².
• Natural spline models (`ns(freq, df=4)`): Include the pure-strain endpoint for the focal strain (the one whose frequency the model predicts) but exclude the non-focal pure strain. For example, when modeling "frequency of ΔrhlI to total population", include pure ΔrhlI (freq=1.0) but exclude pure ΔlasI (freq=0.0). This anchors the spline at the high end where the focal strain dominates.

3. R vs Python splines: R's ns() (from the splines package) and Python's patsy.cr() or scipy BSpline produce DIFFERENT knot placements and boundary conditions. If the question references R's lm(y ~ ns(x, df=4)), use the bundled wrapper which runs R via Rscript:

   python skills/tooluniverse-statistical-modeling/scripts/r_natural_spline_regression.py \
     --csv data.csv --y-col Area --x-col Frequency \
     --df 4 --workdir /tmp/spline_run

That emits R², adjusted R², overall F-test p-value, the coefficient table, and the prediction grid with peak (x, y) and 95% confidence interval at the peak.

Do NOT substitute Python's patsy.dmatrix("cr(x, df=4)") — it will give different R², F-statistics, and predictions.

4. Prediction at peak: The bundled R wrapper already predicts on a fine 1000-point grid and reports which.max(pred[,"fit"]) plus its CI from interval="confidence". If you must do this by hand, follow the same recipe; do not use a coarse grid.

5. Best of {quadratic, cubic, natural-spline}: When the question asks for the "best fitting model among quadratic, cubic and natural spline" or the "maximum colony area at the optimal frequency", use the comparison wrapper which fits all three in R and ranks them:

   python skills/tooluniverse-statistical-modeling/scripts/spline_model_compare.py \
     --csv data.csv --y-col Area --x-col Frequency \
     --ns-df 4 --workdir /tmp/spline_cmp

The output's BEST_BY_ADJ_R2 row tells you which model wins, and BEST_PEAK_Y is the peak y to report.

Support

For detailed examples and troubleshooting:

• Logistic regression: references/logistic_regression.md
• Ordinal models: references/ordinal_logistic.md
• Survival analysis: references/cox_regression.md
• Linear/mixed models: references/linear_models.md
• Common patterns: references/common_patterns.md
• ANOVA and tests: anova_and_tests.md
• Diagnostics: references/troubleshooting.md

ANOVA and Statistical Tests Reference

One-way ANOVA (Single Feature)

from scipy import stats

group1 = df[df['celltype'] == 'CD4']['expression']
group2 = df[df['celltype'] == 'CD8']['expression']
group3 = df[df['celltype'] == 'CD14']['expression']

f_stat, p_value = stats.f_oneway(group1, group2, group3)
print(f"F-statistic: {f_stat:.4f}, p-value: {p_value:.6f}")

Multi-Feature ANOVA Decision Tree

When data has multiple features (genes, miRNAs, metabolites), there are TWO approaches:

Question: "What is the F-statistic comparing [feature] expression across groups?"

DECISION TREE:
|
+-- Does question specify "the F-statistic" (singular)?
|   |
|   +-- YES, singular -> Likely asking for SPECIFIC FEATURE(S) F-statistic
|   |   |
|   |   +-- Are there thousands of features (genes, miRNAs)?
|   |   |   YES -> Per-feature approach (Method B below)
|   |   |
|   |   +-- Is there one feature of interest?
|   |      YES -> Single feature ANOVA (Method A below)
|   |
|   +-- NO, asks about "all features" or "genes" (plural)?
|      YES -> Aggregate approach or per-feature summary
|
+-- When unsure: Calculate PER-FEATURE and report summary statistics

Method A: Aggregate ANOVA (all features combined)

Use when: Testing overall expression differences across all features. Result: Single F-statistic representing global effect.

groups_agg = []
for celltype in ['CD4', 'CD8', 'CD14']:
    samples = df[df['celltype'] == celltype]
    all_values = expression_matrix.loc[:, samples.index].values.flatten()
    groups_agg.append(all_values)

f_stat_agg, p_value = stats.f_oneway(*groups_agg)
print(f"Aggregate F-statistic: {f_stat_agg:.4f}")
# Result: Very large F-statistic (e.g., 153.8)

Method B: Per-Feature ANOVA (RECOMMENDED for gene expression)

Use when: Testing EACH feature individually (most common in genomics). Result: Distribution of F-statistics (one per feature).

import numpy as np
from scipy import stats

per_feature_f_stats = []

for feature in expression_matrix.index:
    groups = []
    for celltype in ['CD4', 'CD8', 'CD14']:
        samples = df[df['celltype'] == celltype]
        values = expression_matrix.loc[feature, samples.index].values
        groups.append(values)

    f_stat, _ = stats.f_oneway(*groups)
    if not np.isnan(f_stat):
        per_feature_f_stats.append((feature, f_stat))

# Summary statistics
f_values = [f for _, f in per_feature_f_stats]
print(f"Per-feature F-statistics:")
print(f"  Median: {np.median(f_values):.4f}")
print(f"  Mean: {np.mean(f_values):.4f}")
print(f"  Range: [{np.min(f_values):.4f}, {np.max(f_values):.4f}]")

# Find features in specific range (e.g., 0.76-0.78)
target_features = [(name, f) for name, f in per_feature_f_stats
                   if 0.76 <= f <= 0.78]
if target_features:
    print(f"Features with F in [0.76, 0.78]: {len(target_features)}")
    for name, f_val in target_features:
        print(f"  {name}: F = {f_val:.6f}")

Key Differences

Real-World Example

• Question: "What is the F-statistic comparing miRNA expression across immune cell types?"
• Method A (aggregate ANOVA): Very large F-statistic (e.g., 153.8) -- WRONG
• Method B (per-miRNA ANOVA): Individual F-statistics per gene -- CORRECT

Default assumption for gene expression data: Use Method B (per-feature).

Other Statistical Tests

t-test (two groups)

t_stat, p_value = stats.ttest_ind(group1, group2)

Chi-square (categorical)

contingency = pd.crosstab(df['exposure'], df['outcome'])
chi2, p_value, dof, expected = stats.chi2_contingency(contingency)

Mann-Whitney U (non-parametric, two groups)

u_stat, p_value = stats.mannwhitneyu(group1, group2, alternative='two-sided')

Kruskal-Wallis (non-parametric, 3+ groups)

h_stat, p_value = stats.kruskal(group1, group2, group3)

Examples: Statistical Modeling Skill

Example 1: Ordinal Logistic Regression - BCG Vaccination and COVID-19 Severity

Pattern: "What is the odds ratio of COVID-19 severity associated with BCG vaccination in ordinal logistic regression?"

Setup

import pandas as pd
import numpy as np
from statsmodels.miscmodels.ordinal_model import OrderedModel

# Example dataset: COVID-19 patients with severity and vaccination status
data = {
    'severity': ['Mild']*40 + ['Moderate']*35 + ['Severe']*25,
    'bcg_vaccinated': [1]*25 + [0]*15 + [1]*10 + [0]*25 + [1]*5 + [0]*20,
    'age': np.random.normal(55, 15, 100).astype(int),
    'male': np.random.binomial(1, 0.55, 100),
}
df = pd.DataFrame(data)

Analysis

# Step 1: Define ordinal outcome
severity_order = ['Mild', 'Moderate', 'Severe']
df['severity'] = pd.Categorical(df['severity'], categories=severity_order, ordered=True)
y = df['severity'].cat.codes

# Step 2: Prepare predictors
X = df[['bcg_vaccinated', 'age', 'male']].astype(float)

# Step 3: Fit ordinal logistic regression
model = OrderedModel(y, X, distr='logit')
fit = model.fit(method='bfgs', disp=0)

# Step 4: Extract odds ratio
bcg_coef = fit.params['bcg_vaccinated']
bcg_or = np.exp(bcg_coef)
bcg_ci = np.exp(fit.conf_int().loc['bcg_vaccinated'])

print(f"BCG Vaccination Odds Ratio: {bcg_or:.4f}")
print(f"95% CI: ({bcg_ci.iloc[0]:.4f}, {bcg_ci.iloc[1]:.4f})")
print(f"P-value: {fit.pvalues['bcg_vaccinated']:.6f}")

Interpretation

• OR < 1: BCG vaccination is associated with lower odds of being in a higher severity category
• OR > 1: BCG vaccination is associated with higher odds of being in a higher severity category
• The proportional odds assumption means this OR applies uniformly across all severity cut points

---

Example 2: Binary Logistic Regression - Treatment Response

Pattern: "What is the odds ratio of treatment response associated with biomarker positivity?"

import statsmodels.formula.api as smf
import numpy as np

# Fit logistic regression
model = smf.logit('response ~ biomarker_positive + age + stage', data=df).fit(disp=0)

# Odds ratios with confidence intervals
or_table = np.exp(model.params)
ci = np.exp(model.conf_int())
ci.columns = ['OR_lower', 'OR_upper']
results = pd.DataFrame({
    'OR': or_table,
    'CI_lower': ci['OR_lower'],
    'CI_upper': ci['OR_upper'],
    'p_value': model.pvalues
})
print(results.round(4))

---

Example 3: Percentage Reduction in Odds Ratio (Confounding Assessment)

Pattern: "What is the percentage reduction in odds ratio for higher severity after adjusting for confounders?"

import statsmodels.formula.api as smf
import numpy as np

# Unadjusted (crude) model
model_crude = smf.logit('outcome ~ exposure', data=df).fit(disp=0)
or_crude = np.exp(model_crude.params['exposure'])

# Adjusted model (with confounders)
model_adj = smf.logit('outcome ~ exposure + age + sex + comorbidity', data=df).fit(disp=0)
or_adj = np.exp(model_adj.params['exposure'])

# Calculate percentage reduction
pct_reduction = ((or_crude - or_adj) / or_crude) * 100

print(f"Crude OR: {or_crude:.4f}")
print(f"Adjusted OR: {or_adj:.4f}")
print(f"Percentage reduction in OR: {pct_reduction:.1f}%")
print(f"Interpretation: Adjusting for confounders reduced the OR by {pct_reduction:.1f}%,")
print(f"suggesting {'substantial' if abs(pct_reduction) > 10 else 'minimal'} confounding.")

---

Example 4: Interaction Effect in Ordered Logit

Pattern: "What is the odds ratio associated with patient interaction using ordered logit model?"

import pandas as pd
import numpy as np
from statsmodels.miscmodels.ordinal_model import OrderedModel

# When "patient interaction" refers to interaction between patient-level variables:
# e.g., interaction between treatment and comorbidity

# Prepare data with interaction term
df['treatment_x_comorbidity'] = df['treatment'] * df['comorbidity']

# Predictors including interaction
X = df[['treatment', 'comorbidity', 'treatment_x_comorbidity', 'age']].astype(float)

# Ordinal outcome
y = df['severity'].cat.codes

model = OrderedModel(y, X, distr='logit')
fit = model.fit(method='bfgs', disp=0)

# Interaction OR
interaction_or = np.exp(fit.params['treatment_x_comorbidity'])
print(f"Interaction OR: {interaction_or:.4f}")
print(f"P-value: {fit.pvalues['treatment_x_comorbidity']:.6f}")

---

Example 5: Cox Proportional Hazards Survival Analysis

Pattern: "What is the hazard ratio for drug treatment in a Cox regression model?"

import pandas as pd
from lifelines import CoxPHFitter

# Load survival data
df = pd.read_csv('survival_data.csv')
# Columns: time, event, treatment, age, stage, biomarker

# Fit Cox model
cph = CoxPHFitter()
cph.fit(df[['time', 'event', 'treatment', 'age', 'stage']],
        duration_col='time', event_col='event')

# Extract results
summary = cph.summary
for covar in summary.index:
    hr = summary.loc[covar, 'exp(coef)']
    ci_low = summary.loc[covar, 'exp(coef) lower 95%']
    ci_up = summary.loc[covar, 'exp(coef) upper 95%']
    p = summary.loc[covar, 'p']
    print(f"{covar}: HR={hr:.4f} (95% CI: {ci_low:.4f}-{ci_up:.4f}), p={p:.6f}")

print(f"\nConcordance index: {cph.concordance_index_:.4f}")

# Check proportional hazards assumption
cph.check_assumptions(df[['time', 'event', 'treatment', 'age', 'stage']],
                      p_value_threshold=0.05, show_plots=False)

---

Example 6: Kaplan-Meier with Group Comparison

Pattern: "What is the median survival time for treatment vs control group?"

import pandas as pd
from lifelines import KaplanMeierFitter
from lifelines.statistics import logrank_test

df = pd.read_csv('survival_data.csv')

kmf = KaplanMeierFitter()

# Fit for each group
for group_name in ['Treatment', 'Control']:
    mask = df['group'] == group_name
    kmf.fit(df.loc[mask, 'time'], df.loc[mask, 'event'], label=group_name)

    median_s = kmf.median_survival_time_
    s_12mo = kmf.predict(12) if 12 <= df['time'].max() else None

    print(f"{group_name}:")
    print(f"  Median survival: {median_s:.1f} months")
    if s_12mo is not None:
        print(f"  12-month survival: {s_12mo:.1%}")

# Log-rank test
g_treat = df['group'] == 'Treatment'
g_ctrl = df['group'] == 'Control'
lr = logrank_test(
    df.loc[g_treat, 'time'], df.loc[g_ctrl, 'time'],
    event_observed_A=df.loc[g_treat, 'event'],
    event_observed_B=df.loc[g_ctrl, 'event']
)
print(f"\nLog-rank test: chi2={lr.test_statistic:.4f}, p={lr.p_value:.6f}")

---

Example 7: Linear Mixed-Effects Model

Pattern: "What is the treatment effect in a mixed-effects model with random intercepts per site?"

import statsmodels.formula.api as smf

# Longitudinal data with repeated measures
model = smf.mixedlm('outcome ~ treatment + time + treatment:time',
                     data=df, groups=df['site_id'])
fit = model.fit(reml=True)

# Fixed effects
print("Fixed Effects:")
for name in fit.fe_params.index:
    coef = fit.fe_params[name]
    p = fit.pvalues[name]
    print(f"  {name}: {coef:.4f} (p={p:.6f})")

# Random effects
group_var = float(fit.cov_re.iloc[0, 0])
resid_var = float(fit.scale)
icc = group_var / (group_var + resid_var)
print(f"\nRandom intercept variance: {group_var:.4f}")
print(f"Residual variance: {resid_var:.4f}")
print(f"ICC: {icc:.4f}")

---

Example 8: Complete Analysis Pipeline

Pattern: Full statistical analysis from data loading to reporting.

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.miscmodels.ordinal_model import OrderedModel
from scipy import stats as scipy_stats

# Step 1: Load and inspect data
df = pd.read_csv('study_data.csv')
print(f"Dataset: {df.shape[0]} observations, {df.shape[1]} variables")
print(f"Outcome: {df['outcome'].value_counts().to_dict()}")

# Step 2: Determine outcome type and select model
outcome_type = 'ordinal'  # Based on inspection: Mild < Moderate < Severe
levels = ['Mild', 'Moderate', 'Severe']

# Step 3: Fit appropriate model
df['outcome_cat'] = pd.Categorical(df['outcome'], categories=levels, ordered=True)
y = df['outcome_cat'].cat.codes
X = df[['treatment', 'age', 'sex']].astype(float)

model = OrderedModel(y, X, distr='logit')
fit = model.fit(method='bfgs', disp=0)

# Step 4: Extract key results
treatment_coef = fit.params['treatment']
treatment_or = np.exp(treatment_coef)
treatment_p = fit.pvalues['treatment']
treatment_ci = np.exp(fit.conf_int().loc['treatment'])

# Step 5: Report
print(f"\n=== RESULTS ===")
print(f"Model: Ordinal Logistic Regression (Proportional Odds)")
print(f"Outcome: {' < '.join(levels)}")
print(f"N = {len(y)}")
print(f"\nTreatment effect:")
print(f"  Odds Ratio: {treatment_or:.4f}")
print(f"  95% CI: ({treatment_ci.iloc[0]:.4f}, {treatment_ci.iloc[1]:.4f})")
print(f"  P-value: {treatment_p:.6f}")
print(f"  Significant: {'Yes' if treatment_p < 0.05 else 'No'}")

if treatment_or > 1:
    print(f"  Interpretation: Treatment increases odds of higher severity by {(treatment_or-1)*100:.1f}%")
else:
    print(f"  Interpretation: Treatment decreases odds of higher severity by {(1-treatment_or)*100:.1f}%")

Quick Start: Statistical Modeling Skill

Example 1: Binary Logistic Regression - Odds Ratios

Question: "What is the odds ratio of disease associated with exposure, adjusting for age and sex?"

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

# Load data
df = pd.read_csv('clinical_data.csv')

# Fit logistic regression
model = smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0)

# Extract odds ratios
odds_ratios = np.exp(model.params)
conf_int = np.exp(model.conf_int())

print(f"Odds Ratio for exposure: {odds_ratios['exposure']:.4f}")
print(f"95% CI: ({conf_int.loc['exposure', 0]:.4f}, {conf_int.loc['exposure', 1]:.4f})")
print(f"P-value: {model.pvalues['exposure']:.6f}")

Example 2: Ordinal Logistic Regression

Question: "What is the odds ratio of COVID-19 severity associated with BCG vaccination?"

import pandas as pd
import numpy as np
from statsmodels.miscmodels.ordinal_model import OrderedModel

# Load data with ordinal outcome
df = pd.read_csv('covid_data.csv')

# Set up ordinal outcome
severity_order = ['Mild', 'Moderate', 'Severe']
df['severity'] = pd.Categorical(df['severity'], categories=severity_order, ordered=True)
y = df['severity'].cat.codes

# Predictors
X = df[['bcg_vaccination', 'age', 'sex']].copy()
X = pd.get_dummies(X, drop_first=True, dtype=float)

# Fit ordered logit
model = OrderedModel(y, X, distr='logit')
fit = model.fit(method='bfgs', disp=0)

# Extract odds ratio for BCG vaccination
bcg_coef = fit.params['bcg_vaccination']
bcg_or = np.exp(bcg_coef)
print(f"Odds Ratio (BCG): {bcg_or:.4f}")
print(f"P-value: {fit.pvalues['bcg_vaccination']:.6f}")

Example 3: Cox Proportional Hazards - Hazard Ratios

Question: "What is the hazard ratio for treatment arm in a Cox model?"

import pandas as pd
from lifelines import CoxPHFitter

# Load survival data
df = pd.read_csv('survival_data.csv')

# Fit Cox PH model
cph = CoxPHFitter()
cph.fit(df[['time', 'event', 'treatment', 'age', 'stage']],
        duration_col='time', event_col='event')

# Print summary with hazard ratios
cph.print_summary()
print(f"\nHR for treatment: {cph.hazard_ratios_['treatment']:.4f}")
print(f"Concordance index: {cph.concordance_index_:.4f}")

Example 4: Kaplan-Meier with Log-Rank Test

Question: "Is there a significant survival difference between treatment groups?"

import pandas as pd
from lifelines import KaplanMeierFitter
from lifelines.statistics import logrank_test

df = pd.read_csv('survival_data.csv')

# Fit KM for each group
kmf = KaplanMeierFitter()
for group in df['treatment'].unique():
    mask = df['treatment'] == group
    kmf.fit(df.loc[mask, 'time'], df.loc[mask, 'event'], label=group)
    print(f"Group {group}: median survival = {kmf.median_survival_time_:.1f}")

# Log-rank test
g1 = df['treatment'] == 'Control'
g2 = df['treatment'] == 'Treatment'
result = logrank_test(
    df.loc[g1, 'time'], df.loc[g2, 'time'],
    event_observed_A=df.loc[g1, 'event'],
    event_observed_B=df.loc[g2, 'event']
)
print(f"Log-rank p-value: {result.p_value:.6f}")

Example 5: Mixed-Effects Model

Question: "What is the treatment effect accounting for repeated measures per patient?"

import pandas as pd
import statsmodels.formula.api as smf

df = pd.read_csv('longitudinal_data.csv')

# Fit linear mixed model with random intercepts for patient
model = smf.mixedlm('outcome ~ treatment + time + treatment:time',
                     data=df, groups=df['patient_id'])
fit = model.fit(reml=True)
print(fit.summary())

# ICC
group_var = float(fit.cov_re.iloc[0, 0])
resid_var = float(fit.scale)
icc = group_var / (group_var + resid_var)
print(f"ICC: {icc:.4f}")

Example 6: Percentage Reduction in Odds Ratio (Confounding)

Question: "What is the percentage reduction in odds ratio for severity after adjusting for age?"

import statsmodels.formula.api as smf
import numpy as np

# Unadjusted model
model_crude = smf.logit('outcome ~ exposure', data=df).fit(disp=0)
or_crude = np.exp(model_crude.params['exposure'])

# Adjusted model
model_adj = smf.logit('outcome ~ exposure + age + sex', data=df).fit(disp=0)
or_adj = np.exp(model_adj.params['exposure'])

# Percentage reduction
pct_reduction = (or_crude - or_adj) / or_crude * 100
print(f"Crude OR: {or_crude:.4f}")
print(f"Adjusted OR: {or_adj:.4f}")
print(f"Percentage reduction: {pct_reduction:.1f}%")

Example 7: Interaction Terms

Question: "What is the interaction effect between treatment and biomarker?"

import statsmodels.formula.api as smf
import numpy as np

# Model with interaction
model = smf.logit('outcome ~ treatment * biomarker + age', data=df).fit(disp=0)

# The interaction term
interaction_coef = model.params['treatment:biomarker']
interaction_or = np.exp(interaction_coef)
interaction_p = model.pvalues['treatment:biomarker']

print(f"Interaction OR: {interaction_or:.4f}")
print(f"P-value: {interaction_p:.6f}")

Example 8: Model Comparison

Question: "Which model fits the data better?"

import statsmodels.formula.api as smf

# Fit multiple models
m1 = smf.logit('outcome ~ exposure', data=df).fit(disp=0)
m2 = smf.logit('outcome ~ exposure + age', data=df).fit(disp=0)
m3 = smf.logit('outcome ~ exposure + age + sex + bmi', data=df).fit(disp=0)

# Compare
print(f"Model 1: AIC={m1.aic:.1f}, BIC={m1.bic:.1f}")
print(f"Model 2: AIC={m2.aic:.1f}, BIC={m2.bic:.1f}")
print(f"Model 3: AIC={m3.aic:.1f}, BIC={m3.bic:.1f}")

# Likelihood ratio test (m1 vs m2)
from scipy import stats
lr_stat = -2 * (m1.llf - m2.llf)
p_value = stats.chi2.sf(lr_stat, m2.df_model - m1.df_model)
print(f"LR test p-value (m1 vs m2): {p_value:.6f}")

Statistical Question Patterns

Common statistical question patterns with solutions.

Pattern 1: Odds Ratio from Binary Logistic Regression

Question format: "What is the odds ratio of [outcome] associated with [exposure]?"

Example: "What is the odds ratio of disease associated with exposure to chemical X?"

Solution:

import statsmodels.formula.api as smf
import numpy as np

# Fit logistic regression
model = smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0)

# Extract odds ratio
or_exposure = np.exp(model.params['exposure'])
ci = np.exp(model.conf_int())
ci_lower = ci.loc['exposure', 0]
ci_upper = ci.loc['exposure', 1]
p_val = model.pvalues['exposure']

# Answer
print(f"Odds Ratio: {or_exposure:.4f}")
print(f"95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
print(f"p-value: {p_val:.6f}")

---

Pattern 2: Odds Ratio from Ordinal Logistic Regression

Question format: "What is the odds ratio of [ordinal outcome] associated with [exposure] in ordinal logistic regression?"

Example: "What is the odds ratio of COVID-19 severity associated with BCG vaccination?"

Solution:

from statsmodels.miscmodels.ordinal_model import OrderedModel
import pandas as pd
import numpy as np

# Set up ordered outcome
severity_order = ['Mild', 'Moderate', 'Severe']
df['severity'] = pd.Categorical(df['severity'], categories=severity_order, ordered=True)
y = df['severity'].cat.codes

# Prepare predictors
X = pd.get_dummies(df[['bcg_vaccination', 'age', 'sex']], drop_first=True, dtype=float)

# Fit ordinal logit
model = OrderedModel(y, X, distr='logit').fit(method='bfgs', disp=0)

# Extract OR for BCG vaccination (first predictor)
or_bcg = np.exp(model.params[0])
ci = np.exp(model.conf_int())
ci_lower = ci.iloc[0, 0]
ci_upper = ci.iloc[0, 1]
p_val = model.pvalues[0]

# Answer
print(f"Odds Ratio (BCG): {or_bcg:.4f}")
print(f"95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
print(f"p-value: {p_val:.6f}")

---

Pattern 3: Percentage Reduction in Odds Ratio

Question format: "What is the percentage reduction in odds ratio for [outcome] after adjusting for [confounders]?"

Example: "What is the percentage reduction in OR for disease after adjusting for age and sex?"

Solution:

# Unadjusted model
model_crude = smf.logit('disease ~ exposure', data=df).fit(disp=0)
or_crude = np.exp(model_crude.params['exposure'])

# Adjusted model
model_adj = smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0)
or_adj = np.exp(model_adj.params['exposure'])

# Percentage reduction
pct_reduction = (or_crude - or_adj) / or_crude * 100

# Answer
print(f"Crude OR: {or_crude:.4f}")
print(f"Adjusted OR: {or_adj:.4f}")
print(f"Percentage reduction: {pct_reduction:.1f}%")

---

Pattern 4: Hazard Ratio from Cox Regression

Question format: "What is the hazard ratio for [exposure] in a Cox proportional hazards model?"

Example: "What is the hazard ratio for treatment in a Cox model adjusting for age and stage?"

Solution:

from lifelines import CoxPHFitter

# Fit Cox model
cph = CoxPHFitter()
cph.fit(df[['time', 'event', 'treatment', 'age', 'stage']],
        duration_col='time', event_col='event')

# Extract HR
hr_treatment = cph.hazard_ratios_['treatment']
summary = cph.summary
ci_lower = summary.loc['treatment', 'exp(coef) lower 95%']
ci_upper = summary.loc['treatment', 'exp(coef) upper 95%']
p_val = summary.loc['treatment', 'p']

# Answer
print(f"Hazard Ratio: {hr_treatment:.4f}")
print(f"95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
print(f"p-value: {p_val:.6f}")
print(f"Concordance: {cph.concordance_index_:.4f}")

---

Pattern 5: Kaplan-Meier Survival Estimate

Question format: "What is the Kaplan-Meier survival estimate at time T?"

Example: "What is the 5-year survival probability in the treatment group?"

Solution:

from lifelines import KaplanMeierFitter

# Subset treatment group
treatment_df = df[df['treatment'] == 1]

# Fit KM
kmf = KaplanMeierFitter()
kmf.fit(treatment_df['time'], treatment_df['event'])

# Survival at 5 years (60 months)
survival_5yr = kmf.predict(60)

# Answer
print(f"5-year survival probability: {survival_5yr:.4f}")
print(f"Median survival time: {kmf.median_survival_time_:.1f} months")

---

Pattern 6: Interaction Effect

Question format: "What is the odds ratio associated with the interaction between [A] and [B]?"

Example: "What is the OR for the interaction between treatment and biomarker status?"

Solution:

# Fit model with interaction
model = smf.logit('outcome ~ treatment * biomarker + age', data=df).fit(disp=0)

# Interaction term
interaction_coef = model.params['treatment:biomarker']
interaction_or = np.exp(interaction_coef)
interaction_p = model.pvalues['treatment:biomarker']

# Answer
print(f"Interaction OR: {interaction_or:.4f}")
print(f"p-value: {interaction_p:.6f}")

# Interpretation
if interaction_p < 0.05:
    print("Significant interaction: effect of treatment varies by biomarker status")

---

Pattern 7: Linear Regression Coefficient

Question format: "What is the coefficient for [predictor] in a linear regression model?"

Example: "What is the coefficient for BMI in a linear regression of blood pressure?"

Solution:

import statsmodels.formula.api as smf

# Fit OLS
model = smf.ols('blood_pressure ~ bmi + age + sex', data=df).fit()

# Extract coefficient
coef_bmi = model.params['bmi']
ci = model.conf_int()
ci_lower = ci.loc['bmi', 0]
ci_upper = ci.loc['bmi', 1]
p_val = model.pvalues['bmi']

# Answer
print(f"Coefficient (BMI): {coef_bmi:.4f}")
print(f"95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
print(f"p-value: {p_val:.6f}")
print(f"R-squared: {model.rsquared:.4f}")

# Interpretation
print(f"Interpretation: Each 1-unit increase in BMI is associated with {coef_bmi:.2f} mmHg change in blood pressure")

---

Pattern 8: Mixed-Effects Model Coefficient

Question format: "What is the coefficient for [predictor] in a mixed-effects model with random intercepts for [grouping]?"

Example: "What is the treatment effect in a mixed model with random intercepts for patient?"

Solution:

import statsmodels.formula.api as smf

# Fit LMM
model = smf.mixedlm('outcome ~ treatment + time', data=df, groups=df['patient_id']).fit(reml=True)

# Extract fixed effect
coef_treatment = model.fe_params['treatment']
se = model.bse_fe['treatment']
p_val = model.pvalues['treatment']
ci = model.conf_int()
ci_lower = ci.loc['treatment', 0]
ci_upper = ci.loc['treatment', 1]

# Answer
print(f"Coefficient (treatment): {coef_treatment:.4f}")
print(f"95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
print(f"p-value: {p_val:.6f}")

# ICC
group_var = float(model.cov_re.iloc[0, 0])
resid_var = float(model.scale)
icc = group_var / (group_var + resid_var)
print(f"ICC: {icc:.4f}")

---

Pattern 9: Model Comparison with AIC/BIC

Question format: "Which model has better fit according to AIC?"

Example: "Compare models with and without interaction term using AIC"

Solution:

# Fit models
model1 = smf.logit('outcome ~ A + B', data=df).fit(disp=0)
model2 = smf.logit('outcome ~ A * B', data=df).fit(disp=0)

# Compare
print("Model 1 (no interaction):")
print(f"  AIC: {model1.aic:.2f}")
print(f"  BIC: {model1.bic:.2f}")

print("\nModel 2 (with interaction):")
print(f"  AIC: {model2.aic:.2f}")
print(f"  BIC: {model2.bic:.2f}")

# Answer
if model2.aic < model1.aic:
    print("\nModel 2 preferred (lower AIC)")
else:
    print("\nModel 1 preferred (lower AIC)")

---

Pattern 10: Proportional Odds Assumption Test

Question format: "Is the proportional odds assumption met?"

Example: "Test if proportional odds assumption holds for ordinal severity model"

Solution:

import statsmodels.api as sm

# Fit binary logits at each cutpoint
severity_order = ['Mild', 'Moderate', 'Severe']
df['severity'] = pd.Categorical(df['severity'], categories=severity_order, ordered=True)
y_codes = df['severity'].cat.codes

X = pd.get_dummies(df[['exposure', 'age']], drop_first=True, dtype=float)
X_const = sm.add_constant(X)

# Fit at each cutpoint
coef_by_cutpoint = {}
for k in range(len(severity_order) - 1):
    y_binary = (y_codes > k).astype(int)
    model = sm.Logit(y_binary, X_const).fit(disp=0)
    coef_by_cutpoint[k] = model.params['exposure']

# Check if coefficients similar
coefs = list(coef_by_cutpoint.values())
coef_range = max(coefs) - min(coefs)

print(f"Coefficients by cutpoint: {coefs}")
print(f"Range: {coef_range:.4f}")

# Answer
if coef_range < 0.5:
    print("Proportional odds assumption likely satisfied")
else:
    print("Proportional odds assumption may be violated")

---

Pattern 11: Log-Rank Test

Question format: "Is there a significant difference in survival between groups?"

Example: "Test if survival differs between treatment and control groups"

Solution:

from lifelines.statistics import logrank_test

# Split by group
treatment_group = df['treatment'] == 1
control_group = df['treatment'] == 0

# Log-rank test
result = logrank_test(
    df.loc[treatment_group, 'time'],
    df.loc[control_group, 'time'],
    df.loc[treatment_group, 'event'],
    df.loc[control_group, 'event']
)

# Answer
print(f"Log-rank test statistic: {result.test_statistic:.4f}")
print(f"p-value: {result.p_value:.6f}")

if result.p_value < 0.05:
    print("Survival curves are significantly different")
else:
    print("No significant difference in survival")

---

Pattern 12: R-squared Interpretation

Question format: "What is the R-squared of the model?"

Example: "What proportion of variance is explained by the linear model?"

Solution:

model = smf.ols('outcome ~ predictor1 + predictor2 + age', data=df).fit()

# R-squared
r2 = model.rsquared
adj_r2 = model.rsquared_adj

# Answer
print(f"R-squared: {r2:.4f}")
print(f"Adjusted R-squared: {adj_r2:.4f}")
print(f"Interpretation: {r2*100:.1f}% of variance in outcome is explained by predictors")

---

Pattern 13: Concordance Index Interpretation

Question format: "What is the concordance index of the Cox model?"

Example: "How well does the Cox model discriminate between patients?"

Solution:

# Already fitted cph model
c_index = cph.concordance_index_

# Answer
print(f"Concordance index: {c_index:.4f}")

# Interpretation
if c_index > 0.7:
    print("Good discrimination (C > 0.7)")
elif c_index > 0.6:
    print("Acceptable discrimination (0.6 < C < 0.7)")
else:
    print("Poor discrimination (C < 0.6)")

---

Pattern 14: Coefficient Change After Adjustment

Question format: "How does the coefficient change after adjusting for confounders?"

Example: "Compare exposure coefficient before and after adjusting for age"

Solution:

# Unadjusted
model_crude = smf.ols('outcome ~ exposure', data=df).fit()
coef_crude = model_crude.params['exposure']

# Adjusted
model_adj = smf.ols('outcome ~ exposure + age + sex', data=df).fit()
coef_adj = model_adj.params['exposure']

# Change
absolute_change = coef_adj - coef_crude
pct_change = (coef_adj - coef_crude) / coef_crude * 100

# Answer
print(f"Crude coefficient: {coef_crude:.4f}")
print(f"Adjusted coefficient: {coef_adj:.4f}")
print(f"Absolute change: {absolute_change:.4f}")
print(f"Percentage change: {pct_change:.1f}%")

---

Pattern 15: Stratified Analysis

Question format: "What is the odds ratio stratified by [variable]?"

Example: "What is the OR for exposure separately in men and women?"

Solution:

# Stratify by sex
results_by_sex = {}

for sex in ['M', 'F']:
    df_subset = df[df['sex'] == sex]
    model = smf.logit('outcome ~ exposure + age', data=df_subset).fit(disp=0)

    or_exposure = np.exp(model.params['exposure'])
    ci = np.exp(model.conf_int())

    results_by_sex[sex] = {
        'OR': or_exposure,
        'CI_lower': ci.loc['exposure', 0],
        'CI_upper': ci.loc['exposure', 1],
        'p_value': model.pvalues['exposure']
    }

# Answer
for sex, result in results_by_sex.items():
    print(f"\n{sex}:")
    print(f"  OR: {result['OR']:.4f}")
    print(f"  95% CI: ({result['CI_lower']:.4f}, {result['CI_upper']:.4f})")
    print(f"  p-value: {result['p_value']:.6f}")

---

Quick Reference Table

---

Common Mistakes to Avoid

1. Forgetting to exponentiate: Logistic/Cox coefficients are log-odds/log-hazards. Must use np.exp() for ORs/HRs.

2. Wrong variable type: Ordinal outcomes need OrderedModel, not Logit.

3. Missing confounders: Always check if question specifies "adjusting for" variables.

4. Interpretation direction: HR > 1 = worse outcome (higher hazard), OR > 1 = higher odds.

5. Precision: Round to requested decimal places (typically 4 for ORs/HRs, 6 for p-values).

6. CI extraction: Use conf_int() method, not conf_int attribute.

7. Formula syntax: Interactions use *, not +. E.g., 'y ~ A * B' includes A, B, and A:B.

8. Duration/event cols: Cox models require explicit duration_col and event_col parameters.

Cox Proportional Hazards and Survival Analysis Reference

Complete guide to survival analysis using Cox regression and Kaplan-Meier estimation.

When to Use Survival Analysis

Use survival analysis when:

• Outcome is time-to-event (time until death, disease progression, recovery)
• Censoring present (some participants didn't experience event by study end)
• Want to model hazard (instantaneous rate of event occurrence)

Examples:

• Time to death, progression-free survival, overall survival
• Time to hospital readmission, disease recurrence
• Duration of remission, time to treatment failure

Cox Proportional Hazards Model

Basic Cox Regression

import pandas as pd
import numpy as np
from lifelines import CoxPHFitter

# Load survival data
# Required columns: duration (time), event (1=event, 0=censored), covariates
df = pd.read_csv('survival_data.csv')

# Initialize and fit Cox model
cph = CoxPHFitter()
cph.fit(df[['time', 'event', 'treatment', 'age', 'stage']],
        duration_col='time',
        event_col='event')

# Print summary with hazard ratios
cph.print_summary()

Extracting Hazard Ratios

# Get summary table
summary = cph.summary

# Hazard ratios (HR)
hrs = cph.hazard_ratios_

# Print results
print("\n=== Hazard Ratios ===")
for var in summary.index:
    hr = hrs[var]
    ci_lower = summary.loc[var, 'exp(coef) lower 95%']
    ci_upper = summary.loc[var, 'exp(coef) upper 95%']
    p_val = summary.loc[var, 'p']
    coef = summary.loc[var, 'coef']
    se = summary.loc[var, 'se(coef)']

    print(f"\n{var}:")
    print(f"  HR: {hr:.4f}")
    print(f"  95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
    print(f"  p-value: {p_val:.6f}")
    print(f"  Coefficient: {coef:.4f} (SE: {se:.4f})")

Interpreting Hazard Ratios

HR > 1: Increased hazard (worse prognosis, faster time to event) HR < 1: Decreased hazard (better prognosis, slower time to event) HR = 1: No effect on hazard

def interpret_hazard_ratio(hr, ci_lower, ci_upper, p_val, var_name):
    """Interpret hazard ratio."""
    if hr > 1:
        pct_increase = (hr - 1) * 100
        interp = f"{var_name} is associated with {pct_increase:.1f}% increase in hazard (worse prognosis)"
    elif hr < 1:
        pct_decrease = (1 - hr) * 100
        interp = f"{var_name} is associated with {pct_decrease:.1f}% decrease in hazard (better prognosis)"
    else:
        interp = f"{var_name} has no effect on hazard"

    sig = "statistically significant" if p_val < 0.05 else "not statistically significant"

    return f"{interp} (HR={hr:.4f}, 95% CI [{ci_lower:.4f}, {ci_upper:.4f}], p={p_val:.6f}, {sig})"

# Example
print(interpret_hazard_ratio(0.65, 0.45, 0.92, 0.015, "Treatment"))
# Output: "Treatment is associated with 35.0% decrease in hazard (better prognosis)
#          (HR=0.6500, 95% CI [0.4500, 0.9200], p=0.015000, statistically significant)"

Model Fit Statistics

# Concordance index (C-index)
# Measures discrimination ability (similar to AUC)
# 0.5 = random, 1.0 = perfect
print(f"Concordance index: {cph.concordance_index_:.4f}")

# Partial log-likelihood
print(f"Partial log-likelihood: {cph.log_likelihood_:.4f}")

# AIC
if hasattr(cph, 'AIC_partial_'):
    print(f"AIC: {cph.AIC_partial_:.4f}")

# Number of events
print(f"Events: {cph.event_observed.sum()}/{len(cph.event_observed)}")

Testing Proportional Hazards Assumption

Critical assumption: The hazard ratio is constant over time.

Schoenfeld Residuals Test

# Test PH assumption
results = cph.check_assumptions(df, p_value_threshold=0.05, show_plots=False)

if len(results) == 0:
    print("✅ Proportional hazards assumption met for all covariates")
else:
    print(f"⚠️  Proportional hazards assumption violated for: {results}")

What to Do if Assumption Violated

1. Stratify by problematic variable

# Stratify by treatment (doesn't estimate HR for treatment)
cph_strat = CoxPHFitter()
cph_strat.fit(df, duration_col='time', event_col='event', strata=['treatment'])

2. Time-varying coefficient

# Allow coefficient to change over time (advanced)
# Use cph.fit(..., formula='...') with time interactions

3. Parametric survival models

from lifelines import WeibullAFTFitter

# Accelerated failure time model (doesn't assume PH)
wf = WeibullAFTFitter()
wf.fit(df, duration_col='time', event_col='event')

Kaplan-Meier Survival Curves

Single Group

from lifelines import KaplanMeierFitter

# Fit Kaplan-Meier estimator
kmf = KaplanMeierFitter()
kmf.fit(df['time'], df['event'], label='All patients')

# Median survival time
median_survival = kmf.median_survival_time_
print(f"Median survival: {median_survival:.1f} months")

# Survival probability at specific times
for t in [12, 24, 36, 60]:
    survival_prob = kmf.predict(t)
    print(f"Survival at {t} months: {survival_prob:.4f}")

# Plot survival curve
kmf.plot_survival_function()

Comparing Groups (Log-Rank Test)

from lifelines.statistics import logrank_test

# Split by treatment group
treatment_group = df['treatment'] == 1
control_group = df['treatment'] == 0

# Fit KM for each group
kmf_treatment = KaplanMeierFitter()
kmf_treatment.fit(df.loc[treatment_group, 'time'],
                  df.loc[treatment_group, 'event'],
                  label='Treatment')

kmf_control = KaplanMeierFitter()
kmf_control.fit(df.loc[control_group, 'time'],
                df.loc[control_group, 'event'],
                label='Control')

# Median survival times
print(f"Median survival (treatment): {kmf_treatment.median_survival_time_:.1f}")
print(f"Median survival (control): {kmf_control.median_survival_time_:.1f}")

# Log-rank test
result = logrank_test(
    df.loc[treatment_group, 'time'],
    df.loc[control_group, 'time'],
    df.loc[treatment_group, 'event'],
    df.loc[control_group, 'event']
)

print(f"\nLog-rank test:")
print(f"  Test statistic: {result.test_statistic:.4f}")
print(f"  p-value: {result.p_value:.6f}")
if result.p_value < 0.05:
    print("  Conclusion: Survival curves significantly different")
else:
    print("  Conclusion: No significant difference in survival")

# Plot both curves
import matplotlib.pyplot as plt
kmf_treatment.plot_survival_function()
kmf_control.plot_survival_function()
plt.title('Survival Curves by Treatment Group')
plt.xlabel('Time (months)')
plt.ylabel('Survival Probability')
plt.show()

Advanced Features

Robust Standard Errors (Clustered Data)

# Use cluster_col for repeated measures or matched data
cph = CoxPHFitter()
cph.fit(df, duration_col='time', event_col='event',
        cluster_col='patient_id')  # Adjusts SEs for clustering

Stratified Cox Model

# Don't estimate HR for strata variable, but adjust for it
cph = CoxPHFitter()
cph.fit(df, duration_col='time', event_col='event',
        strata=['center', 'sex'])  # Stratify by center and sex

Weighted Cox Model

# Weight observations (e.g., for propensity score weighting)
cph = CoxPHFitter()
cph.fit(df, duration_col='time', event_col='event',
        weights_col='propensity_weight')

Categorical Variables

Encoding Categorical Predictors

# One-hot encode (drop first level as reference)
df_encoded = pd.get_dummies(df, columns=['stage'], drop_first=True, dtype=int)

# Fit model
cph = CoxPHFitter()
cph.fit(df_encoded[['time', 'event', 'stage_II', 'stage_III', 'stage_IV', 'age']],
        duration_col='time', event_col='event')

# Interpret: HRs are relative to Stage I (reference)
print(f"HR (Stage II vs I): {cph.hazard_ratios_['stage_II']:.4f}")
print(f"HR (Stage III vs I): {cph.hazard_ratios_['stage_III']:.4f}")
print(f"HR (Stage IV vs I): {cph.hazard_ratios_['stage_IV']:.4f}")

Prediction

Individual Risk Scores

# Predict partial hazard (relative risk)
# Higher score = higher risk
risk_scores = cph.predict_partial_hazard(df)
print(f"Risk score range: {risk_scores.min():.4f} to {risk_scores.max():.4f}")

# Add to dataframe
df['risk_score'] = risk_scores

# Identify high-risk patients
high_risk = df[risk_scores > risk_scores.quantile(0.75)]
print(f"High-risk patients (top 25%): {len(high_risk)}")

Survival Curves for Individuals

# Predict survival function for specific patient profiles
new_patient = pd.DataFrame({
    'treatment': [1],
    'age': [55],
    'stage_II': [0],
    'stage_III': [1],
    'stage_IV': [0]
})

# Get survival function
survival_func = cph.predict_survival_function(new_patient)

# Survival probability at specific times
print(f"Survival at 12 months: {survival_func.loc[12].values[0]:.4f}")
print(f"Survival at 24 months: {survival_func.loc[24].values[0]:.4f}")

Model Comparison

Likelihood Ratio Test

# Compare nested models
cph_reduced = CoxPHFitter()
cph_reduced.fit(df[['time', 'event', 'treatment']], duration_col='time', event_col='event')

cph_full = CoxPHFitter()
cph_full.fit(df[['time', 'event', 'treatment', 'age', 'stage']], duration_col='time', event_col='event')

# LR test
from scipy import stats
lr_stat = -2 * (cph_reduced.log_likelihood_ - cph_full.log_likelihood_)
df_diff = cph_full.params_.shape[0] - cph_reduced.params_.shape[0]
p_value = stats.chi2.sf(lr_stat, df_diff)

print(f"LR statistic: {lr_stat:.4f}")
print(f"p-value: {p_value:.6f}")
if p_value < 0.05:
    print("Full model significantly better")

AIC Comparison

print(f"Reduced model AIC: {cph_reduced.AIC_partial_:.2f}")
print(f"Full model AIC: {cph_full.AIC_partial_:.2f}")

if cph_full.AIC_partial_ < cph_reduced.AIC_partial_:
    print("Full model preferred (lower AIC)")

Complete Example: Cancer Clinical Trial

import pandas as pd
import numpy as np
from lifelines import CoxPHFitter, KaplanMeierFitter
from lifelines.statistics import logrank_test

# Simulated cancer trial data
np.random.seed(42)
n = 200

df = pd.DataFrame({
    'patient_id': range(1, n+1),
    'time': np.random.exponential(scale=24, size=n),
    'event': np.random.binomial(1, 0.6, n),
    'treatment': np.random.binomial(1, 0.5, n),
    'age': np.random.normal(60, 10, n),
    'stage': np.random.choice(['I', 'II', 'III', 'IV'], n, p=[0.1, 0.3, 0.4, 0.2])
})

# Encode stage
df = pd.get_dummies(df, columns=['stage'], drop_first=True, dtype=int)

print("=== Cancer Clinical Trial Survival Analysis ===\n")

# 1. Cox Proportional Hazards
print("1. Cox Proportional Hazards Model\n")
cph = CoxPHFitter()
cph.fit(df[['time', 'event', 'treatment', 'age', 'stage_II', 'stage_III', 'stage_IV']],
        duration_col='time', event_col='event')

print(f"Concordance index: {cph.concordance_index_:.4f}\n")

# Hazard ratios
hrs = cph.hazard_ratios_
summary = cph.summary

for var in ['treatment', 'age', 'stage_II', 'stage_III', 'stage_IV']:
    hr = hrs[var]
    ci_lower = summary.loc[var, 'exp(coef) lower 95%']
    ci_upper = summary.loc[var, 'exp(coef) upper 95%']
    p_val = summary.loc[var, 'p']

    print(f"{var}:")
    print(f"  HR: {hr:.4f} (95% CI: {ci_lower:.4f}-{ci_upper:.4f})")
    print(f"  p-value: {p_val:.6f}")
    print()

# 2. Test proportional hazards assumption
print("\n2. Proportional Hazards Assumption Test\n")
ph_results = cph.check_assumptions(df, p_value_threshold=0.05, show_plots=False)
if len(ph_results) == 0:
    print("✅ All covariates meet PH assumption\n")
else:
    print(f"⚠️  PH violated for: {ph_results}\n")

# 3. Kaplan-Meier by treatment group
print("\n3. Kaplan-Meier Analysis by Treatment\n")

treatment_group = df['treatment'] == 1
control_group = df['treatment'] == 0

kmf_tx = KaplanMeierFitter()
kmf_tx.fit(df.loc[treatment_group, 'time'],
           df.loc[treatment_group, 'event'],
           label='Treatment')

kmf_ctrl = KaplanMeierFitter()
kmf_ctrl.fit(df.loc[control_group, 'time'],
             df.loc[control_group, 'event'],
             label='Control')

print(f"Median survival (treatment): {kmf_tx.median_survival_time_:.1f} months")
print(f"Median survival (control): {kmf_ctrl.median_survival_time_:.1f} months\n")

# Log-rank test
lr = logrank_test(df.loc[treatment_group, 'time'],
                  df.loc[control_group, 'time'],
                  df.loc[treatment_group, 'event'],
                  df.loc[control_group, 'event'])

print(f"Log-rank test:")
print(f"  Test statistic: {lr.test_statistic:.4f}")
print(f"  p-value: {lr.p_value:.6f}")

Common Question Pattern

Typical question: "What is the hazard ratio for treatment in a Cox proportional hazards model?"

Solution: 1. Load survival data (time, event, covariates) 2. Fit Cox model with duration_col and event_col 3. Extract HR from cph.hazard_ratios_['treatment'] 4. Report with CI, p-value, and concordance index

# Answer format
hr_treatment = cph.hazard_ratios_['treatment']
ci_lower = cph.summary.loc['treatment', 'exp(coef) lower 95%']
ci_upper = cph.summary.loc['treatment', 'exp(coef) upper 95%']
p_val = cph.summary.loc['treatment', 'p']

print(f"Answer: {hr_treatment:.4f}")
print(f"95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
print(f"p-value: {p_val:.6f}")
print(f"Concordance: {cph.concordance_index_:.4f}")

Reporting Template

def report_cox_model(cph, duration_col='time', event_col='event'):
    """Generate publication-ready Cox model report."""
    summary = cph.summary
    hrs = cph.hazard_ratios_

    report = []
    report.append("=== Cox Proportional Hazards Model ===\n")
    report.append(f"Events: {cph.event_observed.sum()}/{len(cph.event_observed)}")
    report.append(f"Concordance index: {cph.concordance_index_:.4f}")
    if hasattr(cph, 'AIC_partial_'):
        report.append(f"AIC: {cph.AIC_partial_:.2f}\n")

    report.append("Hazard Ratios (95% CI):\n")
    for var in summary.index:
        hr = hrs[var]
        ci_lower = summary.loc[var, 'exp(coef) lower 95%']
        ci_upper = summary.loc[var, 'exp(coef) upper 95%']
        p_val = summary.loc[var, 'p']
        sig = "*" if p_val < 0.05 else ""

        report.append(f"  {var}: HR={hr:.4f} ({ci_lower:.4f}-{ci_upper:.4f}), p={p_val:.6f}{sig}")

    return "\n".join(report)

print(report_cox_model(cph))

Linear Models and Mixed-Effects Reference

Complete guide to linear regression and mixed-effects models.

Ordinary Least Squares (OLS) Regression

Basic Linear Regression

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

# Method 1: Formula API (recommended)
model = smf.ols('outcome ~ predictor1 + predictor2 + age', data=df).fit()

# Method 2: Matrix API
X = sm.add_constant(df[['predictor1', 'predictor2', 'age']])
y = df['outcome']
model = sm.OLS(y, X).fit()

# Print summary
print(model.summary())

Interpreting Coefficients

# Extract coefficients
coefs = model.params
se = model.bse
t_vals = model.tvalues
p_vals = model.pvalues
ci = model.conf_int()

for var in coefs.index:
    print(f"\n{var}:")
    print(f"  Coefficient: {coefs[var]:.4f}")
    print(f"  Std Error: {se[var]:.4f}")
    print(f"  t-value: {t_vals[var]:.4f}")
    print(f"  p-value: {p_vals[var]:.6f}")
    print(f"  95% CI: ({ci.loc[var, 0]:.4f}, {ci.loc[var, 1]:.4f})")

Model Fit Statistics

# R-squared
print(f"R-squared: {model.rsquared:.4f}")
print(f"Adjusted R-squared: {model.rsquared_adj:.4f}")

# F-statistic
print(f"F-statistic: {model.fvalue:.4f}")
print(f"F-test p-value: {model.f_pvalue:.6f}")

# AIC/BIC
print(f"AIC: {model.aic:.2f}")
print(f"BIC: {model.bic:.2f}")

# Root mean squared error
rmse = np.sqrt(model.mse_resid)
print(f"RMSE: {rmse:.4f}")

Diagnostics

Residual Normality

from scipy import stats as scipy_stats

# Shapiro-Wilk test
residuals = model.resid
sw_stat, sw_p = scipy_stats.shapiro(residuals)

print(f"Shapiro-Wilk test:")
print(f"  Statistic: {sw_stat:.4f}")
print(f"  p-value: {sw_p:.6f}")
if sw_p > 0.05:
    print("  ✅ Residuals appear normally distributed")
else:
    print("  ⚠️  Residuals may not be normally distributed")

Homoscedasticity (Equal Variance)

from statsmodels.stats.diagnostic import het_breuschpagan

# Breusch-Pagan test
bp_stat, bp_p, _, _ = het_breuschpagan(residuals, model.model.exog)

print(f"\nBreusch-Pagan test:")
print(f"  Statistic: {bp_stat:.4f}")
print(f"  p-value: {bp_p:.6f}")
if bp_p > 0.05:
    print("  ✅ Homoscedasticity assumption met")
else:
    print("  ⚠️  Heteroscedasticity detected")
    print("  Consider: robust standard errors, log transformation, or WLS")

Autocorrelation

from statsmodels.stats.stattools import durbin_watson

# Durbin-Watson test
dw = durbin_watson(residuals)

print(f"\nDurbin-Watson statistic: {dw:.4f}")
if 1.5 < dw < 2.5:
    print("  ✅ No autocorrelation detected")
else:
    print("  ⚠️  Possible autocorrelation")

Multicollinearity (VIF)

from statsmodels.stats.outliers_influence import variance_inflation_factor

X = model.model.exog
vif_data = []

for i in range(X.shape[1]):
    var_name = model.model.exog_names[i]
    if var_name != 'Intercept':
        vif = variance_inflation_factor(X, i)
        vif_data.append({'variable': var_name, 'VIF': vif})

vif_df = pd.DataFrame(vif_data)
print("\nVariance Inflation Factors:")
print(vif_df)

# Rule of thumb: VIF > 10 indicates multicollinearity
if (vif_df['VIF'] > 10).any():
    print("\n⚠️  High multicollinearity detected (VIF > 10)")
    print("Consider: removing correlated predictors or PCA")

Linear Mixed-Effects Models (LMM)

When to Use LMM

Use LMM when:

• Repeated measures (same subject measured multiple times)
• Nested/clustered data (patients within hospitals)
• Hierarchical structure (students within schools)
• Need to model between-subject and within-subject variation

Random Intercept Model

import statsmodels.formula.api as smf

# Random intercepts for subjects
model = smf.mixedlm('outcome ~ treatment + time',
                     data=df,
                     groups=df['subject_id'])
fit = model.fit(reml=True)

print(fit.summary())

Random Slope Model

# Random slopes for time (subjects have different time trends)
model = smf.mixedlm('outcome ~ treatment + time',
                     data=df,
                     groups=df['subject_id'],
                     re_formula='~time')
fit = model.fit(reml=True)

Fixed Effects Interpretation

# Extract fixed effects
fe_params = fit.fe_params
fe_pvalues = fit.pvalues
fe_ci = fit.conf_int()

print("Fixed Effects:")
for var in fe_params.index:
    coef = fe_params[var]
    p_val = fe_pvalues[var]
    ci_lower = fe_ci.loc[var, 0]
    ci_upper = fe_ci.loc[var, 1]

    print(f"\n{var}:")
    print(f"  Coefficient: {coef:.4f}")
    print(f"  95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
    print(f"  p-value: {p_val:.6f}")

Random Effects Variance

# Group (subject) variance
cov_re = fit.cov_re
if hasattr(cov_re, 'values'):
    group_var = cov_re.iloc[0, 0]
else:
    group_var = float(cov_re)

# Residual variance
resid_var = fit.scale

print(f"\nRandom Effects Variance:")
print(f"  Between-subject (Group): {group_var:.4f}")
print(f"  Within-subject (Residual): {resid_var:.4f}")

Intraclass Correlation Coefficient (ICC)

# ICC: proportion of variance due to grouping
icc = group_var / (group_var + resid_var)

print(f"\nICC: {icc:.4f}")
print(f"Interpretation: {icc*100:.1f}% of variance is between subjects")

if icc > 0.1:
    print("  ✅ Substantial clustering - LMM appropriate")
else:
    print("  ⚠️  Low clustering - OLS may be sufficient")

Complete Example: Longitudinal Study

import pandas as pd
import numpy as np
import statsmodels.formula.api as smf

# Simulated longitudinal data
np.random.seed(42)
n_subjects = 50
n_timepoints = 4

data = []
for subject_id in range(n_subjects):
    treatment = np.random.binomial(1, 0.5)
    baseline = np.random.normal(100, 10)

    for time in range(n_timepoints):
        # Treatment effect + time trend + random noise
        outcome = baseline + treatment * 5 + time * 2 + np.random.normal(0, 3)
        data.append({
            'subject_id': subject_id,
            'treatment': treatment,
            'time': time,
            'outcome': outcome
        })

df = pd.DataFrame(data)

print("=== Longitudinal Analysis ===\n")

# 1. OLS (ignoring repeated measures - WRONG)
print("1. OLS (incorrect for repeated measures):\n")
ols_model = smf.ols('outcome ~ treatment + time', data=df).fit()
print(f"Treatment effect: {ols_model.params['treatment']:.4f}")
print(f"p-value: {ols_model.pvalues['treatment']:.6f}")
print(f"R-squared: {ols_model.rsquared:.4f}\n")

# 2. LMM with random intercepts (correct)
print("2. LMM with random intercepts (correct):\n")
lmm_model = smf.mixedlm('outcome ~ treatment + time',
                         data=df,
                         groups=df['subject_id'])
lmm_fit = lmm_model.fit(reml=True)

print(f"Treatment effect: {lmm_fit.fe_params['treatment']:.4f}")
print(f"p-value: {lmm_fit.pvalues['treatment']:.6f}")

# ICC
group_var = float(lmm_fit.cov_re.iloc[0, 0])
resid_var = float(lmm_fit.scale)
icc = group_var / (group_var + resid_var)
print(f"ICC: {icc:.4f}\n")

# 3. LMM with random slopes (allow different time trends)
print("3. LMM with random slopes:\n")
lmm_slopes = smf.mixedlm('outcome ~ treatment + time',
                          data=df,
                          groups=df['subject_id'],
                          re_formula='~time')
lmm_slopes_fit = lmm_slopes.fit(reml=True)

print(f"Treatment effect: {lmm_slopes_fit.fe_params['treatment']:.4f}")
print(f"p-value: {lmm_slopes_fit.pvalues['treatment']:.6f}\n")

# Model comparison
print("4. Model Comparison:")
print(f"Random intercepts AIC: {lmm_fit.aic:.2f}")
print(f"Random slopes AIC: {lmm_slopes_fit.aic:.2f}")

if lmm_slopes_fit.aic < lmm_fit.aic:
    print("Random slopes model preferred (lower AIC)")
else:
    print("Random intercepts model sufficient")

Weighted Least Squares (WLS)

When to Use WLS

Use WLS when:

• Heteroscedasticity detected
• Known different precision across observations
• Want to downweight outliers

# Estimate weights from variance model
# Common approach: inverse variance weighting
residuals_sq = model.resid ** 2
weights = 1 / residuals_sq

# Fit WLS
wls_model = sm.WLS(y, X, weights=weights).fit()
print(wls_model.summary())

Robust Standard Errors

Heteroscedasticity-Consistent Standard Errors

# Use robust standard errors (HC3)
model_robust = model.get_robustcov_results(cov_type='HC3')

print("Robust Standard Errors:")
print(model_robust.summary())

Generalized Estimating Equations (GEE)

Alternative to LMM for Clustered Data

import statsmodels.api as sm
from statsmodels.genmod.generalized_estimating_equations import GEE
from statsmodels.genmod.families import Gaussian
from statsmodels.genmod.cov_struct import Exchangeable

# GEE for continuous outcome
gee_model = GEE.from_formula('outcome ~ treatment + time',
                              groups='subject_id',
                              data=df,
                              family=Gaussian(),
                              cov_struct=Exchangeable())
gee_fit = gee_model.fit()

print(gee_fit.summary())

GEE vs LMM:

• GEE: Population-averaged effects, requires weaker assumptions
• LMM: Subject-specific effects, models random effects explicitly
• Use GEE when interested in marginal effects, LMM when interested in individual trajectories

Polynomial and Spline Models

Polynomial Regression

# Add polynomial terms
df['age_sq'] = df['age'] ** 2
df['age_cube'] = df['age'] ** 3

model = smf.ols('outcome ~ age + age_sq + age_cube', data=df).fit()

Natural Cubic Splines

from patsy import dmatrix

# Create spline basis
spline_basis = dmatrix("bs(age, df=4, degree=3)", df, return_type='dataframe')
X_spline = pd.concat([spline_basis, df[['treatment']]], axis=1)
y = df['outcome']

model_spline = sm.OLS(y, X_spline).fit()

Reporting Template

def report_linear_model(model):
    """Generate publication-ready linear model report."""
    report = []
    report.append("=== Linear Regression Results ===\n")
    report.append(f"N = {int(model.nobs)}")
    report.append(f"R² = {model.rsquared:.4f}")
    report.append(f"Adjusted R² = {model.rsquared_adj:.4f}")
    report.append(f"F({model.df_model:.0f}, {model.df_resid:.0f}) = {model.fvalue:.4f}, p = {model.f_pvalue:.6f}\n")

    report.append("Coefficients:\n")
    report.append("Variable | Coef | SE | t | p | 95% CI")
    report.append("---------|------|----|----|---|-------")

    for var in model.params.index:
        coef = model.params[var]
        se = model.bse[var]
        t_val = model.tvalues[var]
        p_val = model.pvalues[var]
        ci = model.conf_int().loc[var]
        sig = "*" if p_val < 0.05 else ""

        report.append(f"{var} | {coef:.4f} | {se:.4f} | {t_val:.4f} | {p_val:.6f}{sig} | ({ci[0]:.4f}, {ci[1]:.4f})")

    return "\n".join(report)

print(report_linear_model(model))

Logistic Regression Reference

Complete guide to binary logistic regression for biomedical data analysis.

Binary Logistic Regression

Basic Model

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

# Load data
df = pd.read_csv('clinical_data.csv')

# Method 1: Formula API (recommended)
model = smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0)

# Method 2: Matrix API
X = sm.add_constant(df[['exposure', 'age', 'sex']])
y = df['disease']
model = sm.Logit(y, X).fit(disp=0)

# Print summary
print(model.summary())

Extracting Odds Ratios

# Coefficients (log odds)
coefs = model.params
print("Log odds (coefficients):")
print(coefs)

# Odds ratios (exponentiate coefficients)
odds_ratios = np.exp(model.params)
print("\nOdds ratios:")
print(odds_ratios)

# 95% Confidence intervals
conf_int = model.conf_int()
conf_int_exp = np.exp(conf_int)

# Pretty print
for var in model.params.index:
    or_val = odds_ratios[var]
    ci_lower = conf_int_exp.loc[var, 0]
    ci_upper = conf_int_exp.loc[var, 1]
    p_val = model.pvalues[var]

    print(f"\n{var}:")
    print(f"  OR: {or_val:.4f}")
    print(f"  95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
    print(f"  p-value: {p_val:.6f}")
    print(f"  Significant: {'Yes' if p_val < 0.05 else 'No'}")

Model Fit Statistics

# Pseudo R-squared
print(f"McFadden's R²: {model.prsquared:.4f}")

# Log-likelihood
print(f"Log-likelihood: {model.llf:.4f}")

# AIC/BIC
print(f"AIC: {model.aic:.2f}")
print(f"BIC: {model.bic:.2f}")

# Number of observations
print(f"N: {int(model.nobs)}")

Categorical Predictors

Manual Dummy Coding

# Create dummy variables
df_encoded = pd.get_dummies(df, columns=['treatment_group'], drop_first=True, dtype=int)

# Fit model
model = smf.logit('disease ~ treatment_group_B + treatment_group_C + age',
                  data=df_encoded).fit(disp=0)

Formula with Categorical Variables

# statsmodels handles categorical automatically with C()
model = smf.logit('disease ~ C(treatment_group) + age', data=df).fit(disp=0)

# Set reference level
model = smf.logit('disease ~ C(treatment_group, Treatment("Control")) + age',
                  data=df).fit(disp=0)

Interaction Terms

Two-Way Interactions

# Interaction between continuous and binary
model = smf.logit('disease ~ exposure * age + sex', data=df).fit(disp=0)

# The interaction term is 'exposure:age'
interaction_coef = model.params['exposure:age']
interaction_or = np.exp(interaction_coef)
print(f"Interaction OR: {interaction_or:.4f}")

Interpreting Interactions

# Main effects + interaction
# Model: logit(p) = β0 + β1*exposure + β2*age + β3*exposure*age

# OR for exposure depends on age:
# OR(exposure) = exp(β1 + β3*age)

# Example: OR at age=30 vs age=50
beta_exposure = model.params['exposure']
beta_interaction = model.params['exposure:age']

or_age30 = np.exp(beta_exposure + beta_interaction * 30)
or_age50 = np.exp(beta_exposure + beta_interaction * 50)

print(f"OR (exposure) at age 30: {or_age30:.4f}")
print(f"OR (exposure) at age 50: {or_age50:.4f}")

Adjusted vs Unadjusted Analysis

Percentage Reduction in OR

# Unadjusted (crude) model
model_crude = smf.logit('disease ~ exposure', data=df).fit(disp=0)
or_crude = np.exp(model_crude.params['exposure'])

# Adjusted model
model_adj = smf.logit('disease ~ exposure + age + sex + bmi', data=df).fit(disp=0)
or_adj = np.exp(model_adj.params['exposure'])

# Calculate percentage reduction
pct_reduction = (or_crude - or_adj) / or_crude * 100

print(f"Crude OR: {or_crude:.4f}")
print(f"Adjusted OR: {or_adj:.4f}")
print(f"Percentage reduction: {pct_reduction:.1f}%")

# Interpretation
if pct_reduction > 10:
    print("Strong confounding detected")
elif pct_reduction > 5:
    print("Moderate confounding")
else:
    print("Minimal confounding")

Model Comparison

Likelihood Ratio Test

from scipy import stats

# Nested models
model_reduced = smf.logit('disease ~ exposure', data=df).fit(disp=0)
model_full = smf.logit('disease ~ exposure + age + sex + bmi', data=df).fit(disp=0)

# LR test statistic
lr_stat = -2 * (model_reduced.llf - model_full.llf)
df_diff = model_full.df_model - model_reduced.df_model
p_value = stats.chi2.sf(lr_stat, df_diff)

print(f"LR statistic: {lr_stat:.4f}")
print(f"df: {df_diff}")
print(f"p-value: {p_value:.6f}")

if p_value < 0.05:
    print("Full model provides significantly better fit")
else:
    print("Additional predictors not significant")

AIC/BIC Comparison

# Fit multiple models
models = {
    'Model 1': smf.logit('disease ~ exposure', data=df).fit(disp=0),
    'Model 2': smf.logit('disease ~ exposure + age', data=df).fit(disp=0),
    'Model 3': smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0),
}

# Compare
print("Model Comparison:")
for name, m in models.items():
    print(f"\n{name}:")
    print(f"  AIC: {m.aic:.2f}")
    print(f"  BIC: {m.bic:.2f}")
    print(f"  Pseudo R²: {m.prsquared:.4f}")

# Best model (lowest AIC)
best_model = min(models.items(), key=lambda x: x[1].aic)
print(f"\nBest model (by AIC): {best_model[0]}")

Prediction

Predicted Probabilities

# Get predicted probabilities for existing data
df['predicted_prob'] = model.predict(df)

# Predict for new data
new_data = pd.DataFrame({
    'exposure': [1],
    'age': [45],
    'sex': ['M']
})
pred_prob = model.predict(new_data)
print(f"Predicted probability: {pred_prob[0]:.4f}")

Classification

# Binary classification with 0.5 threshold
df['predicted_class'] = (model.predict(df) > 0.5).astype(int)

# Confusion matrix
from sklearn.metrics import confusion_matrix, classification_report
cm = confusion_matrix(df['disease'], df['predicted_class'])
print("Confusion Matrix:")
print(cm)

# Accuracy, precision, recall
print("\nClassification Report:")
print(classification_report(df['disease'], df['predicted_class']))

Diagnostics

Influential Observations

# Cook's distance
from statsmodels.stats.outliers_influence import OLSInfluence

# Get influence measures
influence = model.get_influence()

# Standardized residuals
std_resid = influence.resid_studentized

# Plot influential points
import matplotlib.pyplot as plt
plt.scatter(range(len(std_resid)), std_resid)
plt.axhline(y=2, color='r', linestyle='--')
plt.axhline(y=-2, color='r', linestyle='--')
plt.xlabel('Observation')
plt.ylabel('Studentized Residual')
plt.title('Influential Observations')
plt.show()

Hosmer-Lemeshow Test

# Goodness of fit test
from statsmodels.stats.diagnostic import _diagnostic_hl

# Not directly available in statsmodels
# Use manual implementation
def hosmer_lemeshow_test(y_true, y_pred, g=10):
    """Hosmer-Lemeshow goodness of fit test."""
    data = pd.DataFrame({'y': y_true, 'pred': y_pred})
    data['decile'] = pd.qcut(data['pred'], g, duplicates='drop')

    obs = data.groupby('decile')['y'].agg(['sum', 'count'])
    exp = data.groupby('decile')['pred'].agg(['sum', 'count'])

    hl_stat = ((obs['sum'] - exp['sum'])**2 / (exp['sum'] * (1 - exp['sum']/exp['count']))).sum()
    p_value = stats.chi2.sf(hl_stat, g-2)

    return hl_stat, p_value

hl_stat, hl_p = hosmer_lemeshow_test(df['disease'], model.predict(df))
print(f"Hosmer-Lemeshow: χ²={hl_stat:.4f}, p={hl_p:.6f}")

Common Issues

Separation (Perfect Prediction)

Problem: One predictor perfectly predicts the outcome.

Solution: Use Firth logistic regression (penalized likelihood):

# This requires logistf package (not in standard statsmodels)
# Alternative: Remove the problematic predictor or use regularization

# Check for separation
print("Value counts by predictor:")
print(pd.crosstab(df['exposure'], df['disease']))

# If separation detected, try Ridge logistic
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(penalty='l2', C=1.0)
lr.fit(df[['exposure', 'age', 'sex']], df['disease'])

Convergence Failure

Problem: Model doesn't converge (common with small samples or collinearity).

Solution: Increase max iterations or check for collinearity:

# Increase iterations
model = smf.logit('disease ~ exposure + age + sex', data=df).fit(disp=0, maxiter=200)

# Check for collinearity
from statsmodels.stats.outliers_influence import variance_inflation_factor
X = df[['exposure', 'age', 'sex']].copy()
X = pd.get_dummies(X, drop_first=True, dtype=float)
X['const'] = 1

for i, col in enumerate(X.columns[:-1]):
    vif = variance_inflation_factor(X.values, i)
    print(f"{col}: VIF={vif:.2f}")

Quasi-Complete Separation

Problem: Predictor almost perfectly separates outcomes.

Symptoms: Very large coefficients (>10), very large standard errors.

Solution: Use regularization or remove problematic predictor.

Reporting Template

def report_logistic_regression(model):
    """Generate publication-quality report."""
    report = []
    report.append("=== Logistic Regression Results ===\n")
    report.append(f"N = {int(model.nobs)}")
    report.append(f"Pseudo R² = {model.prsquared:.4f}")
    report.append(f"AIC = {model.aic:.2f}\n")

    report.append("Odds Ratios (95% CI):\n")
    ors = np.exp(model.params)
    ci = np.exp(model.conf_int())

    for var in model.params.index:
        if var == 'Intercept':
            continue
        or_val = ors[var]
        ci_lower = ci.loc[var, 0]
        ci_upper = ci.loc[var, 1]
        p_val = model.pvalues[var]
        sig = "*" if p_val < 0.05 else ""

        report.append(f"  {var}: OR={or_val:.4f} ({ci_lower:.4f}-{ci_upper:.4f}), p={p_val:.4f}{sig}")

    return "\n".join(report)

print(report_logistic_regression(model))

Ordinal Logistic Regression Reference

Complete guide to ordinal logistic regression (proportional odds model) for ordered categorical outcomes.

When to Use Ordinal Logistic Regression

Use ordinal logit when your outcome has:

• 3 or more levels (if 2 levels, use binary logistic)
• Natural ordering (mild < moderate < severe)
• No numeric interpretation (can't assume equal spacing between levels)

Examples:

• Disease severity: mild, moderate, severe, critical
• Cancer stage: I, II, III, IV
• Pain score: none, mild, moderate, severe
• Likert scale: strongly disagree, disagree, neutral, agree, strongly agree
• Performance status: 0, 1, 2, 3, 4

Basic Ordinal Logit Model

import pandas as pd
import numpy as np
from statsmodels.miscmodels.ordinal_model import OrderedModel

# Load data
df = pd.read_csv('data.csv')

# Define order of outcome levels
severity_order = ['Mild', 'Moderate', 'Severe', 'Critical']

# Convert to ordered categorical
df['severity'] = pd.Categorical(df['severity'],
                                 categories=severity_order,
                                 ordered=True)

# Encode as integer codes (0, 1, 2, 3)
y = df['severity'].cat.codes

# Prepare predictors (handle categorical variables)
X = df[['exposure', 'age', 'sex']].copy()
X = pd.get_dummies(X, drop_first=True, dtype=float)

# Fit proportional odds model
model = OrderedModel(y, X, distr='logit')
fit = model.fit(method='bfgs', disp=0, maxiter=200)

# Print summary
print(fit.summary())

Extracting Odds Ratios

# Number of outcome levels
n_levels = len(df['severity'].cat.categories)
n_thresholds = n_levels - 1

# Number of predictors
n_predictors = len(X.columns)

# Parameters are: [predictors, thresholds]
# Extract predictor coefficients only
predictor_params = fit.params[:n_predictors]
predictor_names = X.columns.tolist()

# Odds ratios
odds_ratios = np.exp(predictor_params)

# Confidence intervals
conf_int = fit.conf_int()
conf_int_exp = np.exp(conf_int.iloc[:n_predictors, :])

# Print results
print("\n=== Odds Ratios ===")
for i, name in enumerate(predictor_names):
    or_val = odds_ratios[i]
    ci_lower = conf_int_exp.iloc[i, 0]
    ci_upper = conf_int_exp.iloc[i, 1]
    p_val = fit.pvalues[i]

    print(f"\n{name}:")
    print(f"  OR: {or_val:.4f}")
    print(f"  95% CI: ({ci_lower:.4f}, {ci_upper:.4f})")
    print(f"  p-value: {p_val:.6f}")

Interpreting Odds Ratios

Proportional odds assumption: The odds ratio is constant across all levels of the outcome.

# Example: OR = 2.5 for exposure
# Interpretation:
# - Exposure increases odds of being in higher severity category by factor of 2.5
# - This applies to ALL cutpoints:
#   * Odds of moderate vs mild
#   * Odds of severe vs (moderate or mild)
#   * Odds of critical vs (severe or moderate or mild)

Interpretation Function

def interpret_ordinal_or(or_val, ci_lower, ci_upper, p_val, var_name, outcome_name):
    """Generate interpretation for ordinal OR."""
    if or_val > 1:
        direction = "higher"
        magnitude = (or_val - 1) * 100
        interp = f"{var_name} is associated with {magnitude:.1f}% increased odds of being in a higher {outcome_name} category"
    elif or_val < 1:
        direction = "lower"
        magnitude = (1 - or_val) * 100
        interp = f"{var_name} is associated with {magnitude:.1f}% decreased odds of being in a higher {outcome_name} category"
    else:
        interp = f"{var_name} has no association with {outcome_name}"

    sig = "statistically significant" if p_val < 0.05 else "not statistically significant"

    return f"{interp} (OR={or_val:.4f}, 95% CI [{ci_lower:.4f}, {ci_upper:.4f}], p={p_val:.6f}, {sig})"

# Example usage
print(interpret_ordinal_or(2.5, 1.8, 3.5, 0.001, "BCG vaccination", "COVID-19 severity"))

Thresholds (Cut Points)

# Extract thresholds
threshold_params = fit.params[n_predictors:]
threshold_names = [f"Threshold_{i}" for i in range(n_thresholds)]

print("\n=== Thresholds ===")
for i, name in enumerate(threshold_names):
    threshold = threshold_params[i]
    print(f"{name} ({severity_order[i]}|{severity_order[i+1]}): {threshold:.4f}")

Interpretation: Thresholds represent the log-odds cutpoints between adjacent categories when all predictors are 0.

Testing Proportional Odds Assumption

Critical assumption: The odds ratio is the same for all cutpoints.

Brant Test (Approximation)

import statsmodels.api as sm

def test_proportional_odds(df, outcome_col, predictors, order):
    """Test proportional odds assumption."""
    df_test = df.copy()
    df_test[outcome_col] = pd.Categorical(df_test[outcome_col],
                                          categories=order,
                                          ordered=True)
    y_codes = df_test[outcome_col].cat.codes
    n_levels = len(order)

    # Prepare predictors
    X = df_test[predictors].copy()
    X = pd.get_dummies(X, drop_first=True, dtype=float)
    X_const = sm.add_constant(X)

    # Fit binary logit at each cutpoint
    results = {}
    for k in range(n_levels - 1):
        y_binary = (y_codes > k).astype(int)
        try:
            binary_model = sm.Logit(y_binary, X_const).fit(disp=0)
            results[k] = {
                'cutpoint': f"{order[k]}|{order[k+1]}",
                'coefs': binary_model.params[1:].to_dict()  # Skip intercept
            }
        except Exception as e:
            print(f"Failed at cutpoint {k}: {e}")

    # Compare coefficients across cutpoints
    print("\n=== Proportional Odds Test ===")
    print("Coefficients should be similar across cutpoints:\n")

    for pred in X.columns:
        coefs = [results[k]['coefs'][pred] for k in results if pred in results[k]['coefs']]
        if len(coefs) > 1:
            coef_range = max(coefs) - min(coefs)
            print(f"{pred}:")
            for k in results:
                print(f"  {results[k]['cutpoint']}: {results[k]['coefs'].get(pred, 'N/A'):.4f}")
            print(f"  Range: {coef_range:.4f}")
            if coef_range > 0.5:
                print(f"  ⚠️  Warning: Large variation suggests violation of proportional odds")
            else:
                print(f"  ✅ Proportional odds likely satisfied")
            print()

    return results

# Run test
test_proportional_odds(df, 'severity', ['exposure', 'age', 'sex'], severity_order)

What to Do if Assumption Violated

1. Partial proportional odds model - Allow some predictors to vary across cutpoints 2. Multinomial logistic regression - Treat outcome as nominal (loses ordering info) 3. Alternative link functions - Try probit instead of logit 4. Transform outcome - Consider different categorization

# Multinomial logit as alternative (if PO violated)
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import LabelEncoder

le = LabelEncoder()
y_encoded = le.fit_transform(df['severity'])

X_array = pd.get_dummies(df[['exposure', 'age', 'sex']], drop_first=True, dtype=float).values

model_multinom = LogisticRegression(multi_class='multinomial', solver='lbfgs', max_iter=1000)
model_multinom.fit(X_array, y_encoded)

print("Multinomial model fitted (proportional odds not assumed)")

Model Fit Statistics

# Log-likelihood
print(f"Log-likelihood: {fit.llf:.4f}")

# AIC/BIC
print(f"AIC: {fit.aic:.2f}")
print(f"BIC: {fit.bic:.2f}")

# Number of observations
print(f"N: {len(y)}")

# Pseudo R-squared (not directly available, compute manually)
# McFadden's R²
ll_null = OrderedModel(y, np.ones((len(y), 1)), distr='logit').fit(disp=0).llf
pseudo_r2 = 1 - (fit.llf / ll_null)
print(f"Pseudo R²: {pseudo_r2:.4f}")

Prediction

# Predicted probabilities for each outcome level
pred_probs = fit.model.predict(fit.params, exog=X)

# pred_probs is an array of shape (n_obs, n_levels)
# Each row sums to 1

# Add to dataframe
for i, level in enumerate(df['severity'].cat.categories):
    df[f'prob_{level}'] = pred_probs[:, i]

# Predicted category (highest probability)
df['predicted_severity'] = df['severity'].cat.categories[pred_probs.argmax(axis=1)]

# Show first few predictions
print(df[['severity', 'predicted_severity', 'prob_Mild', 'prob_Moderate', 'prob_Severe']].head())

Complete Example: COVID-19 Severity

import pandas as pd
import numpy as np
from statsmodels.miscmodels.ordinal_model import OrderedModel

# Simulated COVID-19 severity data
np.random.seed(42)
n = 500

df = pd.DataFrame({
    'age': np.random.normal(55, 15, n),
    'bcg_vaccination': np.random.binomial(1, 0.6, n),
    'comorbidities': np.random.binomial(1, 0.3, n),
    'severity': np.random.choice(['Mild', 'Moderate', 'Severe'], n, p=[0.5, 0.3, 0.2])
})

# Define ordered outcome
severity_order = ['Mild', 'Moderate', 'Severe']
df['severity'] = pd.Categorical(df['severity'], categories=severity_order, ordered=True)
y = df['severity'].cat.codes

# Prepare predictors
X = df[['age', 'bcg_vaccination', 'comorbidities']].copy()

# Fit model
model = OrderedModel(y, X, distr='logit')
fit = model.fit(method='bfgs', disp=0)

# Extract results
n_predictors = len(X.columns)
ors = np.exp(fit.params[:n_predictors])
ci = np.exp(fit.conf_int().iloc[:n_predictors, :])

print("=== COVID-19 Severity Analysis ===\n")
print("Outcome: Mild → Moderate → Severe\n")

for i, var in enumerate(X.columns):
    or_val = ors[i]
    ci_lower = ci.iloc[i, 0]
    ci_upper = ci.iloc[i, 1]
    p_val = fit.pvalues[i]

    print(f"{var}:")
    print(f"  OR: {or_val:.4f} (95% CI: {ci_lower:.4f}-{ci_upper:.4f})")
    print(f"  p-value: {p_val:.6f}")

    if var == 'bcg_vaccination' and or_val < 1 and p_val < 0.05:
        pct_reduction = (1 - or_val) * 100
        print(f"  Interpretation: BCG vaccination reduces odds of higher severity by {pct_reduction:.1f}%")
    print()

Common Question Pattern

Typical question: "What is the odds ratio of COVID-19 severity associated with BCG vaccination in ordinal logistic regression?"

Solution: 1. Identify ordinal outcome (severity levels) 2. Define ordering (mild → moderate → severe) 3. Fit ordinal logistic regression 4. Extract OR for BCG vaccination predictor 5. Report OR with CI and p-value

# Answer format
or_bcg = ors[X.columns.get_loc('bcg_vaccination')]
ci_lower_bcg = ci.iloc[X.columns.get_loc('bcg_vaccination'), 0]
ci_upper_bcg = ci.iloc[X.columns.get_loc('bcg_vaccination'), 1]

print(f"Answer: {or_bcg:.4f}")
print(f"95% CI: ({ci_lower_bcg:.4f}, {ci_upper_bcg:.4f})")

Common Issues

Issue 1: Convergence Failure

Problem: Model doesn't converge.

Solutions:

# Solution 1: Increase max iterations
fit = model.fit(method='bfgs', disp=0, maxiter=500)

# Solution 2: Try different optimizer
fit = model.fit(method='nm', disp=0)  # Nelder-Mead

# Solution 3: Scale predictors
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = pd.DataFrame(scaler.fit_transform(X), columns=X.columns)
fit = OrderedModel(y, X_scaled, distr='logit').fit(method='bfgs', disp=0)

Issue 2: Separation

Problem: Predictor perfectly separates some categories.

Check:

# Check crosstabs
for pred in ['bcg_vaccination', 'comorbidities']:
    print(f"\n{pred} by severity:")
    print(pd.crosstab(df[pred], df['severity']))

Solution: Remove problematic predictor or use regularization.

Issue 3: Small Sample Size

Problem: Too few observations per category.

Rule of thumb: Need at least 10 events per predictor per category.

# Check sample size requirements
n_predictors = len(X.columns)
n_per_level = df['severity'].value_counts()

print("\nSample size check:")
for level, count in n_per_level.items():
    ratio = count / n_predictors
    print(f"{level}: {count} obs ({ratio:.1f} per predictor)")
    if ratio < 10:
        print(f"  ⚠️  Warning: Small sample size for {level}")

Reporting Template

def report_ordinal_logit(fit, X, outcome_levels):
    """Generate publication-ready report."""
    n_predictors = len(X.columns)

    report = []
    report.append("=== Ordinal Logistic Regression Results ===\n")
    report.append(f"Outcome levels: {' → '.join(outcome_levels)}")
    report.append(f"N = {len(X)}")
    report.append(f"AIC = {fit.aic:.2f}\n")

    report.append("Odds Ratios (95% CI):\n")
    ors = np.exp(fit.params[:n_predictors])
    ci = np.exp(fit.conf_int().iloc[:n_predictors, :])

    for i, var in enumerate(X.columns):
        or_val = ors[i]
        ci_lower = ci.iloc[i, 0]
        ci_upper = ci.iloc[i, 1]
        p_val = fit.pvalues[i]
        sig = "*" if p_val < 0.05 else ""

        report.append(f"  {var}: OR={or_val:.4f} ({ci_lower:.4f}-{ci_upper:.4f}), p={p_val:.6f}{sig}")

    return "\n".join(report)

print(report_ordinal_logit(fit, X, severity_order))

Troubleshooting Guide

Common statistical modeling issues and solutions.

Convergence Issues

Problem: Model doesn't converge

Symptoms:

Warning: Maximum iterations reached
ConvergenceWarning: Maximum Likelihood optimization failed to converge

Solutions:

1. Increase max iterations:

# Logistic regression
model = smf.logit('y ~ x + z', data=df).fit(disp=0, maxiter=500)

# Ordinal logit
model = OrderedModel(y, X, distr='logit').fit(method='bfgs', maxiter=500)

# Cox model (lifelines handles this automatically)

2. Scale predictors:

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = pd.DataFrame(scaler.fit_transform(X), columns=X.columns)
model = OrderedModel(y, X_scaled, distr='logit').fit(method='bfgs')

3. Try different optimization method:

# Try Nelder-Mead instead of BFGS
model = OrderedModel(y, X, distr='logit').fit(method='nm')

4. Check for separation (see below)

---

Separation Problems

Problem: Perfect or quasi-complete separation

Symptoms:

• Very large coefficients (>10)
• Very large standard errors
• Warning: "Perfect separation detected"

Check for separation:

# Crosstab of predictor vs outcome
for pred in ['exposure', 'treatment']:
    print(f"\n{pred} by outcome:")
    print(pd.crosstab(df[pred], df['outcome']))

# Look for cells with 0 counts - that's separation

Solutions:

1. Remove problematic predictor:

# If one predictor causes separation, exclude it
model = smf.logit('outcome ~ age + sex', data=df).fit(disp=0)  # Exclude 'exposure'

2. Use Firth logistic regression (penalized likelihood):

# Requires logistf package (not standard)
# Alternative: use Ridge penalty in sklearn

from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(penalty='l2', C=1.0, solver='lbfgs')
lr.fit(X, y)

3. Collapse categories:

# If categorical predictor has sparse levels, combine them
df['stage_collapsed'] = df['stage'].replace({'I': 'Early', 'II': 'Early',
                                               'III': 'Late', 'IV': 'Late'})

4. Increase sample size (if possible)

---

Multicollinearity

Problem: Predictors highly correlated

Symptoms:

• Large standard errors
• VIF > 10
• Coefficients change dramatically when adding/removing predictors

Check for multicollinearity:

from statsmodels.stats.outliers_influence import variance_inflation_factor

X = df[['age', 'bmi', 'weight', 'height']].copy()
X['const'] = 1

for i, col in enumerate(X.columns[:-1]):
    vif = variance_inflation_factor(X.values, i)
    print(f"{col}: VIF = {vif:.2f}")
    if vif > 10:
        print(f"  ⚠️  High multicollinearity")

Solutions:

1. Remove correlated predictors:

# Check correlation matrix
corr_matrix = df[['age', 'bmi', 'weight', 'height']].corr()
print(corr_matrix)

# Remove one of highly correlated pairs (r > 0.8)
# E.g., remove weight if weight and BMI are r=0.9
model = smf.ols('outcome ~ age + bmi', data=df).fit()  # Exclude weight

2. Use Ridge regression (L2 regularization):

from sklearn.linear_model import Ridge

ridge = Ridge(alpha=1.0)
ridge.fit(X, y)

3. Principal Component Analysis:

from sklearn.decomposition import PCA

pca = PCA(n_components=3)
X_pca = pca.fit_transform(X)
# Use principal components as predictors

---

Heteroscedasticity

Problem: Non-constant variance of residuals

Symptoms:

• Breusch-Pagan test p < 0.05
• Residual plot shows funnel shape

Check:

from statsmodels.stats.diagnostic import het_breuschpagan

residuals = model.resid
bp_stat, bp_p, _, _ = het_breuschpagan(residuals, model.model.exog)

if bp_p < 0.05:
    print("⚠️  Heteroscedasticity detected")

Solutions:

1. Robust standard errors:

model_robust = model.get_robustcov_results(cov_type='HC3')
print(model_robust.summary())

2. Log transformation (if outcome is right-skewed):

df['outcome_log'] = np.log(df['outcome'] + 1)
model = smf.ols('outcome_log ~ x + z', data=df).fit()

3. Weighted Least Squares:

# Weight by inverse variance
residuals_sq = model.resid ** 2
weights = 1 / residuals_sq
wls_model = sm.WLS(y, X, weights=weights).fit()

---

Non-Normality of Residuals

Problem: Residuals not normally distributed

Symptoms:

• Shapiro-Wilk test p < 0.05
• Q-Q plot deviates from line

Check:

from scipy import stats as scipy_stats

residuals = model.resid
sw_stat, sw_p = scipy_stats.shapiro(residuals)

if sw_p < 0.05:
    print("⚠️  Residuals not normally distributed")

Solutions:

1. Transform outcome:

# Log transformation
df['outcome_log'] = np.log(df['outcome'] + 1)

# Square root transformation
df['outcome_sqrt'] = np.sqrt(df['outcome'])

# Box-Cox transformation
from scipy.stats import boxcox
df['outcome_bc'], lambda_param = boxcox(df['outcome'] + 1)

2. Use robust regression:

from statsmodels.robust.robust_linear_model import RLM

rlm_model = RLM.from_formula('outcome ~ x + z', data=df).fit()

3. Use non-parametric methods:

# Bootstrap confidence intervals instead of t-tests

Note: For large samples (n > 30), non-normality is less critical due to Central Limit Theorem.

---

Missing Data

Problem: Missing values in predictors or outcome

Check:

# Count missing values
missing = df.isnull().sum()
print("\nMissing values:")
print(missing[missing > 0])

# Percentage missing
pct_missing = (missing / len(df) * 100)
print("\nPercentage missing:")
print(pct_missing[pct_missing > 0])

Solutions:

1. Complete case analysis (delete rows with missing):

df_complete = df.dropna(subset=['outcome', 'x', 'z'])
model = smf.ols('outcome ~ x + z', data=df_complete).fit()

2. Mean/median imputation (simple):

df['age'].fillna(df['age'].mean(), inplace=True)

3. Multiple imputation (best practice):

from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

imputer = IterativeImputer(random_state=42)
df_imputed = pd.DataFrame(imputer.fit_transform(df),
                          columns=df.columns)

4. Missing indicator method:

# Create indicator for missingness
df['age_missing'] = df['age'].isnull().astype(int)
df['age'].fillna(df['age'].mean(), inplace=True)

# Include indicator in model
model = smf.ols('outcome ~ age + age_missing + x', data=df).fit()

---

Proportional Hazards Violation

Problem: PH assumption violated in Cox model

Check:

# Test PH assumption
results = cph.check_assumptions(df, p_value_threshold=0.05, show_plots=False)

if len(results) > 0:
    print(f"⚠️  PH violated for: {results}")

Solutions:

1. Stratify by problematic variable:

# Don't estimate HR for treatment, but adjust for it
cph_strat = CoxPHFitter()
cph_strat.fit(df, duration_col='time', event_col='event',
              strata=['treatment'])

2. Time-varying coefficients:

# Allow coefficient to change over time (advanced)
# Interact predictor with time
df['treatment_time'] = df['treatment'] * df['time']
cph.fit(df[['time', 'event', 'treatment', 'treatment_time', 'age']],
        duration_col='time', event_col='event')

3. Use parametric survival model:

from lifelines import WeibullAFTFitter

# Accelerated failure time model (no PH assumption)
wf = WeibullAFTFitter()
wf.fit(df, duration_col='time', event_col='event')

---

Proportional Odds Violation

Problem: PO assumption violated in ordinal logit

Check:

# Fit binary logits at each cutpoint, compare coefficients
# See ordinal_logistic.md for full test

Solutions:

1. Partial proportional odds model:

# Allow some predictors to vary across cutpoints (requires mord package)

2. Multinomial logistic regression:

from sklearn.linear_model import LogisticRegression

# Treat outcome as nominal (loses ordering information)
model = LogisticRegression(multi_class='multinomial', solver='lbfgs')
model.fit(X, y)

3. Adjacent category logit (alternative ordinal model)

---

Small Sample Size

Problem: Too few observations per predictor

Rule of thumb:

• Linear regression: ≥20 observations per predictor
• Logistic regression: ≥10 events per predictor
• Ordinal logit: ≥10 observations per outcome level per predictor
• Cox regression: ≥10 events per predictor

Check:

n = len(df)
n_predictors = len(['x', 'z', 'age', 'sex'])  # Your predictors

print(f"Observations per predictor: {n / n_predictors:.1f}")

# For logistic
n_events = df['outcome'].sum()
print(f"Events per predictor: {n_events / n_predictors:.1f}")

if n_events / n_predictors < 10:
    print("⚠️  Small sample size - results may be unreliable")

Solutions:

1. Reduce number of predictors:

# Only include most important predictors
model = smf.logit('outcome ~ exposure + age', data=df).fit(disp=0)

2. Use penalized regression:

from sklearn.linear_model import LogisticRegression

# Ridge penalty helps with small samples
lr = LogisticRegression(penalty='l2', C=1.0)
lr.fit(X, y)

3. Exact logistic regression (for very small samples):

# Requires R or specialized packages

---

Outliers and Influential Points

Problem: Outliers affecting model fit

Check:

from statsmodels.stats.outliers_influence import OLSInfluence

influence = model.get_influence()

# Cook's distance
cooks_d = influence.cooks_distance[0]
influential = cooks_d > 4 / len(df)

print(f"Influential points: {influential.sum()}")
print(f"Indices: {df.index[influential].tolist()}")

Solutions:

1. Remove outliers (if justified):

# Remove points with Cook's distance > 4/n
df_clean = df[~influential]
model_clean = smf.ols('outcome ~ x + z', data=df_clean).fit()

2. Robust regression:

from statsmodels.robust.robust_linear_model import RLM

# Downweights outliers automatically
rlm_model = RLM.from_formula('outcome ~ x + z', data=df).fit()

3. Winsorize extreme values:

from scipy.stats.mstats import winsorize

# Cap extreme values at 5th and 95th percentiles
df['outcome_wins'] = winsorize(df['outcome'], limits=[0.05, 0.05])

---

Model Selection Uncertainty

Problem: Unsure which predictors to include

Solutions:

1. Forward selection:

# Start with null model, add predictors one by one
# Keep if p < 0.05 or AIC improves

2. Backward elimination:

# Start with full model, remove predictors one by one
# Remove if p > 0.10 or AIC improves

3. LASSO for variable selection:

from sklearn.linear_model import LogisticRegressionCV

# LASSO automatically selects variables
lasso = LogisticRegressionCV(penalty='l1', solver='saga', cv=5)
lasso.fit(X, y)

# Non-zero coefficients are selected
selected = X.columns[lasso.coef_[0] != 0]
print(f"Selected variables: {selected.tolist()}")

4. Use domain knowledge:

# Always include clinically important confounders
# Age, sex are usually important in biomedical studies

---

Package-Specific Issues

statsmodels singular matrix error

Problem:

LinAlgError: singular matrix

Cause: Perfect multicollinearity (one predictor is linear combination of others)

Solution:

# Check correlation matrix
corr = X.corr()
print(corr)

# Remove one of perfectly correlated predictors
# Or use pd.get_dummies(..., drop_first=True)

lifelines convergence warning

Problem:

ConvergenceWarning: Newton-Raphson failed to converge

Solutions:

# 1. Check for separation
# 2. Scale predictors
# 3. Use robust=True
cph.fit(df, duration_col='time', event_col='event', robust=True)

---

Quick Diagnostic Checklist

Before finalizing analysis:

• [ ] Check for missing data
• [ ] Check variable distributions (outliers, skewness)
• [ ] Check for multicollinearity (VIF < 10)
• [ ] Check model convergence
• [ ] Check sample size adequacy
• [ ] Run residual diagnostics
• [ ] Test model assumptions
• [ ] Compare alternative models
• [ ] Perform sensitivity analyses
• [ ] Interpret results in context