R/BetaNP.R

In vthorrf/CIRM: Procedures for IRT modeling

BetaNP <- function(x,l){
  require(R2jags)
  ### Assemble data into list for JAGS====
  y     = x
  Nitem = ncol(y)
  Nsubj = nrow(y)
  L     = (l - 1)
  dataList = list( y=y, L=L, Nsubj=Nsubj, Nitem=Nitem )

  # Define the model====
  modelString = "
  model {
    for ( i in 1:Nsubj ) {
      for ( j in 1:Nitem ) {
        y[i,j] ~ dbin( Pr[i,j] , L)
        Pr[i,j] ~ dbeta(Abil[i] + 1, Diff[j] + 1)
      }
    }

  ### Theta distribution
    for ( i in 1:Nsubj ) {
      Abil[i] ~ dgamma( aT , bT )
    }

  ### Delta distribution
    for ( j in 1:Nitem ) {
      Diff[j] ~ dgamma( aD , bD )
    }

  ### Priors
    aT ~ dgamma(1e-2,1e-2) # Ability shape
    bT ~ dgamma(1e-2,1e-2) # Ability rate
    aD ~ dgamma(1e-2,1e-2) # Difficulty shape
    bD ~ dgamma(1e-2,1e-2) # Difficulty rate
  }
  " # close quote for modelString
  model = textConnection(modelString)

  # Run the chains====
  # Name the parameters to be monitored
  params <- c("Pr","Abil","Diff")
  # Random initial values
  inits <- function(){list("Abil" = stats::rgamma(Nsubj,1e-2,1e-2),
                           "Diff"   = stats::rgamma(Nitem,1e-2,1e-2))}
  # Define some MCMC parameters for JAGS
  nthin    = 1    # 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; methods::show(elapsedTime)

  ### Gathering====
  REs <- colMeans(samples$BUGSoutput$sims.list$Pr[,,])
  abil <- colMeans(samples$BUGSoutput$sims.list$Abil)
  diff <- colMeans(samples$BUGSoutput$sims.list$Diff)
  dic <- samples$BUGSoutput$DIC
  full <- samples
  matrix <- ordering(REs,abil,diff)$matrix
  Result <- list("matrix"=matrix,"abil"=abil,"diff"=diff,"dic"=dic,"full"=full)
  return(Result)
}