Simulation-Study-with-cIRT.R
In cIRT: Choice Item Response Theory

## ----sim_setup, eval = F------------------------------------------------------
# ### Variables
# # Y = trial matix
# # C = KN vector of binary choices
# # N = #of subjects
# # J = # of items
# # K = # of choices
# # atrue = true item discriminations
# # btrue = true item locations
# # thetatrue = true thetas/latent performance
# # gamma = fixed effects coefficients
# # Sig = random-effects variance-covariance
# # subid = id variable for subjects
# 
# # install mvtnorm is necessary
# # install.packages("mvtnorm")
# 
# # Load the Library
# library(cIRT)
# 
# # Simulate 2PNO Data
# N = 1000  # Students
# J = 20    # Total numbers of possible items per SA
# 
# # Set a seed for the random generation of a's and b's.
# set.seed(1337)
# 
# # Randomly pick a's and b's
# 
# # Generate as, bs
# atrue=runif(J)+1
# btrue=2*runif(J)-1
# 
# save(atrue, btrue, file="a_b_true.rda")
# 
# 
# # 2 Level Probit Data
# K = 30
# 
# gam_notheta = c(.5,1)
# gam_theta   = c(3,.25)
# gamma = c(gam_notheta,gam_theta)
# 
# Sig = matrix(c(.25,0,0,.125),2,2)
# 
# # Number of replications
# B = 100
# 
# # Begin the Simulation Study
# for(b in 1:B){
# 
#   # Provide user with state information
#   cat(paste0("On Iteration:",b,"\n"))
#   set.seed(b + 1234)
# 
#   # True theta and etay
#   thetatrue = rnorm(N)
#   etay = outer(rep(1,N),atrue) * thetatrue - outer(rep(1,N),btrue)
# 
#   # Generate Y for 2PNO model
#   p.correct = pnorm(etay)
#   Y = matrix(rbinom(N*J, 1, p.correct),N,J)
# 
#   #################################################
#   # Simulating 2 level probit data
#   #################################################
# 
#   subid = expand.grid(cid = 1:K,sid = 1:N)[,2]
# 
#   pred = rnorm(K*N,0,1) # Pred
# 
#   center_pred = center_matrix(as.matrix(pred))
# 
#   Xnotheta = cbind(1,center_pred)
# 
#   Xtheta = rep(thetatrue,each=K)*Xnotheta
#   X = cbind(Xnotheta,Xtheta)
# 
# 
#   zetas = mvtnorm::rmvnorm(N,mean=c(0,0),sigma=Sig) # mvtnorm environment accessed
#   W_veczeta = apply(Xnotheta*zetas[rep(1:N,each=K),],1,sum)
# 
#   etac = X%*%gamma + W_veczeta
#   Zc = rnorm(N*K,mean=etac,sd=1)
#   C = 1*(Zc>0)
# 
#   # Run the Choice Item Response Model
#   out1 = cIRT(subid,
#               Xnotheta,
#               c(1,2),
#               Xnotheta,
#               Y,
#               C,
#               5000)
# 
#   mname = paste0("model_",b)
# 
#   # Assign the data set name
#   assign(mname, out1)
# 
#   # Save the out object
#   save(list=mname, file=paste0(mname,".rda"))
# 
#   # Clean up Export
#   rm(list = c(mname,"out1"))
# }

