一个混合效应模型的“显著”随机效应的毛毛虫图

Question

一个混合效应模型的“显著”随机效应的毛毛虫图

rggplot2lme4

8

我之前在这里得到了很好的帮助经验，现在又希望能得到一些帮助。

我正在估计一个相当大的混合效应模型，其中一个随机效应有超过150个不同的水平。这将使得标准的毛虫图变得难以阅读。

如果可能的话，我想要一个毛虫图，只显示随机效应中“显著”的水平（缺乏更好的术语）。也就是说：我需要一个毛虫图，其中随机拦截或变化系数的随机斜率具有不包括零的“置信区间”（我知道这不完全正确）。

考虑这个标准模型来自lme4的sleepstudy数据。

library(lme4)
fit <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
ggCaterpillar(ranef(fit,condVar=TRUE), QQ=FALSE, likeDotplot=TRUE, reorder=FALSE)[["Subject"]]

我想要得到这个毛毛虫图。

我所使用的毛毛虫图来自于这段代码。请注意，我倾向于使用较少保守的间隔范围（即1.645 * se而不是1.96 * se）。

基本上，我想要一个毛毛虫图，其中只包括308、309、310、330、331、335、337、349、350、352和370这些水平，因为这些水平要么具有截距，要么具有斜率，其区间不包括零。我这样做是因为我的毛毛虫图包含了150多个不同的水平，有点难以阅读，我认为这可能是一个有价值的解决方案。

下面是可重现的代码。非常感谢您的帮助。

# https://dev59.com/EJHea4cB1Zd3GeqPt8gN
ggCaterpillar <- function(re, QQ=TRUE, likeDotplot=TRUE, reorder=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, ])))
if (reorder) {
  ord  <- unlist(lapply(x, order)) + rep((0:(ncol(x) - 1)) * nrow(x), each=nrow(x))
  pDf  <- data.frame(y=unlist(x)[ord],
                     ci=1.645*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)))
} else {
  pDf  <- data.frame(y=unlist(x),
                     ci=1.645*se,
                     nQQ=rep(qnorm(ppoints(nrow(x))), ncol(x)),
                     ID=factor(rep(rownames(x), ncol(x)), levels=rownames(x)),
                     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 of the Random Effect") + ylab("Random Effect")
}

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


library(lme4)
fit <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
ggCaterpillar(ranef(fit,condVar=TRUE), QQ=FALSE, likeDotplot=TRUE, reorder=FALSE)[["Subject"]] 
ggsave(file="sleepstudy.png")

- steve

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

首先，感谢将“significant”用引号括起来...每个阅读此文的人都应该记住，在这个背景下，“significant”没有任何统计学意义（也许最好使用Z统计量（值/ std.error）标准，例如|Z|>1.5或|Z|>1.75，只是为了强调这不是一个推断性的阈值...）

我最终有些过度...我决定将事物重新组织/模块化一下，所以我编写了一个方法（旨在与软件包一起使用），从ranef.mer对象构建有用的数据框...完成后，所需的操作变得非常容易。

我将augment.ranef.mer代码放在我的答案末尾-它有点长（您需要在此处运行代码之前进行源代码处理）。更新：这个方法现在已经是broom.mixed包的一部分了...

library(broom)
library(reshape2)
library(plyr)

将 augment 方法应用于 RE 对象：

rr <- ranef(fit,condVar=TRUE)
aa <- augment(rr)

names(aa)
## [1] "grp"       "variable"  "level"     "estimate"  "qq"        "std.error"
## [7] "p"         "lb"        "ub"

现在ggplot的代码非常基础。我使用geom_errorbarh(height=0)而不是geom_pointrange()+coord_flip()，因为ggplot2无法与facet_wrap(...,scales="free")结合使用coord_flip...

## Q-Q plot:
g0 <- ggplot(aa,aes(estimate,qq,xmin=lb,xmax=ub))+
    geom_errorbarh(height=0)+
    geom_point()+facet_wrap(~variable,scale="free_x")

## regular caterpillar plot:
g1 <- ggplot(aa,aes(estimate,level,xmin=lb,xmax=ub))+
    geom_errorbarh(height=0)+
    geom_vline(xintercept=0,lty=2)+
    geom_point()+facet_wrap(~variable,scale="free_x")

现在找到您想要保留的级别：

aa2 <- ddply(aa,c("grp","level"),
             transform,
             keep=any(p<0.05))
aa3 <- subset(aa2,keep)

更新毛毡图，仅显示具有“显著”斜率或截距的水平：

g1 %+% aa3

如果你只想突出显示“显著”的水平，而不是完全删除“非显著”的水平

ggplot(aa2,aes(estimate,level,xmin=lb,xmax=ub,colour=factor(keep)))+
    geom_errorbarh(height=0)+
    geom_vline(xintercept=0,lty=2)+
    geom_point()+facet_wrap(~variable,scale="free_x")+
    scale_colour_manual(values=c("black","red"),guide=FALSE)

##' @importFrom reshape2 melt
##' @importFrom plyr ldply name_rows 
augment.ranef.mer <- function(x,
                                 ci.level=0.9,
                                 reorder=TRUE,
                                 order.var=1) {
    tmpf <- function(z) {
        if (is.character(order.var) && !order.var %in% names(z)) {
            order.var <- 1
            warning("order.var not found, resetting to 1")
        }
        ## would use plyr::name_rows, but want levels first
        zz <- data.frame(level=rownames(z),z,check.names=FALSE)
        if (reorder) {
            ## if numeric order var, add 1 to account for level column
            ov <- if (is.numeric(order.var)) order.var+1 else order.var
            zz$level <- reorder(zz$level, zz[,order.var+1], FUN=identity)
        }
        ## Q-Q values, for each column separately
        qq <- c(apply(z,2,function(y) {
                  qnorm(ppoints(nrow(z)))[order(order(y))]
              }))
        rownames(zz) <- NULL
        pv   <- attr(z, "postVar")
        cols <- 1:(dim(pv)[1])
        se   <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ])))
        ## n.b.: depends on explicit column-major ordering of se/melt
        zzz <- cbind(melt(zz,id.vars="level",value.name="estimate"),
                     qq=qq,std.error=se)
        ## reorder columns:
        subset(zzz,select=c(variable, level, estimate, qq, std.error))
    }
    dd <- ldply(x,tmpf,.id="grp")
    ci.val <- -qnorm((1-ci.level)/2)
    transform(dd,
              p=2*pnorm(-abs(estimate/std.error)), ## 2-tailed p-val
              lb=estimate-ci.val*std.error,
              ub=estimate+ci.val*std.error)
}