Generalized Linear Models (GLM)

In plain language

A Generalized Linear Model (GLM) is a flexible regression framework that lets analysts pick both the shape of the outcome distribution (called the "family") and the scale on which the treatment effect is measured (called the "link function"), so the model fits binary events, count data, and skewed healthcare costs with equal ease. The family choice answers "how does variance grow with the mean?" while the link choice answers "should the effect be measured as a difference, a ratio of means, or an odds ratio?" — these are two separate decisions even though most textbooks bundle them together. Getting both right is the single most consequential statistical modelling choice in an HEOR analysis, because the same dataset will give a risk difference, a risk ratio, or an odds ratio depending solely on which link you specify.

What a GLM is: three components that click together

Ordinary least squares (OLS) regression assumes that the outcome Y is continuous, has constant variance, and that residuals are approximately normally distributed. Most outcomes in real-world evidence (RWE) and health economics and outcomes research (HEOR) violate at least one of those assumptions: binary endpoints (hospitalization, response, mortality) are Bernoulli-distributed; healthcare costs are right-skewed with a spike at zero and variance that grows with the mean; prescription fill counts are non-negative integers. Forcing these outcomes into an OLS framework produces systematically biased predictions (fitted probabilities outside [0,1] for binary outcomes; negative predicted costs) and invalid standard errors.

A Generalized Linear Model (GLM) solves this by specifying three components explicitly and linking them together:

1. The random component (distributional family): The assumed conditional distribution of Y given the covariates — Binomial, Poisson, Negative Binomial, Gamma, Gaussian, and inverse Gaussian are the families used most commonly in HEOR. The family determines the mean-variance relationship. Gaussian assumes Var(Y) = σ² (constant). Poisson assumes Var(Y) = μ (variance equals mean). Negative Binomial adds an overdispersion parameter k so Var(Y) = μ + μ²/k. Gamma assumes Var(Y) = φμ² (variance proportional to squared mean, appropriate for right-skewed positive costs). Binomial assumes Var(Y) = μ(1-μ)/n.

2. The systematic component (linear predictor): η = β₀ + β₁X₁ + β₂X₂ + … + βₚXₚ. The same ordinary weighted sum of predictors, including continuous variables, indicators, interactions, and offset terms (log of person-time for rate models), that appears in OLS.

3. The link function: g(μ) = η. The link function maps the expected value of Y (μ) to the linear predictor (η). The choice of link determines the scale on which the additive linear predictor operates and therefore the estimand the model's coefficients represent. The model is estimated by iteratively reweighted least squares (IRLS), the algorithm Nelder and Wedderburn (1972) showed solves the maximum likelihood score equations for the entire GLM family in a single unified pass.

Canonical vs chosen links — and why they are independent decisions

Each distributional family has a canonical link — the link for which the IRLS algorithm simplifies and the score equations have closed-form sufficient statistics. Logit is canonical for Binomial; log is canonical for Poisson and Gamma; identity is canonical for Gaussian. Canonical links generally give the best numerical stability and the fastest convergence.

However, the link is a separate decision from the family. An analyst fitting a binomial GLM can choose logit (canonical; produces ORs), log (produces RRs but can produce predicted probabilities above 1 for common exposures), or identity (produces risk differences; requires constrained estimation to keep predictions in [0,1]). Choosing a non-canonical link is legitimate and often scientifically preferable when the effect measure on the canonical scale (an OR) is less meaningful for decision-making than the alternative (an RR or RD).

The two decisions — family and link — are made for different reasons and analysts too often conflate them. The family is a statement about variance structure: choosing Gamma asserts that the variance of the outcome grows with the square of the mean, which is the right structural assumption for healthcare costs. The link is a statement about effect-measure scale: choosing a log link asserts that the covariate acts multiplicatively on the mean (a one-unit change multiplies the mean by exp(β)), which may be the clinically or economically natural scale regardless of the family. You can run a Gamma family with log link (the standard for costs), a Gamma family with identity link (less common), or a Gamma family with inverse link (the canonical; typically least useful for HEOR). The family and link choices are orthogonal.

Interpreting the output

The link function is the master key for coefficient interpretation. Consider a single binary covariate (exposed vs unexposed) with estimated coefficient β = 0.693 (approximately ln 2) under three different links for the same binomial data:

