R/RCOcompare.R

In vthorrf/CIRM: Procedures for IRT modeling

RCOcompare <- function(x,l,p=.5){
  require(R2jags)
  ### Assemble data into list for JAGS====
  y = x
  p = p
  Nitem = ncol(y)
  Nsubj = nrow(y)
  v = ncol(y)
  n = nrow(y)
  L = (l - 1)
  delta <- rep(1,3)
  prec  = tryCatch(prec, error=function(e) apply(y,1,sd)/sqrt(Nitem))
  prec[prec == 0] <- mean(prec[prec != 0])
  thetain  = tryCatch(thetain, error=function(e) rowMeans(y) )
  knot <- tryCatch(knot, error=function(e)
    unique(c(0,quantile(thetain, seq(0,1,len=length(thetain)), type=5),1)))
  nk    = length(knot)
  dataList = list( y=y, L=L, Nsubj=Nsubj, Nitem=Nitem, prec=prec,
                   thetain=thetain, p=p, nk=nk, knots=knot, delta=delta )

  # 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(M == 1, ilogit( AbilR[i] - DiffR[j] ),
         ifelse(M == 2, ((p * AbilC[i]) /( (p * AbilC[i]) + ((1-p) * DiffC[j]) )),
          ilogit( W[i,j] ) ))
      W[i,j]  <- sum(Beta[1:nk,j] * Z[i,1:nk])
      #pi[i,j] ~ dcat(pp[1:Nsubj])
    }
  }

  ### Bayes Factor
  BFrasch <- alpha[1]/alpha[1]
  BFcirm  <- alpha[2]/alpha[1]
  BFosirm <- alpha[3]/alpha[1]

  ### Model selection
  M ~ dcat(alpha[1:3])
  for(d in 1:3) {
    c[d]      ~ dexp( kappa )
    xi_raw[d] ~ dgamma(c[d] + 1, 1)
    alpha[d]  <- xi_raw[d]/sum(xi_raw[1:3])
  }

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

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

  ### Ability estimates
  for ( i in 1:Nsubj ) {
    AbilR[i] ~ dnorm( muAbil , sigmaAbil )
    AbilC[i] ~ dbeta( (ACIRM*(AsCIRM-2)) + 1 ,
                  ((1 - ACIRM)*(AsCIRM-2)) + 1 )
    Abil[i] ~ dbeta( (thetain[i]*(1/prec[i])) + 1 ,
                  ((1 - thetain[i])*(1/prec[i])) + 1 )
  }

  ### Difficulty estimates
  for ( j in 1:Nitem ) {
    DiffR[j] ~ dnorm( muDiff , sigmaDiff )
    DiffC[j] ~ dbeta( (DCIRM*(DsCIRM-2)) + 1 , ((1 - DCIRM)*(DsCIRM-2)) + 1 )
  }

  ### Endless priors
  ACIRM  ~ dbeta(1,1)
  AsCIRM <- AsmTwo + 2
  AsmTwo ~ dgamma(.01,.01)
  DCIRM  ~ dbeta(1,1)
  DsCIRM <- DsmTwo + 2
  DsmTwo ~ dgamma(.01,.01)

  muAbil    ~ dnorm(0,.001)
  sigmaAbil ~ dgamma(.01,.01)
  muDiff    ~ dnorm(0,.001)
  sigmaDiff ~ dgamma(.01,.01)

  #delta ~ dunif( 1e-2, 1e2 )
  kappa ~ dunif( 1e-2, 1e2 )
  betamu ~ dunif( 0, 1 )
  betasd  ~ dunif( 1e-2, 1e2 )

  }
  " # close quote for modelString
  model = textConnection(modelString)

  # Run the chains====
  # Name the parameters to be monitored
  params <- c("alpha","M","BFrasch","BFcirm","BFosirm")
  # Random initial values
  inits <- function(){list("omega"=stats::rbeta(Nsubj,1,1),
                           "kappaMinusTwo"=stats::rgamma(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)

  ### Inspect and diagnose the run====
  PM <- matrix(samples$BUGSoutput$median$alpha,ncol=3)
  colnames(PM) <- c("Rasch","CIRM","BayesOS")
  Result <- list("M"=samples$BUGSoutput$median$M,
                 "PM"=PM,
                 "dic"=samples$BUGSoutput$DIC)
  return(Result)
}