[R] stepAIC() that can use new extractAIC() function implementing AICc

Marc Girondot marc_grt at yahoo.fr
Thu Jun 8 06:54:26 CEST 2017


I would like test AICc as a criteria for model selection for a glm using 
stepAIC() from MASS package.

Based on various information available in WEB, stepAIC() use 
extractAIC() to get the criteria used for model selection.

I have created a new extractAIC() function (and extractAIC.glm() and 
extractAIC.lm() ones) that use a new parameter criteria that can be AIC, 
BIC or AICc.

It works as expected using extractAIC() but when I run stepAIC(), the 
first AIC shown in the result is correct, but after it still shows the 
original AIC:

For example (the full code is below) the "Start:  AIC=70.06" is indeed 
the AICc but after, "<none>      47.548 67.874" is the AIC.

 > stepAIC(G1, criteria="AICc")
Start:  AIC=70.06
x ~ A + B

        Df Deviance    AIC
- A     1   48.506 66.173
<none>      47.548 67.874
- B     1   57.350 68.685

Thanks if you can help me that stepAIC() use always the new extractAIC() 
function.

Marc




library("MASS")
set.seed(1)

df <- data.frame(x=rnorm(15, 15, 2))
for (i in 1:10) {
   df <- cbind(df, matrix(data = rnorm(15, 15, 2), ncol=1, 
dimnames=list(NULL, LETTERS[i])))
}

G1 <- glm(formula = x ~ A + B,
           data=df, family = gaussian(link = "identity"))

extractAIC(G1)
stepAIC(G1)

extractAIC.glm <- function(fit, scale, k = 2, criteria = c("AIC", 
"AICc", "BIC"), ...) {
   criteria <- match.arg(arg=criteria, choice=c("AIC", "AICc", "BIC"))
   AIC <- fit$aic
   edf <- length(fit$coefficients)
   n <- nobs(fit, use.fallback = TRUE)
   if (criteria == "AICc") return(c(edf, AIC + (2*edf*(edf+1))/(n - edf 
-1)))
   if (criteria == "AIC")  return(c(edf, AIC-2*edf + k*edf))
   if (criteria == "BIC")  return(c(edf, AIC -2*edf + log(n)*edf))
}

extractAIC <- extractAIC.lm <- extractAIC.glm

extractAIC(G1, criteria="AIC")
extractAIC(G1, k=log(15))
extractAIC(G1, criteria="BIC")
stats:::extractAIC.glm(G1, k=log(15))
extractAIC(G1, criteria="AICc")

stepAIC(G1)
stepAIC(G1, criteria="AICc")



More information about the R-help mailing list