Identity link (linear predictor = μ directly): - Formal interpretation: β = 0.693 is the conditional additive difference in the probability of the event for exposed vs unexposed patients, holding all other covariates fixed. The 95% CI on β is a CI on the risk difference. The scale is the original probability scale. - Practical interpretation: "Exposed patients are 0.693 percentage points — approximately 69 percentage points — more likely to experience the event than unexposed patients." (If both groups have plausible baseline probabilities, the identity-scale estimate is the most direct input to number-needed-to-treat or budget-impact calculations.)

Log link (linear predictor = log μ): - Formal interpretation: exp(0.693) ≈ 2.0. The exponentiated coefficient is the conditional ratio of expected outcomes (for Binomial: the risk ratio; for Poisson/NB: the rate ratio; for Gamma: the cost ratio or mean utilization ratio). The CI on β maps, after exponentiation, to a CI on the multiplicative ratio. - Practical interpretation: "Exposed patients have approximately twice the risk (rate, mean cost) of the outcome compared with unexposed patients, holding covariates fixed." The log link is the most common choice in HEOR because multiplication is a natural model for how disease exposure, drug effects, and comorbidity burden combine.

Logit link (linear predictor = log(μ/(1-μ))): - Formal interpretation: exp(0.693) ≈ 2.0. The exponentiated coefficient is the conditional odds ratio — the ratio of the odds of the event in exposed vs unexposed patients, holding covariates fixed. Critically, this is NOT a risk ratio; it approximates the RR only when the baseline outcome probability is below approximately 10%. - Practical interpretation: "The odds of the event are approximately twice as high for exposed patients as for unexposed patients." Odds ratios are non-collapsible: the marginal OR (across the whole population) will differ from the conditional OR (within strata), even without confounding, which creates communication challenges when presenting to payers and HTA bodies who think in terms of risks and rates rather than odds.

The choice of link should therefore be driven by the effect measure required for decision-making, not by estimation convenience alone. Risk differences feed into number-needed-to-treat and budget-impact models. Risk ratios and rate ratios are natural inputs to cost-effectiveness models. Odds ratios are required when the study design only identifies odds (case-control) or when a meta-analytic literature reports ORs.

The RWE family/link decision table

• Binary endpoint (hospitalization, mortality, treatment response): Binomial family is
• Count and rate outcomes (hospitalizations per year, prescriptions, office visits):
• Healthcare costs and positive continuous skewed outcomes: Gamma family with log link is
• Zero-heavy semicontinuous costs (many patients with $0 spend): Neither OLS nor a

Quasi-likelihood and robust standard errors

When the chosen distributional family does not fully match the data's variance structure — for example, when counts are overdispersed relative to Poisson — quasi-likelihood (Wedderburn 1974) extends the GLM framework by specifying only the mean-variance relationship (V(μ)) rather than the full distribution. Quasi-Poisson and quasi-binomial keep the mean model intact but scale the standard errors by a dispersion factor estimated from the Pearson chi-square statistic, inflating SE when data are overdispersed. This is simpler to implement than switching to Negative Binomial but does not support AIC-based model comparison (AIC is undefined for quasi-likelihoods). The alternative for robust inference without distributional assumptions is the sandwich (Huber-White) variance estimator, which is available in all three major statistical languages and is widely used in HEOR when the variance assumption cannot be fully verified.

Deviance, AIC, and model checking

Deviance is the GLM analogue of residual sum of squares: it measures the discrepancy between the fitted model and the saturated model. Residual deviance on (n - p) degrees of freedom should be approximately 1.0 for a well-specified family; a residual deviance/df ratio substantially above 1.0 signals overdispersion. AIC (Akaike Information Criterion) allows comparison of GLMs with the same family but different linear predictors or link functions; it cannot compare models with different families (e.g., Poisson vs Gamma). For residual diagnostics, the concepts in regression-diagnostics provide the Pearson and deviance residual checks that confirm family-link adequacy.

Pros, cons, and trade-offs

Pros: - Unified framework: one fitting algorithm (IRLS), one software call (glm/PROC GENMOD), handles binary, count, and continuous skewed outcomes with the same code pattern, making the framework accessible and results comparable across outcome types. - Effect-measure flexibility: the analyst controls the estimand (OR, RR, RD, rate ratio, cost ratio) by choosing the link, rather than being forced into the effect measure implied by a transformation. - Correct variance structure: by choosing the right family, predictions respect the support of the outcome (non-negative costs, bounded probabilities) and SEs are calibrated to the true variance shape. - Natural accommodation of offsets and weights for propensity-score analyses, standardization, and rate models.

