Builds a model using MCMC

cIRT(
  subject_ids,
  fixed_effects,
  B_elem_plus1,
  rv_effects,
  trial_matrix,
  choices_nk,
  burnit,
  chain_length = 10000L
)

Arguments

subject_ids

A vector that contains subject IDs for each line of data in the choice vector (e.g. For 1 subject that made 5 choices, we would have the number 1 appear five times consecutively.)

fixed_effects

A matrix with NK x P1 dimensions that acts as the design matrix for terms WITHOUT theta.

B_elem_plus1

A V[[1]] dimensional column vector indicating which zeta_i relate to theta_i.

rv_effects

A matrix with NK x V dimensions for random effects design matrix.

trial_matrix

A matrix with N x J dimensions, where J denotes the number of items presented. The matrix MUST contain only 1's and 0's.

choices_nk

A vector with NK length that contains the choice value e.g. 0 or 1.

burnit

An int that describes how many MCMC draws should be discarded.

chain_length

An int that controls how many MCMC draws there are. (> 0)

Value

A list that contains:

as

A matrix of dimension chain_length x J

bs

A matrix of dimension chain_length x J

gs

A matrix of dimension chain_length x P_1

Sigma_zeta_inv

An array of dimension V x V x chain_length

betas

A matrix of dimension chain_length x P_2

Author

Steven Andrew Culpepper and James Joseph Balamuta

Examples

if (FALSE) {
# 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

# Load the Package
library(cIRT)

# Load the Data
data(trial_matrix)
data(choice_matrix)

# Thurstone design matrices
all_nopractice = subset(all_data_trials, experiment_loop.thisN > -1)
hard_items = choice_matrix$hard_q_id
easy_items = choice_matrix$easy_q_id

D_easy = model.matrix( ~ -1 + factor(easy_items))
D_hard = -1 * model.matrix( ~ -1 + factor(hard_items))[, -c(5, 10, 15)]

# Defining effect-coded contrasts
high_contrasts = rbind(-1, diag(4))
rownames(high_contrasts) = 12:16
low_contrasts = rbind(-1, diag(2))
rownames(low_contrasts) = 4:6

# Creating high & low factors
high = factor(choice_matrix[, 'high_value'])
low = factor(choice_matrix[, 'low_value'])
contrasts(high) = high_contrasts
contrasts(low) = low_contrasts

fixed_effects = model.matrix( ~ high + low)
fixed_effects_base = fixed_effects[, 1]
fixed_effects_int = model.matrix( ~ high * low)


# Model with Thurstone D Matrix
system.time({
 out_model_thurstone = cIRT(
   choice_matrix[, 'subject_id'],
   cbind(fixed_effects[, -1], D_easy, D_hard),
   c(1:ncol(fixed_effects)),
   as.matrix(fixed_effects),
   as.matrix(trial_matrix),
   choice_matrix[, 'choose_hard_q'],
   20000,
   25000
 )
})


vlabels_thurstone = colnames(cbind(fixed_effects[, -1], D_easy, D_hard))
G_thurstone = t(apply(
 out_model_thurstone$gs0,
 2,
 FUN = quantile,
 probs = c(.5, .025, .975)
))

rownames(G_thurstone) = vlabels_thurstone
B_thurstone = t(apply(
 out_model_thurstone$beta,
 2,
 FUN = quantile,
 probs = c(.5, 0.025, .975)
))

rownames(B_thurstone) = colnames(fixed_effects)

S_thurstone = solve(
  apply(out_model_thurstone$Sigma_zeta_inv, c(1, 2), FUN = mean)
)

inv_sd = diag(1 / sqrt(diag(solve(
 apply(out_model_thurstone$Sigma_zeta_inv, c(1, 2),
       FUN = mean)
))))

inv_sd %*% S_thurstone %*% inv_sd
apply(out_model_thurstone$as, 2, FUN = mean)
apply(out_model_thurstone$bs, 2, FUN = mean)
}