R/coherence.R

In changxiulee/BayesianRasch: Bayesian Rasch Modeling

#' coherence
#'
#' For response category analysis, assess how well the category predict the measure and how the measure
#' predict the category
#'
#' Reference
#' Linacre, J. M. (2002). Optimizing rating scale category effectiveness. Journal of applied measurement, 3(1), 85-106.
#' Linacre, J. M., & Wright, B. D. (2000). Winsteps. URL: http://www. winsteps. com/index. html
#'
#' @param obj An BRSM model
#' @export coherence
#'

coherence = function (obj) {
  if(obj$DIF == TRUE){
    show("This is not available for DIF object")
  }else{
    codaSamples = obj$mcmc
    data = obj$data
    if(min(apply(data, 2, min,na.rm = TRUE)) == 0){data = data + 1}
    item = obj$item
    n.item = length(item)
    data = data[, item]
    P = ncol(data)
    N = nrow(data)
    K = apply(data, 2, max,na.rm = TRUE)
    m = matrix(NA, 2, max(K))

    if(class(obj) == "BPCM"){beta.matrix = beta.specific(codaSamples, data, item)}
    else{
      beta.matrix = beta(obj)[,3]
      thresh = obj$thresh
    }
    EXM.all = NULL
    for(i in 1:P){
      if(class(obj) == "BPCM"){
        if(max(K) == 2){
          itemdiff = beta.matrix[i]
        }else{
          itemdiff = beta.matrix[i,]
        }
      }
      else{
        itemdiff = beta.matrix[i]
      }
      EWX = matrix(NA, N, 3)
      EWX[,3] = data[,i] - 1
      for(j in 1:N){
        E = 0
        w = 0
        person = mean(codaSamples[,paste0("theta[",1,",",j,"]")])
        prob = probgenerator(itemdiff, person, obj)
        for(k in 1:length(prob)){E = E + prob[k] * (k-1)}
        for(k in 1:length(prob)){w = w + (((k-1) - E) ^ 2 ) * prob[k] }
        EWX[j,1:2] = c(E,w)
      }
      EXM = round(EWX[,c(1,3)])
      EXM = cbind(EXM, EXM[,1] == EXM[,2])
      EXM.all = rbind(EXM.all, EXM)
    }

    for(f in 1:max(K)){

      mc = subset(EXM.all, EXM.all[,1] == (f-1))
      m[f,1] = mean(mc[,3]) #M -> C
      cm = subset(EXM.all, EXM.all[,2] == (f-1))
      m[f,2] = mean(cm[,3]) #C -> M
    }
    colnames(m) = c("M --> C", "C --> M")
    rownames(m) = c(1:max(K))
    m
  }
}