Cons and limitations: - Family misspecification is silent: a GLM fits even when the family is wrong; only diagnostic plots and the Park test reveal poor fit. Analysts who never check residuals produce plausible-looking but miscalibrated results. - Non-canonical links can fail to converge, especially logit-vs-log for binary outcomes with common events. Convergence issues are a practical barrier in production pipelines. - The conditional mean is the target estimand; for decision-making that requires marginal (population-average) effects, recycled predictions or average marginal effects must be computed post-estimation. - A single GLM is not the right tool for semicontinuous cost distributions (mass at zero plus right-skewed positive values); two-part models or Tweedie distributions are needed. - Interpretation of the logit-link OR for common outcomes is routinely misunderstood as a risk ratio, which inflates perceived effect sizes and misleads decision-makers.

When to use

• Binary outcomes: whenever logistic regression (OR) or a modified Poisson approach (RR)
• Healthcare cost outcomes: Gamma GLM with log link as the primary analysis whenever the
• Count and utilization outcomes: Poisson with log link and offset as the starting model;
• Covariate-adjusted analyses: GLMs accommodate propensity-score weighting (via the
• When an interpretable effect measure on a specific scale (RR, cost ratio) is required:

When NOT to use

• Censored time-to-event outcomes: survival models (Cox proportional hazards, parametric
• *Repeated measures, clustered, or longitudinal outcomes without accounting for
• Zero-heavy semicontinuous cost distributions: when a substantial fraction of patients
• *When an OLS linear predictor on the outcome scale is adequate and the CLT protects
• As a substitute for causal identification: a GLM adjusts for measured confounders but

Worked example

Scenario

A claims-based cohort study compares a binary outcome (30-day hospital readmission) between 100 patients who received a new care coordination program (exposed) and 100 patients who received standard care (unexposed). The analyst wants to estimate three different effect measures from the same 2x2 event table — the risk difference (RD), the risk ratio (RR), and the odds ratio (OR) — by fitting the same Binomial GLM with three different link functions: identity, log, and logit. This illustrates that the choice of link function, not the data, determines which effect measure falls out of the regression coefficient.

Dataset

Summary 2x2 table: readmission events out of 100 patients per group. The individual-level dataset has 200 rows (one per patient) with outcome = 1 (readmitted) or 0 (not readmitted).

Steps

• Compute group proportions. Exposed: 30/100 = 0.30. Unexposed: 15/100 = 0.15.

• Identity link (linear predictor acts directly on the probability μ). The regression coefficient for the exposed indicator is the additive difference in proportions: RD = 0.30 - 0.15 = 0.15. Interpretation: exposed patients have a 15 percentage-point higher absolute probability of readmission.

• Log link (linear predictor acts on log μ). The regression coefficient for exposed is log(RR), so the exponentiated coefficient is the risk ratio: RR = 0.30/0.15 = 2.0. Interpretation: exposed patients are twice as likely to be readmitted.

• Logit link (linear predictor acts on log(μ/(1-μ))). Compute odds in each group. Odds exposed = 30/70. Odds unexposed = 15/85. OR = (30/70) / (15/85) = (3085) / (7015) = 2550/1050 = 2.43. Interpretation: the odds of readmission are 2.43 times higher for exposed patients. Note that 2.43 ≠ 2.0 (the RR); odds ratios overstate risk ratios when the baseline rate is not rare (here 15% is not rare).

• All three models use the same 200-observation dataset and the same linear predictor (intercept + exposed indicator). The link function alone determines whether the coefficient is a RD, log(RR), or log(OR). The decision should be made based on which effect measure the decision-maker needs, not on which link converges more easily.

Result

Identity link: RD = 0.30 - 0.15 = 0.15 (15 pp higher absolute risk for exposed group). Log link: RR = 0.30/0.15 = 2.0 (exposed have twice the readmission rate). Logit link: OR = 3085/(7015) = 2550/1050 = 2.43 (odds 2.43x higher for exposed). The three models share the same data; the coefficient scale changes with the link.

Runnable example

import numpy as np
import statsmodels.api as sm

# ── Dataset: 200 patients, binary readmission outcome ──────────────────────────
# Exposed group: 30 events in 100 patients
# Unexposed group: 15 events in 100 patients
outcome  = np.array([1]*30 + [0]*70 + [1]*15 + [0]*85)
exposed  = np.array([1]*100 + [0]*100)
X        = sm.add_constant(exposed.astype(float))

# ── 1. Binomial + identity link → Risk Difference ──────────────────────────────
try:
    m_id = sm.GLM(outcome, X,
                  family=sm.families.Binomial(link=sm.families.links.Identity())
                  ).fit(method="newton")
    rd = m_id.params["x1"]
    ci = m_id.conf_int().loc["x1"]
    print(f"Identity link  RD  = {rd:.4f}  95% CI [{ci[0]:.4f}, {ci[1]:.4f}]")
    # Expected: RD ≈ 0.15 (= 0.30 - 0.15)
except Exception as e:
    print(f"Identity link did not converge: {e}  -- use constrained estimation")

# ── 2. Binomial + log link → Risk Ratio ────────────────────────────────────────
# Log-link binomials can fail for high event rates; add starting values or use
# modified Poisson (Poisson + log + robust SE) as a pragmatic RR estimator.
try:
    m_log = sm.GLM(outcome, X,
                   family=sm.families.Binomial(link=sm.families.links.Log())
                   ).fit(method="bfgs", start_params=[np.log(0.15), np.log(2)])
    rr = np.exp(m_log.params["x1"])
    ci_rr = np.exp(m_log.conf_int().loc["x1"])
    print(f"Log link       RR  = {rr:.4f}  95% CI [{ci_rr[0]:.4f}, {ci_rr[1]:.4f}]")
    # Expected: RR = 2.0
except Exception as e:
    print(f"Log-binomial did not converge: {e}")
    # Fallback: modified Poisson (Poisson + log + robust SE) — widely used for RR
    m_mp = sm.GLM(outcome, X,
                  family=sm.families.Poisson(link=sm.families.links.Log())
                  ).fit(cov_type="HC3")
    rr_mp = np.exp(m_mp.params["x1"])
    print(f"Modified Poisson RR = {rr_mp:.4f}  (Poisson + log + robust SE)")

# ── 3. Binomial + logit link → Odds Ratio (canonical, always converges) ────────
m_logit = sm.GLM(outcome, X,
                 family=sm.families.Binomial(link=sm.families.links.Logit())
                 ).fit()
log_or = m_logit.params["x1"]
or_est = np.exp(log_or)
ci_or  = np.exp(m_logit.conf_int().loc["x1"])
print(f"Logit link     OR  = {or_est:.4f}  95% CI [{ci_or[0]:.4f}, {ci_or[1]:.4f}]")
# Expected: OR = (30*85)/(70*15) = 2550/1050 ≈ 2.43

# ── 4. Poisson + log link → Rate Ratio (count/rate outcomes) ───────────────────
# Simulate monthly utilization counts (overdispersed; use quasi or NB in practice)
rng    = np.random.default_rng(42)
counts = rng.poisson(lam=np.where(exposed == 1, 3.0, 1.5))  # rate ratio = 2.0
offset = np.log(np.ones(200))                                # log(person-time)
X_off  = sm.add_constant(exposed.astype(float))
m_pois = sm.GLM(counts, X_off,
                family=sm.families.Poisson(link=sm.families.links.Log()),
                offset=offset).fit()
print(f"Poisson+log    Rate ratio = {np.exp(m_pois.params['x1']):.4f}")

# ── 5. Gamma + log link → Cost Ratio (healthcare costs) ────────────────────────
costs  = rng.gamma(shape=2, scale=np.where(exposed == 1, 6000, 4000))
m_gam  = sm.GLM(costs, sm.add_constant(exposed.astype(float)),
                family=sm.families.Gamma(link=sm.families.links.Log())).fit()
print(f"Gamma+log      Cost ratio = {np.exp(m_gam.params['x1']):.4f}")
print(f"  (Gamma residual deviance / df = {m_gam.deviance / m_gam.df_resid:.3f})")
# A residual deviance/df close to 1.0 suggests the Gamma family is well-specified.

# ── 6. Quasi-Poisson (overdispersion correction) ────────────────────────────────
# Scale SEs by sqrt(Pearson χ²/df); does not change point estimates
m_quasi = sm.GLM(counts, X_off,
                 family=sm.families.Poisson(link=sm.families.links.Log()),
                 offset=offset).fit(scale="X2")   # X2 = Pearson chi-square dispersion
print(f"Quasi-Poisson SE scale factor: {np.sqrt(m_quasi.scale):.3f}")

# ── Dataset: 200 patients, binary readmission outcome ──────────────────────────
outcome <- c(rep(1L, 30), rep(0L, 70), rep(1L, 15), rep(0L, 85))
exposed <- c(rep(1L, 100), rep(0L, 100))
df      <- data.frame(outcome = outcome, exposed = exposed)

# ── 1. Binomial + identity link → Risk Difference ──────────────────────────────
m_id <- tryCatch(
  glm(outcome ~ exposed, data = df, family = binomial(link = "identity")),
  error = function(e) NULL
)
if (!is.null(m_id)) {
  cat("Identity link  RD =", coef(m_id)["exposed"], "\n")
  cat("               95% CI:", confint(m_id)["exposed", ], "\n")
  # Expected: ≈ 0.15
} else {
  cat("Identity-binomial did not converge; use Gaussian GLM for approximate RD.\n")
}

# ── 2. Binomial + log link → Risk Ratio ────────────────────────────────────────
m_log <- tryCatch(
  glm(outcome ~ exposed, data = df, family = binomial(link = "log"),
      start = c(log(0.15), log(2))),   # starting values improve convergence
  error = function(e) NULL
)
if (!is.null(m_log)) {
  cat("Log link       RR =", exp(coef(m_log)["exposed"]), "\n")
} else {
  # Modified Poisson: Poisson + log + sandwich SE — pragmatic RR estimator
  m_mp <- glm(outcome ~ exposed, data = df, family = poisson(link = "log"))
  library(sandwich)
  se_mp <- sqrt(diag(sandwich(m_mp)))["exposed"]
  rr_mp <- exp(coef(m_mp)["exposed"])
  cat("Modified Poisson  RR =", rr_mp, " robust SE =", se_mp, "\n")
}

# ── 3. Binomial + logit link → Odds Ratio (canonical; always converges) ────────
m_logit <- glm(outcome ~ exposed, data = df, family = binomial(link = "logit"))
cat("Logit link     OR =", exp(coef(m_logit)["exposed"]), "\n")
cat("               95% CI:", exp(confint(m_logit)["exposed", ]), "\n")
# Expected: OR = (30*85)/(70*15) = 2550/1050 = 2.43

# ── 4. Poisson + log link → Rate Ratio (counts with offset) ───────────────────
set.seed(42)
counts <- rpois(200, lambda = ifelse(exposed == 1, 3.0, 1.5))
df$counts <- counts
df$lpt    <- 0   # log(person-time) = log(1) = 0 per patient
m_pois <- glm(counts ~ exposed + offset(lpt), data = df,
              family = poisson(link = "log"))
cat("Poisson+log    Rate ratio =", exp(coef(m_pois)["exposed"]), "\n")

# Overdispersion check: residual deviance / df
cat("  Residual deviance/df =", deviance(m_pois) / df.residual(m_pois), "\n")
# If >> 1, switch to quasi-Poisson or Negative Binomial

# ── 5. Quasi-Poisson (overdispersion-corrected SEs) ────────────────────────────
m_quasi <- glm(counts ~ exposed + offset(lpt), data = df,
               family = quasipoisson(link = "log"))
cat("Quasi-Poisson  Rate ratio =", exp(coef(m_quasi)["exposed"]),
    " (SE inflation factor:", sqrt(summary(m_quasi)$dispersion), ")\n")

# ── 6. Gamma + log link → Cost Ratio (healthcare costs) ────────────────────────
set.seed(1)
costs    <- rgamma(200, shape = 2, rate = 2 / ifelse(exposed == 1, 12000, 8000))
df$costs <- costs
m_gam    <- glm(costs ~ exposed, data = df, family = Gamma(link = "log"))
cat("Gamma+log      Cost ratio =", exp(coef(m_gam)["exposed"]), "\n")
cat("  Residual deviance/df =", deviance(m_gam) / df.residual(m_gam), "\n")

# ── 7. AIC comparison across links (same family, same data) ────────────────────
# AIC is valid for comparing models with the SAME family but different links/predictors.
# It CANNOT compare Binomial vs Poisson (different likelihoods).
cat("AIC logit:", AIC(m_logit), " AIC log (if converged):",
    if (!is.null(m_log)) AIC(m_log) else NA, "\n")