# In a GLM, is the log likelihood of the saturated model always zero?

As part of the output of a generalised linear model, the null and residual deviance are used to evaluate the model. I often see the formulas for these quantities expressed in terms of the log likelihood of the saturated model, for example: https://stats.stackexchange.com/a/113022/22199, Logistic Regression : How to obtain a saturated model

The saturated model, as far as I understand it, is the model that perfectly fits the observed response. Thus, in most places I have seen, the log-likelihood of the saturated model is always given as zero.

Yet, the way the formula for deviance is given suggests that sometimes this quantity is non zero. (As if it is zero always, why bother including it?)

In what cases can it be non zero? If it is never non-zero, why include it in the formula for deviance?

If you really meant log-likelihood, then the answer is: it's not always zero.

For example, consider Poisson data: $$y_i \sim \text{Poisson}(\mu_i), i = 1, \ldots, n$$. The log-likelihood for $$Y = (y_1, \ldots, y_n)$$ is given by: $$\ell(\mu; Y) = -\sum_{i = 1}^n \mu_i + \sum_{i = 1}^n y_i \log \mu_i - \sum_{i = 1}^n \log(y_i!). \tag{*}$$

Differentiate $$\ell(\mu; Y)$$ in $$(*)$$ with respect to $$\mu_i$$ and set it to $$0$$ (this is how we obtain the MLE for saturated model): $$-1 + \frac{y_i}{\mu_i} = 0.$$ Solve this for $$\mu_i$$ to get $$\hat{\mu}_i = y_i$$, substituting $$\hat{\mu}_i$$ back into $$(*)$$ for $$\mu_i$$ gives that the log-likelihood of the saturated model is: $$\ell(\hat{\mu}; Y) = \sum_{i = 1}^n y_i(\log y_i - 1) -\sum_{i = 1}^n \log(y_i!) \neq 0$$ unless $$y_i$$ take very special values.

In the help page of the R function glm, under the item deviance, the document explains this issue as follows:

deviance up to a constant, minus twice the maximized log-likelihood. Where sensible, the constant is chosen so that a saturated model has deviance zero.

Notice that it mentioned that the deviance, instead of the log-likelihood of the saturated model is chosen to be zero.

Probably, what you really wanted to confirm is that "the deviance of the saturated model is always given as zero", which is true, since the deviance, by definition (see Section 4.5.1 of Categorical Data Analysis (2nd Edition) by Alan Agresti) is the likelihood ratio statistic of a specified GLM to the saturated model. The constant aforementioned in the R documentation is actually twice the maximized log-likelihood of the saturated model.

Regarding your statement "Yet, the way the formula for deviance is given suggests that sometimes this quantity is non zero.", it is probably due to the abuse of usage of the term deviance. For instance, in R, the likelihood ratio statistic of comparing two arbitrary (nested) models $$M_1$$ and $$M_2$$ is also referred to as deviance, which would be more precisely termed as the difference between the deviance of $$M_1$$ and the deviance of $$M_2$$, if we closely followed the definition as given in Agresti's book.

Conclusion

1. The log-likelihood of the saturated model is in general non-zero.

2. The deviance (in its original definition) of the saturated model is zero.

3. The deviance output from softwares (such as R) is in general non-zero as it actually means something else (the difference between deviances).

The following are the derivation for the general exponential-family case and another concrete example. Suppose that data come from exponential family (see Modern Applied Statistics with S, Chapter $$7$$): $$f(y_i; \theta_i, \varphi) = \exp[A_i(y_i\theta_i - \gamma(\theta_i))/\varphi + \tau(y_i, \varphi/A_i)]. \tag{1}$$ where $$A_i$$ are known prior weights and $$\varphi$$ are dispersion/scale parameter (for many cases such as binomial and Poisson, this parameter is known, while for other cases such as normal and Gamma, this parameter is unknown). Then the log-likelihood is given by: $$\ell(\theta, \varphi; Y) = \sum_{i = 1}^n A_i(y_i \theta_i - \gamma(\theta_i))/\varphi + \sum_{i = 1}^n \tau(y_i, \varphi/A_i).$$ As in the Poisson example, the saturated model's parameters can be estimated by solving the following score function: $$0 = U(\theta_i) = \frac{\partial \ell(\theta, \varphi; Y)}{\partial \theta_i} = \frac{A_i(y_i - \gamma'(\theta_i))}{\varphi}$$

