5
$\begingroup$

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'))

enter image description here

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.

$\endgroup$
0

1 Answer 1

2
$\begingroup$

I have found the answer in Modern Applied Statistics with S by W.N. Venables and B.D. Ripley, p.189:

The deviance residuals $d_i$ are defined as the signed square roots of the summands of the deviance (7.5) taking the same sign as $y_i - \hat\mu_i$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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