Multivariable Regression in Epidemiology: Logistic, Conditional Logistic, Log-Binomial & Poisson Models

In binary logistic regression, the outcome variable $Y_i \in \{0, 1\}$ represents the presence ($1$) or absence ($0$) of disease for individual $i$. Let $\mathbf{x}_i = (1, x_{i1}, x_{i2}, \dots, x_{ip})^T$ denote the vector of explanatory covariates, where $x_{i1}$ represents the primary exposure of interest and $x_{i2}, \dots, x_{ip}$ are potential confounders. The conditional probability of disease is modeled through the logistic (sigmoid) cumulative distribution function:

Applying the logit link function yields the linear predictor:

Interpretation of Regression Coefficients as Adjusted Odds Ratios

Consider the primary exposure variable $x_1$, evaluated for a dichotomous exposure ($x_1 = 1$ vs $x_1 = 0$) while holding all other adjustment covariates $\mathbf{x}_{-1} = (x_2, \dots, x_p)$ strictly constant:

• Log-odds among exposed: $\ln(\text{Odds}_1) = \beta_0 + \beta_1(1) + \sum_{j=2}^p \beta_j x_j$
• Log-odds among unexposed: $\ln(\text{Odds}_0) = \beta_0 + \beta_1(0) + \sum_{j=2}^p \beta_j x_j$

Subtracting the two equations:

Exponentiating both sides:

Thus, the exponentiated coefficient $\exp(\beta_1)$ represents the adjusted odds ratio for disease associated with exposure, holding all other covariates in the model mathematically constant. For a continuous covariate (such as age in years), $\exp(\beta_j)$ represents the multiplicative change in the odds of disease per 1-unit increase in that covariate. For a $c$-unit increase, the adjusted odds ratio is $\exp(c \cdot \beta_j)$.

Maximum Likelihood Estimation via the Newton-Raphson Algorithm

Because $Y_i$ takes values in $\{0, 1\}$, its probability mass function under the Bernoulli distribution is $P(Y_i = y_i) = \pi(\mathbf{x}_i)^{y_i} [1 - \pi(\mathbf{x}_i)]^{1 - y_i}$. Assuming independent observations across $n$ individuals, the likelihood function is:

Taking the natural logarithm yields the log-likelihood function:

Substituting $\ln(1 - \pi(\mathbf{x}_i)) = -\ln(1 + e^{\mathbf{x}_i^T \boldsymbol{\beta}})$ and $\ln \pi(\mathbf{x}_i) = \mathbf{x}_i^T \boldsymbol{\beta} - \ln(1 + e^{\mathbf{x}_i^T \boldsymbol{\beta}})$:

To maximize $\ell(\boldsymbol{\beta})$, we evaluate the first-order partial derivatives with respect to $\boldsymbol{\beta}$ (the Score Vector $\mathbf{U}(\boldsymbol{\beta})$):

Setting $\mathbf{U}(\boldsymbol{\beta}) = \mathbf{0}$ yields a system of non-linear equations with no closed-form algebraic solution. We compute the Maximum Likelihood Estimate $\widehat{\boldsymbol{\beta}}$ iteratively using the Newton-Raphson algorithm (Iteratively Reweighted Least Squares). The second partial derivative matrix is the negative Hessian Matrix (observed Fisher Information $\mathbf{I}(\boldsymbol{\beta})$):

where $\mathbf{W} = \text{diag}(\pi(\mathbf{x}_i)(1 - \pi(\mathbf{x}_i)))$ is an $n \times n$ diagonal weight matrix. The parameter vector is updated at iteration $(t + 1)$ via:

where $\mathbf{z}^{(t)} = \mathbf{X} \boldsymbol{\beta}^{(t)} + [\mathbf{W}^{(t)}]^{-1} (\mathbf{y} - \boldsymbol{\pi}^{(t)})$ is the adjusted working response vector. Upon convergence, the asymptotic covariance matrix of $\widehat{\boldsymbol{\beta}}$ is obtained directly as the inverse Fisher information:

In individually matched case-control studies (e.g., each case matched to $M$ controls on age, sex, and clinic), standard unconditional logistic regression encounters a fatal theoretical obstacle known as the Neyman-Scott Problem.

The Neyman-Scott Paradox

Suppose a study enrolls $K$ individually matched 1:1 pairs ($N = 2K$ total participants). To control for the matching strata in an unconditional logistic regression model, the investigator must include a separate indicator variable (intercept $\alpha_k$) for each matched pair:

