Statistical Background

This page provides background on the statistical models and methods used in LEAP.

Generalized Linear Models

A generalized linear model (GLM) is a generalization of ordinary linear regression that allows for response variables with non-normal error distributions.

Background: Ordinary Linear Regression

In ordinary linear regression, we assume that the response variable \(y\) can be described by the equation:

\[y_i = \beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i} + \ldots + \beta_n x_{n,i}\]

where \(i\) is the index of the sample, \(n\) is the number of features, \(\beta_0, \beta_1, \ldots, \beta_n\) are the model coefficients, and \(x_{1,i}, \ldots, x_{n,i}\) are the features of the \(i\)-th sample.

GLM: Linear Predictor

For a GLM, we assume that the relationship between different terms is linear, but the terms themselves are not necessarily linear. The linear predictor \(\eta_i\) can include polynomial terms, interaction terms, or transformed variables. The key constraint is that the \(\beta\) parameters must enter linearly — exponents on \(\beta\) themselves are not permitted.

Valid examples of linear predictors:

\[\begin{split}\eta_i &= \beta_0 + \beta_{11} x_{1,i} + \beta_{12} x_{1,i}^2 + \beta_{21} x_{2,i} + \beta_{22} x_{2,i}^2 \quad \text{(polynomial)} \\[4pt] \eta_i &= \beta_0 + \beta_{12} x_{1,i} x_{2,i} + \beta_3 x_{3,i} \quad \text{(interaction)} \\[4pt] \eta_i &= \beta_0 + \beta_1 e^{x_{1,i}} + \beta_2 x_{2,i} \quad \text{(transformed covariate)}\end{split}\]

Invalid — the exponent on a \(\beta\) makes this non-linear:

\[\eta_i = \beta_0 + \beta_1 x_{1,i}^{\beta_5} + \beta_2 x_{2,i}\]

GLM: Random Component

The random component assumes that the observed data \(Y\) comes from an exponential family of distributions — such as Gaussian, Poisson, binomial, or gamma:

\[p(y_i;\, \theta_i) = a(\theta_i)\, b(y_i)\, e^{y_i Q(\theta_i)}\]

where \(\theta_i\) is a function of the explanatory variables and \(Q(\theta)\) is the natural parameter.

The link function \(g\) connects the mean \(\mu_i = E(Y \mid X = x_i)\) to the linear predictor:

\[g(\mu_i) = \eta_i\]

The goal of a GLM is to estimate \(\mu\), not \(\eta\) directly. The link function that maps \(\mu\) to \(Q(\theta)\) is called the canonical link function.

For ordinary linear regression, the link function is the identity:

\[g(\mu_i) = \mu_i = \eta_i\]

The random component is the Gaussian distribution:

\[p(y_i;\, \mu_i) = \frac{1}{\sigma\sqrt{2\pi}}\, e^{-\frac{(y_i - \mu_i)^2}{2\sigma^2}}\]

When the response variable represents counts of events (non-negative integers), the Poisson distribution is a natural choice. Incidence and prevalence rates — counts of individuals diagnosed with or living with a condition in a given time interval — follow this pattern.

The log link is the canonical link for the Poisson distribution. The mean \(\mu_i\) is always positive, but the linear predictor \(\eta_i\) can be any real number. The log link satisfies this constraint:

\[g(\mu_i) = \log(\mu_i) = \eta_i\]

The random component is the Poisson distribution:

\[P(Y = y) = \frac{(\mu_i)^y\, e^{-\mu_i}}{y!}\]

This is the distribution family used in LEAP’s occurrence Model 1 to predict population-level asthma incidence and prevalence rates.

The Poisson distribution assumes that the mean and variance are equal:

\[\mu = \sigma^2\]

When the variance exceeds the mean — a common problem in count data known as overdispersion — the Negative Binomial distribution is a better choice. It introduces an extra parameter \(\theta\) that controls the degree of overdispersion:

\[\sigma^2 = \mu + \frac{\mu^2}{\theta}\]

As \(\theta \to \infty\), the Negative Binomial distribution converges to the Poisson distribution.

The standard form of the Negative Binomial probability mass function is:

\[P(Y = k;\, r, p) := \binom{k+r-1}{k}(1-p)^k p^r\]

where \(k\) is the number of failures before \(r\) successes, and \(p\) is the probability of success. We reparametrize using \(\mu\) and \(\theta\) via:

\[\begin{split}p &= \frac{\mu}{\sigma^2} \\ r &= \frac{\mu^2}{\sigma^2 - \mu} \\ \sigma^2 &= \mu + \frac{\mu^2}{\theta}\end{split}\]

Substituting \(\sigma^2\) and simplifying:

\[\begin{split}p &= \frac{\theta}{\theta + \mu} \\ r &= \theta\end{split}\]
Math: \(p\) and \(r\)
\[\begin{split}p &= \dfrac{\mu}{\sigma^2} = \dfrac{\mu}{\mu + \dfrac{\mu^2}{\theta}} = \dfrac{\theta}{\theta + \mu} \\ r &= \dfrac{\mu^2}{\sigma^2 - \mu} = \dfrac{\mu^2}{\mu + \dfrac{\mu^2}{\theta} - \mu} = \dfrac{\mu^2}{\dfrac{\mu^2}{\theta}} = \theta\end{split}\]

Letting \(y = k\), the probability mass function in terms of \(\mu\) and \(\theta\) is:

\[P(Y = y;\, \mu, \theta) = \binom{y + \theta - 1}{y} \frac{\mu^y\, \theta^{\theta}}{(\theta + \mu)^{y+\theta}}\]
Math: \(P(Y = y;\, \mu, \theta)\)
\[\begin{split}P(Y = y; \mu, \theta) &= \binom{y + \theta - 1}{y} \left(1-\dfrac{\theta}{\theta + \mu}\right)^y \left(\dfrac{\theta}{\theta + \mu}\right)^{\theta} \\ &= \binom{y + \theta - 1}{y} \left(\dfrac{\theta + \mu - \theta}{\theta + \mu}\right)^y \left(\dfrac{\theta}{\theta + \mu}\right)^{\theta} \\ &= \binom{y + \theta - 1}{y} \dfrac{\mu^y \theta^{\theta}}{(\theta + \mu)^{y+\theta}}\end{split}\]

The log link is the natural choice, since the mean is always positive but the linear predictor \(\eta^{(i)}\) can be any real number:

\[g(\mu^{(i)}) = \log(\mu^{(i)}) = \eta^{(i)}\]

This is the distribution family used in LEAP’s Antibiotic Exposure Model to predict per-capita antibiotic exposure rates.

Ordinal regression is used when the response variable is ordered but the intervals between levels are arbitrary. Rather than modelling a single mean, the model predicts the cumulative probability of being at or below each level \(k\) using the logistic (sigmoid) function as the link:

\[P(y^{(i)} \leq k) = \sigma(\theta_k + \eta^{(i)})\]

where \(\theta_k\) is the threshold parameter for level \(k\), \(\eta^{(i)}\) is the linear predictor, and \(\sigma(x) = \dfrac{1}{1 + e^{-x}}\) is the logistic function. The probability of being in exactly level \(k\) follows from the cumulative probabilities:

\[P(y^{(i)} = k) = P(y^{(i)} \leq k) - P(y^{(i)} \leq k - 1)\]

A patient-specific random effect \(\beta_0^{(i)} \sim \mathcal{N}(0, \sigma^2)\) can be added to the linear predictor to account for within-subject correlation across repeated measurements, giving:

\[P(y^{(i)} \leq k) = \sigma\!\left(\theta_k + \eta^{(i)} + \beta_0^{(i)}\right)\]

This is the model used in LEAP’s Asthma Control Model to predict asthma control level.

Contingency Tables

A contingency table (also called a cross-tabulation or crosstab) displays the joint frequency distribution of two categorical variables. They are commonly used in statistics to examine the relationship between two variables, and to calculate odds ratios and other measures of association.

Two-by-Two Contingency Table

The simplest form is a 2×2 table, which has two rows and two columns. Consider two binary variables — an exposure and an outcome — for a population of \(n\) individuals:

outcome + outcome - total
exposure + a b n₁
exposure - c d n₀
total n₂ n

where the marginal totals are:

\[\begin{split}n_1 &= a + b \quad \text{(total with exposure)} \\ n_0 &= c + d \quad \text{(total without exposure)} \\ n_2 &= a + c \quad \text{(total with outcome)} \\ n &= a + b + c + d \quad \text{(total population)}\end{split}\]

The odds ratio measures the strength of association between the exposure and the outcome:

\[\omega = \frac{a \cdot d}{b \cdot c}\]

An odds ratio of 1 indicates no association; values greater than 1 indicate that exposure is associated with increased odds of the outcome.

Example: Smoking and Lung Cancer

Suppose we observe n = 300 patients and record whether they smoke and whether they have lung cancer:

lung cancer no lung cancer total
smoker a = 100 b = 50 n₁ = 150
non-smoker c = 25 d = 125 n₀ = 150
total n₂ = 125 n = 300

The odds ratio for this table is:

\[\omega = \frac{100 \times 125}{50 \times 25} = 10\]

This means that smokers have 10 times the odds of developing lung cancer compared to non-smokers in this population.

Usage in LEAP

Contingency tables are used in LEAP to optimize the age-dependent beta parameters for the incidence equation in Occurrence Model 2: Risk Factors. See the Calibrating Age-Dependent Odds Ratios section of the Asthma Occurrence Model for details.

Python Examples

For a hands-on walkthrough of GLMs using statsmodels and plotly — including code you can run locally — see the notebook below.