Bio Clinical Biostatistics Logistic Regression

Version Compatibility

Reference examples tested with: statsmodels 0.14+, scipy 1.12+, numpy 1.26+, pandas 2.1+, firthmodels 0.3+, marginaleffects 0.0.13+ (Python) / 0.20+ (R). R packages cited: RobinCar, marginaleffects, brant, MASS, VGAM.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Logistic Regression for Clinical Outcomes

"Model clinical outcomes with logistic regression" -> Estimate the marginal or conditional treatment effect on a binary or ordinal endpoint using a model that respects randomisation stratification, declares its estimand, and survives covariate misspecification.

Conditional vs marginal (the non-collapsibility subtlety): the OR is non-collapsible. The conditional OR from logistic regression is a different parameter than the marginal OR, even when there is NO confounding and randomisation is perfect. This is mathematical, not statistical bias. FDA 2023 favours marginal RD (via g-computation) for primary reporting to avoid parameter ambiguity. See clinical-biostatistics/effect-measures and Permutt 2020.

Algorithmic Taxonomy

Postdoc reading list: Permutt 2020 Stat Biopharm Res 12:45 (conditional vs marginal estimand); FDA May 2023 Final Guidance "Adjusting for Covariates in RCTs" (the regulatory rulebook); Tsiatis et al 2008 Stat Med 27:4658 (robust ANCOVA framework); Senn 2013 Stat Med 32:1439 (precision from adjustment even under perfect balance); Kahan-Morris 2012 Stat Med 31:328 (must adjust for stratification factors).

Decision Tree by Scenario

Standard Workflow

Goal: Fit a covariate-adjusted logistic regression for a binary clinical endpoint with explicit reference category and pre-specified covariates.

Approach: Use the formula API for automatic categorical handling and explicit reference; report both the conditional log-OR and the marginal RD via g-computation.

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

# CRITICAL: set reference category explicitly. statsmodels default is alphabetical,
# so 'Active' sorts before 'Placebo' and the OR direction silently flips.
model = smf.logit(
    'outcome ~ C(ARM, Treatment(reference="Placebo")) + age + C(sex) + baseline_score',
    data=df
).fit()

# Conditional ORs (Wald)
or_table = pd.DataFrame({
    'OR': np.exp(model.params),
    'Lower_CI': np.exp(model.conf_int()[0]),
    'Upper_CI': np.exp(model.conf_int()[1]),
    'p_value': model.pvalues
})

# Marginal RD via g-computation (the FDA 2023-recommended primary estimand):
df_active = df.assign(ARM='Active')
df_placebo = df.assign(ARM='Placebo')
risk_active = model.predict(df_active).mean()
risk_placebo = model.predict(df_placebo).mean()
marginal_rd = risk_active - risk_placebo
# For SE, use influence-function or bootstrap; see RobinCar / marginaleffects in R

The reference-category trap is the single most common silent bug. smf.logit('y ~ C(ARM)') uses alphabetical ordering, so ARM='Active' becomes 1 and ARM='Placebo' becomes 0 -- but with ARM=['Active','Placebo','Control'], 'Active' is the reference and 'Placebo' is the comparison. Always pass Treatment(reference="Placebo") or equivalent. R glm(y ~ relevel(ARM, ref='Placebo')).

Marginal vs Conditional Estimand -- The FDA 2023 Pivot

The single most important methodological shift in clinical biostatistics 2020-2025. Permutt 2020 established and FDA 2023 codified: the maximum-likelihood coefficient on Z in glm(Y ~ Z + X, family=binomial) is the conditional log-OR -- the OR comparing Z=1 to Z=0 holding X fixed. Due to OR non-collapsibility, this is a different parameter than the marginal log-OR (the OR comparing the entire treated population to the entire control population), even under perfect randomisation.