Denote the solution of the above equation by $$\hat{\theta}_i$$, then the general form of the log-likelihood of the saturated model (treat the scale parameter as constant) is: $$\ell(\hat{\theta}, \varphi; Y) = \sum_{i = 1}^n A_i(y_i \hat{\theta}_i - \gamma(\hat{\theta}_i))/\varphi + \sum_{i = 1}^n \tau(y_i, \varphi/A_i). \tag{**}$$

In my previous answer, I incorrectly stated that the first term on the right side of $$(**)$$ is always zero, the above Poisson data example proves it is wrong. For a more complicated example, consider the Gamma distribution $$\Gamma(\alpha, \beta)$$ given in the appendix.

Proof of the first term in the log-likelihood of saturated Gamma model is non-zero: Given $$f(y; \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)}e^{-\beta y}y^{\alpha - 1}, \quad y > 0, \alpha > 0, \beta > 0,$$ we must do reparameterization first so that $$f$$ has the exponential family form $$(1)$$. It can be verified if letting $$\varphi = \frac{1}{\alpha},\, \theta = -\frac{\beta}{\alpha},$$ then $$f$$ has the representation: $$f(y; \theta, \varphi) = \exp\left[\frac{\theta y - (-\log(-\theta))}{\varphi}+ \tau(y, \varphi)\right],$$ where $$\tau(y, \varphi) = -\frac{\log \varphi}{\varphi} + \left(\frac{1}{\varphi} - 1\right)\log y - \log\Gamma(\varphi^{-1}).$$ Therefore, the MLEs of the saturated model are $$\hat{\theta}_i = -\frac{1}{y_i}$$. Hence $$\sum_{i = 1}^n \frac{1}{\varphi}[\hat{\theta}_iy_i - (-\log(-\hat{\theta}_i))] = \sum_{i = 1}^n \frac{1}{\varphi}[-1 - \log(y_i)] \neq 0,$$ unless $$y_i$$ take very special values.

• Is the loglikelihood zero if and only if the model can assign 100% probability to each of the possible outcomes?
– Alex
Dec 3, 2015 at 2:57
• I don't quite understand what you meant. But from my derivation you might conclude that it is $0$ if and only if the $\tau$ is identically $0$ and there is no dispersion parameter. Dec 3, 2015 at 2:59
• What I meant is that is that if $Y$ was a Poisson random variable, then $P(Y = y_i) < 1$ for any $i$ or Poisson mean, thus it is impossible to find any model parameter that gives a log likelihood of zero for the observed. Maybe I am completely misunderstanding the concept of a saturated model.
– Alex
Dec 3, 2015 at 3:25
• Maybe. I don't think your understanding is $100\%$ correct, for your argument also applies to binary logistic model, but in that case the log-likelihood for saturated model is indeed $0$. You may review the concept of saturated model first then you can understand how my argument works. Dec 3, 2015 at 3:39
• This is my thought for the binary logistic model, for each observation $y_i$, the likelihood is maximised when $P(Y = y_i) = 1$, which is found by setting the binomial mean to be $y_i$. This is what I interpreted to be the 'perfect' fit for a saturated model.
– Alex
Dec 3, 2015 at 3:50

Zhanxiong's answer is already great (+1), but here's a quick demonstration that the log-likelihood of the saturated model is $$0$$ for a logistic regression. I figured I would post because I haven't seen this TeX'd up on this site, and because I just wrote these up for a lecture.

The likelihood is $$L(\mathbf{y} ; \mathbf{X}, \boldsymbol{\beta}) = \prod_{i=1}^n f(y_i ; \mathbf{x}_i, \boldsymbol{\beta}) = \prod_{i=1}^n \pi_i^{y_i}(1-\pi_i)^{1-y_i} = \prod_{i=1}^n\left( \frac{\pi_i}{1-\pi_i}\right)^{y_i} (1 - \pi_i) \tag{1}$$ where $$\pi_i = \text{invlogit}(\mathbf{x}_i^\intercal \boldsymbol{\beta} )$$.

The log-likelihood is \begin{align*} \log L(\mathbf{y} ; \mathbf{X}, \boldsymbol{\beta}) &= \sum_{i=1}^n y_i \log \left( \frac{\pi_i}{1-\pi_i}\right) + \log(1-\pi_i) \\ &= \sum_{i=1}^n y_i \text{logit} \left( \pi_i \right) + \log(1-\pi_i) \\ &= \sum_{i=1}^n y_i \mathbf{x}_i^\intercal \boldsymbol{\beta} + \log( 1 - \text{invlogit}(\mathbf{x}_i^\intercal \boldsymbol{\beta} )) \\ &= \sum_{i=1}^n y_i \mathbf{x}_i^\intercal \boldsymbol{\beta} + \log( \text{invlogit}( - \mathbf{x}_i^\intercal \boldsymbol{\beta} )) \\ &= \sum_{i=1}^n y_i \mathbf{x}_i^\intercal \boldsymbol{\beta} - \log( 1 + \exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}] )) \end{align*}

