How can I do model selection by AIC with a Gamma GLM in R?

Question

As the documentation for glm() explains, the aic component of the value returned by glm() is not a valid AIC:

For gaussian, Gamma and inverse gaussian families the dispersion is estimated from the residual deviance, and the number of parameters is the number of coefficients plus one. For a gaussian family the MLE of the dispersion is used so this is a valid value of AIC, but for Gamma and inverse gaussian families it is not.

Thus a valid AIC needs to obtained in some other way.

e3bo · Answer 1 · 2012-11-16T00:00:13.133

If you want to use the step() or MASS::stepAIC() model selection functions, you could first ensure that the AIC is calculated properly by doing something like this:

GammaAIC <- function(fit){
  disp <- MASS::gamma.dispersion(fit)
  mu <- fit$fitted.values
  p <- fit$rank
  y <- fit$y
  -2 * sum(dgamma(y, 1/disp, scale = mu * disp, log = TRUE)) + 2 * p
}
GammaAICc <- function(fit){
  val <- logLik(fit)
  p <- attributes(val)$df
  n <- attributes(val)$nobs
  GammaAIC(fit) + 2 * p * (p + 1) / (n - p - 1)      
}

my_extractAIC <- function(fit, scale=0, k=2, ...){
  n <- length(fit$residuals)
  edf <- n - fit$df.residual  
  if (fit$family$family == "Gamma"){
    aic <- GammaAIC(fit)
  } else {
    aic <- fit$aic
  }
  c(edf, aic + (k - 2) * edf)
}
assignInNamespace("extractAIC.glm", my_extractAIC, ns="stats")

If you use the glmulti package, you can simply specify the use of the above GammaAIC() or GammaAICc() functions with the crit parameter of glmulti().

How can I do model selection by AIC with a Gamma GLM in R?

1 Answers1