R/BayesOS.R

In vthorrf/CIRM: Procedures for IRT modeling

BayesOS <- function(x,l,thetain,prec,knot){
  require(R2jags)
  ### Assemble data into list for JAGS====
  y     = x
  Nitem = ncol(y)
  Nsubj = nrow(y)
  L     = (l - 1)
  prec  = tryCatch(prec, error=function(e) apply(y,1,sd)/sqrt(Nitem))
  prec[prec == 0] <- mean(prec[prec != 0])
  TP    <- Nitem * L
  thetain  = tryCatch(thetain, error=function(e) rowMeans(y) )
  if (length(unique(thetain)) >= 10){
  knot <- tryCatch(knot, error=function(e)
           unique(c(0,quantile(thetain, seq(0,1,len=10), type=5),1))) } else {
      knot <- tryCatch(knot, error=function(e)
        unique(c(0,quantile(thetain, seq(0,1,len=3), type=5),1)))
    }
  nk    = length(knot)
  dataList = list( y=y, L=L, Nsubj=Nsubj, Nitem=Nitem,
                   prec=prec, thetain=thetain, TP=TP,
                   nk=nk, knots=knot )

  # 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] <- ilogit( W[i,j] )
        W[i,j] <- sum(Beta[1:nk,j] * Z[i,1:nk])
      }
    }

  ### Theta distribution
    for ( i in 1:Nsubj ) {
      Abil[i] ~ dbeta( (thetain[i]*(1/prec[i])) + 1 ,
                       ((1 - thetain[i])*(1/prec[i])) + 1 )
    }
    psi[1:Nsubj] <- Abil[1:Nsubj] * TP

  ### Rademacher basis
    for (j in 1:Nitem) {
      for (l in 1:nk) {
        Beta[l,j] ~ ddexp( mu , sigma )T(0,1)
      }
    }
    for (l in 1:nk) {
      for (i in 1:Nsubj) {
        Z[i,l] <- ifelse( Abil[i] < knots[l], -1, 1 )
      }
    }

  ### Priors
    mu    ~ dunif( 0, 1 ) # Betas mu
    sigma ~ dunif( 1e-2, 1e2 ) # Betas sd
  }
  " # close quote for modelString
  model = textConnection(modelString)

  # Run the chains====
  # Name the parameters to be monitored
  params <- c("Pr","W","Beta","psi")
  # Random initial values
  inits <- function(){list("Abil" = stats::rbeta(Nsubj,1,1),
                           "Pr"   = stats::rbeta(Nitem,1,1))}
  # 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$psi)
  W <- colMeans(samples$BUGSoutput$sims.list$W[,,])
  diff <- 1 - colMeans(y)
  dic <- samples$BUGSoutput$DIC
  full <- samples
  matrix <- ordering(REs,abil,diff)$matrix
  Result <- list("matrix"=matrix,"abil"=abil,"W"=W,"dic"=dic,"full"=full)
  return(Result)
}