As sample size grows by enrolling more pairs ($K \to \infty$), the number of incidental nuisance parameters $\alpha_1, \dots, \alpha_K$ increases at the exact same rate as the sample size! In 1948, Jerzy Neyman and Elizabeth Scott proved that when the number of nuisance parameters grows proportionally with $n$, the unconditional Maximum Likelihood Estimator of $\beta_1$ is asymptotically inconsistent. For 1:1 matched pairs, the unconditional estimator converges to exactly twice the true parameter value:

The Conditional Likelihood Formulation

To eliminate the infinite nuisance intercepts $\alpha_k$, David Cox (1958) formulated Conditional Logistic Regression. Rather than maximizing the unconditional likelihood, we condition on the sufficient statistic for $\alpha_k$: the total number of cases observed in stratum $k$ (which is fixed at 1 for 1:$M$ matched designs).

Consider matched stratum $k$ containing 1 case ($i = 0$) and $M$ controls ($i = 1, \dots, M$). The conditional probability that individual 0 is the case, given that exactly one person in the set developed the disease, is:

Dividing numerator and denominator by $\prod_{i=0}^M [1 - \pi(\mathbf{x}_{ki})]$:

Notice that the stratum intercept $\exp(\alpha_k)$ factors out identically from both the numerator and every term in the denominator:

The nuisance parameters $\alpha_k$ vanish completely! The conditional log-likelihood across all $K$ matched strata is:

Maximizing $\ell_{\text{cond}}(\boldsymbol{\beta})$ yields estimators that are fully consistent, asymptotically normal, and semiparametrically efficient. For 1:1 matched pairs with a single binary exposure, conditional logistic regression reduces exactly to McNemar’s Test and the paired ratio $b / c$.

In prospective cohort studies, randomized clinical trials, and cross-sectional studies where the outcome is common ($P > 10\%$), the odds ratio severely exaggerates the relative risk. Clinicians and policymakers naturally interpret adjusted effect estimates as relative risks. When analysts mistakenly interpret an adjusted odds ratio of $3.5$ from a logistic model as a "3.5-fold increase in risk," they mislead clinical audiences if the true risk ratio is only $1.8$.

The Log-Binomial Regression Model

To estimate the adjusted Relative Risk directly, we specify a Generalized Linear Model with a log link function applied to the binomial probability $\pi(\mathbf{x}_i) = P(Y_i = 1 \mid \mathbf{x}_i)$:

In this model, the parameter for a binary exposure $x_1$ has a direct relative risk interpretation:

However, the log-binomial model suffers from a severe numerical limitation: convergence failure. Because probabilities must satisfy $\pi(\mathbf{x}_i) \le 1.0$, the parameter space is constrained by the boundary condition $\mathbf{x}_i^T \boldsymbol{\beta} \le 0$ for all observed covariate vectors $\mathbf{x}_i$. In multivariable settings with continuous covariates, standard Newton-Raphson optimization algorithms frequently step outside this bounded parameter space, causing standard statistical software to crash with convergence errors.

Modified Poisson Regression with Robust Sandwich Variance

To circumvent log-binomial non-convergence, George Zou (2004) proved that an investigator can fit a Poisson regression model to binary data while applying a Huber-White robust sandwich covariance estimator.

In a Poisson model, the log link is standard: $\ln(\mu_i) = \mathbf{x}_i^T \boldsymbol{\beta}$. When applied to binary $Y_i \in \{0, 1\}$, the expectation is mathematically identical to the probability: $\mu_i = \mathbb{E}[Y_i \mid \mathbf{x}_i] = P(Y_i = 1 \mid \mathbf{x}_i) = \pi(\mathbf{x}_i)$. Thus, the exponentiated coefficients $\exp(\beta_j)$ yield consistent, unbiased estimates of the true adjusted Relative Risks!

However, fitting a Poisson model to binary data violates the mean-variance equality assumption ($\text{Var}(Y) = \mu$), because binary data have variance $\text{Var}(Y) = \pi(1 - \pi) < \pi$. The naive Poisson variance severely overestimates the true sampling variability. Zou demonstrated that applying the Huber-White Sandwich Estimator corrects this variance distortion:

Modified Poisson regression with robust standard errors is guaranteed to converge, handles continuous covariates effortlessly, and provides mathematically sound, asymptotically valid confidence intervals for adjusted Relative Risks.

When the epidemiological study tracks an open dynamic population or cohort where participants accumulate varying person-time at risk ($PT_i$), the outcome of interest is the count of incident events $Y_i \in \{0, 1, 2, \dots\}$.

Poisson Regression with Person-Time Offset

