Насколько достоверны доверительные интервалы для lmer объектов через пакет эффектов?

EffectsПакет предоставляет очень быстрый и удобный способ для построения результатов линейной модели смешанного эффекта, полученных с помощью lme4пакета . В effectфункции вычисляет доверительные интервалы (ДИ) очень быстро, но , как заслуживающие доверия этих доверительные интервалы?

Например:

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()

Согласно КИ, рассчитанным с использованием effectsпакета, партия "Е" не перекрывается с партией "А".

Если я попробую то же самое, используя confint.merModфункцию и метод по умолчанию:

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()

Я вижу, что все CI перекрываются. Я также получаю предупреждения, указывающие на то, что функции не удалось рассчитать достоверные CI. Этот пример и мой фактический набор данных заставляют меня подозревать, что effectsпакет использует ярлыки в расчете CI, которые не могут быть полностью утверждены статистиками. Насколько достоверны CI, возвращаемые effectфункцией из effectsпакета для lmerобъектов?

Что я пробовал: изучая исходный код, я заметил, что effectфункция опирается на Effect.merModфункцию, которая, в свою очередь, указывает на Effect.merфункцию, которая выглядит следующим образом:

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.glmКажется, функция вычисляет матрицу дисперсии-ковариации из lmerобъекта:

effects:::mer.to.glm

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

Это, в свою очередь, вероятно, используется в Effect.defaultфункции для расчета КИ (возможно, я неправильно понял эту часть):

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
...

Я не знаю достаточно о LMM, чтобы судить, является ли это правильным подходом, но, учитывая дискуссию о расчете доверительного интервала для LMM, этот подход кажется подозрительно простым.

Все результаты по сути одинаковы ( для данного конкретного примера ). Некоторые теоретические различия:

• как указывает @rvl, ваша реконструкция КИ без учета ковариации между параметрами просто неверна (извините)
• доверительные интервалы для параметров могут быть основаны на доверительные интервалах Wald (при условии , квадратная логарифмическое правдоподобие поверхности): lsmeans, effects, confint(.,method="Wald"); за исключением того lsmeans, что эти методы игнорируют эффекты конечного размера («степени свободы»), но в этом случае они практически не имеют никакого значения ( df=40практически неотличимы от бесконечногоdf )
• ... или на доверительных интервалах профиля (метод по умолчанию; игнорирует эффекты конечного размера, но допускает неквадратичные поверхности)
• ... или при параметрической начальной загрузке (золотой стандарт - предполагается, что модель верна [ответы нормальны, случайные эффекты распределяются нормально, данные условно независимы и т. д.], но в противном случае делается несколько предположений)

Я думаю, что все эти подходы разумны (некоторые из них более приблизительны, чем другие), но в этом случае едва ли имеет значение, какой вы используете. Если вы обеспокоены, попробуйте несколько контрастных методов на ваших данных или на смоделированных данных, которые похожи на ваши собственные, и посмотрите, что происходит ...

(PS: Я бы не слишком большой вес на том , что доверительные интервалы Aи Eне перекрывают друг друга Вы должны были бы сделать надлежащую процедуру сравнения парное сделать достоверные выводы о различиях между этим. Определенной пары оценок. ..)

95% ДИ:

Код сравнения:

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)

Похоже, что вы сделали во втором методе, чтобы вычислить доверительные интервалы для коэффициентов регрессии, а затем преобразовать их, чтобы получить CI для предсказаний. Это игнорирует ковариации между коэффициентами регрессии.

Попробуйте подобрать модель без перехвата, чтобы batchэффекты фактически были предсказаниями и confintвозвращали нужные вам интервалы.

Приложение 1

Я сделал именно то, что я предложил выше:

> 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

Эти интервалы, похоже, соответствуют результатам из effects.

Приложение 2

Другой альтернативой является пакет lsmeans . Он получает степени свободы и скорректированную ковариационную матрицу из пакета 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 

Они даже больше соответствуют effectрезультатам: стандартные ошибки идентичны, но effectиспользуют разные значения df. confintРезультаты в Приложении 1 даже уже, чем асимптотические, основанные на использовании$\pm1.96\times\mbox{se}$, Так что теперь я думаю, что они не очень заслуживают доверия.

Результаты effectи lsmeansаналогичны, но с несбалансированной многофакторной ситуацией, lsmeansпо умолчанию усредняются по неиспользованным факторам с равными весами, тогда как effectвесы по наблюдаемым частотам (доступны в качестве опции в lsmeans).