The FDA 2023 final guidance ("Adjusting for Covariates in Randomized Clinical Trials," May 2023) endorses covariate-adjusted nonlinear-model analysis provided the analyst targets a marginal estimand. Reporting the conditional OR from a multivariable logistic regression as "the treatment effect" silently changes the estimand from what was pre-specified.

G-computation / Standardisation

# Marginal effects on three scales via g-computation:
df_z1 = df.assign(ARM='Active')
df_z0 = df.assign(ARM='Placebo')
p_z1 = model.predict(df_z1)
p_z0 = model.predict(df_z0)

marg_p1 = p_z1.mean()
marg_p0 = p_z0.mean()
marg_rd = marg_p1 - marg_p0
marg_rr = marg_p1 / marg_p0
marg_or = (marg_p1 / (1 - marg_p1)) / (marg_p0 / (1 - marg_p0))

For valid SE: the delta-method/influence-function variance and HC3 sandwich are required. Python marginaleffects v0.0.13+ implements this via marginaleffects.avg_comparisons(model, variables='ARM', vcov='HC3'). R is more mature: marginaleffects::avg_comparisons (Arel-Bundock-Greifer-Heiss 2024 JSS 111:9), RobinCar (purpose-built for FDA 2023), riskCommunicator.

Tsiatis et al 2008 robustness guarantee: under randomisation Z ⊥ X, the g-computation marginal estimator is consistent for the marginal ATE even if the outcome model is misspecified. Efficiency depends on model quality; consistency does not. This is why FDA 2023 accepts the marginal RD via g-computation without requiring proof of correct logistic mean structure.

Reporting template (post-FDA-2023):

"Primary estimand: marginal risk difference of -8.5 percentage points (95% CI -12.7 to -4.2, HC3 SE), computed by g-computation/standardisation from a logistic regression adjusted for age, sex, and baseline severity. Supportive: conditional OR 0.48 (95% CI 0.34-0.67, Wald) from the same model."

Modified Poisson for Common Outcomes -- Direct RR

When prevalence > 10% and the policy quantity is RR (not OR), modified Poisson with sandwich SE is the de facto modern standard (Zou 2004 AJE 159:702):

import statsmodels.api as sm
import numpy as np

X = sm.add_constant(df[['treatment', 'age', 'baseline_score']])
poisson_model = sm.GLM(df['outcome'], X, family=sm.families.Poisson()).fit(cov_type='HC1')
rr = np.exp(poisson_model.params)
rr_ci = np.exp(poisson_model.conf_int())

cov_type='HC1' corrects the over-dispersion that Poisson inherently assumes. Without it the variance is wrong because binary data are NOT Poisson. For n < 250, HC3 is preferred (Long-Ervin 2000). What postdocs argue about: modified Poisson vs log-binomial -- log-binomial directly models RR but frequently fails to converge when predicted risk approaches 1; modified Poisson always converges but is marginally less efficient when log-binomial works.

Proportional Odds and Brant Test

Proportional odds (PO) assumption: the effect of each predictor is constant across all cut-points of the ordinal outcome. Must be tested or the model parameters are invalid.

from statsmodels.miscmodels.ordinal_model import OrderedModel
from statsmodels.api import MNLogit
import pandas as pd

df['severity'] = pd.Categorical(df['severity'], categories=['mild', 'moderate', 'severe'], ordered=True)

po_model = OrderedModel.from_formula('severity ~ treatment + age', data=df, distr='logit').fit(method='bfgs', disp=0)
mn_model = MNLogit.from_formula('severity ~ treatment + age', data=df).fit(disp=0)

# LR test for PO vs MN (omnibus)
lr_stat = 2 * (mn_model.llf - po_model.llf)
lr_df = mn_model.df_model - po_model.df_model
from scipy.stats import chi2
lr_p = 1 - chi2.cdf(lr_stat, lr_df)

# Per-coefficient Brant test in R: brant::brant(po_model_object)
# The Brant test localises which predictor(s) violate PO -- essential before deciding on partial PO

