R/item.fit.R

In changxiulee/BayesianRasch: Bayesian Rasch Modeling

#' item.fit
#'
#' Compute the item fit statistics for BRPM model
#'
#' Reference
#' Wright, B. D., & Masters, G. N. (1982). Rating scale analysis. MESA press.
#'
#' @param obj An BPCM model or BRSM model
#' @param g If BPCM / BRSM has DIF, specify group of interest
#' @export item.fit

item.fit = function(obj, g = 1) {
  item = obj$item
  data = obj$data
  DIF = obj$DIF
  if(DIF == TRUE){
    DIF.var = obj$DIF.var
    codaSamples = obj$mcmc[[g]]
    data = data[data[,DIF.var] == g,]
  }else{
    codaSamples = obj$mcmc
  }
  data = data[, item]
  if(min(apply(data, 2, min,na.rm = TRUE)) == 0){data = data + 1}
  N = obj$N
  n.item = length(item)
  P = ncol(data)
  K = apply(data, 2, max,na.rm = TRUE)
  item.infit = c()
  item.outfit = c()

  if(class(obj) == "BPCM"){
    beta.matrix = beta.specific(codaSamples, data, item)
  }else{
    beta.matrix = beta(obj)
  }

  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,3]
    }

    EWX = matrix(NA, N[g], 3)
    EWX[,3] = data[,i] - 1
    for(j in 1:N[g]){
      E = 0
      w = 0
      person = mean(codaSamples[,paste0("theta[",g,",",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)
    }
    z = (EWX[,3] - EWX[,1]) / sqrt(EWX[,2])

    #Outfit
    z.sum = sum(z ^ 2, na.rm = TRUE)
    outfit = z.sum / N[g]

    #Infit
    w.z = EWX[,2] * (z ^ 2)
    w.z.sum = sum(w.z, na.rm = TRUE)
    w.sum = sum(EWX[,2], na.rm = TRUE)
    infit = w.z.sum / w.sum

    item.infit = c(item.infit, infit)
    item.outfit = c(item.outfit, outfit)
  }
  matrix = rbind(item.infit, item.outfit)
  rownames(matrix) = c("infit", "outfit")
  colnames(matrix) = item
  matrix
}