Как рассчитывается AIC в stepAIC

Question

Вот очень простая модель lm от? 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)

Если я использую stepAIC для lm.D9, в самой первой строке будет указано AIC = -12.58

require(MASS)
stepAIC(lm.D9)

Если я использую AIC непосредственно на lm.D9, это дает другое значение 46.17648

AIC(lm.D9)

У меня вопрос, почему значения 2 AIC разные. Спасибо!

Answer 1

Это меня раздражало, поэтому я решил выработать это из первых принципов.

Переустановите модель:

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)

Значения, возвращаемые стандартными средствами доступа:

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

Реконструкция из первых принципов:

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

Это крошечный бит, не нашли сбой.

Количество параметров:

npar <- length(coef(lm.D9))+1 


(AIC2 <- -2*LL2+2*npar)
> 46.1756

Еще больше, чем числовой пух, но только одна часть на миллион.

Теперь посмотрим, что делает stepAIC:

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

Вы можете найти формулу, использованную выше, из sigma-hat = RSS / n - или посмотреть Venables и Ripley MASS для получения.

Добавить пропущенные термины: несчетный параметр дисперсии плюс нормализационная постоянная

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

Это точно так же (до 1e-14), как AIC1 выше

Answer 2

AIC определяется только с точностью до произвольной константы. Если при сравнении AIC для разных моделей используется одно и то же значение константы, это не имеет значения. Если вы посмотрите на ?extractAIC и ?AIC, вы найдете формулы, используемые обоими методами.

В основном, используйте extractAIC или AIC, но не оба одновременно.

Answer 3

Спасибо @benbolker за подробный ответ. Вы упомянули:

Это крошечный бит, не нашли сбой.

Я посмотрел на это и обнаружил, что если вы измените эту строку:

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

к этому:

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

тогда результат будет точно таким же.