The omnibus LR test is necessary but not sufficient. Brant test (Brant 1990 Biometrics 46:1171; R brant::brant) provides per-coefficient PO tests, identifying which predictor violates PO. If only one or two predictors violate PO, fit a partial proportional odds model (R VGAM::vglm(..., cumulative(parallel=FALSE~X1+X2))) rather than abandoning the model entirely.

OrderedModel intercept gotcha: do NOT add an intercept term. Threshold parameters (cut-points between ordinal levels) replace the intercept. An explicit constant causes non-identifiability and optimiser failure.

Separation and Firth Penalty

Detection:

from firthmodels import FirthLogisticRegression, detect_separation
import numpy as np

sep_result = detect_separation(X, y)
if sep_result.separation:
    print(sep_result.summary())
# Manual signs: coefficient > 10, SE > 100, convergence warnings

Firth penalty (Firth 1993 Biometrika 80:27):

firth = FirthLogisticRegression()
firth.fit(X, y)
or_firth = np.exp(firth.coef_)

# IMPORTANT: Wald p-values from Firth are LIBERAL (anti-conservative).
# Prefer the penalised likelihood-ratio test (PLRT).
# Some Python `firthlogist` releases expose a PLRT attribute (e.g. `pvalues_lrt_`)
# but its presence and name vary by release -- check `dir(firth)` against the
# installed version, otherwise compute the PLRT manually from the penalised
# log-likelihoods of nested models.

Heinze-Schemper 2002 Stat Med 21:2409 showed Wald inference from Firth is liberal (anti-conservative); the penalised likelihood-ratio test (PLRT) is the recommended inference. PLRT attributes (e.g. pvalues_lrt_) appear in some Python Firth packages but the exact attribute name varies by release -- introspect the installed package; compute the PLRT manually if no attribute is exposed:

def penalised_lrt(firth_full, firth_reduced):
    # Compute 2 * (penalised log-lik full - penalised log-lik reduced) ~ chi-square_df
    pass  # implementation depends on package version; see Heinze-Schemper 2002

Firth's method was originally designed for finite-sample bias reduction, not separation per se -- it adds the Jeffreys prior penalty to the likelihood, keeping coefficients finite under separation as a side effect. Also recommended for rare events (<5% prevalence) where ML bias is non-negligible.

Hauck-Donner Effect Detection

The Hauck-Donner effect (1977 JASA 72:851; revived by Yee 2022 JASA 117:1763): the Wald test statistic is non-monotonic in the parameter estimate near the boundary. A large log-OR can produce a small Wald chi-square -- so the Wald test fails to reject when LR/profile-likelihood would. Common in small samples with strong predictors.

# In R, detect with VGAM::hdeff() and replace Wald with profile likelihood:
# MASS::confint.glm(model) returns profile-likelihood CIs as default
# Python equivalent: bootstrap or manual profile likelihood

When to suspect Hauck-Donner: large coefficient magnitude with non-significant Wald p; large SE with finite estimate; switching from Wald to LR test changes significance. The fix is always: switch to profile likelihood or LR inference.

Reconciliation: When Methods Disagree

Per-Method Failure Modes

Reference-category silent reversal

• Trigger: Treatment variable coded as character with Active and Placebo; reference not explicitly set.
• Mechanism: statsmodels defaults to alphabetical ordering; 'Active' becomes the reference.
• Symptom: OR for "ARM[T.Placebo]" appears in output; direction is the opposite of intended.
• Fix: Always pass C(ARM, Treatment(reference="Placebo")) or numeric encoding with explicit comment.

Adjusting for a mediator

• Trigger: Including a post-randomisation variable that lies on the causal pathway.
• Mechanism: Adjustment blocks the causal effect, attenuating treatment effect estimate toward null.
• Symptom: Significant unadjusted effect becomes non-significant after adjustment for "mechanism" variable (e.g., inflammation marker).
• Fix: Use causal DAG to distinguish confounder from mediator; never adjust for post-treatment variables in the primary analysis. See ICH E9(R1).