Hello Dr. Valpine,
This is an example code of what I want to accomplish. The code will fail because of the "which" function inside nimble code. The below code is similar to a multidimensional item response theory model but with the latent skills being discrete variable.
library(CDM)
library(nimble)
data( data.fraction1 )
X = data.fraction1$data
Q = data.fraction1$q.matrix
fullModel <- nimbleCode({
#coeff. for the skills
for(k in 1:K){
m_A[k] ~ dbeta(1,1)
}
# Prior for item difficulty
for(j in 1:J){
D[j] ~ dbeta(1,1)
}
for (i in 1:N){
for(k in 1:K){
lambda[i,k] ~ dbeta(1,1)
skills[i,k] ~ dbern(lambda[i,k])
}
}
for(i in 1:N){
for(j in 1:J){
# Modeling item response probability ========
# Logic to compute effective ability ========
q_j[1:K] <- which(Q[j,1:K] == 1)
eff_ability[i,j] <- pmin(m_A[k]*skills[i,q_j])
# =========================================
theta[i,j] <- eff_ability[i,j] - D[j]
pi[i,j] <- ilogit(1.7*theta[i,j])
# Item-Response Likelihood
X[i,j] ~ dbern(pi[i,j])
# ===============================================
}
}
})
modelConstants <- list(
N = dim(X)[1],
J = dim(X)[2],
K = dim(Q)[2]
)
data <- list(
X = X,
Q = Q
)
initial_values <- NULL
M0 <- nimbleModel(
code = fullModel,
data = data,
inits = initial_values,
constants = modelConstants,
name = "baseline model"
)
dataNodes <- M0$getNodeNames(dataOnly = TRUE)
parentNodes <- M0$getParents(dataNodes, stochOnly = TRUE)
simNodes <- M0$getDependencies(parentNodes, self = FALSE)
n.iter <- 15000
n.burnun <- 2000
n.chain <- 3
params.to.save <- parentNodes
mcmc_fit<- nimbleMCMC( model = cmodel,
monitors = params.to.save,
niter = n.iter,
nburnin = n.burnun,
nchains = n.chain,
samplesAsCodaMCMC=TRUE, WAIC=TRUE)
I am also attaching the code as an attachment.