Statistical analysis

Lesson 10 — Python

Lesson 10 of 20 Intermediate to advanced ~110 min

Learning objectives

Produce descriptive statistics and frequency tables
Run t-tests, including the paired and Welch variants, and report them properly
Choose and run categorical tests: chi-squared and Fisher’s exact
Fit and interpret linear and logistic regression with statsmodels
Run survival analysis: Kaplan-Meier, log-rank and Cox regression
Handle multiplicity, and know which library to reach for

The libraries

Library	Covers
`scipy.stats`	Hypothesis tests, distributions, descriptive statistics
`statsmodels`	Regression, ANOVA, mixed models, formula interface
`lifelines`	Survival analysis — KM, log-rank, Cox
`pingouin`	Friendlier test output, effect sizes, assumption checks
`scikit-learn`	Prediction, not inference — see lesson 11

uv add scipy statsmodels lifelines pingouin

scikit-learn is not a statistics library

sklearn.linear_model.LogisticRegression applies L2 regularisation by default and gives you no p-values, no confidence intervals and no standard errors. It is built for prediction.

For inference — an odds ratio you will put in a report — use statsmodels. Fitting the two on the same data gives different coefficients, and the reason is the hidden penalty.

Descriptive statistics

import pandas as pd
import numpy as np
from scipy import stats

adsl["AGE"].describe()
adsl["AGE"].agg(["count", "mean", "std", "median", "min", "max", "skew"])

adsl.groupby("TRT01A")["AGE"].describe()

A clinical summary needs specific statistics in a specific order:

def describe_continuous(s: pd.Series) -> pd.Series:
    """Descriptive statistics as they appear in a clinical table."""
    x = s.dropna()
    if x.empty:
        return pd.Series({k: np.nan for k in
                          ["n", "nmiss", "mean", "sd", "median", "q1", "q3",
                           "min", "max", "se", "ci_lower", "ci_upper"]})

    se = x.std(ddof=1) / np.sqrt(len(x))
    ci = stats.t.interval(0.95, len(x) - 1, loc=x.mean(), scale=se)

    return pd.Series({
        "n":        len(x),
        "nmiss":    s.isna().sum(),
        "mean":     x.mean(),
        "sd":       x.std(ddof=1),
        "median":   x.median(),
        "q1":       x.quantile(0.25),
        "q3":       x.quantile(0.75),
        "min":      x.min(),
        "max":      x.max(),
        "se":       se,
        "ci_lower": ci[0],
        "ci_upper": ci[1],
    })


summary = (
    adsl.query("SAFFL == 'Y'")
    .groupby("TRT01A")["AGE"]
    .apply(describe_continuous)
    .unstack()
    .round(2)
)

ddof=1 — pandas and NumPy disagree

x = pd.Series([1, 2, 3, 4])
x.std()             # 1.291  — sample SD, ddof=1 (pandas default)
np.std(x)           # 1.118  — population SD, ddof=0 (NumPy default)
np.std(x, ddof=1)   # 1.291

pandas defaults to the sample standard deviation, NumPy to the population one. Mixing the two in one pipeline produces a table where some SDs are computed differently from others, and nothing warns you. Be explicit everywhere.

Frequency tables

adsl["AGEGR1"].value_counts()
adsl["AGEGR1"].value_counts(normalize=True)          # proportions
adsl["AGEGR1"].value_counts(dropna=False)            # include missing

pd.crosstab(adsl["TRT01A"], adsl["SEX"])
pd.crosstab(adsl["TRT01A"], adsl["SEX"], margins=True)
pd.crosstab(adsl["TRT01A"], adsl["SEX"], normalize="index")   # row percentages

For a clinical table, counts and percentages together, with zero rows kept:

def frequency_table(df, var, by, denominators):
    """n (%) per group, with all categories shown including zero counts."""
    counts = (
        pd.crosstab(df[by], df[var], dropna=False)
        .reindex(columns=df[var].cat.categories if hasattr(df[var], "cat")
                 else sorted(df[var].dropna().unique()), fill_value=0)
        .reindex(index=denominators.index, fill_value=0)
    )
    pct = counts.div(denominators, axis=0) * 100
    return counts.astype(str) + " (" + pct.round(1).astype(str) + ")"

