R/npCIRM.R

In vthorrf/CIRM: Procedures for IRT modeling

npCIRM <- function(x,l,p=.5){
  require(R2jags)

  ### Assemble data into list for JAGS====
  y = x
  p = p
  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, p=p, Nsubj=Nsubj, Nitem=Nitem, L=L)

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

  ### Dirichlet distribution prior
    for(i in 1:Nsubj) {
      c[i]      ~ dexp( delta )
      xi_raw[i] ~ dgamma(c[i] + 1, 1)
      pp[i]     <- xi_raw[i]/sum(xi_raw[1:Nsubj])

      g[i]      ~ dexp( alpha )
      gi_raw[i] ~ dgamma(g[i] + 1, 1)
      pg[i]     <- gi_raw[i]/sum(gi_raw[1:Nsubj])
    }

  ### Abil estimates
    for ( i in 1:Nsubj ) {
      subjAbil[i] <- Abil[gi[i]]
      subjAbil[i] ~ dbeta( (mu*(sigma-2)) + 1 , ((1 - mu)*(sigma-2)) + 1  )
      gi[i] ~ dcat(pg[1:Nsubj])
    }

  ### Diff estimates
    for ( j in 1:Nitem ) {
      itemDiff[j] ~ dbeta( (omega*(kappa-2)) + 1 , ((1 - omega)*(kappa-2)) + 1 )
    }

  ### Priors
    mu            ~ dbeta(1,1)
    sigma         <- sigmaMinusTwo + 2
    sigmaMinusTwo ~ dgamma(.01,.01)
    omega         ~ dbeta(1,1)
    kappa         <- kappaMinusTwo + 2
    kappaMinusTwo ~ dgamma(.01,.01)
    delta         ~ dunif( 1e-2, 1e2 )
    alpha         ~ dunif( 1e-2, 1e2 )
  }
  " # close quote for modelString
  model = textConnection(modelString)

  # Run the chains====
  # Name the parameters to be monitored
  params <- c("pCorr","subjAbil","itemDiff","basebeta")
  # Random initial values
  inits <- function(){list("mu"=stats::rbeta(Nsubj,1,1),
                           "sigmaMinusTwo"=stats::rgamma(Nsubj,
                           shape=.01,rate=.01),
                           "omega"=stats::rbeta(Nitem,1,1),
                           "kappaMinusTwo"=stats::rgamma(Nitem,
                           shape=.01,rate=.01),
                           "alpha"=stats::runif(1e-4,1e4),
                           "wrinkleMinusTwo"=stats::rgamma((Nitem * Nsubj),
                              shape=.01,rate=.01) )}
  # Define some MCMC parameters for JAGS
  nthin    = 1    # How Much Thinning?
  nchains  = 3    # How Many Chains?
  nburnin  = 100  # How Many Burn-in Samples?
  nsamples = 1000 # 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====
  ZZZ <- invisible(print(samples))
  REs <- matrix(ZZZ$mean$pCorr,n,v,byrow = T)
  abil <- ZZZ$mean$subjAbil
  easi <- ZZZ$mean$itemDiff*-1
  dic <- ZZZ$DIC
  full <- ZZZ
  matrix <- ordering(REs,abil,easi)$matrix
  Result <- list("matrix"=matrix,"abil"=abil,"easi"=easi,"full"=full)
  return(Result)
}