R/Rasch.PCA.Bayes.R

In changxiulee/BayesianRasch: Bayesian Rasch Modeling

#' Rasch Residual Principal Component Analysis
#'
#' This is for unidimensionality assumption check. If the data is unidimensional,
#' the eigenvalue for the first constrast should be less than 2. However, when
#' the eigenvalue is more than 2, it can either indicate a local change in intensity
#' or the multidimensionality. To know this, use the item.selection method. Rasch's parameter in this function
#' is estimated using Bayesian.
#'
#' @param obj A BPCM or BRSM obj
#' @export Rasch.PCA.Bayes

Rasch.PCA.Bayes = function(obj) {

  if(obj$DIF == TRUE){
    stop("This is not available for DIF object")
  }

  codaSamples = obj$mcmc
  data = obj$data
  item = obj$item

  data = data[,item]
  P = length(item)
  N = nrow(data)
  K = max(apply(na.omit(data),2, max))
  g = 1

  std.res = matrix(NA, nrow(data), ncol(data))
  raw.res = matrix(NA, nrow(data), ncol(data))
  colnames(std.res) = item

  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, 3)
    EWX[,3] = data[,i] - 1
    for(j in 1:nrow(data)){
      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)
    }
    raw.res[,i] = EWX[,3] - EWX[,1]
    std.res[,i] = (EWX[,3] - EWX[,1]) / sqrt(EWX[,2])
  }

  OV = sum(apply(na.omit(data), 2, var))
  UV = sum(apply(na.omit(raw.res), 2, function(x){mean(x ^ 2)}))
  EV = OV - UV
  UV.percent = UV * 100 / OV
  OEigen = 100/ UV.percent * length(item)
  EEigen = OEigen - length(item)
  eigenValues = eigen(cor(na.omit(std.res)))$values
  for(i in 1:length(eigenValues)){
    if(eigenValues[i] < 1) {
      eigenValuesIndex = i - 1
      break
    }
  }
  out = matrix(NA, 3 + eigenValuesIndex, 3)
  out[,1] = c("Total Variance", "Rasch Variance", "Residual Variable", paste0("Constrast ", 1:eigenValuesIndex))
  out[1, 2:3] = c(OEigen, "100%")
  out[2, 2:3] = c(EEigen, paste0(round(100 - UV.percent), "%"))
  out[3, 2:3] = c(OEigen - EEigen, paste0(round(UV.percent), "%"))
  out[4:(4 + (eigenValuesIndex) - 1), 2] = eigenValues[1:eigenValuesIndex]
  out[,2] = round(as.numeric(out[,2]), 2)
  colnames(out) = c("", "Eigenvalues", "Percentage of Variance")
  out
}