I'll use an example so that you can reproduce the results
# 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)
Now, fit a model for mortality data
fit = lm(mort ~ trend + temp + temp2 + part, na.action=NULL)
What I want now is to reproduce the result of the AIC command
AIC(fit)
[1] 3332.282
According to R's help file for AIC, AIC = -2 * log.likelihood + 2 * npar. If I'm correct I think that log.likelihood is given using the following formula:
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
This is approximately equal to
logLik(fit)
'log Lik.' -1660.141 (df=6)
As far as I can tell, the number of parameters in the model are 5 (how can I get this number programmatically ??). So AIC should be given by:
-2 * log.likelihood + 2 * 5
[1] 3330.271
Ooops, it seems like I should have used 6 instead of 5 as the number of parameters. What is wrong with those calculations?
