R/item.selection.R

In changxiulee/BayesianRasch: Bayesian Rasch Modeling

#'Item Selection
#'
#'Given that a set of item is not unidimensional, this function helps to determine which item should be removed.
#'To do this, user need to first select a small set of item (core.item) that is known for sure to be unidimenional. This small set
#'of item will be subject to Rasch PCA to verify this belief. Then, for each item outside core.item (peripheral item),
#'a Rasch Analysis will be conducted together with the core item and the item fit is computed. If an item does not share
#'a common dimensional with the core item, the item fit of the peripheral item is worse (> 1.3). Please note that, in this
#'function, all Rasch model is estimate using ltm / MMLE, not Bayes because 1) estimating using Bayes is too slow and not
#'suitable when we want to estimate a model repeatedly 2) We do not utilized the uncertainty information
#'in this function, 3) LTM's MMLE is identical to Bayesian's Maximum a Posterior most of the time and it is faster.
#'
#'@param data A data frame containing the data
#'@param core.item A set of item that is obviously measuring the dimension in question
#'@param peripheral.item Item outside of core item
#'@export item.selection

item.selection = function(data, core.item, peripheral.item) {
  obj = list()
  class(obj) = "BPCM"

  #Core item fitness check
  core.data = data[, core.item]
  fit = ltm::gpcm(core.data, constraint = "rasch")
  theta.pattern = ltm::factor.scores.gpcm(fit)$score.dat
  check = Rasch.PCA.ltm(data, core.item)
  core.data = theta.matching(data, core.item, theta.pattern)
  beta.matrix = NULL
  for(i in 1:length(core.item)){
    beta.matrix = rbind(beta.matrix, fit$coefficients[[i]])
  }
  beta.matrix = beta.matrix[, -ncol(beta.matrix)]
  K = max(apply(na.omit(core.data),2, max))
  N = nrow(core.data)
  item.infit = c()
  item.outfit = c()
  for(i in 1:length(core.item)){
    if(max(K) == 2){
      itemdiff = beta.matrix[i]
    }else{
      itemdiff = beta.matrix[i,]
    }
    EWX = matrix(NA, N, 3)
    EWX[,3] = core.data[,i] - 1
    for(j in 1:N){
      E = 0
      w = 0
      person = core.data[,"theta"][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

    #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) = core.item

  show(check)
  show(matrix)

  f = length(core.item) + 1
  N = nrow(data)
  item.infit = c()
  item.outfit = c()
  K = max(apply(na.omit(data), 2, max))
  for(i in 1:length(peripheral.item)){
    temp.data = data[, c(core.item, peripheral.item[i])]
    fit = ltm::gpcm(temp.data, constraint = "rasch")
    beta.matrix = NULL
    P = ncol(temp.data)
    for(i in 1:P){
      beta.matrix = rbind(beta.matrix, fit$coefficients[[i]])
    }
    beta.matrix = beta.matrix[, -ncol(beta.matrix)]
    if(max(K) == 2){
      itemdiff = beta.matrix[f]
    }else{
      itemdiff = beta.matrix[f,]
    }
    EWX = matrix(NA, N, 3)
    EWX[,3] = data[,f] - 1
    for(j in 1:N){
      E = 0
      w = 0
      person = core.data[,"theta"][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

    #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) = peripheral.item
  matrix
}