Effects paquete ofrece una forma muy rápida y cómoda de trazado de los resultados del modelo lineal de efectos mixtos obtenido a través de lme4 paquete . El effect calcula los intervalos de confianza (IC) muy rápidamente, pero cómo de confianza ¿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 IC calculados con effects paquete, el lote "E" no se solapa con el lote "A".
Si intento lo mismo utilizando confint.merMod y el método por defecto:
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 todas las IC se solapan. También obtengo advertencias que indican que la función no pudo calcular CIs fiables. Este ejemplo, y mi conjunto de datos real, me hace sospechar que effects toma atajos en el cálculo del IC que podrían no ser del todo aprobados por los estadísticos. ¿Qué grado de fiabilidad tienen los IC devueltos por effect función de effects paquete para lmer ¿Objetos?
Qué he probado: Buscando en el código fuente, me di cuenta de que effect se basa en Effect.merMod que a su vez dirige a Effect.mer que tiene el siguiente aspecto:
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 parece calcular la Matriz de Varianza-Covariante a partir de la lmer objeto:
effects:::mer.to.glm
function (mod)
{
...
mod2$vcov <- as.matrix(vcov(mod))
...
mod2
}Esto, a su vez, se utiliza probablemente en Effect.default para calcular los IC (puede que haya 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 los 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 los LMM, este enfoque parece sospechosamente simple.
2 votos
Cuando tengas líneas largas de código, te agradecería mucho que las dividieras en varias líneas para que no tengamos que desplazarnos para verlo todo.
1 votos
@rvl El código debería ser más legible ahora.