d-prime reliability simulation

dprime_sim.R

library(plyr)
library(data.table) #for rbindlist
library(MASS) #for mvrnorm
N = 110
P = .25
K = 440
intercept = -0.7793
effect = 1.0772
intercept_sd = 0.3757
effect_sd = 0.2081
ie_corr = .22

generate_data = function(
	n # number of units
	, k # number of trials within each condition within each unit
	, 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
	, means = NULL
){
	if(is.null(means)){
		require(MASS)
		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] = rbinom(k*2,1,pnorm(as.numeric(mm %*% means[i,])))
		data$true_intercept[data$unit==i] = means[i,1]
		data$true_effect[data$unit==i] = means[i,2]
	}
	data$unit = factor(data$unit)
	data$A = factor(data$A)
	contrasts(data$A) = contr.sum
	# return(data)
	return(list(data=data,means=means))
}



out = generate_data(
	n = N
	, k  = K
	, I = intercept
	, vI = intercept_sd
	, A = effect
	, vA = effect_sd
	, rIA = ie_corr
)
data = out$data
data = data[!(data$A=='a1' & data$trial>(K*P)),]
mean(data$value[data$A=='a1'])
mean(data$value[data$A=='a2'])

fit = glmer(
	data = data
	, formula = value ~ ( 1 + A | unit ) + A
	, family = binomial(link='probit')
)
fit
plot(means[,2],ranef(fit)$unit$A1)

out = llply(
	.data = 1:1e2
	, .fun = function(iteration){
		out1 = generate_data(
			n = N
			, k  = K
			, I = intercept
			, vI = intercept_sd
			, A = effect
			, vA = effect_sd
			, rIA = ie_corr
		)
		data1 = out1$data
		data1 = data1[!(data1$A=='a1' & data1$trial>(K*P)),]
		out2 = generate_data(
			n = N
			, k  = K
			, I = intercept
			, vI = intercept_sd
			, A = effect
			, vA = effect_sd
			, rIA = ie_corr
			, means = out1$means
		)
		data2 = out2$data
		data2 = data2[!(data2$A=='a1' & data2$trial>(K*P)),]
		data1$session = 1
		data2$session = 2
		both = rbind(data1,data2)
		session_scores = dlply(
			.data = both
			, .variables = .(unit,session)
			, .fun = function(x){
				this_fit = glm(
					data = x
					, formula = value ~ A
					, family = binomial(link='probit')
				)
				phit = mean(x$value[(x$A=='a1')])
				pfa = mean(x$value[(x$A=='a2')])
				needs_correction = ifelse(
					(phit==0)|(phit==1)|(pfa==0)|(pfa==1)
					, TRUE
					, FALSE
				)
				if(phit==0){
					phit = 1/(nrow(x[(x$A=='a1'),])*2)
				}
				if(phit==1){
					phit = 1 - 1/(nrow(x[(x$A=='a1'),])*2)
				}
				if(pfa==0){
					pfa = 1/(nrow(x[(x$A=='a2'),])*2)
				}
				if(pfa==1){
					pfa = 1 - 1/(nrow(x[(x$A=='a2'),])*2)
				}
				to_return = data.frame(
					unit = x$unit[1]
					, session = x$session[1]
					, true_intercept = x$true_intercept[1]
					, true_effect = x$true_effect[1]
					, bias = coef(this_fit)[1]-plogis(P)
					, discrim = coef(this_fit)[2]
					, converged = this_fit$converged
					, needs_correction = needs_correction
					, dprime = qnorm(phit)-qnorm(pfa)
				)
				return(to_return)
			}
		)
		session_scores = rbindlist(session_scores)
		half_scores = dlply(
			.data = data1
			, .variables = .(unit)
			, .fun = function(z){
				z = dlply(
					.data = z
					, .variables = .(A)
					, .fun = function(x){
						x$half = ((1:nrow(x))%%2)[order(rnorm(nrow(x)))]+1
						return(x)
					}
				)
				z = rbindlist(z)
				to_return = dlply(
					.data = z
					, .variables = .(half)
					, .fun = function(x){
						this_fit = glm(
							data = x
							, formula = value ~ A
							, family = binomial(link='probit')
						)
						phit = mean(x$value[(x$A=='a1')])
						pfa = mean(x$value[(x$A=='a2')])
						needs_correction = ifelse(
							(phit==0)|(phit==1)|(pfa==0)|(pfa==1)
							, TRUE
							, FALSE
						)
						if(phit==0){
							phit = 1/(nrow(x[(x$A=='a1'),])*2)
						}
						if(phit==1){
							phit = 1 - 1/(nrow(x[(x$A=='a1'),])*2)
						}
						if(pfa==0){
							pfa = 1/(nrow(x[(x$A=='a2'),])*2)
						}
						if(pfa==1){
							pfa = 1 - 1/(nrow(x[(x$A=='a2'),])*2)
						}
						to_return = data.frame(
							unit = x$unit[1]
							, half = x$half[1]
							, true_intercept = x$true_intercept[1]
							, true_effect = x$true_effect[1]
							, bias = coef(this_fit)[1]-plogis(P)
							, discrim = coef(this_fit)[2]
							, converged = this_fit$converged
							, needs_correction = needs_correction
							, dprime = qnorm(phit)-qnorm(pfa)
						)
						return(to_return)				
					}
				)
				return(rbindlist(to_return))
			}
		)
		half_scores = rbindlist(half_scores)
		to_return = data.frame(
			iteration = iteration
			, dprime_session1 = with(
				session_scores
				, cor(dprime[session==1],dprime[session==2])
			)
			, dprime_session2 = with(
				session_scores[!(session_scores$unit%in% session_scores$unit[session_scores$needs_correction]),]
				, cor(dprime[session==1],dprime[session==2])
			)
			, discrim_session = with(
				session_scores[!(session_scores$unit%in% session_scores$unit[session_scores$needs_correction]),]
				, cor(discrim[session==1],discrim[session==2])
			)
			, dprime_half1 = with(
				half_scores
				, cor(dprime[half==1],dprime[half==2])
			)
			, dprime_half2 = with(
				half_scores[!(half_scores$unit%in% half_scores$unit[half_scores$needs_correction]),]
				, cor(dprime[half==1],dprime[half==2])
			)
			, discrim_half = with(
				half_scores[!(half_scores$unit%in% half_scores$unit[half_scores$needs_correction]),]
				, cor(discrim[half==1],discrim[half==2])
			)
		)
		return(to_return)
	}
	, .progress = 'time'
)
out = rbindlist(out)

summary(out)