If you take the derivatives with respect to all of the coefficients you get $$\nabla \ell(\boldsymbol{\beta}) = \sum_{i=1}^n y_i \mathbf{x}_i - \frac{\exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}]}{( 1 + \exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}] ) }\mathbf{x}_i \tag{2}.$$

Setting this expression equal to $$\mathbf{0}$$ and solving for $$\boldsymbol{\beta}$$ will give you your answer. Usually this can't be done analytically, which explains the popularity/necessity of using iterative algorithms to fit this model, but in the case of a saturated model, it is possible.

To find the saturated model, we give each row it's own coefficent. So $$\boldsymbol{\beta} \in \mathbb{R}^n$$ and the design matrix times the coefficient vector is $$\mathbf{X}\boldsymbol{\beta} = \begin{bmatrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1\\ \end{bmatrix} \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{bmatrix}.$$

Note that in particular, $$\mathbf{x}_i^\intercal \boldsymbol{\beta} = \beta_i$$.

So taking the $$j$$th row of equation (2) gives us $$\sum_{i=1}^n y_i x_{i,j} = \sum_{i=1}^n\frac{\exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}]}{( 1 + \exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}] ) }x_{i,j}$$

which can only be true if for each observation $$i$$:

$$y_i = \text{invlogit}(\beta_i )$$ or in other words each $$\beta_i$$ is plus or minus infinity (if $$y_i$$ is $$1$$ or $$0$$, respectively). We can plug these parameters back into (1) to get the maximized likelihood: $$\prod_{i=1}^n \hat{\pi}_i^{y_i}(1-\hat{\pi}_i)^{1-y_i} = 1^n = 1.$$ Clearly the log of this is $$0$$.

• But this assumes ungrouped data. If you have groups with $n_i>1$ (and the same covariate values) (in R, forexample using the form glm( cbind(k, n-k) ~ x + ... ) then the saturated model do not have loglikelihood zero. Jun 27, 2019 at 7:18
• @kjetilbhalvorsen oh good point. I never tried that let me check Jun 27, 2019 at 13:52
• Taylor: Hi, i have a question concerning your post please. In the part in which there is a multiplication of 2 matrix, how could you be sure that the design matrix X contains 1 on its diagonal and all 0 elsewhere ? I think it is not true. Thank you for your response! Jun 8, 2021 at 22:04
• @InTheSearchForKnowledge are you talking about this equation: $\mathbf{X}\boldsymbol{\beta} = \begin{bmatrix} 1 & 0 & \cdots & 0\\ 0 & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1\\ \end{bmatrix} \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{bmatrix}$? It's only true sometimes--when the model is saturated or in other words, when there is a $\beta$ parameter for each individual observation/row/datapoint. Jun 9, 2021 at 0:40
• @Taylor: Hi, thank you very much for your response. What I mean is that the matrix $\mathbb{X}$ you mention is not the same as the original design matrix. Could you please tell me how to construct the "new" design matrix X from the original design matrix ? And moreover, I see that you estimate the paramter $\beta$ instead of estimate the parameter of distribution (i.e. the $\pi_i$ like in the post of Zhanxiong), could you tell me are the 2 method of estimation (via $\beta$ et via $\mu_i$) equivalent, and why ? Thank you so so much for your help! Jun 14, 2021 at 18:41

@Alex: yes, thats right. at least for discrete distributions. for continuous distributions, it would come down to letting the density be equal 1, which is not necessarily meaningful and therefore not a sensible thing to try and achieve. slightly more generally, the log-likelihood of the saturated model gives you an upper bound for the performance of any model that follows your assumption of the underlying distribution family. In other words, the log-likelihood of a saturated binomial model it is "as good as it gets" for the given data set (X,Y) assuming Y is binomial. It makes sense to compare your glm model to this upper bound as opposed to, say, 100% (or similar), since your model is inherently constrained by your assumption on the response distribution. The deviance as defined by @Zhanxiong therefore gives you a good idea how well your model performs w.r.t to its inherit limitations that come from assuming a certain response type.