¿Cómo se calcula la AIC en stepAIC

Question

Aquí es un modelo muy simple de 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) 

Si uso stepAIC a lm.D9, en la primera línea que dice AIC = -12,58

require(MASS) 
stepAIC(lm.D9) 

Si uso AIC directamente en lm.D9, da un valor diferente 46,17648

AIC(lm.D9) 

Mi pregunta es por qué los 2 valores AIC son diferentes. ¡Gracias!

Answer 1

AIC solo se define hasta como una constante arbitraria. Siempre que se use el mismo valor de la constante al comparar los AIC para diferentes modelos, no importa. Si mira ?extractAIC y ?AIC, encontrará las fórmulas utilizadas por ambos métodos.

Básicamente, utilice extractAIC o AIC, pero no ambos al mismo tiempo.

Answer 2

Esto me molestaba, así que decidí trabajar desde los primeros principios.

Vuelva a ajustar el modelo:

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) 

valores devueltos por métodos de acceso estándar:

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

reconstruir a partir de primeros principios:

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 

Esta es una pequeña poco fuera , no he encontrado el fallo.

Número de parámetros:

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


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

Aún fuera por más de pelusa numérica, pero sólo alrededor de una parte en un millón.

Ahora vamos a ver lo que está haciendo 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 

se puede trabajar la fórmula utilizada por encima de Sigma-hat = RSS/n - o ver Venables y Ripley MASAS para la derivación.

añadir términos que faltan: parámetro de varianza no se cuentan, además de constantes

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

normalización Este es exactamente el mismo (a 1e-14) como AIC1 anterior

Answer 3

Gracias @benbolker de respuesta detallada. Usted ha mencionado:

Ésta es una pequeña poco fuera, no han encontrado el problema técnico.

miré en él y encontraron que si se modifica esta línea:

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

a esto:

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

entonces el resultado sería exactamente el mismo.