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

Question

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?

Answer 1

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.

Answer 2

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$.

Answer 3

@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.