使用qqmath或dotplot从lmer（lme4包）绘制随机效应：如何使其看起来漂亮？

Question

qqmath函数使用lmer软件包的输出来生成随机效应的优秀毛虫图。也就是说，qqmath非常适合绘制分层模型中拦截的估计值及其误差的图表。下面是使用lme4软件包中名为Dyestuff的内置数据的lmer和qqmath函数示例。该代码将生成分层模型以及使用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组合，则可以从这些问题中受益：

1. 是否有一种方法可以采取上面的qqmath函数并添加一些选项，例如使某些点为空心而不是填充的，或者为不同的点使用不同的颜色？例如，您可以使批次变量的A、B和C点填充，但然后使其余点为空？
2. 是否可能为每个点添加轴标签（例如，沿着顶部或右侧y轴）？
3. 我的数据更接近于45个截距，因此可能需要在标签之间添加间距，以便它们不互相重叠？主要是，我对区分/标记图中的点感兴趣，这似乎在ggmath函数中很麻烦/不可能。

到目前为止，在qqmath函数中添加任何其他选项都会产生错误，如果是标准图则不会出现这些错误，因此我感到迷茫。

此外，如果您觉得有更好的绘制lmer输出的截距的软件包/函数，我很乐意听取！ （例如，是否可以使用dotplot完成1-3点？）

编辑： 如果可以合理地格式化，我也可以接受另一种dotplot。 我只是喜欢ggmath图表的外观，因此从一个问题开始。

- Captain Murphy

Answer 1

Didzis的回答非常好！为了让它更加简洁明了，我将它放到了一个自己的函数中，这个函数与qqmath.ranef.mer()和dotplot.ranef.mer()有很多相似之处。除了Didzis的回答外，它还处理具有多个相关随机效应的模型（就像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)
}

- caracal

---

这个运行非常出色。但是如果要生成一个输出表，比如LaTeX，怎么办？ - bshor

@caracal 当你计算1.96*se[ord]时，为什么不需要考虑每个组中的观测数量？ - user3022875

很棒的功能，但同时会抛出一个警告。不过根据这个答案，我们只需要稍微改变对ggCaterpillar(ranef(fit, condVar=TRUE), QQ=FALSE, likeDotplot=FALSE)的调用即可。 - jay.sf

@jaySf 谢谢提醒！已修复。 - caracal

Answer 2

一个可能的选择是使用库 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)])

这段代码会生成一个图表。根据因子水平，现在点的形状将有所不同。

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)

- Didzis Elferts

---

非常感谢！看起来很不错。但在我给出奖励之前，我遇到了两个错误，分别是：无法找到您绘图代码中的“guides”函数和“theme”函数。我已经加载了ggplot2和scales库，但仍然出现错误。有任何想法吗？这些是不同的包吗？我仍然可以打印出一个图，但由于这些错误，它与原图不完全相同。此外，是否可能翻转轴，使水平线为Y轴（误差条将是水平的）？ - Captain Murphy

你应该更新你的ggplot（和scales）版本。最近的版本已经有了重大变化，包括使用theme（而不是opts）。 - mnel

@CaptainMurphy，sessionInfo（）显示你使用的是哪个版本的ggplot2？上述代码应该可以在最新版本的ggplot2中运行。 - MattBagg

@CaptainMurphy 我更新了我的解决方案以翻转轴。这个图是使用ggplot2版本0.9.3生成的。要使用这个版本的ggplot2，你的R版本至少应该是2.14。 - Didzis Elferts

@mohvd：根据?ranef，postVar已经过时了，但请注意，您可以使用condVar，即“condVar的（已弃用）同义词”；另一个令人困惑的点是，尽管使用了新的condVar名称，但相关属性仍然被命名为postVar。 - davedgd

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

生成的图表如下：

一个很好的特性是，那些置信区间不与零重叠的值会用黑色突出显示。您可以使用plotREsim中的level参数来调整区间的宽度，根据需要缩小或扩大置信区间。

- jknowles

Answer 4

通过集成在 sjPlot 包中的 plot_model() 命令，可以另一种获得所需的绘图方式。优点是该命令返回一个 ggplot 对象，因此有许多选项可以根据需要调整图像。我将示例保持简单，因为有许多选项可以个性化可视化 - 只需检查 ?plot_model 获取所有选项。

library(lme4)
library(sjPlot)
#?plot_model

data(Dyestuff, package = "lme4")
summary(Dyestuff)             

fit1 <- lmer(Yield ~ 1 + (1|Batch), Dyestuff) 
summary(fit1)

plot_model(fit1, type="re",
           vline.color="#A9A9A9", dot.size=1.5,
           show.values=T, value.offset=.2)

- Julian