Estimating deviance from the saturated model

Question

I am trying to consider how one would actually calculate deviance of a linear model in R. I am trying to implement: $$D = -2 \times (L_p - L_s)$$

Where $L_p$ and $L_s$ are the log-likelihood of the proposed and saturated model respectively.

I have seen a few things that have confused me slightly:

1. Gordon Smyth mentioned that you would never actually fit a saturated model in a comment on this question. If that's the case, how do we know $L_s$?
2. Mike Love calculated deviance by comparing $-2\times L_p$ with $-2 \times$a simulated negative-binomial distribution on Biostars. I expect the answer to (1) might be in here somewhere, but I am unsure what is being calculated here.

I also understand that the deviance of the saturated model is zero in a logistic regression, so we can take $D = -2 \times L_p$ as the deviance.

How about for linear regression? How would I estimate $L_s$ without fitting a saturated model in R? I'm definitely missing some understanding here, so any reading material you can recommend to would be appreciated. Thanks.

Answer 1

Normally when we compare likelihoods it is only meaningful to compare ratios of likelihoods (or differences of log-likelihoods). Therefore, if we want some kind of absolute measure (starting from zero) of the goodness of fit of a model, we need to compare it either to the saturated model (how good is this model relative to the best possible fit?) or the null model (how good is it relative to the worst possible/most naive fit?). (See also here)

For models with a discrete response, the saturated model is not perfect. Consider the trivial case of fitting a Poisson model to a single data point $y$. The MLE is $\lambda =y$, but the probability of getting an outcome of $y$ when $\lambda=y$ is not 1! Instead,

$$ L_s = -\lambda + y \log \lambda - \log(y!) = - y + y \log y - \log(y!) $$

and

$$ L_p = -\hat y + y \log \hat y - \log(y!) $$

so $$ L_p -L_s = y (\log y - \log \hat y) - (y - \hat y) = y \log(y/\hat y) - (y-\hat y) $$

If you look at the R code for poisson()$dev.resid (the deviance-residual function for the Poisson family), you'll see the same thing, up to some complications about dealing with weights and zero-valued responses (for which the log-likelihood of the saturated model is exactly 0, and where the expression given above will blow up)

function (y, mu, wt) 
{
    r <- mu * wt
    p <- which(y > 0)
    r[p] <- (wt * (y * log(y/mu) - (y - mu)))[p]
    2 * r
}

For continuous-valued responses, the saturated model is perfect (in some sense), and we set $L_s=0$ (to be honest I'm not quite sure how this works formally; if we were to fit the saturated model to a Gaussian, for example, we would get an estimated residual variance of zero and hence an infinite probability density ...)