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

36

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

ingrese la descripción de la imagen aquí

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

ingrese la descripción de la imagen aquí

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.

Mikko
fuente
1
Cuando tenga largas líneas de código, le agradecería mucho si las divide en varias líneas para que no tengamos que desplazarnos para verlo todo.
rvl
1
@rvl El código debería ser más legible ahora.
Mikko

Respuestas:

52

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%:

enter image description here

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)
Ben Bolker
fuente
2
Acepto esta respuesta, ya que es correcta y da una buena comparación entre diferentes métodos. Sin embargo, consulte la excelente respuesta de rlv para obtener más información.
Mikko
Gracias por una excelente y muy útil respuesta. ¿Entiendo correctamente que uno no puede usar CI para comparar grupos / lotes, pero es posible comparar efectos? Digamos que tuve dos tratamientos, varios individuos y varias mediciones dentro de los individuos. Usaría individuos como efecto aleatorio ya que cada uno de ellos contendría x medidas. Entonces quise saber si estos dos tratamientos dieron como resultado una respuesta diferente. ¿Podría usar el effectspaquete y la superposición de CI en este caso?
Mikko
55
Esta es una pregunta más general que es relevante para cualquier enfoque estándar basado en modelos. Podría valer una pregunta por separado. (1) En general, la forma en que uno responde preguntas sobre las diferencias entre los tratamientos es configurar el modelo de manera que la diferencia entre los tratamientos focales sea un contraste (es decir, un parámetro estimado) en el modelo, y luego calcular un valor p o verifique si los intervalos de confianza en un nivel alfa particular incluyen cero. (continuación)
Ben Bolker
44
(2) la superposición de IC es, en el mejor de los casos, un criterio conservador y aproximado para las diferencias entre los parámetros (hay varios documentos publicados sobre este tema). (3) Hay un problema separado / ortogonal con las comparaciones por pares, que es que uno tiene que controlar adecuadamente la multiplicidad y la no independencia de las comparaciones (esto se puede hacer, por ejemplo, por los métodos en el multcomppaquete, pero requiere al menos un un poco de cuidado)
Ben Bolker
1
¿Para qué? Es posible que desee hacer una nueva pregunta.
Ben Bolker,
20

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±1,96×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).

rvl
fuente
Gracias por esta solución Los intervalos son ahora más similares, aunque no exactamente iguales. Su respuesta aún no responde a la pregunta de si effectsse puede confiar en los CI del paquete para los lmerobjetos. Estoy considerando usar los resultados en una publicación y quiero estar seguro de que los IC se calculan utilizando un método aprobado para LMM.
Mikko
¿Podría decir: en su Anexo 1, los dos primeros parámetros .sig01y .sigmaproducidos por confint, son esos intervalos de confianza para la varianza ? o intervalo de confianza de la desviación estándar ?
ABC
Son elementos de configuración para cualquier parámetro etiquetado de esa manera en el modelo. Debe consultar la documentación para lmerobtener una respuesta definitiva. Sin embargo, las personas generalmente usan anotaciones como sigmapara referirse a desviaciones estándar y / sigma.squareo sigma^2para referirse a variaciones.
rvl
Is it better to use lmertest, lsmeans or mertools?
skan