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

36

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, этот подход кажется подозрительно простым.

Микко
источник
1
Когда у вас есть длинные строки кода, я был бы очень признателен, если бы вы разбили их на несколько строк, чтобы нам не пришлось прокручивать все это.
RVL
1
@rvl Код должен быть более читабельным.
Микко

Ответы:

52

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

  • как указывает @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)
Бен Болкер
источник
2
Я принимаю этот ответ, поскольку он подходит к сути и дает хорошее сравнение между различными методами. Тем не менее, проверьте превосходный ответ RLV для получения дополнительной информации.
Микко
Спасибо за отличный и очень полезный ответ. Правильно ли я понимаю, что нельзя использовать CI для сравнения групп / партий, но можно сравнивать эффекты. Скажем, у меня было два курса лечения, несколько человек и несколько измерений у отдельных людей. Я использовал бы людей как случайный эффект, поскольку каждый из них содержал бы x измерений. Затем я хотел узнать, привели ли эти два метода к разным реакциям. Могу ли я использовать effectsпакет и CI перекрываются в этом случае?
Микко
5
Это более общий вопрос, который относится к любому стандартному модельно-ориентированному подходу. Может быть стоит отдельного вопроса. (1) В общем случае ответ на вопросы о различиях между обработками заключается в том, чтобы настроить модель таким образом, чтобы разница между фокусными обработками представляла собой контраст (т. Е. Оценочный параметр) в модели, а затем рассчитать значение p или проверьте, включают ли нулевые доверительные интервалы на конкретном альфа-уровне. (продолжение)
Бен Болкер
4
(2) перекрывающиеся CI являются в лучшем случае консервативным и приблизительным критерием различий между параметрами (существует несколько опубликованных работ на эту тему). (3) Существует отдельная / ортогональная проблема с парными сравнениями, которая заключается в том, что нужно надлежащим образом контролировать множественность и независимость сравнений (это можно сделать, например, с помощью методов в multcompпакете, но для этого требуется как минимум немного заботы)
Бен Болкер
1
Для чего? Вы можете задать новый вопрос.
Бен Болкер,
20

Похоже, что вы сделали во втором методе, чтобы вычислить доверительные интервалы для коэффициентов регрессии, а затем преобразовать их, чтобы получить 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 даже уже, чем асимптотические, основанные на использовании±1,96×се, Так что теперь я думаю, что они не очень заслуживают доверия.

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

RVL
источник
Спасибо за это решение. Интервалы теперь больше похожи, хотя и не совсем одинаковые. Ваш ответ по-прежнему не отвечает на вопрос, effectsможно ли доверять элементам конфигурации из пакета для lmerобъектов. Я рассматриваю возможность использования результатов в публикации и хочу убедиться, что CI рассчитываются с использованием утвержденного метода для LMM.
Микко
Не могли бы вы сказать: в вашем Приложении 1 первые два параметра .sig01и .sigmaпроизводные confint, являются ли эти доверительные интервалы для дисперсии ? или доверительный интервал стандартного отклонения ?
ABC
Они являются КЭ для любых параметров, помеченных таким образом в модели. Вы должны посмотреть на документацию для lmerокончательного ответа. Тем не менее, люди обычно используют обозначения любят sigmaдля обозначения стандартных отклонений, и sigma.squareили sigma^2для обозначения отклонений.
rvl
Лучше использовать lmertest, lsmeans или mertools?
Скан