R/BCM.R

In vthorrf/CIRM: Procedures for IRT modeling

BCM <- function(x,l){
  require(R2jags)

  ### Assemble data into list for JAGS====
  y = x/(l-1)
  x = x
  itemID = colnames(y)
  subjID = rownames(y)
  Nitem = ncol(y)
  Nsubj = nrow(y)
  v = ncol(y)
  n = nrow(y)
  L = (l - 1)
  dataList = list( y=y, x=x, Nsubj=Nsubj, Nitem=Nitem, L=L)

  # Define the model====
  modelString = "
  model {
  for ( i in 1:Nsubj ) {
    for ( j in 1:Nitem ) {
      ypred[i,j] <- yM[i,j] * L
      x[i,j] ~ dbinom( y[i,j], L)
      y[i,j] ~ dbeta( (yM[i,j]*(yS-2)) + 1 , ((1 - yM[i,j])*(yS-2)) + 1  )
      yM[i,j] <- exp( log(subjAbil[i]) + log(itemDiff[j]) )
    }
  }
  for ( subjIdx in 1:Nsubj ) {
    subjAbil[subjIdx] ~ dbeta( (mu*(sigma-2)) + 1 , ((1 - mu)*(sigma-2)) + 1  )
  }
  for ( itemIdx in 1:Nitem ) {
    itemDiff[itemIdx] ~ dbeta( (omega*(kappa-2)) + 1 , ((1 - omega)*(kappa-2)) + 1 )
  }

  yS            <- ySMinusTwo + 2
  ySMinusTwo    ~ dunif(1e-4,1e4)
  mu            ~ dbeta(1,1)
  sigma         <- sigmaMinusTwo + 2
  sigmaMinusTwo ~ dunif(1e-4,1e4)
  omega         ~ dbeta(1,1)
  kappa         <- kappaMinusTwo + 2
  kappaMinusTwo ~ dunif(1e-4,1e4)
  }
  " # close quote for modelString
  model = textConnection(modelString)

  # Run the chains====
  # Name the parameters to be monitored
  params <- c("subjAbil","itemDiff","yM","yS","ypred")
  # Random initial values
  inits <- function(){list("subjAbil"=stats::rbeta(Nsubj,1,1),
                           "itemDiff"=stats::rbeta(Nitem,1,1))}
  # Define some MCMC parameters for JAGS
  nthin    = 1    # How Much Thinning?
  nchains  = 3    # How Many Chains?
  nburnin  = 10  # How Many Burn-in Samples?
  nsamples = 110 # How Many Recorded Samples?
  ### Calling JAGS to sample
  startTime = proc.time()
  samples <- R2jags::jags(dataList, NULL, params, model.file =model,
                          n.chains=nchains, n.iter=nsamples, n.burnin=nburnin,
                          n.thin=nthin, DIC=T, jags.seed=666)
  stopTime = proc.time(); elapsedTime = stopTime - startTime; methods::show(elapsedTime)

  ### Inspect and diagnose the run====
  REs <- colMeans(samples$BUGSoutput$sims.list$pCorr[,,])
  abil <- colMeans(samples$BUGSoutput$sims.list$subjAbil)
  diff <- colMeans(samples$BUGSoutput$sims.list$itemDiff)
  dic <- samples$BUGSoutput$DIC
  rmsd <- mean(samples$BUGSoutput$sims.list$MAE)
  full <- samples
  matrix <- ordering(REs,abil,diff)$matrix
  Result <- list("matrix"=matrix,"abil"=abil,"diff"=diff,"dic"=dic,"full"=full)
  return(Result)
}