total - modelos de efectos aleatorios ejercicios resueltos

En R, trazando los efectos aleatorios de lmer(paquete lme4) usando qqmath o dotplot: ¿cómo hacer que se vea elegante?

Question

La función qqmath hace grandes tramas de oruga de efectos aleatorios usando la salida del paquete lmer. Es decir, qqmath es excelente para trazar las intersecciones de un modelo jerárquico con sus errores en torno a la estimación puntual. Un ejemplo de las funciones lmer y qqmath están debajo utilizando los datos incorporados en el paquete lme4 llamado colorante. El código producirá el modelo jerárquico y una buena trama usando la función ggmath.

library("lme4")
data(package = "lme4")

# Dyestuff 
# a balanced one-way classiï¬cation of Yield 
# from samples produced from six Batches

summary(Dyestuff)             

# Batch is an example of a random effect
# Fit 1-way random effects linear model
fit1 <- lmer(Yield ~ 1 + (1|Batch), Dyestuff) 
summary(fit1)
coef(fit1) #intercept for each level in Batch 

# qqplot of the random effects with their variances
qqmath(ranef(fit1, postVar = TRUE), strip = FALSE)$Batch

La última línea de código produce una trama realmente bonita de cada intersección con el error alrededor de cada estimación. Pero formatear la función qqmath parece ser muy difícil, y he estado luchando por formatear la trama. He llegado a algunas preguntas que no puedo responder, y creo que otros también podrían beneficiarse si usan la combinación lmer / qqmath:

1. ¿Hay alguna manera de llevar a cabo la función qqmath anterior y agregar algunas opciones, como por ejemplo, hacer que ciertos puntos estén vacíos o no rellenos, o diferentes colores para diferentes puntos? Por ejemplo, ¿puede completar los puntos para A, B y C de la variable Lote, pero luego el resto de los puntos está vacío?
2. ¿Es posible agregar etiquetas de eje para cada punto (tal vez a lo largo del eje y superior o derecho, por ejemplo)?
3. Mis datos tienen más de 45 interceptaciones, por lo que es posible agregar espacio entre las etiquetas para que no se encuentren entre sí. PRINCIPALMENTE, estoy interesado en distinguir / etiquetar entre puntos en el gráfico, que parece ser engorroso / imposible en la función de ggmath.

Hasta ahora, agregar cualquier opción adicional en la función qqmath produce errores en los que no obtendría errores si se tratara de un diagrama estándar, por lo que estoy perdido.

TAMBIÉN, si sientes que hay un paquete / función mejor para trazar intersecciones de la salida de lmer, ¡me encantaría escucharlo! (por ejemplo, ¿puedes hacer los puntos 1-3 usando dotplot?)

Gracias.

EDITAR: También estoy abierto a un gráfico de puntos alternativo si se puede formatear razonablemente. Simplemente me gusta el aspecto de una trama de ggmath, así que estoy empezando con una pregunta sobre eso.

Answer 1

¡La respuesta de Didzis es genial! Solo para resumirlo un poco, lo puse en su propia función que se comporta muy parecido a qqmath.ranef.mer() y dotplot.ranef.mer() . Además de la respuesta de Didzis, también maneja modelos con múltiples efectos aleatorios correlacionados (como qqmath() y dotplot() do). Comparación con qqmath() :

require(lme4)                            ## for lmer(), sleepstudy
require(lattice)                         ## for dotplot()
fit <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
ggCaterpillar(ranef(fit, condVar=TRUE))  ## using ggplot2
qqmath(ranef(fit, condVar=TRUE))         ## for comparison

Comparación con dotplot() :

ggCaterpillar(ranef(fit, condVar=TRUE), QQ=FALSE)
dotplot(ranef(fit, condVar=TRUE))

A veces, puede ser útil tener diferentes escalas para los efectos aleatorios, algo que dotplot() impone. Cuando traté de relajar esto, tuve que cambiar la faceta (ver esta answer ).

ggCaterpillar(ranef(fit, condVar=TRUE), QQ=FALSE, likeDotplot=FALSE)

