Recalculer la log-vraisemblance à partir d'un simple modèle R lm

10

J'essaie simplement de recalculer avec dnorm () la log-vraisemblance fournie par la fonction logLik à partir d'un modèle lm (dans R).

Cela fonctionne (presque parfaitement) pour un grand nombre de données (par exemple 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

mais pour les petits ensembles de données, il existe des différences claires:

> 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

En raison de l'effet de petit ensemble de données, je pensais que cela pourrait être dû aux différences d'estimation de la variance résiduelle entre lm et glm, mais l'utilisation de lm fournit le même résultat que 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

Où ai-je tort?

Gilles
la source
2
Avec lm(), vous utilisez au lieu de . σ^σ^
Stéphane Laurent
Merci Stéphane pour la correction mais ça ne semble toujours pas fonctionner
Gilles
essayez de regarder le code source:stats:::logLik.glm
supposé normal
Je l'ai fait, mais cette fonction inverse simplement la fente aic de l'objet glm pour retrouver la probabilité de journalisation. Et je ne vois rien d'aic dans la fonction glm ...
Gilles
Je soupçonne que cela a quelque chose à voir avec LogLik et AIC (qui sont liés ensemble à la hanche) en supposant que trois paramètres sont estimés (la pente, l'ordonnée à l'origine et la dispersion / l'erreur standard résiduelle) tandis que la dispersion / l'erreur standard résiduelle est calculée en supposant deux paramètres sont estimés (pente et interception).
Tom

Réponses:

12

logLik()βjXβσϵ^je2nσ^=ϵ^je2n-2σ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
la source
Soit dit en passant, vous devez également être prudent avec l'option REML / ML pour les modèles lme / lmer.
Stéphane Laurent
(1) Est - il n-1 ou est - ce en effet n-2 dans le dénominateur σσ^ ?
Patrick Coulombe
@PatrickCoulombe No: interception + pente
Stéphane Laurent
Ok, parfaitement clair maintenant. Merci beaucoup ! Mais que voulez-vous dire avec REML / ML (quelque chose à voir avec mon dernier post sur GuR je suppose)? Veuillez expliquer (peut-être). Je veux apprendre !
Gilles
Les estimations REML des composantes de la variance dans un modèle mixte sont similaires aux estimations ML «corrigées du biais». Je n'ai pas encore vu votre post sur GuR :)
Stéphane Laurent