Firstly, thanks for quoting the "significant" ones ... anyone reading this should remember that the value has no statistical value in this context (it would be better to use the Z statistics (value / std. Error ), for example, | Z |> 1.5 or | Z |> 1.75, just to emphasize that this is not a withdrawal threshold ...)
In the end, I got a little carried away ... I decided that it would be better to refactor / modulate things a bit, so I wrote the augment method (designed to work with the broom package) that creates useful data frames from ranef.mer objects ... as soon as this will be done, the required manipulations are quite simple.
I put the augment.ranef.mer code at the end of my answer - it is a bit long (you will need to specify it before you can run the code here).
library(broom) library(reshape2) library(plyr)
Apply the augment method to the RE object:
rr <- ranef(fit,condVar=TRUE) aa <- augment(rr) names(aa) ## [1] "grp" "variable" "level" "estimate" "qq" "std.error" ## [7] "p" "lb" "ub"
Now the ggplot code ggplot pretty simple. I use geom_errorbarh(height=0) rather than geom_pointrange()+coord_flip() because ggplot2 cannot use coord_flip with facet_wrap(...,scales="free") ...
## QQ 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")
Now find the levels you want to keep:
aa2 <- ddply(aa,c("grp","level"), transform, keep=any(p<0.05)) aa3 <- subset(aa2,keep)
Refresh the track pattern only with levels of “significant” slopes or intercepts:
g1 %+% aa3
If you want to highlight only "significant" levels, and not completely remove the "non-essential" levels
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) } ## QQ 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, ]))) ## nb: 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) }