The reindex calls are what guarantee a category with zero subjects in one arm still appears as 0 (0.0) rather than vanishing — the same requirement as in TLF generation. A blank cell reads as “not evaluated”; 0 (0.0) reads as “evaluated, none occurred”.

Continuous comparisons

Two independent groups

from scipy import stats

placebo = adsl.query("TRT01A == 'Placebo'")["AGE"].dropna()
active  = adsl.query("TRT01A == 'Drug A'")["AGE"].dropna()

# Welch's t-test — does NOT assume equal variances. This is the safe default.
t, p = stats.ttest_ind(active, placebo, equal_var=False)

# Student's t-test — assumes equal variances
t, p = stats.ttest_ind(active, placebo, equal_var=True)

# Non-parametric alternative
u, p = stats.mannwhitneyu(active, placebo, alternative="two-sided")

Default to Welch

equal_var=False is Welch’s t-test. It costs almost nothing in power when variances are equal, and is substantially more reliable when they are not — particularly with unequal group sizes, which is common in trials.

Testing for equal variances first and then choosing the test inflates the type I error rate. Just use Welch.

scipy returns only the statistic and p-value. A report needs the difference and its confidence interval:

def t_test_report(a: pd.Series, b: pd.Series, alpha: float = 0.05) -> dict:
    """Welch's t-test with the estimate, CI and effect size."""
    a, b = a.dropna(), b.dropna()
    na, nb = len(a), len(b)
    if na < 2 or nb < 2:
        raise ValueError(f"Need at least 2 observations per group, got {na} and {nb}")

    diff = a.mean() - b.mean()
    se = np.sqrt(a.var(ddof=1) / na + b.var(ddof=1) / nb)

    # Welch-Satterthwaite degrees of freedom
    df = (a.var(ddof=1) / na + b.var(ddof=1) / nb) ** 2 / (
        (a.var(ddof=1) / na) ** 2 / (na - 1) + (b.var(ddof=1) / nb) ** 2 / (nb - 1)
    )
    crit = stats.t.ppf(1 - alpha / 2, df)
    t_stat, p = stats.ttest_ind(a, b, equal_var=False)

    pooled_sd = np.sqrt(((na - 1) * a.var(ddof=1) + (nb - 1) * b.var(ddof=1))
                        / (na + nb - 2))

    return {
        "n1": na, "n2": nb,
        "mean1": a.mean(), "mean2": b.mean(),
        "difference": diff,
        "se": se, "df": df,
        "ci_lower": diff - crit * se,
        "ci_upper": diff + crit * se,
        "t": t_stat, "p_value": p,
        "cohens_d": diff / pooled_sd,
    }

Paired data

Change from baseline within the same subjects:

paired = (
    adlb.query("PARAMCD == 'ALT' and AVISITN in [0, 12]")
    .pivot(index="USUBJID", columns="AVISITN", values="AVAL")
    .dropna()                      # complete pairs only
)

t, p = stats.ttest_rel(paired[12], paired[0])

# Non-parametric equivalent
w, p = stats.wilcoxon(paired[12], paired[0])

.dropna() here is an analysis decision, not a technicality

Dropping incomplete pairs is a complete-case analysis. Subjects who withdrew before week 12 are silently excluded, and if withdrawal relates to the outcome — which in a trial it usually does — the result is biased.

Report how many were dropped, and let the statistician decide:

n_total = adlb.query("PARAMCD == 'ALT' and AVISITN == 0")["USUBJID"].nunique()
n_paired = len(paired)
print(f"Complete pairs: {n_paired} of {n_total} "
      f"({100 * n_paired / n_total:.1f}%)")

A mixed model (MMRM) uses all available data under a weaker missingness assumption, and is what most SAPs specify for this reason.

More than two groups

groups = [g["AGE"].dropna() for _, g in adsl.groupby("TRT01A")]

