Last modified: 2026-10-27 23:24:14 (PDT)
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].
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].
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.
# 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.693The 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.
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.
Definition 2 (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
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.
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}}\).
Solution 1. Here the rows of \(X\) are \((1, -1)\) for the three patients in the group coded \(-1\) and \((1, 1)\) for the three treated patients, and \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{y}}\vphantom{y}}= {(0, 0, 1, 0, 1, 1)}^{\top}\).
1. At \(\mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{\beta}}\vphantom{\beta}}^* = \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{0}}\vphantom{0}}\) every fitted probability is \({\operatorname{expit}\mathopen{}\left(0\right)\mathclose{}} = 1/2\), so every weight is \(w_{i} = (1/2)(1/2) = 1/4\). The working responses are \(z_{i} = 0 + \mathopen{}\left(y_{i} - 1/2\right)\mathclose{} / (1/4) = 4 y_{i} - 2\):
\[ \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}} = {(-2, -2, 2, -2, 2, 2)}^{\top}. \]
2. With \(W = \tfrac{1}{4} \mathbf{I}_6\),
\[ {X}^{\top} W X = \tfrac{1}{4} \begin{bmatrix} 6 & 0 \\ 0 & 6 \end{bmatrix} = \begin{bmatrix} 3/2 & 0 \\ 0 & 3/2 \end{bmatrix}, \]
because the entries of the first column of \(X\) sum to \(6\), the squares in the second column sum to \(6\), and the two columns multiply to a sum of \(-3 + 3 = 0\). The first entry of \({X}^{\top} W \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}}\) is \(\tfrac{1}{4}(-2 - 2 + 2 - 2 + 2 + 2) = 0\), and the second is \(\tfrac{1}{4}\mathopen{}\left((-1)(-2 - 2 + 2) + (1)(-2 + 2 + 2)\right)\mathclose{} = \tfrac{1}{4}(2 + 2) = 1\), so \({X}^{\top} W \mathord{\smash[b]{\underset{\raise.3ex\hbox{$\scriptstyle\sim$}}{z}}\vphantom{z}} = {(0, 1)}^{\top}\).
3. The updated coefficients are
\[ \begin{bmatrix} 2/3 & 0 \\ 0 & 2/3 \end{bmatrix} \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 0 \\ 2/3 \end{bmatrix}. \]
The intercept stays at \(0\) and the slope moves from \(0\) to \(2/3\).
Theorem 1 (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
\[ \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}\]
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 5). 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} \]
Example 2 (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\).
Definition 3 (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].
Source
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.
Example 3 (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\).
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.
Definition 4 (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\).
Example 4 (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\).
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.
Theorem 2 (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}. \]
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}\).
Corollary 1 (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.
Proof. By Theorem 2 the drop in deviance is the likelihood ratio statistic of \(q\) constraints, and Wilks’ theorem gives its approximate distribution.
Example 6 (Testing the treatment effect) For the statistic of Example 5, the p-value comes from the chi-square distribution with \(q = 1\) degree of freedom.
The p-value is \(0.41\), so six patients give no evidence of a treatment effect.
Definition 5 (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\).
Definition 6 (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.
Definition 7 (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\).
Example 7 (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.
Definition 8 (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.
Theorem 3 (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}. \]
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\).
Definition 9 (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}\]
Theorem 4 (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.
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 3, 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.
Corollary 2 (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}}. \]
Proof. Definition 3 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.
Remark 7 (Why the shift does not change the drop in deviance). The constant \(2 \ell_{\mathrm{sat}}\) in Corollary 2 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 2, 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}\).
Example 8 (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 7. R’s deviance() is the deviance of Definition 3 for the first fit and the deviance relative to the saturated model of Definition 9 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.
Example 9 (The shift between the two deviances) For the two fits of Example 8, Corollary 2 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 3.
# 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.818Theorem 5 (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}\]
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 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.
Corollary 3 (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{}}. \]
Example 10 (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 6.
Corollary 4 (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 3, 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\).
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}\).
Example 11 (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.