R/BRPM.R

In changxiulee/BayesianRasch: Bayesian Rasch Modeling

#' Bayesian Rasch Polytomous Model
#'
#' This function analyse your data using a polytomous Rasch model estimated with Bayes. Please note:
#' 1) When your data has only two categories, the polytomous Rasch model will fall back to the
#'    dichotomous Rasch model. Therefore, both choice of model is equivalent.
#' 2) To ensure the PCM and RSM model is identifiable, ability is assumed to be normally distributed (theta ~ N(0, 1)).
#' 3) In RSM model, my intention is to make the andrich threshold follow the sum-to-zero constraint. However, I could not find a
#'    way to do it in JAGS. Therefore, the first andrich threshold is assumed to be 0 in JAGS. Then, the andrich threshold and
#'    beta is transformed to allow the andrich threshold to follow the sum-to-zero constraint
#'
#' Wright, B. D., & Masters, G. N. (1982). Rating scale analysis. MESA press.
#' Kruschke, J. (2014). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press.
#'
#' @param data A data frame with your data
#' @param item Item in the data to be included in the model
#' @param model Choose between 'BRSM' (Bayesian Rating Scale Model) and 'BPCM' (Bayesian Partial Credit Model). Default is set to BRSM.
#' @param n.chains Number of chains for Markov Chain Monte Carlo. If you have multicores computer, you may increase the chain
#' @export BRPM

BRPM = function(data, item, n.chains, model = "BRSM"){

  #Initialization
  item.variable = data[, item]
  if(min(apply(item.variable, 2, min,na.rm = TRUE)) == 0){item.variable = item.variable + 1}
  p = ncol(item.variable)
  K = apply(item.variable, 2, max,na.rm = TRUE)
  G = 1
  N  = c(nrow(item.variable))
  group = cbind(rep(1, N), seq(1, N, 1))
  model.name = model
  #Analysis
  if(model.name == "BPCM"){
    model = runjags::run.jags(model = system.file("JAGS","PCM.txt", package = "BayesianRasch" ), monitor = c("beta", "theta"),
                     data = list(Y = as.matrix(item.variable), N = N, group = group,
                                 G = G, K = K, p = p), n.chains = n.chains,
                     burnin = 0, adapt = 1000, sample = 2000 / n.chains, method = "rjparallel")
  }
  else{
    model = runjags::run.jags(model =  system.file("JAGS","RSM.txt", package = "BayesianRasch" ), monitor = c("beta", "theta", "thres"),
                     data = list(Y = as.matrix(item.variable), N = N, group = group,
                                 G = G, K = max(K), p = p), n.chains = n.chains,
                     burnin = 0, adapt = 1000, sample = floor(2000 / n.chains),  method = "rjparallel")
  }

  #Make BRPM object
  mcmc = coda::as.mcmc.list(model)
  n.chains = length(mcmc)
  mcmc.matrix = NULL
  for(i in 1:n.chains){mcmc.matrix = rbind(mcmc.matrix, mcmc[[i]][,])}
  if(model.name == "BRSM"){
    mcmc.matrix = BRSM.transform(mcmc.matrix, item, max(K))
  }
  obj = list(mcmc = mcmc.matrix, data = data, item = item, N = N, DIF = FALSE)
  if(model.name == "BPCM"){
    class(obj) <- "BPCM"
  }
  else{
    class(obj) <- "BRSM"
  }
  obj
}