an attempt to compute reliability (and contributing variances) using variance estimates from a mixed effects model

get_vars_and_rel.R

#define two useful helper functions
colMean = function(x){
	dimx = dim(x)
	.Internal(colMeans(x,dimx[1],dimx[2],na.rm=TRUE))
}
colVar = function(x){
	dimx = dim(x)
	x.mean = .Internal(colMeans(x,dimx[1],dimx[2],na.rm=TRUE))
	err = t(t(x)-x.mean)
	err.sq = err*err
	sum.err.sq = .Internal(colSums(err.sq,dimx[1],dimx[2],na.rm=TRUE))
	n = .Internal(colSums(!is.na(x),dimx[1],dimx[2],na.rm=TRUE))
	sum.err.sq/(n-1)
}

#a possibly statistically sketchy to compute reliability and its CI
get_vars_and_rel = function(
	fit
	, iterations = 1e4 #1e4 takes a couple seconds for reasonably simple models
){

	# Extract BLUPS and associated variances
	r = ranef(fit, postVar = TRUE)[[1]]
	postVar = attr( r, "postVar")
	effects = names(r)
	r = matrix(unlist(r),dim(r)[1],dim(r)[2])
	postVar_dims = dim(postVar)
	within = matrix(NA,nrow=postVar_dims[3],ncol=postVar_dims[1])
	for(i in 1:postVar_dims[3]){
		within[i,] = diag(postVar[,,i])
	}
	
	# Bootstrap distributions for:
	#    1) between-Ss variance of BLUPs
	#    2) mean BLUP within-Ss variance
	#    3) analytic reliability (a straight arithmetic function of #1 & #2)
	
	within_boots = matrix(NA,nrow=iterations,ncol=ncol(within))
	between_boots = matrix(NA,nrow=iterations,ncol=ncol(within))
	rel_boots = matrix(NA,nrow=iterations,ncol=ncol(within))
	for(i in 1:iterations){
		rows = sample(nrow(within),replace=T)
		within_boots[i,] = colMean(within[rows,])
		between_boots[i,] = colVar(r[rows,])
		rel_boots[i,] = 1/(1+within_boots[i,]/between_boots[i,])
	}
	
	# Compute CIs and organize into data frame	
	within_mean = sqrt(colMean(within_boots))
	within_lo = sqrt(apply(within_boots,2,quantile,.025))
	within_hi = sqrt(apply(within_boots,2,quantile,.975))
	within = as.data.frame(cbind(within_mean,within_lo,within_hi))
	names(within) = c('Mean','lo','hi')
	within$effect = effects
	within$type = 'within'
	
	between_mean = sqrt(colMean(between_boots))
	between_lo = sqrt(apply(between_boots,2,quantile,.025))
	between_hi = sqrt(apply(between_boots,2,quantile,.975))
	between = as.data.frame(cbind(between_mean,between_lo,between_hi))
	names(between) = c('Mean','lo','hi')
	between$effect = effects
	between$type = 'between'
	
	rel_mean = colMean(rel_boots)
	rel_lo = apply(rel_boots,2,quantile,.025)
	rel_hi = apply(rel_boots,2,quantile,.975)
	rel = as.data.frame(cbind(rel_mean,rel_lo,rel_hi))
	names(rel) = c('Mean','lo','hi')
	rel$effect = effects
	rel$type = 'rel'
	
	to_return = rbind(within,between,rel)
	return(to_return)
}

########
#example
########

#use example data from the ez package
library(ez)
data(ANT)

#This approach only makes sense (I think) for 2-level effects (because I use sum contrasts, which I think de-confounds the effect variance from the intercept variance), so toss one level of the flank variable
ANT = ANT[ANT$flank!='Neutral',] #toss the neutral condition 
ANT$flank = factor(ANT$flank)

#do the same for the cue
ANT = ANT[ANT$cue %in% c('None','Center'),] #keep only 2 conditions
ANT$cue = factor(ANT$cue)

#fit the model with sum contrasts
options(contrasts=c('contr.sum','contr.poly'))
fit = lmer(
    formula = rt ~ flank + cue + ( ( flank + cue ) | subnum )
    , data = ANT[ANT$er==0,]
)

get_vars_and_rel( fit )