[R] logLik.lm()

[Achim Zeileis]

On Wednesday 25 June 2003 20:23, Edward Dick wrote:

> Hello,
>
> I'm trying to use AIC to choose between 2 models with
> positive, continuous response variables and different error
> distributions (specifically a Gamma GLM with log link and a
> normal linear model for log(y)). I understand that in some
> cases it may not be possible (or necessary) to discriminate
> between these two distributions. However, for the normal
> linear model I noticed a discrepancy between the output of
> the AIC() function and my calculations done "by hand."
> This is due to the output from the function logLik.lm(),
> which does not match my results using the dnorm() function
> (see simple regression example below).
>
> x <- c(43.22,41.11,76.97,77.67,124.77,110.71,144.46,188.05,171.18,
>       
> 204.92,221.09,178.21,224.61,286.47,249.92,313.19,332.17,374.35) y <-
> c(5.18,12.47,15.65,23.42,27.07,34.84,31.03,30.87,40.07,57.36,
> 47.68,43.40,51.81,55.77,62.59,66.56,74.65,73.54)
> test.lm <- lm(y~x)
> y.hat <- fitted(test.lm)
> sigma <- summary(test.lm)$sigma
> logLik(test.lm)
> # `log Lik.' -57.20699 (df=3)
> sum(dnorm(y, y.hat, sigma, log=T))
> # [1] -57.26704
>
> The difference in this simple example is slight, but
> it is magnified when using my data.

That is because you are not using the ML estimate of the variance.

sigmaML <- sqrt(mean(residuals(test.lm)^2))
sum(dnorm(y, y.hat, sigmaML, log=T))
# [1] -57.20699

Best,
Z

> My understanding is
> that it is necessary to use the complete likelihood functions
> for both the Gamma model and the 'lognormal' model (no constants
> removed, etc.) in order to accurately compare using AIC.
>
> Can someone point out my error, or explain this discrepancy?
>
> Thanks in advance!
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help