Recalcular a probabilidade de log a partir de um modelo Rmm simples

10

Estou simplesmente tentando recalcular com dnorm () a probabilidade de log fornecida pela função logLik de um modelo lm (em R).

Funciona (quase perfeitamente) para um grande número de dados (por exemplo, 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

mas para pequenos conjuntos de dados há diferenças claras:

> 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

Devido ao pequeno efeito do conjunto de dados, pensei que poderia ser devido às diferenças nas estimativas de variação residual entre lm e glm, mas o uso de lm fornece o mesmo resultado que o 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

Onde eu estou errado?

Gilles
fonte
2
Com lm(), você está usando vez de . σσ^σ^
Stéphane Laurent
Graças Stéphane para a correção, mas ainda não parece trabalho
Gilles
tente olhar para o código fonte:stats:::logLik.glm
assumenormal
Eu fiz isso, mas essa função apenas inverte o slot aic do objeto glm para encontrar novamente a probabilidade de log. E eu não vejo nada sobre aic na função glm ...
Gilles
Suspeito que isso tenha algo a ver com o LogLik e o AIC (que estão amarrados no quadril), assumindo que três parâmetros estão sendo estimados (inclinação, interceptação e dispersão / erro padrão residual) enquanto a dispersão / erro padrão residual é calculada assumindo dois parâmetros são estimados (inclinação e interceptação).
Tom

Respostas:

12

A logLik()função fornece a avaliação da probabilidade logarítmica, substituindo as estimativas de ML dos parâmetros pelos valores dos parâmetros desconhecidos. Agora, as estimativas de probabilidade máxima dos parâmetros de regressão (os em ) coincidem com as estimativas de mínimos quadrados, mas a estimativa ML de é , enquanto você está usando , que é a raiz quadrada dos estimativa de . X ββjXβσ σ=ϵ^i2n σ2σ^=ϵ^i2n2σ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
fonte
A propósito, você deve ter o mesmo cuidado com a opção REML / ML para modelos lme / lmer.
Stéphane Laurent
(+1) É n-1 ou é realmente n-2 no denominador de ? σ^
precisa
@PatrickCoulombe No: interceptação + inclinação
Stéphane Laurent
Ok, perfeitamente claro agora. Muito obrigado ! Mas o que você quer dizer com REML / ML (algo a ver com meu último post no GuR, eu acho)? Por favor, explique (talvez haja). Eu quero aprender !
Gilles
As estimativas REML dos componentes de variância em modelos mistos são como as estimativas ML corrigidas para viés. Eu não vi o seu post sobre Gur ainda :)
Stéphane Laurent