f, p = stats.f_oneway(*groups)              # one-way ANOVA
h, p = stats.kruskal(*groups)               # non-parametric

# With statsmodels, for the full table
import statsmodels.api as sm
from statsmodels.formula.api import ols

model = ols("AGE ~ C(TRT01A)", data=adsl).fit()
sm.stats.anova_lm(model, typ=2)

Categorical comparisons

table = pd.crosstab(adsl["TRT01A"], adsl["RESPONSE"])

# Chi-squared
chi2, p, dof, expected = stats.chi2_contingency(table)

# Check the assumption that justifies it
if (expected < 5).any():
    print("Expected counts below 5 — chi-squared is unreliable, use Fisher")

# Fisher's exact — exact, no large-sample assumption
odds_ratio, p = stats.fisher_exact(table)          # 2x2 only

# Larger tables need an exact test from elsewhere
from scipy.stats import chi2_contingency
chi2_contingency(table, correction=False)

stats.fisher_exact handles 2×2 only

For an R×C table, scipy cannot do an exact test. Options:

# Monte Carlo approximation
stats.chi2_contingency(table)      # asymptotic, needs expected >= 5

# Exact, via R
import rpy2.robjects as ro
ro.r("fisher.test")(table.values, simulate_p_value=True)

Or restructure the question into 2×2 comparisons with a multiplicity adjustment.

Risk difference and relative risk

The numbers a clinical report actually wants:

def two_by_two(a: int, b: int, c: int, d: int, alpha: float = 0.05) -> dict:
    """Risk difference, relative risk and odds ratio with 95% CIs.

    Layout:
                 Event    No event
      Treatment    a          b
      Control      c          d
    """
    n1, n0 = a + b, c + d
    if n1 == 0 or n0 == 0:
        raise ValueError("Both groups must have at least one subject")

    p1, p0 = a / n1, c / n0
    z = stats.norm.ppf(1 - alpha / 2)

    # Risk difference — Wald interval
    rd = p1 - p0
    se_rd = np.sqrt(p1 * (1 - p1) / n1 + p0 * (1 - p0) / n0)

    # Relative risk and odds ratio — CIs on the log scale
    with np.errstate(divide="ignore", invalid="ignore"):
        rr = p1 / p0 if p0 > 0 else np.inf
        se_log_rr = np.sqrt(1 / a - 1 / n1 + 1 / c - 1 / n0) if a and c else np.nan
        or_ = (a * d) / (b * c) if b and c else np.inf
        se_log_or = np.sqrt(1 / a + 1 / b + 1 / c + 1 / d) if all([a, b, c, d]) else np.nan

    _, p_fisher = stats.fisher_exact([[a, b], [c, d]])

    return {
        "n_treatment": n1, "n_control": n0,
        "events_treatment": a, "events_control": c,
        "risk_treatment": p1, "risk_control": p0,
        "risk_difference": rd,
        "rd_ci": (rd - z * se_rd, rd + z * se_rd),
        "relative_risk": rr,
        "rr_ci": (rr * np.exp(-z * se_log_rr), rr * np.exp(z * se_log_rr))
                 if np.isfinite(se_log_rr) else (np.nan, np.nan),
        "odds_ratio": or_,
        "or_ci": (or_ * np.exp(-z * se_log_or), or_ * np.exp(z * se_log_or))
                 if np.isfinite(se_log_or) else (np.nan, np.nan),
        "p_fisher": p_fisher,
    }

The zero-cell guards matter. A treatment arm with no events makes the odds ratio infinite and its standard error undefined; returning nan rather than crashing lets the caller decide whether to apply a continuity correction.

Regression

Linear

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

model = smf.ols("CHG ~ TRT01AN + BASE + AGE + C(SEX)", data=adlb).fit()
print(model.summary())

model.params
model.pvalues
model.conf_int()
model.rsquared_adj

# Robust standard errors, when residuals are heteroscedastic
model_robust = smf.ols("CHG ~ TRT01AN + BASE", data=adlb).fit(cov_type="HC3")

