# group coded -1 (control) and +1 (treatment); y = 1 means recovered
wf_x <- rep(c(-1, 1), each = 3)
wf_y <- c(0, 0, 1, 0, 1, 1)
# maximum likelihood fit of the logistic model
wf_fit <- glm(wf_y ~ wf_x, family = binomial)
# the fitted intercept and slope
round(coef(wf_fit), 3)
#> (Intercept) wf_x
#> 0.000 0.693Inference for Logistic Regression
1 Introduction
Logistic regression models a binary outcome, such as whether a patient recovers, from one or more predictors. Its fitted coefficients are maximum likelihood estimates, so the likelihood theory of Introduction to maximum likelihood gives them standard errors, confidence intervals and tests. This page specializes that theory to the logistic model and introduces the deviance, the measure of fit that the likelihood ratio test of a logistic model is built from [1].
DefDefinition 1 (Logistic regression)
Definition 1 (Logistic regression) A logistic regression model for binary outcomes \(y_1, \dots, y_n \in {\mathopen{}\left\{0, 1\right\}\mathclose{}}\) with predictor vectors \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_1, \dots, \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_n\) says that the outcomes are independent and that
\[ {\Pr\mathopen{}\left(Y_i= 1 \mid \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i\right)\mathclose{}} = \pi_i = {\operatorname{expit}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i\right)\mathclose{}}, \qquad {\operatorname{expit}\mathopen{}\left(u\right)\mathclose{}} = \frac{1}{1 + e^{-u}}, \tag{1}\]
where \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i= {\mathopen{}\left(1, x_{i1}, \dots, x_{ip}\right)\mathclose{}}^{\top}\) and \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}= {\mathopen{}\left(\beta_{0}, \beta_{1}, \dots, \beta_{p}\right)\mathclose{}}^{\top}\). The leading \(1\) in \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i\) makes \(\beta_{0}\) the intercept, and \(\beta_{j}\) for \(j \ge 1\) is the slope for the \(j\)th predictor.
Source: Wikipedia contributors [2].
ExmExample 1 (A six-patient trial)
Example 1 (A six-patient trial) Six patients are split into a control group, coded \(x= -1\), and a treatment group, coded \(x= +1\). The outcome is \(y= 1\) if the patient recovered. Fitting the model of Equation 1 by maximum likelihood gives the coefficients below.
The fitted slope is \(\hat{\beta}_{1} = 0.693\). The control group recovered at rate \(1/3\) and the treatment group at \(2/3\), so the model fits both group proportions exactly.
2 Odds and the coefficients
LemLemma 1 (The log-odds in a logistic regression model)
Lemma 1 (The log-odds in a logistic regression model) In the logistic regression model of Equation 1, the log-odds of \(Y= 1\) at predictor vector \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) are
\[ {\eta\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}} = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}. \]
ProofProof
Proof. Write \(u \stackrel{\text{def}}{=}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) and \(\pi\stackrel{\text{def}}{=}{\operatorname{expit}\mathopen{}\left(u\right)\mathclose{}}\), so that
\[ \begin{aligned} 1 - \pi &= 1 - \frac{1}{1 + \text{e}^{-u}} && \text{(the logistic model)} \\ &= \frac{1 + \text{e}^{-u}}{1 + \text{e}^{-u}} - \frac{1}{1 + \text{e}^{-u}} && \text{(write } 1 \text{ over the denominator } 1 + \text{e}^{-u} \text{)} \\ &= \frac{\mathopen{}\left(1 + \text{e}^{-u}\right)\mathclose{} - 1}{1 + \text{e}^{-u}} && \text{(subtract fractions with the same denominator)} \\ &= \frac{\text{e}^{-u}}{1 + \text{e}^{-u}}. && \text{(} \mathopen{}\left(1 + a\right)\mathclose{} - 1 = a \text{)} \end{aligned} \]
The first step substitutes the model of Equation 1. So the odds are
\[ \begin{aligned} \frac{\pi}{1 - \pi} &= \frac{1 / \mathopen{}\left(1 + \text{e}^{-u}\right)\mathclose{}}{\text{e}^{-u} / \mathopen{}\left(1 + \text{e}^{-u}\right)\mathclose{}} && \text{(the logistic model, and the display for } 1 - \pi\text{)} \\ &= \frac{1}{\text{e}^{-u}} && \text{(divide fractions with the same denominator)} \\ &= \text{e}^{u}, && \text{(} 1 / \text{e}^{-u} = \text{e}^{u} \text{)} \end{aligned} \]
and the log-odds are \({\operatorname{log}\mathopen{}\left\{\text{e}^{u}\right\}\mathclose{}} = u = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\).
ExmExample 2 (The log-odds of recovery in the six-patient trial)
Example 2 (The log-odds of recovery in the six-patient trial)
In the six-patient trial, a treated patient has \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}= {\mathopen{}\left(1, 1\right)\mathclose{}}^{\top}\). Under the maximum likelihood fit, the code computes the probability of recovery, its log-odds, and the dot product \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\).
# fitted probability of recovery for a treated patient (x = +1)
wf_treated_data <- tibble::tibble(wf_x = 1)
wf_prob_treated <- predict(wf_fit, wf_treated_data, type = "response")
wf_prob_treated <- unname(wf_prob_treated)
wf_log_odds_treated <- log(wf_prob_treated / (1 - wf_prob_treated))
wf_coef <- coef(wf_fit)
wf_dot_treated <- sum(wf_coef)
c(prob = wf_prob_treated, log_odds = wf_log_odds_treated, dot = wf_dot_treated)
#> prob log_odds dot
#> 0.666667 0.693147 0.693147That probability is \(0.6667\), and both the log-odds and the dot product are \(0.6931\), as Lemma 1 states.
DefDefinition 2 (Odds multiplier)
Definition 2 (Odds multiplier) Let:
- \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) be the coefficient vector of a logistic regression model.
- \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) be a change in the predictor vector, whose first entry is \(0\), because every predictor vector has first entry \(1\).
The odds multiplier of \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) is the exponential of the dot product of \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) with \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\):
\[ {\theta_{\omega}\mathopen{}\left(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}} \stackrel{\text{def}}{=}\text{e}^{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}}. \]
ExmExample 3 (The odds multiplier of moving from control to treatment)
Example 3 (The odds multiplier of moving from control to treatment)
In the six-patient trial, moving a patient from the control group (\(x= -1\)) to the treatment group (\(x= +1\)) changes the predictor vector by \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}= {\mathopen{}\left(0, 2\right)\mathclose{}}^{\top}\). The code computes the odds multiplier of that change under the maximum likelihood fit.
That odds multiplier is \(\text{e}^{2 \hat{\beta}_{x}} = 4\).
ExmExample 4 (The risk ratio is not the odds multiplier)
Example 4 (The risk ratio is not the odds multiplier)
In the six-patient trial, the probability of recovery under the maximum likelihood fit is \(2/3\) in the treatment group and \(1/3\) in the control group. The code computes the risk ratio of treatment against control, and shows beside it the odds multiplier of the change \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}= {\mathopen{}\left(0, 2\right)\mathclose{}}^{\top}\) from Example 3.
# fitted probabilities of recovery: control (x = -1), treatment (x = +1)
wf_group_codes <- c(-1, 1)
wf_groups <- tibble::tibble(wf_x = wf_group_codes)
wf_probs <- predict(wf_fit, wf_groups, type = "response")
wf_probs <- unname(wf_probs)
wf_prob_ratio <- wf_probs[2] / wf_probs[1]
c(risk_ratio = wf_prob_ratio, odds_multiplier = wf_multiplier)
#> risk_ratio odds_multiplier
#> 2 4The risk ratio is \(2\), but the odds multiplier is \(4\). The odds multiplier compares odds, \(2\) against \(1/2\), while the risk ratio compares probabilities.
LemTheorem 1 (A change in the predictors multiplies the odds by its odds multiplier)
Theorem 1 (A change in the predictors multiplies the odds by its odds multiplier) Let:
- \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) be a predictor vector in the logistic regression model of Equation 1.
- \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) be a change in the predictor vector, with first entry \(0\).
- \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}' \stackrel{\text{def}}{=}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}+ \Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\).
Then the odds of \(Y= 1\) at \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'\) are the odds at \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) times the odds multiplier of \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\):
\[ {\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'\right\}\mathclose{}} = {\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right\}\mathclose{}} \, {\theta_{\omega}\mathopen{}\left(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}}. \]
ProofProof
Proof. The log-odds are the natural logarithm of the odds, so the odds at a predictor vector are \(\text{e}^{}\) raised to the log-odds there, which Lemma 1 gives.
\[ \begin{aligned} {\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'\right\}\mathclose{}} &= \text{e}^{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'} && \text{(log-odds of the model at } \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}' \text{)} \\ &= \text{e}^{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}+ \Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}} && \text{(definition of } \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}' \text{)} \\ &= \text{e}^{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}} + \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}} && \text{(the dot product is linear in each slot)} \\ &= \text{e}^{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}} \, \text{e}^{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}} && \text{(} \text{e}^{a + b} = \text{e}^{a} \text{e}^{b} \text{)} \\ &= {\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right\}\mathclose{}} \, \text{e}^{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}} && \text{(log-odds of the model at } \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\text{)} \\ &= {\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right\}\mathclose{}} \, {\theta_{\omega}\mathopen{}\left(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}}. && \text{(definition of the odds multiplier)} \end{aligned} \]
The first and fifth steps use Lemma 1, the third step uses the linearity of the dot product, and the last step uses Definition 2.
ExmExample 5 (The treated odds from the control odds)
Example 5 (The treated odds from the control odds)
In the six-patient trial, take \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}= {\mathopen{}\left(1, -1\right)\mathclose{}}^{\top}\), a control patient, and \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}= {\mathopen{}\left(0, 2\right)\mathclose{}}^{\top}\), so that \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}' = {\mathopen{}\left(1, 1\right)\mathclose{}}^{\top}\) is a treated patient. Under the maximum likelihood fit, the code computes the odds of recovery at both predictor vectors, and the control odds times the odds multiplier of \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\).
# fitted odds of recovery: control (x = -1), treatment (x = +1)
wf_odds <- wf_probs / (1 - wf_probs)
c(
odds_control = wf_odds[1],
odds_treated = wf_odds[2],
control_times_multiplier = wf_odds[1] * wf_multiplier
)
#> odds_control odds_treated control_times_multiplier
#> 0.5 2.0 2.0The control odds are \(0.5\), and \(0.5 \times 4 = 2\) equals the treated odds, as Theorem 1 states.
CorCorollary 1 (The odds ratio between two predictor vectors)
Corollary 1 (The odds ratio between two predictor vectors) Let \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) and \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'\) be predictor vectors in the logistic regression model of Equation 1. The odds ratio of \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'\) compared with \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) is the odds multiplier of their difference:
\[ {\text{ϙ}_{\omega}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}', \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}} = {\theta_{\omega}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}' - \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}}. \]
ProofProof
Proof. Let \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\stackrel{\text{def}}{=}\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}' - \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\), so that \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}' = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}+ \Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\), and the first entry of \(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\) is \(1 - 1 = 0\).
\[ \begin{aligned} {\text{ϙ}_{\omega}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}', \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}} &= \frac{{\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}'\right\}\mathclose{}}}{{\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right\}\mathclose{}}} && \text{(definition of the odds ratio)} \\ &= \frac{{\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right\}\mathclose{}} \, {\theta_{\omega}\mathopen{}\left(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}}}{{\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right\}\mathclose{}}} && \text{(odds multiplier theorem)} \\ &= {\theta_{\omega}\mathopen{}\left(\Delta\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right)\mathclose{}}. && \text{(divide numerator and denominator by } {\operatorname{odds}\mathopen{}\left\{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}\right\}\mathclose{}} > 0 \text{)} \end{aligned} \]
The first step is the definition of the odds ratio, and the second is Theorem 1.
ExmExample 6 (The odds ratio of treatment against control from the group counts)
Example 6 (The odds ratio of treatment against control from the group counts)
In the six-patient trial, \(2\) of the \(3\) treated patients and \(1\) of the \(3\) control patients recovered. The code computes the odds ratio of treatment against control from those counts, and the odds multiplier under the maximum likelihood fit of \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}' - \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}= {\mathopen{}\left(1, 1\right)\mathclose{}}^{\top} - {\mathopen{}\left(1, -1\right)\mathclose{}}^{\top}\).
# odds ratio from the counts: treated 2 of 3, control 1 of 3
count_or <- ((2 / 3) / (1 / 3)) / ((1 / 3) / (2 / 3))
wf_difference <- c(1, 1) - c(1, -1)
wf_difference_terms <- coef(wf_fit) * wf_difference
wf_difference_dot <- sum(wf_difference_terms)
wf_difference_multiplier <- exp(wf_difference_dot)
c(count_or = count_or, odds_multiplier = wf_difference_multiplier)
#> count_or odds_multiplier
#> 4 4Both are \(4\), as Corollary 1 states, because the model fits the two group proportions exactly.
3 Fitting by iteratively reweighted least squares
Note 1: Source and license: iteratively reweighted least squares
This section is adapted from the generalized linear model part of the Newton’s method lecture of Zhou [3], used under the MIT License. The license text is:
Copyright (c) 2025 ucla-biostat-257
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.The source writes the logistic score and information, notes that Newton’s method and Fisher scoring coincide, and rewrites the step as a weighted least squares fit. Here the step is stated as a theorem and proved, with the step length fixed at \(1\). The example, using the six-patient trial, is new. The source’s treatment of general generalized linear models, Gauss-Newton and Levenberg-Marquardt is not part of this section. The general optimization methods behind IRLS are on the Optimization page of the math notes.
DefDefinition 3 (Iteratively reweighted least squares (IRLS))
Definition 3 (Iteratively reweighted least squares (IRLS)) Suppose the outcomes follow the logistic regression model. Let \(X\) be the matrix whose rows are \({\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_{i}}^{\top}\). From a current estimate \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^*\) with fitted probabilities \(\pi_{i}^*\) (collected in \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\pi}}\vphantom{\pi}}^*\)), form
- the weights \(w_{i} = \pi_{i}^* \mathopen{}\left(1 - \pi_{i}^*\right)\mathclose{}\), collected in the diagonal matrix \(W\);
- the working responses \(z_{i} = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^* \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_{i} + \mathopen{}\left(y_{i} - \pi_{i}^*\right)\mathclose{} / w_{i}\), collected in \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}}\).
Iteratively reweighted least squares replaces \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^*\) by the weighted least squares fit of \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}}\) on \(X\), the minimizer of \(\sum_{i=1}^{n} w_{i} \mathopen{}\left(z_{i} - \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_{i}\right)\mathclose{}^2\), which is
\[ {\mathopen{}\left(\mathopen{}\left({X}^{\top} W X\right)\mathclose{}\right)\mathclose{}}^{-1} {X}^{\top} W \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}}, \]
and repeats until the estimate stops changing.
ExrExercise 1 (One IRLS step for the six-patient trial)
Exercise 1 (One IRLS step for the six-patient trial) In the six-patient trial, start from \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^* = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{0}}\vphantom{0}}\).
- Find the fitted probabilities, the weights and the working responses.
- Compute \({X}^{\top} W X\) and \({X}^{\top} W \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}}\).
- Compute the updated coefficients \({\mathopen{}\left(\mathopen{}\left({X}^{\top} W X\right)\mathclose{}\right)\mathclose{}}^{-1} {X}^{\top} W \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}}\).
LemTheorem 2 (An IRLS step is a Newton-Raphson step and a Fisher scoring step)
Theorem 2 (An IRLS step is a Newton-Raphson step and a Fisher scoring step) Suppose the outcomes follow the logistic regression model. One step of iteratively reweighted least squares from \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^*\) equals
- one step of the Newton-Raphson algorithm on the log-likelihood from \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^*\), and
- one step of Fisher scoring from \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^*\):
\[ \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^* + {\mathopen{}\left(\mathopen{}\left({X}^{\top} W X\right)\mathclose{}\right)\mathclose{}}^{-1} {X}^{\top} \mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}- \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\pi}}\vphantom{\pi}}^*\right)\mathclose{}. \tag{2}\]
ProofProof
Proof. The log-likelihood has score \({X}^{\top} \mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}- \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\pi}}\vphantom{\pi}}\right)\mathclose{}\) and Hessian \(-{X}^{\top} W X\) at every \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\) (with \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\pi}}\vphantom{\pi}}\) and \(W\) computed from \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}\)), and the Hessian does not involve the outcomes (proof of Theorem 6). So the observed information \({X}^{\top} W X\) equals the expected information, and the Newton-Raphson and Fisher scoring steps both equal Equation 2.
For IRLS, the working responses satisfy \(W \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}} = W X \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^* + \mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}- \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\pi}}\vphantom{\pi}}^*\right)\mathclose{}\), so
\[ \begin{aligned} {\mathopen{}\left(\mathopen{}\left({X}^{\top} W X\right)\mathclose{}\right)\mathclose{}}^{-1} {X}^{\top} W \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}} &= {\mathopen{}\left(\mathopen{}\left({X}^{\top} W X\right)\mathclose{}\right)\mathclose{}}^{-1} {X}^{\top} W X \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^* + {\mathopen{}\left(\mathopen{}\left({X}^{\top} W X\right)\mathclose{}\right)\mathclose{}}^{-1} {X}^{\top} \mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}- \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\pi}}\vphantom{\pi}}^*\right)\mathclose{} && \text{(substitute } W \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}} \text{)} \\ &= \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^* + {\mathopen{}\left(\mathopen{}\left({X}^{\top} W X\right)\mathclose{}\right)\mathclose{}}^{-1} {X}^{\top} \mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}- \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\pi}}\vphantom{\pi}}^*\right)\mathclose{}. && \text{(cancel } {X}^{\top} W X \text{)} \end{aligned} \]
ExmExample 7 (IRLS on the six-patient trial)
Example 7 (IRLS on the six-patient trial) The code runs IRLS from \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^* = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{0}}\vphantom{0}}\) on the data of the six-patient trial and compares the result with the coefficients from glm().
# design matrix: intercept, group (-1 control, +1 treatment); y = 1: recovered
irls_x <- cbind(1, rep(c(-1, 1), each = 3))
irls_y <- c(0, 0, 1, 0, 1, 1)
# start from zero coefficients and take six IRLS steps
irls_coef <- c(0, 0)
irls_path <- data.frame(step = 0, intercept = 0, slope = 0)
for (step in 1:6) {
prob <- plogis(drop(irls_x %*% irls_coef))
weight <- prob * (1 - prob)
work_resp <- drop(irls_x %*% irls_coef) + (irls_y - prob) / weight
irls_coef <- drop(solve(crossprod(irls_x, weight * irls_x),
crossprod(irls_x, weight * work_resp)))
irls_path[step + 1, ] <- c(step, irls_coef)
}
irls_pathThe first step moves the slope to \(0.6667\), which is \(2/3\) (the hand computation). After six steps the slope is \(0.693147\), and glm() reports \(0.693147\). The slope settles at \(\ln 2 = 0.693147\), half the natural logarithm of the odds ratio between the treatment group (recovery odds \(2\)) and the other group (recovery odds \(1/2\)), since the groups are coded \(\pm 1\).
4 Deviance
DefDefinition 4 (Deviance)
Definition 4 (Deviance) The deviance of a fitted logistic regression with labels \(y_i\in {\mathopen{}\left\{0, 1\right\}\mathclose{}}\) is minus twice its maximized log-likelihood:
\[ \operatorname{dev}\stackrel{\text{def}}{=}-2 \, {\ell\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\right)\mathclose{}} = -2 \sum_{i=1}^n\mathopen{}\left\{y_i\ln \hat\pi_i + \mathopen{}\left(1 - y_i\right)\mathclose{} {\ln\mathopen{}\left(1 - \hat\pi_i\right)\mathclose{}}\right\}\mathclose{}, \tag{3}\]
where \(\hat\pi_i = {\operatorname{expit}\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}} \cdot \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i\right)\mathclose{}}\) is the fitted probability.
Source: Wikipedia contributors [4].
NoteSource
Nelder and Wedderburn [5] introduced the deviance for the whole family of generalized linear models; logistic regression is the case of Bernoulli outcomes with the logit link. Hastie et al. [6, Sec. 4.4.2] defines the residual deviance of a fitted logistic regression as minus twice its log-likelihood.
ExmExample 8 (The deviance of the six-patient trial)
Example 8 (The deviance of the six-patient trial) For the fit of Example 1, R’s glm reports the deviance. We check it against Equation 3 computed by hand.
# fitted probabilities of recovery
wf_p <- fitted(wf_fit)
# minus twice the log-likelihood, from the formula
wf_dev_hand <- -2 * sum(wf_y * log(wf_p) + (1 - wf_y) * log(1 - wf_p))
# the deviance that glm reports
wf_dev_glm <- deviance(wf_fit)
c(by_hand = wf_dev_hand, glm = wf_dev_glm)
#> by_hand glm
#> 7.63817 7.63817Both give \(\operatorname{dev}= 7.638\).
RemarkRemark 2 (Why the deviance measures a shortfall)
Remark 2 (Why the deviance measures a shortfall). For labels \(y_i\in {\mathopen{}\left\{0, 1\right\}\mathclose{}}\) the largest log-likelihood any model can reach is \(0\), by giving each label probability \(1\). So \(\operatorname{dev}\) is twice the amount by which the maximized log-likelihood falls short of that largest value.
DefDefinition 5 (Null deviance)
Definition 5 (Null deviance) The null deviance is the deviance of the model with an intercept and no other predictors:
\[ \operatorname{dev}_0= -2 n \mathopen{}\left\{\bar y\ln \bar y+ \mathopen{}\left(1 - \bar y\right)\mathclose{} {\ln\mathopen{}\left(1 - \bar y\right)\mathclose{}}\right\}\mathclose{}, \tag{4}\]
where \(\bar y\) is the fraction of labels equal to \(1\).
ExmExample 9 (The null deviance of the six-patient trial)
Example 9 (The null deviance of the six-patient trial) Three of the six patients recovered, so \(\bar y= 1/2\) and Equation 4 gives \(\operatorname{dev}_0= 2 \times 6 \ln 2\).
Both give \(\operatorname{dev}_0= 8.318\).
RemarkRemark 3 (Where the null deviance formula comes from)
Remark 3 (Where the null deviance formula comes from). The intercept-only model gives every observation the same fitted probability, and the maximizing value is the overall fraction \(\bar y\) of labels equal to \(1\). Substituting \(\bar y\) for each \(\hat\pi_i\) in Equation 3 gives Equation 4.
LemTheorem 3 (The drop in deviance is the likelihood ratio statistic)
Theorem 3 (The drop in deviance is the likelihood ratio statistic) Let \(\operatorname{dev}_0\) be the null deviance and \(\operatorname{dev}\) the deviance of a logistic regression with \(q\) predictors besides the intercept. Then the likelihood ratio statistic of \(H_0\): all \(q\) slopes are zero is
\[ \Lambda= \operatorname{dev}_0- \operatorname{dev}. \]
ProofProof
Proof. In the notation of the likelihood ratio statistic, the unrestricted MLE \(\hat\theta_{\text{ML}}\) is the fitted coefficient vector \(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\), and \(\hat{\theta}_0\) is the fit of the intercept-only model that \(H_0\) imposes. Then \({\ell\mathopen{}\left(\hat\theta_{\text{ML}}\right)\mathclose{}} = -\operatorname{dev}/ 2\) and \({\ell\mathopen{}\left(\hat{\theta}_0\right)\mathclose{}} = -\operatorname{dev}_0/ 2\). Substituting both into the definition gives \(\Lambda= 2 \mathopen{}\left(-\operatorname{dev}/2 + \operatorname{dev}_0/2\right)\mathclose{} = \operatorname{dev}_0- \operatorname{dev}\).
ExmExample 10 (The likelihood ratio statistic of the treatment effect)
Example 10 (The likelihood ratio statistic of the treatment effect) In the six-patient trial, \(q = 1\).
# drop in deviance, the likelihood ratio statistic
wf_g2 <- wf_fit$null.deviance - deviance(wf_fit)
wf_g2
#> [1] 0.679596The statistic is \(0.68\).
CorCorollary 2 (Distribution of the drop in deviance)
Corollary 2 (Distribution of the drop in deviance) Suppose the regularity conditions of Wilks’ theorem hold for the logistic regression model. Then for a large sample the drop in deviance \(\operatorname{dev}_0- \operatorname{dev}\) has approximately a \(\chi^2_q\) distribution when \(H_0\) holds.
ProofProof
Proof. By Theorem 3 the drop in deviance is the likelihood ratio statistic of \(q\) constraints, and Wilks’ theorem gives its approximate distribution.
ExmExample 11 (Testing the treatment effect)
Example 11 (Testing the treatment effect) For the statistic of Example 10, the p-value comes from the chi-square distribution with \(q = 1\) degree of freedom.
# upper tail of the chi-square distribution with 1 degree of freedom
wf_g2_p <- pchisq(wf_g2, df = 1, lower.tail = FALSE)
wf_g2_p
#> [1] 0.409726The p-value is \(0.41\), so six patients give no evidence of a treatment effect.
5 Covariate patterns and the saturated model
DefDefinition 6 (Covariate pattern)
Definition 6 (Covariate pattern) A covariate pattern is a distinct value of the predictor vector among the observations \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_1, \dots, \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_n\). The \(K\) distinct values are written \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^{(1)}, \dots, \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^{(K)}\), so \(K\le n\).
DefDefinition 7 (Replicates)
Definition 7 (Replicates) The replicates of the covariate pattern \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^{(k)}\) are the observations \(i\) with \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i= \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^{(k)}\). Their number is written \(m_k\), so \(\sum_{k=1}^{K}m_k = n\).
Note 4: Source: covariate patterns and replicates
Dobson and Barnett [7, Sec. 5.5] introduces both terms: observations that share a combination of factor levels and the same values of any continuous predictors have the same covariate pattern, and are called replicates.
DefDefinition 8 (Events at a covariate pattern)
Definition 8 (Events at a covariate pattern) The number of events at the covariate pattern \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^{(k)}\) is the sum of the binary outcomes over its replicates:
\[ s_k = \sum_{i:\, \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i= \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^{(k)}} y_i, \]
so \(0 \le s_k \le m_k\).
ExmExample 12 (Covariate patterns in the low-birth-weight data)
Example 12 (Covariate patterns in the low-birth-weight data) The birthwt data record whether each of 189 babies had low birth weight (low = 1, under 2.5 kg), along with the mother’s smoking status and race, among other predictors [8]. Take smoke (0 or 1) and race (1, 2 or 3) as the predictors. Counting the distinct pairs, and the babies and low-weight births at each, gives one row per covariate pattern.
There are \(K= 6\) covariate patterns among \(n= 189\) babies. The pattern with the most replicates has smoke = 0 and race = 3, and it has 55 of them.
DefDefinition 9 (Saturated model (full model, maximal model))
Definition 9 (Saturated model (full model, maximal model)) The saturated model (also called the full model or maximal model [7, Sec. 5.5]) for binary outcomes gives each pattern its own free event probability:
\[ {\Pr\mathopen{}\left(Y_i= 1 \mid \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i= \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}^{(k)}\right)\mathclose{}} = \pi_k, \qquad k = 1, \dots, K, \tag{5}\]
with the outcomes independent given the predictors.
Source: Wikipedia contributors [9].
Note 5: Source: saturated model
Dobson and Barnett [7, Sec. 5.5] defines the saturated model as the model with the maximum number of parameters that can be estimated, one per distinct covariate pattern, and says it is also called a maximal or full model. In model selection, “full model” often means something else: the model with every candidate predictor [7, Sec. 6.3]. That model is generally smaller than the saturated model.
LemTheorem 4 (The saturated model fits each pattern’s observed proportion)
Theorem 4 (The saturated model fits each pattern’s observed proportion) The maximum likelihood estimate of \(\pi_k\) in the saturated model is the observed proportion of events among the replicates of pattern \(k\):
\[ \hat\pi_k = \frac{s_k}{m_k}. \]
ProofProof
Proof. Independence gives each pattern its own log-likelihood term, and the outcomes at pattern \(k\) contribute \(s_k \ln \pi_k + \mathopen{}\left(m_k - s_k\right)\mathclose{} {\ln\mathopen{}\left(1 - \pi_k\right)\mathclose{}}\). The terms share no parameters, so each is maximized on its own. First suppose \(0 < s_k < m_k\), so that both logarithms are defined for \(0 < \pi_k < 1\). Setting the derivative with respect to \(\pi_k\) to zero gives \(s_k / \pi_k - \mathopen{}\left(m_k - s_k\right)\mathclose{} / \mathopen{}\left(1 - \pi_k\right)\mathclose{} = 0\). Multiplying through by \(\pi_k \mathopen{}\left(1 - \pi_k\right)\mathclose{}\) gives \(s_k \mathopen{}\left(1 - \pi_k\right)\mathclose{} = \mathopen{}\left(m_k - s_k\right)\mathclose{} \pi_k\). Expanding gives \(s_k - s_k \pi_k = m_k \pi_k - s_k \pi_k\). Adding \(s_k \pi_k\) to both sides gives \(s_k = m_k \pi_k\). Dividing both sides by \(m_k\) gives \(\hat\pi_k = s_k / m_k\). The log-likelihood term increases to this point and decreases after it, so the point is a maximum. Now suppose \(s_k = 0\). A count of \(0\) times a logarithm is taken as \(0\), so the term is \(m_k {\ln\mathopen{}\left(1 - \pi_k\right)\mathclose{}}\), which decreases in \(\pi_k\) and so peaks at \(\pi_k = 0\). Next suppose \(s_k = m_k\). By the same convention, the term is \(m_k \ln \pi_k\), which peaks at \(\pi_k = 1\). Both agree with \(s_k / m_k\).
DefDefinition 10 (Deviance relative to the saturated model)
Definition 10 (Deviance relative to the saturated model) The deviance relative to the saturated model of a fitted model is twice the amount by which its maximized log-likelihood \({\ell\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\right)\mathclose{}}\) falls short of the maximized log-likelihood \(\ell_{\mathrm{sat}}\) of the saturated model:
\[ \operatorname{dev}_{\mathrm{sat}}= 2 \mathopen{}\left(\ell_{\mathrm{sat}}- {\ell\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\right)\mathclose{}}\right)\mathclose{}. \tag{6}\]
LemTheorem 5 (The deviance relative to the saturated model, from grouped counts)
Theorem 5 (The deviance relative to the saturated model, from grouped counts) For a fitted logistic regression with fitted event probabilities \(\hat\pi_k\) at the \(K\) covariate patterns,
\[ \operatorname{dev}_{\mathrm{sat}} = 2 \sum_{k=1}^{K}\mathopen{}\left\{s_k {\ln\mathopen{}\left(\frac{s_k}{m_k \hat\pi_k}\right)\mathclose{}} + \mathopen{}\left(m_k - s_k\right)\mathclose{} {\ln\mathopen{}\left(\frac{m_k - s_k}{m_k \mathopen{}\left(1 - \hat\pi_k\right)\mathclose{}}\right)\mathclose{}}\right\}\mathclose{}, \tag{7}\]
where a term with \(s_k = 0\) or \(s_k = m_k\) contributes \(0\) in its zero-count part.
ProofProof
Proof. Grouping the outcomes by pattern, the fitted model’s log-likelihood is \({\ell\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\right)\mathclose{}} = \sum_{k=1}^{K}\mathopen{}\left\{s_k \ln \hat\pi_k + \mathopen{}\left(m_k - s_k\right)\mathclose{} {\ln\mathopen{}\left(1 - \hat\pi_k\right)\mathclose{}}\right\}\mathclose{}\). By Theorem 4, the saturated log-likelihood is the same sum with \(\hat\pi_k\) replaced by \(s_k / m_k\). Subtracting the fitted from the saturated log-likelihood, term by term, replaces each \({\ln\mathopen{}\left(s_k / m_k\right)\mathclose{}} - \ln \hat\pi_k\) by \({\ln\mathopen{}\left(s_k / \mathopen{}\left(m_k \hat\pi_k\right)\mathclose{}\right)\mathclose{}}\), and likewise for the non-events. Multiplying by \(2\), as Equation 6 does, gives Equation 7. A count of \(0\) times \(\ln 0\) is taken as \(0\), its limit.
Note 6: Source: deviance from grouped counts
Dobson and Barnett [7, Sec. 5.6.1] derives the deviance of a binomial model relative to the saturated model, and Section 7.5 writes it in the form of Equation 7. Section 7.5 also notes that when every observation has its own covariate pattern, so each count is \(0\) or \(1\), this deviance is not a useful measure of fit.
CorCorollary 3 (The deviance is the saturated-model deviance shifted by a constant)
Corollary 3 (The deviance is the saturated-model deviance shifted by a constant) The deviance of a fitted logistic regression and its deviance relative to the saturated model differ by a constant that depends on the data and not on the fitted model:
\[ \operatorname{dev}= \operatorname{dev}_{\mathrm{sat}}- 2 \ell_{\mathrm{sat}}. \]
ProofProof
Proof. Definition 4 gives \(\operatorname{dev}= -2 {\ell\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\right)\mathclose{}}\). Distributing the \(2\) in Equation 6 gives \(\operatorname{dev}_{\mathrm{sat}}= 2 \ell_{\mathrm{sat}}- 2 {\ell\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\right)\mathclose{}}\). Substituting \(\operatorname{dev}= -2 {\ell\mathopen{}\left(\hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\right)\mathclose{}}\) into this expression gives \(\operatorname{dev}_{\mathrm{sat}}= 2 \ell_{\mathrm{sat}}+ \operatorname{dev}\). Subtracting \(2 \ell_{\mathrm{sat}}\) from both sides gives the claim.
RemarkRemark 7 (Why the shift does not change the drop in deviance)
Remark 7 (Why the shift does not change the drop in deviance). The constant \(2 \ell_{\mathrm{sat}}\) in Corollary 3 is the same for every model fit to the same data, because the saturated model depends only on the data. So the difference between the deviances of two models, such as the drop in deviance of Theorem 3, is the same whether it is computed from \(\operatorname{dev}\) or from \(\operatorname{dev}_{\mathrm{sat}}\). The two deviances differ in their size, though: when every pattern has one replicate, the saturated model fits every label exactly, \(\ell_{\mathrm{sat}}= 0\), and \(\operatorname{dev}_{\mathrm{sat}}= \operatorname{dev}\).
ExmExample 13 (The deviance of the low-birth-weight model, ungrouped and grouped)
Example 13 (The deviance of the low-birth-weight model, ungrouped and grouped) Fit a logistic regression of low birth weight on smoking and race, once to the 189 individual babies and once to the 6 grouped counts of Example 12. R’s deviance() is the deviance of Definition 4 for the first fit and the deviance relative to the saturated model of Definition 10 for the second. We check the second against Equation 7 computed by hand.
# fit to the individual babies
bw_fit <- glm(low ~ smoke + race, family = binomial, data = bw)
# fit to the grouped counts: events and non-events at each pattern
bw_fit_grouped <- glm(
cbind(nevents, nreps - nevents) ~ smoke + race,
family = binomial,
data = bw_grouped
)
# fitted probability at each pattern
bw_phat <- fitted(bw_fit_grouped)
# one term of the grouped formula; a zero count contributes 0
bw_term <- function(count, expected) {
ifelse(count == 0, 0, count * log(count / expected))
}
bw_sat_dev_hand <- 2 * sum(
bw_term(bw_grouped$nevents, bw_grouped$nreps * bw_phat) +
bw_term(
bw_grouped$nreps - bw_grouped$nevents,
bw_grouped$nreps * (1 - bw_phat)
)
)
c(
by_hand = bw_sat_dev_hand,
glm_grouped = deviance(bw_fit_grouped),
glm_ungrouped = deviance(bw_fit)
)
#> by_hand glm_grouped glm_ungrouped
#> 3.15694 3.15694 219.97471The deviance relative to the saturated model is \(\operatorname{dev}_{\mathrm{sat}}= 3.157\) on the grouped counts, and the deviance is \(\operatorname{dev}= 219.975\) on the individual babies.
ExmExample 14 (The shift between the two deviances)
Example 14 (The shift between the two deviances) For the two fits of Example 13, Corollary 3 says that \(\operatorname{dev}- \operatorname{dev}_{\mathrm{sat}}= -2 \ell_{\mathrm{sat}}\). The left side is 216.818. We compute the right side from the observed proportions of Theorem 4.
# twice the saturated log-likelihood, from the observed proportions
bw_phat_sat <- bw_grouped$nevents / bw_grouped$nreps
bw_two_llik_sat <- 2 * sum(
ifelse(bw_grouped$nevents == 0, 0, bw_grouped$nevents * log(bw_phat_sat)) +
ifelse(
bw_grouped$nreps == bw_grouped$nevents,
0,
(bw_grouped$nreps - bw_grouped$nevents) * log(1 - bw_phat_sat)
)
)
c(
difference = deviance(bw_fit) - deviance(bw_fit_grouped),
minus_two_llik_sat = -bw_two_llik_sat
)
#> difference minus_two_llik_sat
#> 216.818 216.8186 Confidence intervals for the coefficients
LemTheorem 6 (Approximate distribution of the fitted coefficients)
Theorem 6 (Approximate distribution of the fitted coefficients) Suppose the outcomes follow the logistic regression model. Let \(X\) be the matrix whose rows are \({\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{x}}\vphantom{x}}_i}^{\top}\), and let \(W\) be the diagonal matrix with entries \(\pi_i \mathopen{}\left(1 - \pi_i\right)\mathclose{}\). Under regularity conditions, for a large sample the fitted coefficients are approximately Gaussian,
\[ \hat{\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}}\approx \operatorname{N}\mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}, {\mathopen{}\left({X}^{\top} W X\right)\mathclose{}}^{-1}\right)\mathclose{}. \tag{8}\]
ProofProof
Proof. The log-likelihood has score \({X}^{\top} \mathopen{}\left(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}- \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\pi}}\vphantom{\pi}}\right)\mathclose{}\), and its Hessian is \(-{X}^{\top} W X\), which does not involve the outcomes. The observed information and expected information are therefore both \({X}^{\top} W X\). The claim is the Gaussian approximation of the central limit theorem for MLEs, with this information in place of the inverse covariance. The result is stated for logistic regression by Hastie et al. [6, Sec. 4.4.3]. For independent but non-identically distributed outcomes it holds under conditions on the predictors that Fahrmeir and Kaufmann [10] give.
CorCorollary 4 (Wald interval for a logistic coefficient)
Corollary 4 (Wald interval for a logistic coefficient) Let \(\hat W\) be \(W\) with each \(\pi_i\) replaced by its fitted value. The estimated standard error of \(\hat{\beta}_{j}\) is the square root of the diagonal entry of \({\mathopen{}\left({X}^{\top} \hat W X\right)\mathclose{}}^{-1}\) in the row and column of \(\beta_{j}\); because \(\beta_{0}\) comes first, that is the \((j+1)\)th diagonal entry counting from \(1\), and a Wald confidence interval with level \(1 - \alpha\) for \(\beta_{j}\) is
\[ \hat{\beta}_{j} \pm z_{1 - \alpha/2} \, {\mathop{\widehat{\operatorname{se}}}\nolimits\mathopen{}\left(\hat{\beta}_{j}\right)\mathclose{}}. \]
ExmExample 15 (A Wald interval from six patients)
Example 15 (A Wald interval from six patients) For the fit of Example 1, vcov returns \({\mathopen{}\left({X}^{\top} \hat W X\right)\mathclose{}}^{-1}\).
# standard errors: square roots of the diagonal of (X'WX)^-1
wf_se <- sqrt(diag(vcov(wf_fit)))
wf_slope <- coef(wf_fit)[["wf_x"]]
# 95% Wald interval for the slope
wf_ci <- wf_slope + c(-1, 1) * qnorm(0.975) * wf_se[["wf_x"]]
wf_summary <- c(
slope = wf_slope, se = wf_se[["wf_x"]], lower = wf_ci[1], upper = wf_ci[2]
)
round(wf_summary, 3)
#> slope se lower upper
#> 0.693 0.866 -1.004 2.391The interval for \(\beta_{1}\) is \((-1, 2.39)\). It contains \(0\), consistent with the test of Example 11.
CorCorollary 5 (Interval for the odds factor of a logistic coefficient)
Corollary 5 (Interval for the odds factor of a logistic coefficient) If \(\mathopen{}\left(L, U\right)\mathclose{}\) is a confidence interval for \(\beta_{j}\) from Corollary 4, then \(\mathopen{}\left(e^{L}, e^{U}\right)\mathclose{}\) is a confidence interval with the same level for \(e^{\beta_{j}}\), the factor by which a one-unit increase in \(x_j\), holding the other predictors fixed, multiplies the odds of \(y= 1\).
ProofProof
Proof. The exponential function is increasing, so \(L \le \beta_{j} \le U\) holds exactly when \(e^{L} \le e^{\beta_{j}} \le e^{U}\).
ExmExample 16 (How little six patients tell us)
Example 16 (How little six patients tell us) The two groups of Example 1 sit two units apart on the \(x\) scale, so the odds ratio between them is \(e^{2 \beta_{1}}\).
The fitted odds of recovery in the treatment group are \(4\) times those in the control group, with interval \((0.13, 119.2)\). The interval runs from a strong harm to a large benefit. Its width shows how little information six observations carry about the slope.
7 Further reading
- Hua Zhou’s course page for “UCLA Biostat 257 - Statistical Computing” (2025 Spring) has lectures on fitting generalized linear models and on quasi-Newton optimization. Its source is on GitHub under the MIT License.
- David Rocke’s EPI 204 lecture Logistic Regression and Likelihood (slides) fits logistic regression by maximum likelihood and compares models by deviance and AIC [11].
- Josh Starmer’s StatQuest video Logistic Regression Details Pt 2: Maximum Likelihood shows how maximum likelihood estimation fits a logistic regression curve [12].
Definitions on this page
- Definition 1 (Logistic regression)
- Definition 2 (Odds multiplier)
- Definition 3 (Iteratively reweighted least squares (IRLS))
- Definition 4 (Deviance)
- Definition 5 (Null deviance)
- Definition 6 (Covariate pattern)
- Definition 7 (Replicates)
- Definition 8 (Events at a covariate pattern)
- Definition 9 (Saturated model (full model, maximal model))
- Definition 10 (Deviance relative to the saturated model)
References
[1]
Wikipedia contributors, “Likelihood-ratio test — Wikipedia, the free encyclopedia.” 2026. Available: https://en.wikipedia.org/wiki/Likelihood-ratio_test
[2]
Wikipedia contributors, “Logistic regression — Wikipedia, the free encyclopedia.” 2026. Available: https://en.wikipedia.org/wiki/Logistic_regression
[3]
H. Zhou, “Newton’s method and Fisher scoring.” Lecture notes for Biostat M257, Statistical Computing, University of California, Los Angeles, 2025. Accessed: Oct. 26, 2026. [Online]. Available: https://ucla-biostat-257.github.io/2025spring/slides/23-newton/newton.html
[4]
Wikipedia contributors, “Deviance (statistics) — Wikipedia, the free encyclopedia.” 2026. Available: https://en.wikipedia.org/wiki/Deviance_%28statistics%29
[5]
J. A. Nelder and R. W. M. Wedderburn, “Generalized linear models,” Journal of the Royal Statistical Society: Series A (General), vol. 135, no. 3, pp. 370–384, 1972, doi: 10.2307/2344614.
[6]
T. Hastie, R. Tibshirani, and J. Friedman, The elements of statistical learning: Data mining, inference, and prediction, 2nd ed. Springer, 2009. doi: 10.1007/978-0-387-84858-7.
[7]
A. J. Dobson and A. G. Barnett, An introduction to generalized linear models, 4th ed. CRC press, 2018. doi: 10.1201/9781315182780.
[8]
D. W. Hosmer and S. Lemeshow, Applied logistic regression. New York: Wiley, 1989.
[9]
Wikipedia contributors, “Saturated model — Wikipedia, the free encyclopedia.” 2026. Available: https://en.wikipedia.org/wiki/Saturated_model
[10]
L. Fahrmeir and H. Kaufmann, “Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models,” The Annals of Statistics, vol. 13, no. 1, pp. 342–368, 1985, doi: 10.1214/aos/1176346597.
[11]
D. M. Rocke, “EPI 204 lecture notes, Spring 2021.” University of California, Davis, 2021. Available: https://dmrocke.ucdavis.edu/Class/EPI204-Spring-2021/EPI204-Spring-2021.html
[12]
J. Starmer, “Logistic regression details, part 2: Maximum likelihood.” StatQuest, video on YouTube, 2018. Available: https://www.youtube.com/watch?v=BfKanl1aSG0