R/NPT.R

In vthorrf/CIRM: Procedures for IRT modeling

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

  ### Assemble data into list for JAGS====
  y = x
  st = as.vector(scale(rowMeans(y)))
  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 , Nsubj=Nsubj, Nitem=Nitem , L=L, st=st)

  # Define the model====
  modelString = "
  model {
  for ( i in 1:Nsubj ) {
    for ( j in 1:Nitem ) {
      y[i,j] ~ dbin( pCorr[i,j] , L)
      pCorr[i,j] <- ifelse(Pr[pi[i,j],j] >= subjAbil[i], 1e-8, (1-1e-8))
      Pr[i,j] <- itemDisc[j] * (subjAbil[i] - itemDiff[j])
      pi[i,j] ~ dcat(pp[1:Nsubj])
    }
  }

  ### Dirichlet distribution prior
  pp[1:Nsubj] ~ ddirch(alpha[1:Nsubj])

  ### Ability estimates
  for ( subjIdx in 1:Nsubj ) {
    subjAbil[subjIdx] ~ dnorm( TheMean[subjIdx] , sigmaAbil )
    TheMean[subjIdx]  <- (select * st[subjIdx]) + ((1 - select) * muAbil)
    alpha[subjIdx] ~ dgamma(kappa + 1, 1)
  }

  ### Itens estimates
  for ( itemIdx in 1:Nitem ) {
    itemDiff[itemIdx] ~ dnorm( muDiff , sigmaDiff )
    itemDisc[itemIdx] ~ dnorm( muDisc , sigmaDisc )T(0,)
  }

  ### Priors
  select    ~ dbern( sigma )
  sigma     ~ dbeta( 1, 1 )
  muAbil    ~ dnorm( 0 , 1 )
  sigmaAbil ~ dunif( 1e-2, 1e2 )
  muDiff    ~ dnorm( 0 , 1 )
  sigmaDiff ~ dunif( 1e-2, 1e2 )
  muDisc    ~ dnorm( 1 , .1 )T(0,)
  sigmaDisc ~ dunif( 1e-2, 1e2 )
  kappa     ~ dunif( 0, 1e4 )
  }
  " # close quote for modelString
  model = textConnection(modelString)

  # Run the chains====
  # Name the parameters to be monitored
  params <- c("pCorr","subjAbil","itemDiff","itemDisc")
  # Random initial values
  inits <- function(){list("subjAbil"=truncnorm::rtruncnorm(Nsubj,a=-3,b=3,0,1),
                           "itemDiff"=truncnorm::rtruncnorm(Nitem,a=-3,b=3,0,1),
                           "muDisc"=truncnorm::rtruncnorm(Nitem,a=0))}
  # Define some MCMC parameters for JAGS
  nthin    = 2    # How Much Thinning?
  nchains  = 3    # How Many Chains?
  nburnin  = 100  # How Many Burn-in Samples?
  nsamples = 1100 # 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

  ### 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)
  disc <- colMeans(samples$BUGSoutput$sims.list$itemDisc)
  dic <- samples$BUGSoutput$DIC
  full <- samples
  matrix <- ordering(REs,abil,diff)$matrix
  Result <- list("matrix"=matrix,"abil"=abil,"diff"=diff,"disc"=disc,
                 "dic"=dic,"full"=full)
  methods::show(elapsedTime)
  return(Result)
}