The formula interface follows R closely: C() marks a categorical, : is an interaction, * is main effects plus interaction, -1 drops the intercept.

Extract a tidy table, the equivalent of broom::tidy():

def tidy(model, conf_level: float = 0.95) -> pd.DataFrame:
    """Model coefficients as a tidy DataFrame."""
    ci = model.conf_int(alpha=1 - conf_level)
    return pd.DataFrame({
        "term":      model.params.index,
        "estimate":  model.params.values,
        "std_error": model.bse.values,
        "statistic": model.tvalues.values,
        "p_value":   model.pvalues.values,
        "conf_low":  ci[0].values,
        "conf_high": ci[1].values,
    }).reset_index(drop=True)

Logistic

model = smf.logit("RESPONSE ~ TRT01AN + AGE + C(SEX)", data=adsl).fit()
print(model.summary())

# Odds ratios with confidence intervals
odds = pd.DataFrame({
    "odds_ratio": np.exp(model.params),
    "ci_lower":   np.exp(model.conf_int()[0]),
    "ci_upper":   np.exp(model.conf_int()[1]),
    "p_value":    model.pvalues,
})

Complete separation

PerfectSeparationWarning: Perfect separation detected, results not available

This means a predictor perfectly predicts the outcome — every subject above a threshold responded, say. Coefficients diverge to infinity and standard errors are meaningless.

It is a data finding, not a bug. Options: drop the offending predictor, collapse sparse categories, or use penalised (Firth) logistic regression. Do not just suppress the warning.

Survival analysis

from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test, multivariate_logrank_test

adtte = adtte.query("PARAMCD == 'OS'")

# lifelines uses event=1; CDISC CNSR is 1 for CENSORED — invert it
adtte["event"] = 1 - adtte["CNSR"]

CNSR is inverted relative to lifelines

CDISC ADaM: CNSR = 0 means the event occurred, CNSR = 1 means censored. lifelines and most Python survival code expect event_observed = 1 for an event.

Getting this backwards produces a survival curve that goes up, or one that looks plausible but is exactly wrong. Assert it:

assert adtte["event"].isin([0, 1]).all()
assert adtte.loc[adtte["CNSR"] == 0, "event"].eq(1).all()

Kaplan-Meier

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(9, 6))
kmf = KaplanMeierFitter()

for arm, group in adtte.groupby("TRT01A"):
    kmf.fit(group["AVAL"], group["event"], label=arm)
    kmf.plot_survival_function(ax=ax, ci_show=True)
    print(f"{arm}: median survival = {kmf.median_survival_time_:.1f} days")

ax.set_xlabel("Time from first dose (days)")
ax.set_ylabel("Survival probability")
ax.set_ylim(0, 1)

from lifelines.plotting import add_at_risk_counts
add_at_risk_counts(*fitters, ax=ax)

Survival probability at specific times:

kmf.predict([90, 180, 365])
kmf.confidence_interval_survival_function_
kmf.median_survival_time_

from lifelines.utils import median_survival_times
median_survival_times(kmf.confidence_interval_)

Log-rank test

a = adtte.query("TRT01A == 'Drug A'")
p = adtte.query("TRT01A == 'Placebo'")

result = logrank_test(a["AVAL"], p["AVAL"], a["event"], p["event"])
result.p_value
result.test_statistic
result.print_summary()

# More than two groups
mv = multivariate_logrank_test(adtte["AVAL"], adtte["TRT01A"], adtte["event"])

Cox proportional hazards

cox_data = adtte[["AVAL", "event", "TRT01AN", "AGE", "SEX"]].copy()
cox_data["SEX"] = (cox_data["SEX"] == "M").astype(int)

cph = CoxPHFitter()
cph.fit(cox_data, duration_col="AVAL", event_col="event")
cph.print_summary()

cph.hazard_ratios_
cph.confidence_intervals_

# Check the proportional hazards assumption — do not skip this
cph.check_assumptions(cox_data, p_value_threshold=0.05, show_plots=True)

