Utilizaré un ejemplo para que puedas reproducir los resultados
# mortality
mort = ts(scan("http://www.stat.pitt.edu/stoffer/tsa2/data/cmort.dat"),start=1970, frequency=52)
# temperature
temp = ts(scan("http://www.stat.pitt.edu/stoffer/tsa2/data/temp.dat"), start=1970, frequency=52)
#pollutant particulates
part = ts(scan("http://www.stat.pitt.edu/stoffer/tsa2/data/part.dat"), start=1970, frequency=52)
temp = temp-mean(temp)
temp2 = temp^2
trend = time(mort)
Ahora, ajuste un modelo para los datos de mortalidad
fit = lm(mort ~ trend + temp + temp2 + part, na.action=NULL)
Lo que quiero ahora es reproducir el resultado del comando AIC
AIC(fit)
[1] 3332.282
Según el archivo de ayuda de R para AIC, AIC = -2 * log.likelihood + 2 * npar. Si estoy en lo cierto, creo que log.likelihood se da usando la siguiente fórmula:
n = length(mort)
RSS = anova(fit)[length(anova(fit)[,2]),2] # there must be better ways to get this, anyway
(log.likelihood <- -n/2*(log(2*pi)+log(RSS/n)+1))
[1] -1660.135
Esto es aproximadamente igual a
logLik(fit)
'log Lik.' -1660.141 (df=6)
Por lo que veo, el número de parámetros del modelo es de 5 (¿cómo puedo obtener este número programáticamente?). Así que el AIC debe ser dado por:
-2 * log.likelihood + 2 * 5
[1] 3330.271
Ooops, parece que debería haber utilizado 6 en lugar de 5 como número de parámetros. ¿Qué hay de malo en esos cálculos?