Solved – How does extractAIC() calculate the AIC of an lm object

by

I'm having trouble to replicate how extractAIC() is calculating the AIC for a linear model (lm). There are a lot of different formulas for the AIC out there, I've tried them all but none of them yielded the exact same value as extractAIC(). Here is what I've tried (extractAIC(finalModel) is 24823.52):

finalModel = step(lm(formula, data = myData), scope = scopeFormula ,direction="both", k=2)  n = df.residual(finalModel) edf = length(finalModel$coefficients) DVVariance = var(companyData$MarketValueTotalFiscal) RSS = sum(resid(finalModel) ^ 2)  #AIC formula using maximum likelyhood estimation (which is not used by extractAIC() with an lm object) -2*logLik(finalModel)+2*edf AIC(finalModel)  #AIC formula using OLS, which extractAIC() should theoretically be using with lm according to help n*log(RSS/n)-n+n*log(2*pi) +2*edf  #AIC formula Venables with known variance (of the DV of the poulation I guess) RSS/DVVariance + 2*edf  #AIC formula according to Venables "Modern Applied Statistics" 2002 p. 174 with unknown variance n * log(RSS/n) + 2*edf # +constant, but he doesn't explain what constant is  #AIC forumla accordint to Venables with the constant term according to Wikipedia about AIC n * log(RSS/n) + 2*edf + (n/2)*log(2*pi)  #Comparison to the extractAIC() function that is used by the step() function (both are in the MASS package) extractAIC(finalModel)

Here is the output:

> #AIC formula using maximum likelyhood estimation (which is not used by extractAIC() with an lm object) > -2*logLik(finalModel)+2*edf 'log Lik.' 29738.72 (df=50) > AIC(finalModel) [1] 29740.72 >  > #AIC formula using OLS, which extractAIC() should theoretically be using with lm according to help > n*log(RSS/n)-n+n*log(2*pi) +2*edf [1] 25582.46 >  > #AIC formula Venables with known variance (of the DV of the poulation I think) > RSS/DVVariance + 2*edf [1] 274.651 >  > #AIC formula according to Venables "Modern Applied Statistics" 2002 p. 174 with unknown variance > n * log(RSS/n) + 2*edf # +constant, but he doesn't explain what constant is [1] 24172.31 >  > #AIC forumla accordint to Venables with the constant term according to Wikipedia about AIC > n * log(RSS/n) + 2*edf + (n/2)*log(2*pi) [1] 25718.88 >  > #Comparison to the extractAIC() function that is used by the step() function (both are in the MASS package) > extractAIC(finalModel) [1]    49.00 24823.52
Contents hide
1 Best Answer
2 Similar Posts:

Best Answer

You can always check the source code since R is open-source:

> stats:::extractAIC.lm function (fit, scale = 0, k = 2, ...)  {     n <- length(fit$residuals)     edf <- n - fit$df.residual     RSS <- deviance.lm(fit)     dev <- if (scale > 0)          RSS/scale - n     else n * log(RSS/n)     c(edf, dev + k * edf) } <bytecode: 0x000000001f88adb8> <environment: namespace:stats>

This is also described in the documentation ?extractAIC:

The criterion used is

AIC = - 2*log L + k * edf,

where L is the likelihood and edf the equivalent degrees of freedom (i.e., the number of free parameters for usual parametric models) of fit.

For linear models with unknown scale (i.e., for lm and aov), -2 log L is computed from the deviance and uses a different additive constant to logLik and hence AIC. If RSS denotes the (weighted) residual sum of squares then extractAIC uses for -2 log L the formulae RSS/s - n (corresponding to Mallows' Cp) in the case of known scale s and n log (RSS/n) for unknown scale. AIC only handles unknown scale and uses the formula n*log(RSS/n) + n + n*log 2pi - sum(log w) where w are the weights. Further AIC counts the scale estimation as a parameter in the edf and extractAIC does not.

Similar Posts:

Rate this post

Leave a Comment