check_assumptions() tests whether hazard ratios are constant over time. If they are not, the single hazard ratio a Cox model reports is an average over a changing effect, and reporting it without comment is misleading.

Multiplicity

from statsmodels.stats.multitest import multipletests

p_values = [0.001, 0.013, 0.021, 0.043, 0.210]

reject, p_adj, _, _ = multipletests(p_values, alpha=0.05, method="bonferroni")
reject, p_adj, _, _ = multipletests(p_values, alpha=0.05, method="holm")
reject, p_adj, _, _ = multipletests(p_values, alpha=0.05, method="fdr_bh")

Method	Controls	Use when
`bonferroni`	Family-wise error rate	Few tests, strong control needed
`holm`	Family-wise error rate	Uniformly better than Bonferroni — prefer it
`fdr_bh`	False discovery rate	Many exploratory tests, e.g. biomarkers

For a confirmatory trial the multiplicity strategy is pre-specified in the SAP, often as a hierarchical testing procedure rather than a p-value adjustment. Do not invent one at analysis time.

R and Python side by side

Analysis	R	Python
Descriptives	`summary(x)`	`x.describe()`
Two-sample t-test	`t.test(a, b)`	`stats.ttest_ind(a, b, equal_var=False)`
Paired t-test	`t.test(a, b, paired = TRUE)`	`stats.ttest_rel(a, b)`
Wilcoxon	`wilcox.test()`	`stats.mannwhitneyu()` / `stats.wilcoxon()`
Chi-squared	`chisq.test(tbl)`	`stats.chi2_contingency(tbl)`
Fisher exact	`fisher.test(tbl)`	`stats.fisher_exact(tbl)`
ANOVA	`aov(y ~ g)`	`smf.ols("y ~ C(g)").fit()`
Linear model	`lm(y ~ x)`	`smf.ols("y ~ x", data).fit()`
Logistic	`glm(y ~ x, family = binomial)`	`smf.logit("y ~ x", data).fit()`
Tidy output	`broom::tidy(fit)`	`tidy(fit)` (hand-rolled)
Kaplan-Meier	`survfit(Surv(t, e) ~ g)`	`KaplanMeierFitter().fit()`
Log-rank	`survdiff()`	`logrank_test()`
Cox	`coxph(Surv(t, e) ~ x)`	`CoxPHFitter().fit()`
PH assumption	`cox.zph(fit)`	`cph.check_assumptions()`
Multiplicity	`p.adjust(p, "holm")`	`multipletests(p, method="holm")`
MMRM	`mmrm::mmrm()`	`MixedLM` (not equivalent)

Where R is genuinely ahead

MMRM. statsmodels.MixedLM fits mixed models but does not implement Kenward-Roger or Satterthwaite denominator degrees of freedom, which most SAPs specify for a repeated-measures primary endpoint. R’s mmrm package does, and matches SAS PROC MIXED.

If your primary analysis is an MMRM, run it in R or SAS. Use Python for the exploratory work around it.

Common mistakes

Mistake	Consequence	Fix
`sklearn` for inference	No p-values; hidden regularisation	`statsmodels`
Mixing `ddof=0` and `ddof=1`	Inconsistent SDs in one table	Be explicit
Testing variances then choosing a test	Inflated type I error	Always use Welch
`.dropna()` on paired data without reporting	Silent complete-case bias	Report the count; consider MMRM
Chi-squared with expected counts < 5	Unreliable p-value	Check `expected`, use Fisher
`CNSR` passed as the event indicator	Survival curve inverted	`event = 1 - CNSR`
Cox without checking PH	Misleading single hazard ratio	`check_assumptions()`
Unadjusted multiple tests	Inflated false positives	Pre-specify the strategy
Suppressing a separation warning	Meaningless coefficients	Investigate the data

Exercise 10.1 — A complete efficacy analysis

For a continuous endpoint, produce: descriptive statistics by arm, a Welch t-test with the difference and 95% CI, an ANCOVA adjusting for baseline, and a formatted results row. Explain why ANCOVA is preferred to a simple t-test on change.

Show solution

import pandas as pd
import numpy as np
from scipy import stats
import statsmodels.formula.api as smf


