Recalcular la probabilidad de registro a partir de un modelo R lm simple

Question

Simplemente estoy tratando de volver a calcular con dnorm() el logaritmo de la probabilidad proporcionados por el logLik función de un modelo lm (R).

Funciona casi a la perfección) por la gran cantidad de datos (por ejemplo, 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

pero para pequeños conjuntos de datos hay claras diferencias :

> 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

Porque de pequeño conjunto de datos efecto pensé que podría ser debido a las diferencias en las estimaciones de la varianza residual entre la película y el glm pero el uso de lm proporciona el mismo resultado que el 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

Donde estoy equivocado ?

Answer 1

La función logLik() proporciona la evaluación de la probabilidad logarítmica sustituyendo las estimaciones de ML de los parámetros por los valores de los parámetros desconocidos. Ahora, las estimaciones de máxima verosimilitud de los parámetros de regresión (el$\beta_j$ 's en$X{\boldsymbol \beta}$) coinciden con las estimaciones de mínimos cuadrados, pero la estimación de LD de$\sigma$ es$\frac{\sum \hat\epsilon_i^2}{n}$ , mientras que estás usando$\hat\sigma = \sqrt{\frac{\sum \hat\epsilon_i^2}{n-2}}$, esa es la raíz cuadrada de la estimación imparcial de$\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)