Qual é a confiabilidade dos intervalos de confiança para objetos antigos através do pacote de efeitos?

EffectsO pacote fornece uma maneira muito rápida e conveniente de plotar resultados lineares de modelo de efeito misto obtidos através do lme4pacote . A effectfunção calcula intervalos de confiança (ICs) muito rapidamente, mas quão confiáveis são esses intervalos de confiança?

Por exemplo:

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

De acordo com os ICs calculados usando o effectspacote, o lote "E" não se sobrepõe ao lote "A".

Se eu tentar o mesmo usando a confint.merModfunção e o método padrão:

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

Vejo que todos os ICs se sobrepõem. Também recebo avisos indicando que a função falhou ao calcular ICs confiáveis. Este exemplo e meu conjunto de dados real me fazem suspeitar que o effectspacote usa atalhos no cálculo do IC que podem não ser totalmente aprovados pelos estatísticos. Qual é a confiabilidade dos ICs retornados por effectfunção do effectspacote para lmerobjetos?

O que tentei: Examinando o código-fonte, notei que a effectfunção depende da Effect.merModfunção, que por sua vez direciona a Effect.merfunção, que se parece com isso:

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.glmA função parece calcular a matriz covariável de variância a partir do lmerobjeto:

effects:::mer.to.glm

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

Isso, por sua vez, provavelmente é usado em Effect.defaultfunção para calcular ICs (eu poderia ter entendido mal esta parte):

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

Eu não sei o suficiente sobre LMMs para julgar se essa é uma abordagem correta, mas, considerando a discussão sobre o cálculo do intervalo de confiança para LMMs, essa abordagem parece suspeitamente simples.

Todos os resultados são essencialmente os mesmos ( para este exemplo em particular ). Algumas diferenças teóricas são:

• como @rvl aponta, sua reconstrução de ICs sem levar em conta a covariância entre parâmetros está errada (desculpe)
• intervalos de confiança para parâmetros podem ser baseados em intervalos de confiança Wald (assumindo uma superfície de probabilidade logarítmica quadrática): lsmeans, effects, confint(.,method="Wald"); exceto por lsmeans, esses métodos ignoram efeitos de tamanho finito ("graus de liberdade"), mas, neste caso, quase não faz diferença ( df=40é praticamente indistinguível de infinito df)
• ... ou em intervalos de confiança do perfil (o método padrão; ignora efeitos de tamanho finito, mas permite superfícies não quadráticas)
• ... ou no bootstrapping paramétrico (o padrão-ouro - assume que o modelo está correto [respostas são normais, efeitos aleatórios são normalmente distribuídos, dados são condicionalmente independentes, etc.], mas, caso contrário, faz poucas suposições)

Eu acho que todas essas abordagens são razoáveis (algumas são mais aproximadas que outras), mas nesse caso, quase não faz diferença qual delas você usa. Se estiver preocupado, experimente vários métodos contrastantes nos seus dados ou em dados simulados que se assemelhem aos seus e veja o que acontece ...

(PS: Eu não daria muita importância ao fato de os intervalos de confiança Ae Enão se sobreporem. Você teria que fazer um procedimento de comparação pareada adequado para fazer inferências confiáveis ​​sobre as diferenças entre esse par de estimativas em particular . ..)

ICs de 95%:

Código de comparação:

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)

Parece que o que você fez no segundo método é ter intervalos de confiança calculados para os coeficientes de regressão e depois transformá-los para obter ICs para as previsões. Isso ignora as covariâncias entre os coeficientes de regressão.

Tente ajustar o modelo sem interceptar, para que os batchefeitos sejam realmente as previsões e confintretorne os intervalos necessários.

Adenda 1

Fiz exatamente o que sugeri acima:

> 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

Esses intervalos parecem combinar com os resultados de effects.

Adenda 2

Outra alternativa é o pacote lsmeans . Ele obtém graus de liberdade e uma matriz de covariância ajustada do pacote 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}$

Os resultados de effecte lsmeanssão semelhantes, mas com uma situação multifatorial desequilibrada, lsmeanspor padrão, calcula a média dos fatores não utilizados com pesos iguais, enquanto os effectpesos pelas frequências observadas (disponível como opção em lsmeans).