R/RCcompare.R

In vthorrf/CIRM: Procedures for IRT modeling

RCcompare <- 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)
      pCorr[i,j] <- (ilogit( AbilR[i] - DiffR[j] ) * (1-m)) +
      (((p * AbilC[i]) /( (p * AbilC[i]) + ((1-p) * DiffC[j]) )) * m)
    }
  }
  for ( subjIdx in 1:Nsubj ) {
    AbilR[subjIdx] ~ dnorm( muAbil , sigmaAbil )
    AbilC[subjIdx] ~ dbeta( (ACIRM*(AsCIRM-2)) + 1 ,
                          ((1 - ACIRM)*(AsCIRM-2)) + 1 )
    }

  for ( itemIdx in 1:Nitem ) {
    DiffC[itemIdx] ~ dbeta( (DCIRM*(DsCIRM-2)) + 1 , ((1 - DCIRM)*(DsCIRM-2)) + 1 )
    DiffR[itemIdx] ~ dnorm( muDiff , sigmaDiff )
  }
  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)

  BF <- pM/(1-pM)
  m ~ dbern(pM)
  pM ~ dbeta( omega*(kappa-2)+1 ,
             (1-omega)*(kappa-2)+1 ) # prior probability for m=1
  omega         ~ dbeta(1,1)
  kappa         <- kappaMinusTwo + 2
  kappaMinusTwo ~ dgamma(.01,.01)
  }
  " # close quote for modelString
  model = textConnection(modelString)

  # Run the chains====
  # Name the parameters to be monitored
  params <- c("BF")
  # 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
  Test <- samples$BUGSoutput$sims.list$BF
  getmode <- function(v) {
    Dens <- density(v,bw="SJ",from=0)
    return(Dens$x[which.max(Dens$y)])
    }
  Mode <- getmode(Test)
  library(coda)
  HDI <- HPDinterval(as.mcmc(Test), prob=0.95)
  Result <- list("BFs"=Test,"HDI"=HDI,"Mode"=Mode)
  return(Result)
}