Seberapa tepercaya interval kepercayaan untuk objek lmer melalui paket efek?

Effectspaket menyediakan cara yang sangat cepat dan mudah untuk memplot hasil model efek campuran linier yang diperoleh melalui lme4paket . The effectinterval fungsi menghitung kepercayaan (CI) sangat cepat, tapi bagaimana dapat dipercaya adalah interval keyakinan ini?

Sebagai contoh:

library(lme4)
library(effects)
library(ggplot)

data(Pastes)

fm1  <- lmer(strength ~ batch + (1 | cask), Pastes)
effs <- as.data.frame(effect(c("batch"), fm1))
ggplot(effs, aes(x = batch, y = fit, ymin = lower, ymax = upper)) + 
  geom_rect(xmax = Inf, xmin = -Inf, ymin = effs[effs$batch == "A", "lower"],
        ymax = effs[effs$batch == "A", "upper"], alpha = 0.5, fill = "grey") +
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

Menurut CI yang dihitung menggunakan effectspaket, batch "E" tidak tumpang tindih dengan batch "A".

Jika saya mencoba confint.merModfungsi menggunakan yang sama dan metode default:

a <- fixef(fm1)
b <- confint(fm1)
# Computing profile confidence intervals ...
# There were 26 warnings (use warnings() to see them)

b <- data.frame(b)
b <- b[-1:-2,]

b1 <- b[[1]]
b2 <- b[[2]]

dt <- data.frame(fit   = c(a[1],  a[1] + a[2:length(a)]), 
                 lower = c(b1[1],  b1[1] + b1[2:length(b1)]), 
                 upper = c(b2[1],  b2[1] + b2[2:length(b2)]) )
dt$batch <- LETTERS[1:nrow(dt)]

ggplot(dt, aes(x = batch, y = fit, ymin = lower, ymax = upper)) +
  geom_rect(xmax = Inf, xmin = -Inf, ymin = dt[dt$batch == "A", "lower"], 
        ymax = dt[dt$batch == "A", "upper"], alpha = 0.5, fill = "grey") + 
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

Saya melihat bahwa semua CI tumpang tindih. Saya juga mendapatkan peringatan yang menunjukkan bahwa fungsi gagal menghitung CI yang dapat dipercaya. Contoh ini, dan set data aktual saya, membuat saya curiga bahwa effectspaket mengambil pintasan dalam perhitungan CI yang mungkin tidak sepenuhnya disetujui oleh ahli statistik. Seberapa tepercaya CI yang dikembalikan dengan effectfungsi dari effectspaket untuk lmerobjek?

Apa yang saya coba: Melihat ke kode sumber, saya perhatikan bahwa effectfungsi bergantung pada Effect.merModfungsi, yang pada gilirannya mengarahkan ke Effect.merfungsi, yang terlihat seperti ini:

effects:::Effect.mer
function (focal.predictors, mod, ...) 
{
    result <- Effect(focal.predictors, mer.to.glm(mod), ...)
    result$formula <- as.formula(formula(mod))
    result
}
<environment: namespace:effects>

mer.to.glmfungsi tampaknya menghitung Variance-Covariate Matrix dari lmerobjek:

effects:::mer.to.glm

function (mod) 
{
...
mod2$vcov <- as.matrix(vcov(mod))
...
mod2
}

Ini, pada gilirannya, mungkin digunakan dalam Effect.defaultfungsi untuk menghitung CI (saya mungkin salah paham bagian ini):

effects:::Effect.default
...
     z <- qnorm(1 - (1 - confidence.level)/2)
        V <- vcov.(mod)
        eff.vcov <- mod.matrix %*% V %*% t(mod.matrix)
        rownames(eff.vcov) <- colnames(eff.vcov) <- NULL
        var <- diag(eff.vcov)
        result$vcov <- eff.vcov
        result$se <- sqrt(var)
        result$lower <- effect - z * result$se
        result$upper <- effect + z * result$se
...

Saya tidak cukup tahu tentang LMM untuk menilai apakah ini merupakan pendekatan yang tepat, tetapi mempertimbangkan diskusi tentang perhitungan interval kepercayaan untuk LMM, pendekatan ini tampak sederhana dan mencurigakan.

Semua hasilnya pada dasarnya sama ( untuk contoh khusus ini ). Beberapa perbedaan teoretis adalah:

• seperti yang ditunjukkan @rvl, rekonstruksi CI Anda tanpa memperhitungkan kovarian antar parameter adalah salah (maaf)
• interval kepercayaan untuk parameter dapat didasarkan pada interval kepercayaan Wald (dengan asumsi permukaan log-likelihood kuadrat): lsmeans, effects, confint(.,method="Wald"); kecuali untuk lsmeans, metode ini mengabaikan efek ukuran terbatas ("derajat kebebasan"), tetapi dalam kasus ini hampir tidak ada bedanya ( df=40praktis tidak bisa dibedakan dari tak terbatas df)
• ... atau pada interval kepercayaan profil (metode standar; mengabaikan efek ukuran terbatas tetapi memungkinkan untuk permukaan non-kuadrat)
• ... atau pada bootstrap parametrik (standar emas - mengasumsikan modelnya benar [respons normal, efek acak terdistribusi normal, data independen secara kondisional, dll], tetapi sebaliknya membuat beberapa asumsi)

Saya pikir semua pendekatan ini masuk akal (beberapa lebih mendekati daripada yang lain), tetapi dalam hal ini hampir tidak ada bedanya yang mana yang Anda gunakan. Jika Anda khawatir, cobalah beberapa metode kontras pada data Anda, atau pada data simulasi yang menyerupai milik Anda, dan lihat apa yang terjadi ...

(PS: Saya tidak akan terlalu menekankan fakta bahwa interval kepercayaan Adan Etidak tumpang tindih. Anda harus melakukan prosedur perbandingan berpasangan yang tepat untuk membuat kesimpulan yang andal tentang perbedaan antara pasangan perkiraan tertentu ini . ..)

95% CI:

Kode perbandingan:

library(lme4)
fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
c0 <- confint(fm2,method="Wald")
c1 <- confint(fm2)
c2 <- confint(fm2,method="boot")
library(effects)
library(lsmeans)
c3 <- with(effect("batch",fm2),cbind(lower,upper))
c4 <- with(summary(lsmeans(fm2,spec="batch")),cbind(lower.CL,upper.CL))
tmpf <- function(method,val) {
    data.frame(method=method,
               v=LETTERS[1:10],
               setNames(as.data.frame(tail(val,10)),
                        c("lwr","upr")))
}
library(ggplot2); theme_set(theme_bw())
allCI <- rbind(tmpf("lme4_wald",c0),
      tmpf("lme4_prof",c1),
      tmpf("lme4_boot",c2),
      tmpf("effects",c3),
               tmpf("lsmeans",c4))
ggplot(allCI,aes(v,ymin=lwr,ymax=upr,colour=method))+
    geom_linerange(position=position_dodge(width=0.8))

ggsave("pastes_confint.png",width=10)

Sepertinya apa yang telah Anda lakukan dalam metode kedua adalah menghitung interval kepercayaan untuk koefisien regresi, kemudian mentransformasikannya untuk mendapatkan CI untuk prediksi. Ini mengabaikan kovariansi antara koefisien regresi.

Coba pas model tanpa intersep, sehingga batchefeknya benar-benar akan menjadi prediksi, dan confintakan mengembalikan interval yang Anda butuhkan.

Adendum 1

Saya melakukan persis apa yang saya sarankan di atas:

> fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
> confint(fm2)
Computing profile confidence intervals ...
           2.5 %    97.5 %
.sig01  0.000000  1.637468
.sigma  2.086385  3.007380
batchA 60.234772 64.298581
batchB 57.268105 61.331915
batchC 60.018105 64.081915
batchD 57.668105 61.731915
batchE 53.868105 57.931915
batchF 59.001439 63.065248
batchG 57.868105 61.931915
batchH 61.084772 65.148581
batchI 56.651439 60.715248
batchJ 56.551439 60.615248

Interval ini tampaknya cocok dengan hasil dari effects.

Tambahan 2

Alternatif lain adalah paket lsmeans . Ia memperoleh derajat kebebasan dan matriks kovarians yang disesuaikan dari paket pbkrtest .

> library("lsmeans")
> lsmeans(fm1, "batch")
Loading required namespace: pbkrtest
 batch   lsmean       SE    df lower.CL upper.CL
 A     62.26667 1.125709 40.45 59.99232 64.54101
 B     59.30000 1.125709 40.45 57.02565 61.57435
 C     62.05000 1.125709 40.45 59.77565 64.32435
 D     59.70000 1.125709 40.45 57.42565 61.97435
 E     55.90000 1.125709 40.45 53.62565 58.17435
 F     61.03333 1.125709 40.45 58.75899 63.30768
 G     59.90000 1.125709 40.45 57.62565 62.17435
 H     63.11667 1.125709 40.45 60.84232 65.39101
 I     58.68333 1.125709 40.45 56.40899 60.95768
 J     58.58333 1.125709 40.45 56.30899 60.85768

Confidence level used: 0.95 

effecteffectconfint$\pm1.96\times\mbox{se}$

Hasil dari effectdan lsmeansserupa, tetapi dengan situasi multi-faktor yang tidak seimbang, lsmeanssecara default rata-rata atas faktor-faktor yang tidak digunakan dengan bobot yang sama, sedangkan effectbobot dengan frekuensi yang diamati (tersedia sebagai opsi dalam lsmeans).