Przelicz prawdopodobieństwo dziennika na podstawie prostego modelu R lm

Po prostu próbuję ponownie obliczyć za pomocą dnorm () prawdopodobieństwo dziennika podane przez funkcję logLik z modelu lm (w języku R).

Działa (prawie idealnie) dla dużej liczby danych (np. 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

ale w przypadku małych zestawów danych istnieją wyraźne różnice:

> 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

Ze względu na mały efekt zestawu danych pomyślałem, że może to wynikać z różnic w szacunkach wariancji rezydualnej między lm i glm, ale użycie lm daje ten sam wynik co 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

Gdzie się mylę?

— Gilles
źródło

lm()

\sqrt{\hat{σ}}

$\sqrt{\hat\sigma}$

\hat{σ}

$\hat\sigma$

Dzięki Stéphane za korektę, ale nadal nie działa

— Gilles,

spróbuj spojrzeć na kod źródłowy:stats:::logLik.glm

— zakłada się normalny

Zrobiłem to, ale funkcja ta po prostu odwraca szczelinę aic od obiektu glm, aby znaleźć prawdopodobieństwo dziennika. I nie widzę nic o aic w funkcji glm ...

— Gilles,

Podejrzewam, że ma to coś wspólnego z LogLik i AIC (które są powiązane razem na biodrze) przy założeniu, że szacowane są trzy parametry (nachylenie, przecięcie i dyspersja / błąd resztkowy standardowy), podczas gdy dyspersja / błąd resztowy standardowy jest obliczany przy założeniu oszacowano dwa parametry (nachylenie i przecięcie).

— Tom

logLik() $\beta_j$ $X{\boldsymbol \beta}$ $\sigma$ $\sqrt{\frac{\sum \hat\epsilon_i^2}{n}}$ $\hat\sigma = \sqrt{\frac{\sum \hat\epsilon_i^2}{n-2}}$ $\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
źródło

Nawiasem mówiąc, musisz podobnie uważać z opcją REML / ML dla modeli lme / lmer.

— Stéphane Laurent,

\hat{σ}

$\hat\sigma$

@PatrickCoulombe No: przechwytywanie + zbocze

— Stéphane Laurent,

Ok, teraz całkowicie jasne. Wielkie dzięki ! Ale co masz na myśli przez REML / ML (chyba coś wspólnego z moim ostatnim postem na GuR)? Proszę wyjaśnić (może tam). Chcę się uczyć !

— Gilles,

Oszacowania REML składników wariancji w modelach mieszanych są jak oszacowania ML „skorygowane o odchylenie”. Nie widziałem jeszcze twojego postu na GuR :)

— Stéphane Laurent,