A function to obtain predictions from a fitted lmer object.

ezPredict.R

#A function to obtain predictions from a fitted lmer object. 
#If to_predict is null, the function will attempt to build a 
#design matrix for prediction from the information in the fit 
#itself. If a .() object is passed to as to_predict, variables 
#listed in to_return are used to create the design matrix. A
#data frame may also be passed as to_predict.

ezPredict <-
function(
	fit
	, to_predict = NULL
	, numeric_res = 0
){
	data = attr(fit,'frame')
	vars = as.character(attr(attr(data,'terms'),'variables'))
	dv = as.character(vars[2])
	if(!is.data.frame(to_predict)){
		if(length(grep('poly(',vars,fixed=TRUE))>0){
			stop('Cannot auto-create "to_predict" when the fitted model contains poly(). Please provide a non-null value to the "to_return" argument.')
		}
		if(is.null(to_predict)){
			vars = vars[3:length(vars)]
		}else{
			vars = vars[vars%in%to_predict]
			temp = attr(attr(fit,'X'),'contrasts')
			if(!all(temp[names(temp)%in%vars]%in%c('contr.sum','contr.helmert','contr.poly'))){
				#stop('When using ezPredict to obtain predictions for a subset of variables in a model, the model must be fit after setting either "options(contrasts=c(\'contr.sum\',\'contr.poly\')" or "options(contrasts=c(\'contr.helmert\',\'contr.poly\')".')
			}
		}
		temp = list()
		for(i in 1:length(vars)){
			this_fixed_data = data[,names(data)==vars[i]]
			if(is.numeric(this_fixed_data)&(numeric_res>0)){
				temp[[i]] = seq(
					min(this_fixed_data)
					, max(this_fixed_data)
					, length.out=numeric_res
				)
			}else{
				temp[[i]] = unique(this_fixed_data)
			}
		}
		to_return = data.frame(expand.grid(temp))
		names(to_return) = vars
	}else{
		to_return = to_predict
	}
	to_return$ezDV = 0
	names(to_return)[ncol(to_return)] = dv
	requested_terms = terms(eval(parse(text=paste(dv,'~',paste(vars,collapse='*')))))
	mm = model.matrix(requested_terms,to_return)
	f = fixef(fit)
	indicies = (1:length(names(f)))[names(f) %in% dimnames(mm)[[2]]]
	f = f[indicies]
	value = mm %*% f
	to_return$value = as.numeric(value[,1])
	vf = vcov(fit)[indicies,indicies]
	tc = Matrix::tcrossprod(vf,mm)
	to_return$var = Matrix::diag(mm %*% tc)
	to_return = to_return[,names(to_return)!=dv]	
	return(to_return)
}