11
$\begingroup$

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?

$\endgroup$
16
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Is the loglikelihood zero if and only if the model can assign 100% probability to each of the possible outcomes? $\endgroup$ – Alex Dec 3 '15 at 2:57
  • $\begingroup$ 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. $\endgroup$ – Zhanxiong Dec 3 '15 at 2:59
  • $\begingroup$ Your derivation is very good but the the formal proof is a little bit above my head at the moment. Thank you for your example with the Poisson model. What I took out of this example is that the Poisson model cannot assign 100% probability to the observed outcome given any value for the Poisson mean, thus the likelihood cannot be zero. $\endgroup$ – Alex Dec 3 '15 at 3:08
  • $\begingroup$ The statement "model assign $100\%$ probability to the observed outcome" sounds weird to me. Do you mean that given the observations $y_1, \ldots, y_n$, and if $Y$ is a Poisson random variable, $P(Y= y_1) + P(Y = y_2) + \cdots + P(Y = y_n) < 1$? $\endgroup$ – Zhanxiong Dec 3 '15 at 3:13
  • 1
    $\begingroup$ 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. $\endgroup$ – Alex Dec 3 '15 at 3:25
2
$\begingroup$

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{1}{( 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{1}{( 1 + \exp[ \mathbf{x}_i^\intercal \boldsymbol{\beta}] ) }x_{i,j} $$

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

$$ y_i = \frac{1}{( 1 + \exp[ \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$.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – kjetil b halvorsen Jun 27 at 7:18
  • $\begingroup$ @kjetilbhalvorsen oh good point. I never tried that let me check $\endgroup$ – Taylor Jun 27 at 13:52
1
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.