The Poisson probability mass function for event count $Y_i$ is:

The expected count $\mu_i$ is the product of the underlying incidence rate $\lambda(\mathbf{x}_i)$ and the accumulated person-time $PT_i$:

Taking the natural logarithm:

The term $\ln(PT_i)$ is included as a regression covariate whose coefficient is pre-fixed to exactly $1.0$. In generalized linear modeling, this fixed predictor is termed an offset.

The exponentiated regression coefficient $\exp(\beta_1)$ for a binary exposure represents the adjusted Incidence Rate Ratio ($IRR$):

Overdispersion and Negative Binomial Regression

A strict property of the Poisson distribution is equidispersion: the conditional variance must equal the conditional mean ($\text{Var}(Y_i \mid \mathbf{x}_i) = \mathbb{E}[Y_i \mid \mathbf{x}_i] = \mu_i$). In real-world disease surveillance, biological clustering, unobserved heterogeneity, and contagion frequently cause overdispersion, where the observed variance substantially exceeds the mean ($\text{Var}(Y_i) > \mu_i$).

To diagnose overdispersion, we evaluate the Pearson dispersion statistic:

If $\phi \gg 1.0$, standard Poisson standard errors are falsely deflated, leading to inflated Type I error rates. In such settings, epidemiologists employ the Negative Binomial Regression Model, which models unobserved subject heterogeneity through a gamma-distributed latent disturbance $\nu_i \sim \text{Gamma}(1/\alpha, \alpha)$. The resulting negative binomial variance is quadratic:

where $\alpha \ge 0$ is the dispersion parameter. As $\alpha \to 0$, the negative binomial model converges back to the standard Poisson model.

A central danger in multivariable epidemiological modeling is the uncritical use of automated algorithmic variable selection (such as forward stepwise selection or backward elimination based on $p$-values). Stepwise algorithms select variables based entirely on statistical significance in the sample, completely blind to the underlying biological causal structure.

Directed Acyclic Graph (DAG) Guided Covariate Selection

Modern epidemiologic practice mandates that covariate selection be guided by an explicit causal DAG formulated a priori from biological and clinical domain knowledge:

1. Identify all variables that satisfy Pearl's Backdoor Criterion (common causes of exposure and outcome). Include these confounders in the adjustment set.
2. Never adjust for an intermediate variable on the causal pathway between exposure and disease ($E \to M \to D$), unless the explicit scientific goal is to estimate the direct effect via formal mediation analysis. Adjusting for an intermediate blocks the indirect causal mechanism, biasing total effect estimates toward the null (overadjustment bias).
3. Never adjust for a collider ($E \to C \leftarrow D$), as conditioning on a collider opens a non-causal path, inducing artificial selection bias.

The Change-in-Estimate Criterion

In empirical data modeling, an informal criterion for assessing whether a candidate covariate $Z$ acts as an active confounder is the Change-in-Estimate Criterion:

If removing $Z$ alters the primary exposure effect estimate by more than $10\%$ (i.e., $\Delta \beta > 0.10$), $Z$ is retained as an essential empirical confounder regardless of its individual $p$-value.

In healthcare utilization studies and environmental exposure monitoring, count data often exhibit a massive excess of zero observations far beyond what even a negative binomial distribution can accommodate. For example, in an annual survey of emergency department visits or days absent from work due to acute occupational asthma, $85\%$ of healthy workers may report exactly zero events.

Epidemiologists model this phenomenon through two distinct statistical paradigms:

1. Zero-Inflated Poisson (ZIP) and Zero-Inflated Negative Binomial (ZINB) Models:

Postulate that the population is a mixture of two distinct latent subpopulations: an "always-zero" immune group who cannot experience the event (with probability $p_i$), and a susceptible count group whose event frequencies follow a Poisson or negative binomial distribution:

The zero probability $p_i$ is modeled via a logistic link ($\text{logit}(p_i) = \mathbf{z}_i^T \boldsymbol{\gamma}$), while the count process $\mu_i$ is modeled via a log link ($\ln(\mu_i) = \mathbf{x}_i^T \boldsymbol{\beta}$).

2. Two-Part Hurdle Models:

Formulate a sequential hurdle process: first, a binary logistic model determines whether an individual crosses the hurdle from zero to at least one event ($P(Y_i > 0)$); second, conditional on crossing the hurdle ($Y_i \ge 1$), a truncated-at-zero count distribution governs event frequency. Unlike zero-inflated models, hurdle models cleanly separate the initial propensity to use healthcare from subsequent frequency of use.