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")