## ----sim_results, eval = F----------------------------------------------------
# # E[theta] - theta
# bias = function(theta.star, theta.true){
#   matrix(mean(theta.star) - theta.true,ncol=1)
# }
# 
# bias2 = function(theta.star, theta.true){
#   matrix(theta.star - theta.true,ncol=1)
# }
# 
# # sqrt ( 1/n * sum( (y_i - y.hat_i)^2 )
# RMSE = function(y,y.hat){
#   sqrt(  (y-y.hat)^2 )
# }
# 
# # Change true values if needed
# # True Values
# gam_notheta = c(.5,1)
# gam_theta   = c(3,.25)
# gamma = c(gam_notheta,gam_theta)
# # Loads a and b values
# load("a_b_true.rda")
# 
# Sig = as.numeric(matrix(c(.25,0,0,.125),2,2))[c(1,2,4)]
# B = 100
# 
# # Storage to hold bootstrap replications
# a_result = matrix(0, B, 20)
# 
# b_result = matrix(0, B, 20)
# 
# gs0_result = matrix(0, B, 2)
# 
# beta_result = matrix(0, B, 2)
# 
# sig_result = array(NA, dim=c(2,2,B))
# 
# 
# for(b in 1:B){
# 
#   mname = paste0("model_",b)
# 
#   load(paste0(mname, ".rda"))
# 
#   d = get(mname)
# 
#   a_result[i,] = apply(d$as, 2, FUN = mean)
#   b_result[i,] = apply(d$bs, 2, FUN = mean)
# 
#   gs0_result[i,] = apply(d$gs0, 2, FUN = mean)
# 
#   beta_result[i,] = apply(d$betas, 2, FUN = mean)
# 
#   sig_result[,,i] = solve(apply(d$Sigma_zeta_inv, c(1,2), FUN = mean))
# }
# 
# # Obtain an overall mean for each of the following:
# 
# m_a_result = apply(a_result, 2, FUN = mean)
# 
# m_b_result = apply(b_result, 2, FUN = mean)
# 
# m_gs0_result = apply(gs0_result, 2, FUN = mean)
# 
# m_beta_result = apply(beta_result, 2, FUN = mean)
# 
# m_sig_result = as.numeric(apply(sig_result, c(1,2), FUN = mean))[c(1,2,4)]
# 
# 
# # Perform a bias evaluation given the true values:
# 
# a_bias = bias2(m_a_result,atrue)
# 
# b_bias = bias2(m_b_result,btrue)
# 
# gs0_bias = bias2(m_gs0_result,gamma[1:2])
# 
# beta_bias = bias2(m_beta_result,gamma[3:4])
# 
# sig_bias = bias2(m_sig_result, Sig)
# 
# # Perform the RMSE under the supplied results:
# 
# a_RMSE = RMSE(m_a_result,atrue)
# 
# b_RMSE = RMSE(m_b_result,btrue)
# 
# gs0_RMSE = RMSE(m_gs0_result,gamma[1:2])
# 
# beta_RMSE = RMSE(m_beta_result,gamma[3:4])
# 
# sig_RMSE = RMSE(m_sig_result, Sig)

## ----out, eval = F------------------------------------------------------------
# # Make a results export
# results_a = cbind(m_a_result, atrue, a_bias, a_RMSE)
# 
# rownames(results_a) = paste0("Item ",1:length(atrue))
# 
# colnames(results_a) = c("$a$ Estimate", "$a$ True", "$a$ Bias", "$a$ RMSE")
# 
# results_b = cbind(m_b_result,  btrue, b_bias, b_RMSE)
# 
# rownames(results_b) = paste0("Item ",1:length(btrue))
# colnames(results_b) = c("$b$ Estimate", "$b$ True", "$b$ Bias", "$b$ RMSE")
# 
# results_gs0 = cbind(m_gs0_result,gamma[1:2],gs0_bias,gs0_RMSE)
# rownames(results_gs0) = paste0("$\\gamma_",1:2,"$")
# colnames(results_gs0) = c("$\\beta$ Estimate", "$\\beta$ True", "$\\beta$ Bias", "$\\beta$ RMSE")
# 
# results_beta = cbind(m_beta_result,gamma[3:4],beta_bias,beta_RMSE)
# rownames(results_beta) = paste0("$\\beta_",1:2,"$")
# colnames(results_beta) = c("$\\beta$ Estimate", "$\\beta$ True", "$\\beta$ Bias", "$\\beta$ RMSE")
# 
# results_sig_zeta = cbind(m_sig_result,Sig,sig_bias,sig_RMSE)
# rownames(results_sig_zeta) = c("$\\Zeta_{1,1}$","$\\Zeta_{2,1} = $\\Zeta_{1,2}$", "$\\Zeta_{2,2}$")
# colnames(results_sig_zeta) = c("$\\Sigma_{\\Zeta}$ Estimate", "$\\beta$ True", "$\\beta$ Bias", "$\\beta$ RMSE")
# 
# save(results_a, results_b, results_gs0, results_beta, results_sig_zeta, file="sim_results.rda")