def efficacy_analysis(
    data: pd.DataFrame,
    endpoint: str = "CHG",
    baseline: str = "BASE",
    arm: str = "TRT01A",
    reference: str = "Placebo",
) -> dict:
    """Descriptives, unadjusted t-test and baseline-adjusted ANCOVA."""
    d = data.dropna(subset=[endpoint, baseline, arm])
    arms = [a for a in d[arm].unique() if a != reference]
    if len(arms) != 1:
        raise ValueError(f"Expected exactly one active arm, found {arms}")
    active = arms[0]

    a = d.loc[d[arm] == active,    endpoint]
    p = d.loc[d[arm] == reference, endpoint]

    # --- Descriptives -------------------------------------------------------
    desc = (d.groupby(arm)[endpoint]
              .agg(n="count", mean="mean", sd=lambda x: x.std(ddof=1),
                   median="median", min="min", max="max")
              .round(3))

    # --- Unadjusted Welch t-test --------------------------------------------
    diff = a.mean() - p.mean()
    se = np.sqrt(a.var(ddof=1) / len(a) + p.var(ddof=1) / len(p))
    df = (a.var(ddof=1)/len(a) + p.var(ddof=1)/len(p))**2 / (
        (a.var(ddof=1)/len(a))**2/(len(a)-1) + (p.var(ddof=1)/len(p))**2/(len(p)-1))
    crit = stats.t.ppf(0.975, df)
    t_stat, p_unadj = stats.ttest_ind(a, p, equal_var=False)

    unadjusted = {
        "difference": diff,
        "ci": (diff - crit*se, diff + crit*se),
        "p_value": p_unadj,
    }

    # --- ANCOVA -------------------------------------------------------------
    d = d.assign(_arm=pd.Categorical(d[arm],
                                     categories=[reference, active]))
    model = smf.ols(f"{endpoint} ~ _arm + {baseline}", data=d).fit()

    term = [t for t in model.params.index if t.startswith("_arm")][0]
    ci = model.conf_int().loc[term]

    adjusted = {
        "difference": model.params[term],
        "se":         model.bse[term],
        "ci":         (ci[0], ci[1]),
        "p_value":    model.pvalues[term],
        "baseline_coef": model.params[baseline],
        "r_squared":  model.rsquared,
    }

    return {"descriptives": desc, "unadjusted": unadjusted,
            "adjusted": adjusted, "model": model}


res = efficacy_analysis(adlb.query("PARAMCD == 'ALT' and AVISITN == 12"))

print(res["descriptives"])
print(f"\nUnadjusted difference: {res['unadjusted']['difference']:.2f} "
      f"(95% CI {res['unadjusted']['ci'][0]:.2f} to {res['unadjusted']['ci'][1]:.2f}), "
      f"p = {res['unadjusted']['p_value']:.4f}")
print(f"ANCOVA difference:     {res['adjusted']['difference']:.2f} "
      f"(95% CI {res['adjusted']['ci'][0]:.2f} to {res['adjusted']['ci'][1]:.2f}), "
      f"p = {res['adjusted']['p_value']:.4f}")

              n    mean      sd  median     min    max
TRT01A
Drug A      102  -4.213   9.847  -3.900  -31.20  22.40
Placebo      98  -1.104  10.221  -1.050  -28.70  25.10

Unadjusted difference: -3.11 (95% CI -6.01 to -0.21), p = 0.0356
ANCOVA difference:     -3.34 (95% CI -5.42 to -1.26), p = 0.0018

A formatted table row:

def format_row(res: dict, label: str) -> str:
    a = res["adjusted"]
    return (f"{label:<30} "
            f"{a['difference']:>7.2f} "
            f"({a['ci'][0]:>6.2f}, {a['ci'][1]:>6.2f}) "
            f"{a['p_value']:>8.4f}")

print(format_row(res, "ALT change at Week 12"))
#> ALT change at Week 12            -3.34 ( -5.42,  -1.26)   0.0018

Why ANCOVA rather than a t-test on change

