Likelihood and the generalized linear model

When you run glm() in R or fit a model in statsmodels, the software has to pick numbers for the coefficients. It does not guess. It chooses the values that make the data you actually collected as probable as the model allows. That rule is maximum likelihood, and it is the engine under linear, logistic, and Poisson regression alike. This lesson opens the engine so the output of every regression you have fitted in Parts A to D starts to read as one idea instead of four.

The likelihood of the data

Start with a single number to estimate. A tuberculosis index case is found in a Klang household, and you tuberculin-test the 12 close contacts. Three test positive, nine negative, so the sample proportion is . You do not care about these 12 people for their own sake. You want the true within-household transmission risk, .

Pick a candidate value for , say 0.2, and ask: how probable is the result I saw (3 of 12 positive) if the true risk were 0.2? The binomial distribution answers that. The probability of positives out of is proportional to . Read as a function of the data with fixed, that is a probability. Read the other way around, with the data fixed at and varying, the same expression is the likelihood: a score for how compatible each value of is with what you observed.

Maximum likelihood and the log-likelihood

Sweep from 0 to 1 and the likelihood rises, peaks, and falls. The value at the peak is the maximum likelihood estimate (MLE): the parameter value that makes your data most probable. For a single proportion the peak sits exactly at , which here is 0.25, the sample proportion you would have written down anyway. Likelihood does not overturn the simple estimate; it explains why the simple estimate is the right one, and it keeps working when the simple formula runs out.

In practice you maximise the log-likelihood, the natural log of the likelihood, rather than the likelihood itself. Two reasons. First, a likelihood multiplies one term per observation, so with 200 patients it underflows to a number the computer rounds to zero; logs turn that product into a sum that stays well-behaved. Second, sums are easier to differentiate, which is how the maximum is found. The log is a steady increasing transform, so it does not move the location of the peak.

Worked example 1: the MLE of a proportion

Confirm that the peak of the log-likelihood for the TB contacts lands at 3/12, two ways: by scanning a grid of candidate values, and by the closed form . The two must agree.

1. Build a grid of values from 0.01 to 0.99.
2. Evaluate the log-likelihood at each one.
3. Read off the with the highest log-likelihood, and check it equals .

d <- 3
n <- 12
grid <- seq(0.01, 0.99, by = 0.01)

loglik <- # d*log(pi) + (n-d)*log(1-pi) over the grid

loglik     <- d * log(grid) + (n - d) * log(1 - grid)
mle_grid   <- grid[which.max(loglik)]
mle_closed <- d / n
round(c(grid = mle_grid, closed = mle_closed), 4)
#   grid closed
#   0.25   0.25

import numpy as np
d, n = 3, 12
grid = np.arange(0.01, 1.00, 0.01)

loglik = # d*log(pi) + (n-d)*log(1-pi) over the grid

loglik     = d * np.log(grid) + (n - d) * np.log(1 - grid)
mle_grid   = grid[np.argmax(loglik)]
mle_closed = d / n
print(round(mle_grid, 4), round(mle_closed, 4))
# 0.25 0.25

The likelihood ratio and the LRT

One parameter value has the highest likelihood, but a band of nearby values fit almost as well. To compare two values you take their likelihood ratio: the likelihood at one divided by the likelihood at the MLE. The ratio is 1 at the MLE and shrinks as you move away. This becomes a hypothesis test when one of the values is a null. Compare a model that includes a term against a nested model that drops it. (Nested means the smaller model is the larger one with some coefficients forced to zero.) The likelihood ratio test (LRT) uses minus twice the difference in log-likelihoods:

The statistic follows a chi-squared distribution with degrees of freedom equal to the number of parameters dropped. A large value means the extra terms bought a real gain in fit, so the data argue against the reduced model.

The generalized linear model

Linear, logistic, and Poisson regression look like three methods. Likelihood shows them as one template, the generalized linear model (GLM). Every GLM has the same two pieces. A linear predictor, the familiar weighted sum of the exposures, and a link function that connects that sum to the mean of the outcome.

What changes between models is the outcome's distribution and the link. The link keeps the prediction in the range the outcome allows: a probability stays between 0 and 1, a rate stays positive.

Ordinary linear regression is the GLM with a normal outcome and the identity link, so the mean difference you read there is a coefficient on the same footing as a log odds ratio in logistic regression. One fitting routine, maximum likelihood, serves all three.

Wald test versus likelihood ratio test

