0
$\begingroup$

I've been trying to generate the deviance residuals for a model I've made using R. It's a Gompertz regression with a number of covariates in the regression and the data is left-truncated, right-censored with an event being a death (i.e. there is an entry age, exit age and an indicator of censoring).

I've used R's phreg() function to handle the truncation/censoring and was hoping to assess the validity of the model using deviance residuals but I'm not entirely sure how to calcuate them straight from the phreg() output.

A manual attempt I tried is as follows:

A paper I've been following "APPLYING SURVIVAL MODELS TO PENSIONER MORTALITY DATA" Richards (2008) states that the the number of deaths in any age range can be assumed Poisson distributed with parameter equal to the sum of the integrated hazard function over that age range over all individuals at risk.

The actual number of deaths is just the calculated from the data.

So, from what I gather, the deviance residuals are then:

$r_{i}=Sign(D_{i}-\hat{D_{i}})\,sqrt(2\times(D_{i}\times log(D_{i}/\hat{D_{i}})-(D_{i}-\hat{D_{i}}))$

where, $D$ is the actual number of deaths in age $i$ across all individuals and $\hat{D}$ is the estimated number of deaths according to the model.

Have I interpreted this correctly? Apologies if I've left any obvious required information out.

$\endgroup$
1
$\begingroup$

The equation you wrote is almost what you are looking for,

$Deviance =\sum{r_i}$

I don't have my books, so I can't peer-review your equation, but it seems ok. I can't tell if it is appropriate to use a deviance statistic with parametric survival regression. But I can help you about the deviance assuming a Poisson distribution.

Generally speaking, the deviance is

$Deviance = 2(l(y) - l(\hat{y}))$

with $l(.)$ the log-likelihood. $l(y)$ is for a full (saturated/overfit) model, so you have $l(y)=l(D)$. Also $l(\hat{y}) = l(\hat{D})$. What you wrote is a simplification of my deviance equation.

You could try this

logLik.Poisson <- function(Deaths, fit, weight){
  sum(
    ( -lgamma(Deaths + 1) + Deaths*log(fit) - fit ) * weight
  )
}

2*(logLik.Poisson(actual.deaths, actual.deaths, 1) - logLik.Poisson(actual.deaths, fit.deaths, 1))

with actual.deaths and fit.deaths vectors of equal length. If you use all of your data, weight=1 . R will repeat the 1 to equal the vectors length. Otherwise, if you exclude some data, have weight to be a vector of 1 and 0, the 0 at the index of the excluded data.

One last thing, the code use $ln(\Gamma(x+1))$ for the log-likelihood, since $\Gamma(x+1) = x!$. R provides this function so I used it.

Well, here is a (dummy) reproductive example

actual.deaths <- sort(rpois(100, 60))
actual.exposure <- rep(100,100)
age <-1:100
fit <- glm(actual.deaths~age, offset=log(actual.exposure), family=poisson)

logLik.Poisson <- function(Deaths, fit, weight){
  sum(
    ( -lgamma(Deaths + 1) + Deaths*log(fit) - fit ) * weight
  )
}

fit.deaths <- fit$fitted
fit$deviance
2*(logLik.Poisson(actual.deaths, actual.deaths, 1) - logLik.Poisson(actual.deaths, fit.deaths, 1))

The deviance calculated from glm (fit$deviance) is the same as we calculated

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