Сюжет случайных эффектов от lmer (пакет lme4) с использованием qqmath или dotplot: как сделать так, чтобы это выглядело модно?

Функция qqmath создает отличные графики случайных эффектов, используя выходные данные пакета lmer. То есть qqmath отлично подходит для построения графиков пересечений из иерархической модели с их ошибками вокруг точечной оценки. Ниже приведены примеры функций lmer и qqmath с использованием встроенных данных в пакете lme4 под названием Dyestuff. Код создаст иерархическую модель и хороший график с использованием функции 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

Последняя строка кода создает действительно хороший график каждого перехвата с ошибкой вокруг каждой оценки. Но форматирование функции qqmath кажется очень сложным, и я изо всех сил пытался отформатировать сюжет. У меня есть несколько вопросов, на которые я не могу ответить, и я думаю, что другие могут также выиграть, если они используют комбинацию lmer/qqmath:

Есть ли способ взять функцию qqmath выше и добавить несколько такие варианты, как, например, сделать определенные точки пустыми или заполненными, или разные цвета для разных точек? Например, можно ли заполнить точки для A, B и C переменной Batch, но затем оставшиеся точки будут пустыми?
Можно ли добавить метки оси для каждой точки (возможно, вдоль например, верхняя или правая ось у)?
Мои данные ближе к 45 перехватам, поэтому можно добавить расстояние между метками, чтобы они не сталкивались друг с другом? В основном, я заинтересован в различении/маркировке точек на график, который кажется громоздким/невозможным в функции ggmath.

Пока что добавление любой дополнительной опции в функцию qqmath приводит к ошибкам, при которых я не получил бы ошибок, если бы это был стандартный график, поэтому я в растерянности.

Кроме того, если вы чувствуете, что есть лучший пакет/функция для построения перехватов из вывода lmer, я бы хотел это услышать! (например, можете ли вы использовать точки 1-3 с помощью точечного графика?)

ОБНОВЛЕНИЕ: Я также открыт для альтернативного точечного графика, если он может быть разумно отформатирован. Мне просто нравится вид сюжета ggmath, поэтому я начинаю с вопроса об этом.

Ответы

Ответ 1

Одна из возможностей - использовать библиотеку ggplot2 для рисования аналогичного графика, а затем вы можете настроить внешний вид вашего сюжета.

Сначала объект ranef сохраняется как randoms. Затем дисперсии перехватов сохраняются в объекте qq.

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

Объект rand.interc содержит только случайные перехваты с именами уровней.

rand.interc<-randoms$Batch

Все объекты помещаются в один фрейм данных. Для интервалов ошибок sd.interc вычисляется как дисперсия в 2 раза квадратный корень.

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

Если вам нужно, чтобы перехваты упорядочивались в зависимости от значения, тогда lev.names следует переупорядочить. Эта строка может быть пропущена, если перехваты должны быть упорядочены по именам уровней.

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

Этот код создает график. Теперь точки будут отличаться на shape в соответствии с уровнями факторов.

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)

Ответ 2

Дидзис ответ велик! Просто чтобы немного обернуть его, я включил его в собственную функцию, которая во многом похожа на qqmath.ranef.mer() и dotplot.ranef.mer(). В дополнение к ответу Дидзиса, он также обрабатывает модели с несколькими коррелированными случайными эффектами (как это делают qqmath() и dotplot()). Сравнение с 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

Сравнение с dotplot():

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

Иногда может быть полезно иметь разные шкалы для случайных эффектов - то, что обеспечивает dotplot(). Когда я попытался это ослабить, мне пришлось сменить фасетку (см. этот ответ).

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

Ответ 3

Другой способ сделать это - извлечь симулированные значения из распределения каждого из случайных эффектов и построить их. Используя пакет merTools, можно легко получить симуляции от объекта lmer или glmer и нарисовать их.

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

randoms теперь является объектом, который выглядит так:

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

Он предоставляет имя фактора группировки, уровень фактора, который мы получаем для оценки, член в модели и среднее, среднее и стандартное отклонение имитируемых значений. Мы можем использовать это для создания графика гусеницы, аналогичного приведенному выше:

plotREsim(randoms)

Что производит:

Одна приятная особенность заключается в том, что значения, имеющие доверительный интервал, который не перекрывает нуль, выделяются черным цветом. Вы можете изменить ширину интервала, используя параметр level, чтобы plotREsim сделать более широкие или более узкие доверительные интервалы на основе ваших потребностей.