lmer_oneway_sim

lmer_oneway_sim.R

library(plyr)
library(MASS)
library(lme4a)

generate_data = function(
	n # number of units
	, k # number of trials within each condition within each unit
	, noise # measurement noise variance
	, I # population intercept
	, vI # across-units variance of intercepts
	, A # population A effect
	, vA # across-units variance of A effects
	, rIA # across-units correlation between intercepts and A effects
){
	Sigma = c(
		vI , sqrt(vI*vA)*rIA
		, sqrt(vI*vA)*rIA , vA
	)
	Sigma = matrix(Sigma,2,2)
	means = mvrnorm(n,c(I,A),Sigma)
	temp = expand.grid(A=c('a1','a2'),value=0)
	temp$A = factor(temp$A)
	contrasts(temp$A) = contr.sum
	from_terms = terms(value~A)
	mm = model.matrix(from_terms,temp)
	data = expand.grid(A=c('a1','a2'),unit=1:n,trial=1:k)
	for(i in 1:n){
		data$value[data$unit==i] = as.numeric(mm %*% means[i,]) + rnorm(k*2,0,sqrt(noise))
	}
	data$unit = factor(data$unit)
	data$A = factor(data$A)
	# contrasts(data$A) = contr.sum
	data$A = as.numeric(data$A)-1.5
	return(data)
}

grid = expand.grid(
	iterations = 1:1e3
	, N = c(10,100)
	, k = c(10,100)
	, A = c(0,1)
	, vA = c(0,1)
	, rIA = c(0,.9)
)
grid = grid[grid$A!=0 | (grid$vA==0 & grid$rIA==0),]
grid = grid[grid$vA!=0 | (grid$rIA==0),]

grid = ddply(
	.data = grid
	, .variables = .(iterations,N,k,A,vA,rIA)
	, .fun = function(x){
		set.seed(as.numeric(row.names(x)))
		data = generate_data(
			n = x$N # number of units
			, k = x$k # number of trials within each condition within each unit
			, noise = 1 # measrurement noise variance
			, I = 0 # population intercept
			, vI = 1 # across-units variance of intercepts
			, A = x$A # population A effect
			, vA = x$vA # across-units variance of A effects
			, rIA = x$rIA # across-units correlation between intercepts and A effects
		)
		m0 = log2(exp(1))*AIC(lmer(
			data = data
			, formula = value ~ (1|unit)
			, REML = F
		))
		m1 = log2(exp(1))*AIC(lmer(
			data = data
			, formula = value ~ (1|unit) + A
			, REML = F
		))
		m2 = log2(exp(1))*AIC(lmer(
			data = data
			, formula = value ~ (1|unit) + A + (0+A|unit)
			, REML = F
		))
		m3 = log2(exp(1))*AIC(lmer(
			data = data
			, formula = value ~ (A|unit) + A
			, REML = F
		))
		to_return = data.frame(
			m0 = m0
			, type = c('m1','m2','m3')
			, value = c(m0-m1,m0-m2,m0-m3)
		)
	}
	, .progress = 'text'
)