12

I'm simply trying to recalculate with dnorm() the log-likelihood provided by the logLik function from a lm model (in R).

It works (almost perfectly) for high number of data (eg n=1000) : 

> n <- 1000
> x <- 1:n
> set.seed(1)
> y <- 10 + 2*x + rnorm(n, 0, 2)
> mod <- glm(y ~ x, family = gaussian)
> logLik(mod)
'log Lik.' -2145.562 (df=3)
> sigma <- sqrt(summary(mod)$dispersion)
> sum(log(dnorm(x = y, mean = predict(mod), sd = sigma)))
[1] -2145.563
> sum(log(dnorm(x = resid(mod), mean = 0, sd = sigma)))
[1] -2145.563

but for small datasets there are clear differences : 

> n <- 5
> x <- 1:n
> set.seed(1)
> y <- 10 + 2*x + rnorm(n, 0, 2)
> 
> mod <- glm(y ~ x, family = gaussian)
> logLik(mod)
'log Lik.' -8.915768 (df=3)
> sigma <- sqrt(summary(mod)$dispersion)
> sum(log(dnorm(x = y, mean = predict(mod), sd = sigma)))
[1] -9.192832
> sum(log(dnorm(x = resid(mod), mean = 0, sd = sigma)))
[1] -9.192832

Because of small dataset effect I thought it could be due to the differences in residual variance estimates between lm and glm but using lm provides the same result as glm : 

> modlm <- lm(y ~ x)
> logLik(modlm)
'log Lik.' -8.915768 (df=3)
> 
> sigma <- summary(modlm)$sigma
> sum(log(dnorm(x = y, mean = predict(modlm), sd = sigma)))
[1] -9.192832
> sum(log(dnorm(x = resid(modlm), mean = 0, sd = sigma)))
[1] -9.192832

Where am I wrong ?

Gilles
  • 405
  • 1
  • 3
  • 9
  • 2
    With `lm()`, you are using $\sqrt{\hat\sigma}$ instead of $\hat\sigma$. – Stéphane Laurent Oct 18 '13 at 22:40
  • Thanks Stéphane for the correction but it still doesn't seem to work – Gilles Oct 19 '13 at 09:31
  • try looking at the source code: `stats:::logLik.glm` – assumednormal Oct 19 '13 at 09:45
  • I did this but this function just reverse the aic slot from the glm object to find back the log-likelihood. And I don't see anything about aic in the glm function... – Gilles Oct 19 '13 at 10:20
  • I suspect this has something to do with LogLik and AIC (which are tied together at the hip) assuming that three parameters are being estimated (the slope, intercept, and dispersion/residual standard error) whereas the dispersion/residual standard error is calculated assuming two parameters are estimated (slope and intercept). – Tom Oct 19 '13 at 12:04

1 Answers1

15

The logLik() function provides the evaluation of the log-likelihood by substituting the ML estimates of the parameters for the values of the unknown parameters. Now, the maximum likelihood estimates of the regression parameters (the $\beta_j$'s in $X{\boldsymbol \beta}$) coincide with the least-squares estimates, but the ML estimate of $\sigma$ is $\sqrt{\frac{\sum \hat\epsilon_i^2}{n}}$, whereas you are using $\hat\sigma = \sqrt{\frac{\sum \hat\epsilon_i^2}{n-2}}$, that is the square root of the unbiased estimate of $\sigma^2$.

>  n <- 5
>  x <- 1:n
>  set.seed(1)
>  y <- 10 + 2*x + rnorm(n, 0, 2)
>  modlm <- lm(y ~ x)
>  sigma <- summary(modlm)$sigma
> 
>  # value of the likelihood with the "classical" sigma hat
>  sum(log(dnorm(x = y, mean = predict(modlm), sd = sigma)))
[1] -9.192832
> 
>  # value of the likelihood with the ML sigma hat
>  sigma.ML <- sigma*sqrt((n-dim(model.matrix(modlm))[2])/n) 
>  sum(log(dnorm(x = y, mean = predict(modlm), sd = sigma.ML)))
[1] -8.915768
>  logLik(modlm)
'log Lik.' -8.915768 (df=3)
Stéphane Laurent
  • 17,425
  • 5
  • 59
  • 101