How AIC is calculated on stepAIC

Question

Here is a very simple lm model from? lm

ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14) trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69) group <- gl(2,10,20, labels=c("Ctl","Trt")) weight <- c(ctl, trt) lm.D9 <- lm(weight ~ group) 

If I use stepAIC for lm.D9, the very first line says AIC = -12.58

 require(MASS) stepAIC(lm.D9) 

If I use AIC directly on lm.D9, it gives a different value of 46.17648

 AIC(lm.D9) 

My question is why the 2 AIC values ​​are different. Thanks!

Answer 1

AIC is defined only with an arbitrary constant. As long as the same constant value is used when comparing AIC for different models, it does not matter. If you look at ?extractAIC and ?AIC , you will find the formulas used by both methods.

Basically, use extractAIC or AIC , but not both at the same time.

Answer 2

It annoyed me, so I decided to deal with the first principles.

Repeat model installation:

 d <- data.frame(weight= c(ctl=c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14), trt=c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)), group=gl(2,10,20, labels=c("Ctl","Trt"))) lm.D9 <- lm(weight ~ group, d) 

Values ​​returned by standard accessories:

 (AIC1 <- AIC(lm.D9)) > 46.17468 (LL1 <- logLik(lm.D9)) > -20.08824 (df=3) 

Restore from first principles:

 n <- nrow(d) ss0 <- summary(lm.D9)$sigma ss <- ss0*(n-1)/n (LL2 <- sum(dnorm(d$weight,fitted(lm.D9), sd=ss,log=TRUE))) > -20.08828 

This is a tiny bit, did not find a glitch.

Number of parameters:

 npar <- length(coef(lm.D9))+1 (AIC2 <- -2*LL2+2*npar) > 46.1756 

Disabled by more than numerical down, but only about one part in a million.

Now let's see what stepAIC does:

 MASS::stepAIC(lm.D9) ## start: AIC = -12.58 extractAIC(lm.D9) ## same value (see MASS::stepAIC for details) stats:::extractAIC.lm ## examine the code RSS1 <- deviance(lm.D9) ## UNSCALED sum of squares RSS2 <- sum((d$weight-fitted(lm.D9))^2) ## ditto, from first principles AIC3 <- n*log(RSS1/n)+2*2 ## formula used within extractAIC 

You can work out the formula used above from sigma-hat = RSS / n - or see Venables and Ripley MASS for output.

Add missing terms: uncertain variance parameter, plus normalization constant

 (AIC3 + 2 - 2*(-n/2*(log(2*pi)+1))) 

This is the same (up to 1e-14) as AIC1 above

Answer 3

Thanks @benbolker for the detailed answer. You mentioned:

This is a tiny bit, did not find a glitch.

I looked through it and found that if you change this line:

 ss <- ss0*(n-1)/n 

:

 ss <- sqrt( (ss0)^2 * (n - length(coef(lm.D9))) / n ) 

then the result will be exactly the same.