Obliczanie AIC „ręcznie” w R.

16

Próbowałem obliczyć AIC regresji liniowej w R, ale bez użycia AICfunkcji:

lm_mtcars <- lm(mpg ~ drat, mtcars)

nrow(mtcars)*(log((sum(lm_mtcars$residuals^2)/nrow(mtcars))))+(length(lm_mtcars$coefficients)*2)
[1] 97.98786

Jednak AICdaje inną wartość:

AIC(lm_mtcars)
[1] 190.7999

Czy ktoś mógłby mi powiedzieć, co robię źle?

luciano
źródło
5
(bez sprawdzania odpowiedzi): niekoniecznie robisz coś złego, ponieważ prawdopodobieństwo jest zdefiniowane tylko do stałej multiplikatywnej; dwie osoby mogą obliczyć prawdopodobieństwo dziennika i uzyskać różne liczby (ale różnice w prawdopodobieństwie dziennika są takie same).
Glen_b
1
Myślę, że odpowiedź Hong Oois jest związana z tym pytaniem. Formuła AICużywana przez tę funkcję to -2*as.numeric(logLik(lm_mtcars))+2*(length(lm_mtcars$coefficients)+1).
COOLSerdash
luciano: „+1” w tej formule @COOLSerdash wskazuje na wynik parametru wariacji. Zauważ też, że funkcja logLikmówi, że dla lmmodeli zawiera „wszystkie stałe” ... więc log(2*pi)gdzieś tam będzie
Glen_b
1
@Glen_b: Dlaczego powiedzmy, że prawdopodobieństwo jest zdefiniowane tylko do stałej multiplikatywnej? W końcu, porównując nie zagnieżdżone modele z różnych rodzin dystrybucji (np. Z AIC lub testem Coxa), musisz pamiętać o tej stałej.
Scortchi - Przywróć Monikę
@Scortchi definicja nie jest moja! Musisz wziąć to z RAFisher. Myślę, że tak było od samego początku (1921). Że wciąż jest tak zdefiniowane, przynajmniej w przypadku ciągłym, patrz tutaj na przykład zdanie rozpoczynające się „Dokładniej”.
Glen_b

Odpowiedzi:

19

Zauważ, że pomoc dla funkcji logLikw R mówi, że dla lmmodeli zawiera ona „wszystkie stałe” ... więc gdzieś tam będzie log(2*pi), a także inny stały termin na wykładnik prawdopodobieństwa. Nie można również zapomnieć o tym, że jest parametrem.σ2

L(μ^,σ^)=(12πsn2)nexp(12i(ei2/sn2))

2logL=nlog(2π)+nlogsn2+i(ei2/sn2)

=n[log(2π)+logsn2+1]

AIC=2p2logL

ale zauważ, że dla modelu z 1 zmienną niezależną p = 3 (współczynnik x, stała i )σ2

Co oznacza, że ​​tak otrzymujesz ich odpowiedź:

nrow(mtcars)*(log(2*pi)+1+log((sum(lm_mtcars$residuals^2)/nrow(mtcars))))
       +((length(lm_mtcars$coefficients)+1)*2)
Glen_b - Przywróć Monikę
źródło
Dlaczego w obliczeniach dzielisz tylko przez n, a nie n - p ? s2nnp
Luke Thorburn,
1
2logL(θ^)+2pθθ^nσ^2σ2
Glen_b - Przywróć Monikę
2logL=nlog(2π)+nlogsn+i(ei2/sn2)2πsn2
11

AICFunkcja daje2)k-2)logL., gdzie L. jest prawdopodobieństwo i k is the number of estimated parameters (including the intercept, & the variance). You're using nlogSrn+2(k1), where Sr is the residual sum of squares, & n is the sample size. These formulæ differ by an additive constant; so long as you're using the same formula & looking at differences in AIC between different models where the constants cancel, it doesn't matter.

Scortchi - Reinstate Monica
źródło