For a single coefficient you can test the null two ways. The Wald test divides the estimate by its standard error, , and is the z and P-value printed next to each coefficient in a regression table. The likelihood ratio test refits the model without that term and compares log-likelihoods. The Wald test reads one fitted model. The LRT needs two fits, but it does not depend on the scale of the parameter, which makes it the safer default. For a categorical exposure carried by several indicators, that difference matters. When the sample is large and the model well behaved, the two agree closely.

Worked example 2: an LRT between nested GLMs

Take a dengue admissions dataset. The outcome is whether a patient developed severe disease. You already adjust for age, and you want to know whether the admission platelet count adds anything. Fit the reduced model (age only) and the full model (age plus platelet), then compare them with an LRT. Read the platelet Wald test from the full model alongside it.

set.seed(1)
n  <- 200
age      <- runif(n, 5, 70)
platelet <- runif(n, 20, 200)
lp       <- -2 + 0.03 * age - 0.02 * platelet
severe   <- rbinom(n, 1, 1 / (1 + exp(-lp)))

m1 <- # glm(severe ~ age, family = binomial)
m2 <- # add platelet

m1 <- glm(severe ~ age, family = binomial)
m2 <- glm(severe ~ age + platelet, family = binomial)
# Wald test for platelet (from the model summary)
round(summary(m2)$coefficients["platelet", ], 4)
# Likelihood ratio test comparing the nested models
anova(m1, m2, test = "LRT")
# 2*(logLik(m2) - logLik(m1)) is the LR statistic, 1 d.f.

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import chi2
rng = np.random.default_rng(1)
n = 200
age      = rng.uniform(5, 70, n)
platelet = rng.uniform(20, 200, n)
lp       = -2 + 0.03 * age - 0.02 * platelet
severe   = rng.binomial(1, 1 / (1 + np.exp(-lp)))
df = pd.DataFrame({"severe": severe, "age": age, "platelet": platelet})

m1 = # glm severe ~ age
m2 = # glm severe ~ age + platelet

B  = sm.families.Binomial()
m1 = smf.glm("severe ~ age", data=df, family=B).fit()
m2 = smf.glm("severe ~ age + platelet", data=df, family=B).fit()
# Wald test for platelet
print("Wald z, p:", round(m2.tvalues["platelet"], 4),
      round(m2.pvalues["platelet"], 4))
# Likelihood ratio test comparing the nested models
lrs  = 2 * (m2.llf - m1.llf)
pval = chi2.sf(lrs, df=1)
print("LR stat, p:", round(lrs, 4), round(pval, 4))

Both languages report a large LR statistic and a small P-value, so platelet count earns its place in the model. The Wald P-value next to the platelet coefficient lands close to the LRT P-value, as expected when the sample is reasonable and the term is a single parameter.

Common mistakes

• Reading the likelihood as a probability of the parameter. is the probability of the data given , not the probability that takes a value. The likelihood does not integrate to 1 over , so it is not a distribution for the parameter.
• Running an LRT on non-nested models. The chi-squared result needs one model to be a restricted form of the other. Comparing age-only against platelet-only is not a valid LRT, because neither nests inside the other. Use an information criterion such as AIC for non-nested comparisons.
• Comparing log-likelihoods from different samples. Missing data, a changed transform of the outcome, or a different weighting all break the comparison. The two models must use identical rows and the same outcome.
• Trusting a single Wald P-value for a multi-category exposure. A factor with four levels carries three coefficients. The three separate Wald tests do not answer whether the factor matters overall. Drop the whole factor and run one LRT with 3 degrees of freedom.
• Forgetting the link when reading a coefficient. A logistic coefficient is a log odds ratio and a Poisson coefficient is a log rate ratio. Exponentiate before you interpret, or you will report the wrong scale.

Tips

• In R, logLik(model) returns the log-likelihood and its degrees of freedom, and anova(m1, m2, test = "LRT") does the comparison for you. In statsmodels the fitted result exposes .llf.
• Prefer the LRT to the Wald test for any term with more than one parameter, and whenever the estimate sits near a boundary. The LRT does not change when you rescale the predictor.
• Use confint(model) in R for profile-likelihood confidence intervals, which match the LRT. The default confint.default gives Wald intervals instead, which can differ for small samples.
• When a GLM warns about non-convergence, suspect perfect separation or a category with no events. More iterations rarely fixes it; collapse sparse categories or rethink the predictor.
• Pick the link from the outcome, not from habit. Counts and person-time go with the log link and a rate ratio; a yes-or-no outcome goes with the logit and an odds ratio.