require(arm) #you'll need this later to look up the key in a list reverse.lookup<-function(a.list, a.value){ names(a.list)[which(as.numeric(a.list)==a.value)] } #given a fixed effect, which grouping make it vary? get.names.and.lengths.relevant.ranefs<-function(lmer.obj, rans, fixed){ relevant.rans<-list() for (ran in rans) { if(ranef(lmer.obj)[[ran]][[fixed]][1]){relevant.rans[[ran]]<-length(ranef(lmer.obj)[[ran]][[fixed]])} } return(relevant.rans) } #I'm sure this is something done elsewhere more elegantly #but, this takes two vectors of differing lengths, and #looks at all possible additions of them crossadd.vectors<-function(v1, v2){ as.vector(apply(as.array(v1), 1, function (x) apply (as.array(v2), 1, function(y) x+y))) } #This is like the above function, only it combines strings crosspaste.vectors<-function(v1, v2, sepstring=","){ as.vector(apply(as.array(v1), 1, function (x) apply (as.array(v2), 1, function(y) paste(x,y,sep=sepstring)))) } #take a vector of points and some measure around them (assumes it's SE) #and plots them. Also, it can color by certain types, and make divisions #between groupings (e.g., for a model with multiple factors) plot.mean.and.se<-function(point.vector, se.vector, name.vector, fixed, nse=1, top=NA, bottom=NA, plot.types=1, div.lines=0, ...){ if(is.na(top)) top<-max(point.vector+nse*se.vector) if(is.na(bottom)) bottom<-min(point.vector-nse*se.vector) types<-"black" if(plot.types>1) types<-rainbow(plot.types) plot(1:length(point.vector), point.vector, ylim=c(bottom,top), ylab=paste(fixed, "coefficient +/-",nse, "SE", sep=" "), xlab="groupings", col=types, ...) for (x in 1:length(point.vector)){ lines(c(x,x),c(point.vector[x]+nse*se.vector[x],point.vector[x]-nse*se.vector[x])) } if(div.lines>0){ print(div.lines) len<-length(point.vector) for(div in seq(div.lines,len, by=div.lines)){ lines(c(div+0.5,div+0.5), c(top,bottom), lty=2) } } } #takes an lmer object, and plots the coefficient estimates #+/- 2 standard errors, after accounting for groupings ploteffects.lmer <- function(lmer.obj,min.max=F, nse=2, plot.types=F, draw.divs=F,...){ require(arm) #what are the fixed and "random" effects fixef.names<-names(fixef(lmer.obj)) ranef.names<-names(ranef(lmer.obj)) par(ask=T, pch=19, cex=1.6) #for each fixed effect for (fixed in fixef.names) { #figure out which random effects apply relevant.rans<-get.names.and.lengths.relevant.ranefs(lmer.obj, ranef.names, fixed) #now, put together mean and se vectors num.fields<-sort(as.numeric(relevant.rans), decreasing = min.max) point.vector<-fixef(lmer.obj)[[fixed]] se.vector<-0 name.vector<-c("") for (field.num in num.fields) { field<-reverse.lookup(relevant.rans, field.num) point.vector<-crossadd.vectors(point.vector,ranef(lmer.obj)[[field]][[fixed]]) se.vector<-crossadd.vectors(se.vector, as.data.frame(se.ranef(lmer.obj)[[field]])[[fixed]]) name.vector<-crosspaste.vectors(name.vector,rownames(ranef(lmer.obj)[[field]]),sepstring=",") } types<-1 div.lines<-0 if(plot.types) types<-num.fields[length(num.fields)] if(length(num.fields)>1 & draw.divs==T) div.lines<-num.fields[length(num.fields)] plot.mean.and.se(point.vector, se.vector, name.vector, fixed, nse=nse, plot.types=types, div.lines=div.lines, ...) print(name.vector) } }