只是混合效应模型中"重要"随机效应的毛虫图

首先,感谢将"重要"放在引号中...每个阅读此内容的人都应该记住,在这种情况下,重要性没有任何统计意义(使用Z-statistic(value/std.error)标准可能更好例如| Z |> 1.5或| Z |> 1.75而不是强调这不是推理阈值......)

我最终得到了一点点......我决定稍微重构/模块化一些东西会更好,所以我写了一个augment方法(设计用于处理broom包),从ranef.mer对象构造有用的数据帧......一旦完成,您想要的操作非常简单.

我把augment.ranef.mer代码放在我的答案的最后 - 它有点长(你需要在你可以运行代码之前获取它).

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(free")...

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