¿Qué tan confiables son los intervalos de confianza para los objetos lmer a través del paquete de efectos?

EffectsEl paquete proporciona una manera muy rápida y conveniente de graficar los resultados del modelo de efectos mixtos lineales obtenidos a través delme4 paquete . La effectfunción calcula los intervalos de confianza (IC) muy rápidamente, pero ¿qué tan confiables son estos intervalos de confianza?

Por ejemplo:

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

Según los CI calculados usando el effectspaquete, el lote "E" no se superpone con el lote "A".

Si intento lo mismo usando la confint.merModfunción y el método predeterminado:

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

Veo que todos los CI se superponen. También recibo advertencias que indican que la función no pudo calcular los CI confiables. Este ejemplo, y mi conjunto de datos real, me hace sospechar que el effectspaquete toma atajos en el cálculo de CI que podrían no estar completamente aprobados por los estadísticos. ¿Cuán confiables son los CI devueltos por effectfunción del effectspaquete para lmerobjetos?

¿Qué he intentado? Al mirar el código fuente, noté que la effectfunción se basa en la Effect.merModfunción, que a su vez se dirige a la Effect.merfunción, que se ve así:

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.glmLa función parece calcular la matriz de varianza-covariable a partir del lmerobjeto:

effects:::mer.to.glm

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

Esto, a su vez, probablemente se usa en la Effect.defaultfunción para calcular los IC (podría haber 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
...

No sé lo suficiente sobre LMM para juzgar si este es un enfoque correcto, pero teniendo en cuenta la discusión sobre el cálculo del intervalo de confianza para LMM, este enfoque parece sospechosamente simple.

Todos los resultados son esencialmente los mismos ( para este ejemplo en particular ). Algunas diferencias teóricas son:

• Como @rvl señala, su reconstrucción de CI sin tener en cuenta la covarianza entre los parámetros es simplemente incorrecta (lo siento)
• los intervalos de confianza para los parámetros se pueden basar en los intervalos de confianza Wald (suponiendo una superficie de probabilidad logarítmica cuadrática): lsmeans, effects, confint(.,method="Wald"); a excepción de lsmeans, estos métodos ignoran los efectos de tamaño finito ("grados de libertad"), pero en este caso apenas hace ninguna diferencia ( df=40es prácticamente indistinguible de infinitodf )
• ... o en intervalos de confianza de perfil (el método predeterminado; ignora los efectos de tamaño finito pero permite superficies no cuadráticas)
• ... o en bootstrapping paramétrico (el estándar de oro: asume que el modelo es correcto [las respuestas son normales, los efectos aleatorios están normalmente distribuidos, los datos son condicionalmente independientes, etc.], pero por lo demás hacen pocas suposiciones)

Creo que todos estos enfoques son razonables. (algunos son más aproximados que otros), pero en este caso apenas hace ninguna diferencia cuál usa. Si le preocupa, pruebe varios métodos contrastantes en sus datos, o en datos simulados que se parecen a los suyos, y vea qué sucede ...

(PD: no pondría demasiado peso en el hecho de que los intervalos de confianza AyE no se superponen. Tendría que hacer un procedimiento de comparación por pares adecuado para hacer inferencias confiables sobre las diferencias entre este par de estimaciones en particular . ..)

IC 95%:

Código de comparación:

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 lo que ha hecho en el segundo método es haber calculado intervalos de confianza para los coeficientes de regresión, luego los ha transformado para obtener CI para las predicciones. Esto ignora las covarianzas entre los coeficientes de regresión.

Intente ajustar el modelo sin una intercepción, de modo que los batchefectos realmente sean las predicciones, yconfint le devuelvan los intervalos que necesita.

Anexo 1

Hice exactamente lo que sugerí anteriormente:

> 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

Estos intervalos parecen coincidir con los resultados de effects .

Anexo 2

Otra alternativa es el paquete lsmeans . Obtiene grados de libertad y una matriz de covarianza ajustada del paquete 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 

Estos están aún más en línea con los effectresultados: los errores estándar son idénticos, pero effectusan df diferentes. Los confintresultados en el Anexo 1 son aún más estrechos que los asintóticos basados ​​en el uso$\pm1.96\times\mbox{se}$. Así que ahora creo que esos no son muy confiables.

Los resultados de effecty lsmeansson similares, pero con una situación de factores múltiples desequilibrados, lsmeanspor defecto promedios sobre factores no utilizados con pesos iguales, mientras que los effectpesos por las frecuencias observadas (disponibles como una opción en lsmeans).