# plot MCMC-estimated coefficients of linear mixed effects model # "lines" adds some lines to represent the intercept # pass it either a linear m-e model, or the output of pvals.fnc # returns "pvals.fnc" output # based on code by Florian Jaeger and Corey McMillan plot.coeffs <- function(y,lines=FALSE,pub=FALSE,names=FALSE,title=TRUE,ylim='auto',sym=FALSE) { # require(Hmisc, quietly=TRUE) bp <- function(...) { # do the plot return(barplot(c, names = n, ylim= y.limits, main= m, ylab= "coefficient estimates", las = 2, border='black', cex= .9, cex.main=1.2, cex.lab=.9, cex.axis=.7, ... )) } # recent versions of lme4 have class "mer" instead of "glmer" if(class(y)[1] == "mer" || class(y)[1] == "glmer") { require(languageR, quietly=TRUE) cat("running pvals.fnc...\n") pls <- pvals.fnc(y,addPlot=FALSE) } else { pls <- y } # vectors of relevant values (not the first, which is the intercept) c <- as.vector(as.numeric(pls[["fixed"]]$MCMCmean))[-1] lc <- as.vector(as.numeric(pls[["fixed"]]$HPD95lower))[-1] hc <- as.vector(as.numeric(pls[["fixed"]]$HPD95upper))[-1] if (class(names) != 'logical') { n <- names } else { n <- row.names(pls[["fixed"]])[-1] } if (title) { m <- "model estimates with 95% confidence intervals" } else { m <- '' } intercept <- as.numeric(pls[["fixed"]]$MCMCmean)[1] hi=as.numeric(pls[["fixed"]]$HPD95upper)[1] lo=as.numeric(pls[["fixed"]]$HPD95lower)[1] # change lower margin for longer column names # (better to do this ad-hoc rather than within a function) # saved.par <- par(no.readonly=TRUE) # on.exit(par(saved.par)) # par(mar=c(8,4,4,2)+0.1) if (class(ylim) != 'numeric') { if (lines) { y.limits=c(min(c(lc,0,lo-intercept)),max(c(hc,0,hi-intercept))) } else { y.limits=c(min(c(lc,0)), max(c(hc,0))) } if (sym) { if (y.limits[1]<0 & y.limits[2]>0) { if (abs(y.limits[1])>y.limits[2]) { y.limits[2]=abs(y.limits[1]) } else { y.limits[1]=-y.limits[2] } } } } else { y.limits=ylim } if (pub) { x <- bp(col='black',density=9) } else { x <- bp(col=c(rep("#CCCCCC",3))) } # add error bars (x is a vector of column positions) # errbar(x,c,hc,lc,lwd=1.5,add=TRUE,pch=' ') arrows(x,hc,x,lc,length=0.10,angle=90,code=3,lwd=1.5) # add lines if requested; subtitle appropriately to ensure no info lost if (lines) { hi = hi - intercept lo = lo - intercept abline(h=hi,lty="dashed",lwd=.5) abline(h=lo,lty="dashed",lwd=.5) abline(h=0,lty="solid",lwd=1) if (title) { mtext(paste('estimated intercept',intercept,'(dotted lines show 95% ci)'),line=.5) } } else { if (title) { mtext(paste(sep='','estimated intercept ',intercept," (ci ",lo," - ",hi,")"),line=.5) } } # return "pvals.fnc" list return(pls) }