Notice the confidence interval narrowed from ±2.90 to ±2.08 and the p-value dropped by a factor of twenty — on the same data. Three reasons:

Precision. Baseline is correlated with change, so including it explains variance that would otherwise be residual error. The standard error shrinks and power rises, with no cost and no assumption beyond linearity.
Chance baseline imbalance. Randomisation balances baseline in expectation, not in any particular trial. If the active arm happened to start higher, some of the observed change is regression to the mean. ANCOVA adjusts for that; a t-test on change does not.
Regulatory expectation. The EMA guideline on baseline covariates recommends adjusting for baseline when it is prognostic, and most SAPs specify ANCOVA for exactly this endpoint shape.

Two cautions on the code

pd.Categorical(categories=[reference, active]) sets the reference level explicitly. Without it, statsmodels orders alphabetically — so "Drug A" would become the reference and the sign of the treatment effect would flip. That is a silent, plausible-looking error.

.dropna(subset=[endpoint, baseline, arm]) is again a complete-case analysis. Report the excluded count, and if missingness is substantial the SAP should specify a mixed model instead:

n_excluded = len(data) - len(d)
if n_excluded:
    print(f"Note: {n_excluded} records excluded for missing endpoint or baseline")

Exercise 10.2 — Survival analysis with a proper summary

Produce a Kaplan-Meier analysis: median survival with CIs per arm, survival probabilities at fixed timepoints, a log-rank test, a Cox hazard ratio, and a check of the proportional hazards assumption.

Show solution

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


def survival_analysis(
    adtte: pd.DataFrame,
    paramcd: str = "OS",
    arm: str = "TRT01A",
    reference: str = "Placebo",
    timepoints: tuple = (90, 180, 365),
) -> dict:
    """Kaplan-Meier, log-rank and Cox analysis of a time-to-event parameter."""
    d = adtte.query("PARAMCD == @paramcd").copy()

    # CDISC CNSR: 1 = censored. lifelines wants 1 = event.
    if "CNSR" not in d.columns:
        raise ValueError("Expected a CNSR column following CDISC conventions")
    d["event"] = 1 - d["CNSR"]

    if not d["event"].isin([0, 1]).all():
        raise ValueError("CNSR must be 0 or 1")
    n_events = int(d["event"].sum())
    if n_events == 0:
        raise ValueError("No events observed — survival analysis is not possible")

    # --- Per-arm Kaplan-Meier ----------------------------------------------
    per_arm, fitters = [], {}
    for a, g in d.groupby(arm):
        kmf = KaplanMeierFitter(label=a).fit(g["AVAL"], g["event"])
        fitters[a] = kmf

        ci = median_survival_times(kmf.confidence_interval_)
        surv = {f"surv_{t}d": float(kmf.predict(t)) for t in timepoints}

        per_arm.append({
            arm: a,
            "n": len(g),
            "events": int(g["event"].sum()),
            "censored": int((g["event"] == 0).sum()),
            "median": kmf.median_survival_time_,
            "median_ci_lower": ci.iloc[0, 0],
            "median_ci_upper": ci.iloc[0, 1],
            **surv,
        })
    summary = pd.DataFrame(per_arm)

    # --- Log-rank -----------------------------------------------------------
    arms = [a for a in d[arm].unique() if a != reference]
    if len(arms) != 1:
        raise ValueError(f"Expected one active arm besides {reference}, got {arms}")
    active = arms[0]

    ga, gp = d[d[arm] == active], d[d[arm] == reference]
    lr = logrank_test(ga["AVAL"], gp["AVAL"], ga["event"], gp["event"])

    # --- Cox ----------------------------------------------------------------
    cox_df = d[["AVAL", "event", arm]].copy()
    cox_df["treated"] = (cox_df[arm] == active).astype(int)
    cox_df = cox_df.drop(columns=arm)

    cph = CoxPHFitter().fit(cox_df, duration_col="AVAL", event_col="event")
    hr = float(cph.hazard_ratios_["treated"])
    hr_ci = np.exp(cph.confidence_intervals_.loc["treated"]).values

    # --- Proportional hazards check -----------------------------------------
    try:
        ph = cph.check_assumptions(cox_df, p_value_threshold=0.05, show_plots=False)
        ph_ok = True
    except Exception:
        ph_ok = None

    return {
        "summary": summary,
        "n_events": n_events,
        "logrank": {"statistic": lr.test_statistic, "p_value": lr.p_value},
        "cox": {"hazard_ratio": hr, "ci": tuple(hr_ci),
                "p_value": float(cph.summary.loc["treated", "p"])},
        "ph_assumption_checked": ph_ok,
        "fitters": fitters,
        "model": cph,
    }

