在R中,使用qqmath或dotplot绘制来自lmer(lme4包)的随机效果:如何使它看起来很花哨？

Question

qqmath函数使用lmer软件包的输出产生很好的随机效应图.也就是说,qqmath非常适合绘制层次模型中的截距及其在点估计周围的误差.下面使用名为Dyestuff的lme4包中的内置数据,给出了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函数并添加一些选项,例如,使某些点为空而不是填充,或不同点的不同颜色？例如,您是否可以填充Batch变量的A,B和C的点数,但其余的点是否为空？
2. 是否可以为每个点添加轴标签(例如,可能沿顶部或右侧y轴)？
3. 我的数据接近45个截距,因此可以在标签之间添加间距,这样它们就不会相互碰撞？主要是,我感兴趣的是在图表上的点之间进行区分/标记,这在ggmath函数中似乎很麻烦/不可能.

到目前为止,在qqmath函数中添加任何附加选项会产生错误,如果它是标准图,我不会得到错误,所以我很茫然.

另外,如果你觉得有一个更好的包装/功能来绘制来自lmer输出的拦截,我很乐意听到它!(例如,你能用dotplot做点1-3吗？)

谢谢.

编辑:如果可以合理格式化,我也可以使用替代的dotplot.我只是喜欢ggmath情节的外观,所以我开始提出一个问题.

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

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

此代码生成绘图.现在,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)

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)

哪个产生:

一个很好的特性是具有不重叠零的置信区间的值以黑色突出显示.您可以使用level参数plotREsim根据需要更改或更窄的置信区间来修改间隔的宽度.