Hitung ulang log-kemungkinan dari model Rm sederhana

10

Saya hanya mencoba untuk menghitung ulang dengan dnorm () log-kemungkinan yang disediakan oleh fungsi logLik dari model lm (dalam R).

Ini berfungsi (hampir sempurna) untuk jumlah data yang tinggi (mis. 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

tetapi untuk dataset kecil ada perbedaan yang jelas:

> 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

Karena efek dataset kecil saya pikir itu bisa jadi karena perbedaan estimasi varians residual antara lm dan glm tetapi menggunakan lm memberikan hasil yang sama seperti 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

Dimana saya salah

Gilles
sumber
2
Dengan lm(), Anda menggunakan alih-alih . σ^σ^
Stéphane Laurent
Terima kasih Stéphane untuk koreksi tetapi sepertinya masih tidak berhasil
Gilles
coba lihat kode sumber:stats:::logLik.glm
diasumsikan normal
Saya melakukan ini tetapi fungsi ini hanya membalikkan slot aic dari objek glm untuk menemukan kembali kemungkinan log. Dan saya tidak melihat apa pun tentang aic dalam fungsi glm ...
Gilles
Saya menduga ini ada hubungannya dengan LogLik dan AIC (yang diikat bersama di pinggul) dengan asumsi bahwa tiga parameter diperkirakan (kemiringan, intersep, dan dispersi / kesalahan standar residu) sedangkan kesalahan standar dispersi / residu dihitung dengan asumsi dua parameter diperkirakan (kemiringan dan penyadapan).
Tom

Jawaban:

12

The logLik()berfungsi memberikan evaluasi log-kemungkinan oleh menggantikan perkiraan ML parameter untuk nilai-nilai parameter yang tidak diketahui. Sekarang, perkiraan kemungkinan maksimum dari parameter regresi ( dalam ) bertepatan dengan estimasi kuadrat-terkecil, tetapi estimasi ML dari adalah , sedangkan Anda menggunakan , itu adalah akar kuadrat dari peta bias estimasi .βjXβσϵ^i2nσ^=ϵ^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
sumber
Ngomong-ngomong Anda juga harus berhati-hati dengan opsi REML / ML untuk model lme / lmer.
Stéphane Laurent
(+1) Apakah n-1 atau memang n-2 dalam penyebut ? σ^
Patrick Coulombe
@PatrickCoulombe No: intercept + slope
Stéphane Laurent
Ok, sangat jelas sekarang. Terima kasih banyak ! Tapi apa yang Anda maksud dengan REML / ML (ada hubungannya dengan posting terakhir saya di GuR kurasa)? Tolong jelaskan (mungkin). Saya ingin belajar !
Gilles
Perkiraan REML dari komponen varians dalam model campuran seperti perkiraan ML "dikoreksi untuk bias". Saya belum melihat pos Anda di GuR :):
Stéphane Laurent