res = survival_analysis(adtte)

print(res["summary"].round(3).to_string(index=False))
print(f"\nLog-rank: chi2 = {res['logrank']['statistic']:.3f}, "
      f"p = {res['logrank']['p_value']:.4f}")
print(f"Cox HR  = {res['cox']['hazard_ratio']:.3f} "
      f"(95% CI {res['cox']['ci'][0]:.3f} to {res['cox']['ci'][1]:.3f}), "
      f"p = {res['cox']['p_value']:.4f}")

 TRT01A   n  events  censored  median  median_ci_lower  median_ci_upper  surv_90d  surv_180d  surv_365d
 Drug A  84      31        53   412.0            356.0            498.0     0.940      0.833      0.548
Placebo  86      47        39   287.0            241.0            339.0     0.895      0.721      0.395

Log-rank: chi2 = 6.412, p = 0.0113
Cox HR  = 0.612 (95% CI 0.389 to 0.963), p = 0.0337

The plot, with the at-risk table that a reviewer will expect:

import matplotlib.pyplot as plt
from lifelines.plotting import add_at_risk_counts

fig, ax = plt.subplots(figsize=(9, 7))
for a, kmf in res["fitters"].items():
    kmf.plot_survival_function(ax=ax, ci_show=True)

add_at_risk_counts(*res["fitters"].values(), ax=ax)
ax.set_xlabel("Time from first dose (days)")
ax.set_ylabel("Survival probability")
ax.set_ylim(0, 1)
ax.text(0.55, 0.9,
        f"HR = {res['cox']['hazard_ratio']:.2f} "
        f"({res['cox']['ci'][0]:.2f}–{res['cox']['ci'][1]:.2f})\n"
        f"Log-rank p = {res['logrank']['p_value']:.4f}",
        transform=ax.transAxes, fontsize=10,
        bbox=dict(boxstyle="round", facecolor="white", edgecolor="grey"))
plt.tight_layout()

Four things this gets right

event = 1 - CNSR, asserted. The single most common error in Python survival code for clinical data. Passing CNSR directly produces a curve that rises, or — worse — one that looks reasonable and is exactly backwards.

A median with a confidence interval. A median survival of 412 days means little without knowing whether the interval is 400–420 or 356–498. Reporting the point estimate alone is a common omission.

Survival at fixed timepoints. When the median is not reached in one arm — common in oncology with good outcomes — median_survival_time_ returns inf and the fixed-timepoint estimates are all you can report.

Checking proportional hazards. If the assumption fails, the hazard ratio is an average of an effect that changes over time. Reporting HR = 0.61 without that caveat is misleading when, say, the curves cross at nine months. Look at the Schoenfeld residual plots, not just the p-value.

Zero events is guarded explicitly. In an interim analysis with few events lifelines produces confusing errors deep in the fit; failing early with a clear message is better.

Recap

statsmodels for inference, sklearn for prediction — they are not interchangeable
pandas defaults to ddof=1, NumPy to ddof=0; be explicit
Welch’s t-test as the default; never test variances first and then choose
Check expected counts before trusting chi-squared; Fisher for sparse 2×2
CNSR is inverted relative to lifelines: event = 1 - CNSR
Always check proportional hazards before reporting a single hazard ratio
Set the categorical reference level explicitly or the effect sign may flip
MMRM with Kenward-Roger is not available in Python — use R or SAS

Next: Machine learning.