R/BRPM.DIF.R

In changxiulee/BayesianRasch: Bayesian Rasch Modeling

#' BRPM.DIF
#'
#' This function assess the differential item functioning in your data
#' 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 model is identifiable, ability is assumed to be normally distributed. (theta ~ N(0, 1))
#' 3) To ensure the RSM model is identifiable, ability is assumed to be normally distributed. (theta ~ N(0, 1)) and
#'    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.
#' 4) I follow Winsteps anchor theta method. The first analysis is to get the mean theta for each group. The mean is
#'    then used to anchor the theta in the next analysis
#'
#'
#' Reference
#' Linacre, J. M., & Wright, B. D. (2000). Winsteps. URL: http://www. winsteps. com/index.html.
#' Soares, T. M., Gonçalves, F. B., & Gamerman, D. (2009). An integrated Bayesian model for DIF analysis. Journal of Educational and Behavioral Statistics, 34(3), 348-377.
#'
#' @param data A data frame or matrix
#' @param item Item to be included in your model
#' @param DIF.var grouping variable. Group should be in numeric
#' @param n.chains Number of chains to be included in MCMC
#' @param model Choose between 'BRSM' (Bayesian Rating Scale Model) and 'BPCM' (Bayesian Partial Credit Model).
#' @param g group to be include in analysis
#' @param ROPE region of practical equivalence. DIF for each item won't be equal to 0. But not all non-zero DIF is practical significant. Set a region that has practically insignificant DIF
#' @export BRPM.DIF

BRPM.DIF = function(data, item, DIF.var, n.chains, model = "BRSM", g = c(1,2), ROPE = c(-0.5, 0.5)){
  #Initialization
  item.variable = data[, item]
  if(min(apply(item.variable, 2, min,na.rm = TRUE)) == 0){item.variable = item.variable + 1}
  G = max(data[,DIF.var])
  N  = table(data[,DIF.var])
  id = c()
  for(i in 1:G){id = c(id, seq(1, N[i], 1))}
  group = cbind(data[, DIF.var], id)
  K = apply(item.variable, 2, max,na.rm = TRUE)
  p = ncol(item.variable)
  model.name = model

  #Main 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 = 2000 / n.chains, method = "rjparallel")
  }
  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"
  }

  #Find group differences to anchor
  group.mean = c()
  group.sd = c()
  for(i in 1:G){
    group.mean = c(group.mean, mean(theta(obj, g = i)[,3]))
    sd.v = c()
    for(j in 1:5){sd.v = c(sd.v, sd(mcmc.matrix[j, paste0("theta[", i,",", 1:N[i],"]")]))}
    group.sd = c(group.sd, mean(sd.v))
  }

  #DIF Analysis
  if(model.name == "BPCM"){
    model = runjags::run.jags(model = system.file("JAGS","PCM.DIF.txt", package = "BayesianRasch" ), monitor = c("beta", "theta", "d"),
                     data = list(Y = as.matrix(item.variable), N = N, group = group,
                                 G = G, K = K, p = p, mu = group.mean, sd = group.sd), n.chains = n.chains,
                     burnin = 0, adapt = 1000, sample = 2000 / n.chains, method = "rjparallel")
  }else{
    model = runjags::run.jags(model = system.file("JAGS","RSM.DIF.txt", package = "BayesianRasch" ), monitor = c("beta", "theta", "thres", "d"),
                     data = list(Y = as.matrix(item.variable), N = N, group = group,
                                 G = G, K = max(K), p = p, mu = group.mean, sd = group.sd), n.chains = n.chains,
                     burnin = 0, adapt = 1000, sample = 2000 / n.chains, method = "rjparallel")
  }

  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]][,])}
  mat = NULL

  #DIF result
  if(model.name == "BPCM"){
    for(j in 1:p){
      for(k in 2:max(K)){
        mat = rbind(mat, c(HDInterval::hdi(mcmc.matrix[,paste0("d[", 2,",", j,",",k,"]")]), mean(mcmc.matrix[,paste0("d[", 2,",", j,",",k,"]")])))
      }
    }
    success.v = c()
    for(j in 1:p){
      success = 0
      for(i in 1:nrow(mcmc.matrix)){
        if (mcmc.matrix[i, paste0("d[", 2,",", j,",",k,"]")] < ROPE[1] | mcmc.matrix[i, paste0("d[", 2,",", j,",",k,"]")] > ROPE[2]){
          success = success + 1
        }
      }
      success.v = c(success.v, success / nrow(mcmc.matrix))
    }
    mat = cbind(mat, success.v)
  }else {
    for(j in 1:p){mat = rbind(mat, c(HDInterval::hdi(mcmc.matrix[,paste0("d[",2,",", j,"]")]), mean(mcmc.matrix[,paste0("d[", 2,",", j,"]")])))}
    success.v = c()

    for(j in 1:p){
      success = 0
      for(i in 1:nrow(mcmc.matrix)){
        if (mcmc.matrix[i, paste0("d[", 2,",", j,"]")] < ROPE[1] | mcmc.matrix[i, paste0("d[", 2,",", j,"]")] > ROPE[2]){
          success = success + 1
        }
      }
      success.v = c(success.v, success / nrow(mcmc.matrix))
    }
    mat = cbind(mat, success.v)
  }
  colnames(mat) = c("lower", "upper", "mean", "Prob.DIF")

  #Make an BRPM object
  if(model.name == "BPCM"){
    mcmc.new = BPCM.DIF.transform(mcmc.matrix, item, p, N, max(K), G)
  }else{
    mcmc.new = BRSM.DIF.transform(mcmc.matrix, item, p, N, max(K), G)
  }

  obj = list(result = list(mat = mat, group.mean = group.mean),
             mcmc = mcmc.new, data = data, item = item, N = N, DIF = TRUE, DIF.var = DIF.var)
  if(model.name == "BPCM"){
    class(obj) <- "BPCM"
  }else{
    class(obj) <- "BRSM"
  }
  obj
}