I'm trying to understand deviance and deviance residuals using a simple Poisson regression model as an example. Let's say we have a response variable $$ y_i \sim \text{Pois}(\lambda_i)$$ and we assume $$ \lambda_i = \exp(a + b x_i) $$
## generate random data
y <- rpois(25, 2)
x <- rnorm(25) + (log(replace(y, y == 0, 0.25)))
We can estimate $\hat a$ and $\hat b$ by maximising the log-likelihood function $$ \ell(a,b) = \sum_{i=1}^n \left[ y_i (a+bx_i) - \exp(a+bx_i) - \log(y_i!) \right ]$$
mod1 <- optim(c(1, 1), function(p)
with(list(mu=p[1] + p[2] * x), -sum(y * mu - exp(mu))))
mod2 <- glm(y ~ x, family='poisson', data=data.frame(y=y, x=x))
## check that the parameters coincide
mod1$par
coefficients(mod2)
Now, we can rewrite the log-likelihood function in terms of $\lambda$ $$ \ell(\lambda) = \sum_{i=1}^n \left[ y_i \lambda_i - \exp(\lambda_i) - \log(y_i!) \right ]$$ The log-likelihood of the model is $\ell(\hat\lambda)$ where $\hat\lambda=\exp(\hat a + \hat b x)$.
## eta = a + bx
eta <- cbind(1, x) %*% mod1$par
tail(eta)
tail(predict(mod2))
## yhat = exp(eta)
yhat <- exp(eta)
tail(yhat)
tail(fitted.values(mod2))
## log-likelihood
sum(dpois(y, yhat, log=TRUE))
logLik(mod2)[1]
On the other hand, the log-likelihood of the "saturated" model is $\ell(\hat\lambda_s)$ where $\hat\lambda_s=y$. The deviance of the model (the so-called "residual deviance") is $$D=2\left[\ell(\hat\lambda_s)-\ell(\hat\lambda)\right]$$ in other words it's two times the difference between the likelihood of the saturated and the likelihood of the model.
## deviance: D = 2(logLik_fullmodel - logLik_model)
2 * (sum(dpois(y, y, log=TRUE)) - sum(dpois(y, yhat, log=TRUE)))
mod2$deviance
Now is when I get lost. I know that $$d^2=2\left[ \hat y - \hat y_s - y (\log\hat\lambda-\log\hat\lambda_s) \right]$$ are the squared "deviance residuals" and it turns out that $\sum_i d^2_i$ equals the residual deviance.
tail(2 * (yhat - y - y * (eta - log(y))))
tail(resid(mod2)^2)
sum(resid(mod2)^2)
However I'm at a loss where this equation comes from. Can somebody shed light on this?
Edit
Based on the fact that $\sum d^2 = D$ I worked out that $$d^2_i = 2 \left[ \ell_i(\hat\lambda_s) - \ell_i(\hat\lambda) \right] $$ So basically $\ell_i(\hat\lambda_s)$ and $\ell_i(\hat\lambda)$ are two estimates of $\log\Pr(y_i)$. Visually, suppose $y_i=2$, then these 2 estimates are the 2 bars corresponding to "y=2" in the plot below
barplot(rbind(dpois(0:4, y[1], log=TRUE),
dpois(0:4, yhat[1], log=TRUE)),
beside=TRUE, names.arg=paste0('y=', 0:4),
legend.text=c('saturated', 'model'), args.legend=list(x='bottom'))
However given that $\ell_i(\hat\lambda_s) - \ell_i(\hat\lambda) \geq 0$ I still don't know where the sign of $d_i$ comes from.
