首页 文章

在R中，使用qqmath或dotplot绘制来自lmer（lme4包）的随机效果：如何使它看起来很花哨？

提问于
浏览 次
31

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组合，其他人也可以从中受益：

  • 有没有办法获取上面的qqmath函数并添加一些选项，例如，使某些点为空而不是填充，或者为不同的点设置不同的颜色？例如，您是否可以填充Batch变量的A，B和C的点数，但其余的点是否为空？

  • 是否可以为每个点添加轴标签（例如，可能沿顶部或右侧y轴）？

  • 我的数据接近45个截距，因此可以在标签之间添加间距，这样它们就不会相互碰撞？主要是，我感兴趣的是在图表上的点之间进行区分/标记，这在ggmath函数中似乎很麻烦/不可能 .

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

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

谢谢 .

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

r ggplot2 lme4 lmer random-effects

3 回答

  • 14

    一种可能性是使用库 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)

    enter image description here

    回复于
  • 34

    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

    enter image description here

    与 dotplot() 的比较：

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

    enter image description here

    有时，为随机效应设置不同的比例可能是有用的 - 这是 dotplot() 强制执行的 . 当我试图放松这个时，我不得不改变刻面（见这个answer） .

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

    enter image description here

    ## 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)
}
    回复于
  • 36

    另一种方法是从每个随机效应的分布中提取模拟值并绘制这些值 . 使用 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)

    哪个产生：

    A caterpillar plot of random effects

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

    回复于

相关问题