## re = object of class ranef.mer
ggCaterpillar <- function(re, QQ=TRUE, likeDotplot=TRUE) {
    require(ggplot2)
    f <- function(x) {
        pv   <- attr(x, "postVar")
        cols <- 1:(dim(pv)[1])
        se   <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ])))
        ord  <- unlist(lapply(x, order)) + rep((0:(ncol(x) - 1)) * nrow(x), each=nrow(x))
        pDf  <- data.frame(y=unlist(x)[ord],
                           ci=1.96*se[ord],
                           nQQ=rep(qnorm(ppoints(nrow(x))), ncol(x)),
                           ID=factor(rep(rownames(x), ncol(x))[ord], levels=rownames(x)[ord]),
                           ind=gl(ncol(x), nrow(x), labels=names(x)))

        if(QQ) {  ## normal QQ-plot
            p <- ggplot(pDf, aes(nQQ, y))
            p <- p + facet_wrap(~ ind, scales="free")
            p <- p + xlab("Standard normal quantiles") + ylab("Random effect quantiles")
        } else {  ## caterpillar dotplot
            p <- ggplot(pDf, aes(ID, y)) + coord_flip()
            if(likeDotplot) {  ## imitate dotplot() -> same scales for random effects
                p <- p + facet_wrap(~ ind)
            } else {           ## different scales for random effects
                p <- p + facet_grid(ind ~ ., scales="free_y")
            }
            p <- p + xlab("Levels") + ylab("Random effects")
        }

        p <- p + theme(legend.position="none")
        p <- p + geom_hline(yintercept=0)
        p <- p + geom_errorbar(aes(ymin=y-ci, ymax=y+ci), width=0, colour="black")
        p <- p + geom_point(aes(size=1.2), colour="blue") 
        return(p)
    }

    lapply(re, f)
}

Answer 2

Otra forma de hacerlo es extraer valores simulados de la distribución de cada uno de los efectos aleatorios y representarlos. Usando el paquete merTools , es posible obtener fácilmente las simulaciones de un objeto lmer o glmer , y glmer .

library(lme4); library(merTools)       ## for lmer(), sleepstudy
fit <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
randoms <- REsim(fit, n.sims = 500)

randoms ahora es un objeto con lo que se ve así:

head(randoms)
groupFctr groupID        term       mean     median       sd
1   Subject     308 (Intercept)   3.083375   2.214805 14.79050
2   Subject     309 (Intercept) -39.382557 -38.607697 12.68987
3   Subject     310 (Intercept) -37.314979 -38.107747 12.53729
4   Subject     330 (Intercept)  22.234687  21.048882 11.51082
5   Subject     331 (Intercept)  21.418040  21.122913 13.17926
6   Subject     332 (Intercept)  11.371621  12.238580 12.65172

Proporciona el nombre del factor de agrupación, el nivel del factor para el que estamos obteniendo una estimación, el término en el modelo y la media, mediana y desviación estándar de los valores simulados. Podemos usar esto para generar un diagrama de oruga similar a los anteriores:

plotREsim(randoms)

Que produce:

Una buena característica es que los valores que tienen un intervalo de confianza que no se superpone a cero están resaltados en negro. Puede modificar el ancho del intervalo utilizando el parámetro de level para plotREsim haciendo intervalos de confianza más amplios o más estrechos según sus necesidades.

Answer 3

Una posibilidad es usar la biblioteca ggplot2 para dibujar un gráfico similar y luego puede ajustar la apariencia de su trazado.

Primero, el objeto ranef se guarda como randoms . Entonces las varianzas de las intersecciones se guardan en el objeto qq .

randoms<-ranef(fit1, postVar = TRUE)
qq <- attr(ranef(fit1, postVar = TRUE)[[1]], "postVar")

El objeto rand.interc contiene solo intercepciones aleatorias con nombres de nivel.

rand.interc<-randoms$Batch

Todos los objetos puestos en un marco de datos. Para los intervalos de error sd.interc se calcula como 2 veces la raíz cuadrada de la varianza.

df<-data.frame(Intercepts=randoms$Batch[,1],
              sd.interc=2*sqrt(qq[,,1:length(qq)]),
              lev.names=rownames(rand.interc))

Si necesita que las interceptaciones se ordenen en parcela de acuerdo con su valor, entonces los lev.names deben reordenarse. Esta línea se puede omitir si las interceptaciones se deben ordenar por nombres de nivel.

df$lev.names<-factor(df$lev.names,levels=df$lev.names[order(df$Intercepts)])

Este código produce un diagrama. Ahora los puntos difieren según la shape según los niveles de los factores.

library(ggplot2)
p <- ggplot(df,aes(lev.names,Intercepts,shape=lev.names))

#Added horizontal line at y=0, error bars to points and points with size two
p <- p + geom_hline(yintercept=0) +geom_errorbar(aes(ymin=Intercepts-sd.interc, ymax=Intercepts+sd.interc), width=0,color="black") + geom_point(aes(size=2)) 

#Removed legends and with scale_shape_manual point shapes set to 1 and 16
p <- p + guides(size=FALSE,shape=FALSE) + scale_shape_manual(values=c(1,1,1,16,16,16))

#Changed appearance of plot (black and white theme) and x and y axis labels
p <- p + theme_bw() + xlab("Levels") + ylab("")

#Final adjustments of plot
p <- p + theme(axis.text.x=element_text(size=rel(1.2)),
               axis.title.x=element_text(size=rel(1.3)),
               axis.text.y=element_text(size=rel(1.2)),
               panel.grid.minor=element_blank(),
               panel.grid.major.x=element_blank())

#To put levels on y axis you just need to use coord_flip()
p <- p+ coord_flip()
print(p)