R/Rasch.PCA.ltm.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. Please note that
#' this function utlized ltm's MMLE to estimate the parameter. It should be almost identical to
#' the result from Rasch.PCA.Bayes but faster.
#'
#' @param data A data frame containing the data
#' @param item Item to be included in Rasch
#' @export Rasch.PCA.ltm

Rasch.PCA.ltm = function(data, item) {
  data = na.omit(data[,item])
  P = length(item)
  N = nrow(data)
  K = max(apply(na.omit(data),2, max))
  obj = list()
  class(obj) = "BPCM"
  std.res = matrix(NA, nrow(data), ncol(data))
  raw.res = matrix(NA, nrow(data), ncol(data))
  colnames(std.res) = item

  A = sum(data - 1) / (nrow(data) * ncol(data))

  fit = ltm::gpcm(data, constraint = c("rasch"))
  theta.pattern = ltm::factor.scores.gpcm(fit)$score.dat
  beta.matrix = NULL

  for(i in 1:P){
    beta.matrix = rbind(beta.matrix, fit$coefficients[[i]])
  }

  data = theta.matching(data, item, theta.pattern)
  Emat = NULL
  Wmat = NULL
  for(i in 1:P){
    if(max(K) == 2){
      itemdiff = beta.matrix[i]
    }else{
      itemdiff = beta.matrix[i,]
    }
    EWX = matrix(NA, N, 3)
    EWX[,3] = data[,i] - 1
    for(j in 1:nrow(data)){
      E = 0
      w = 0
      person = 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)

    }
    Emat = cbind(Emat,EWX[,1])
    Wmat = cbind(Wmat,EWX[,2])
    raw.res[,i] = EWX[,3] - EWX[,1]
    std.res[,i] = (EWX[,3] - EWX[,1]) / sqrt(EWX[,2])
  }
  colnames(raw.res) = item
  colnames(std.res) = item

  OV = sum(((data[,1:P] - 1) - A) ^ 2)
  UV = sum(raw.res ^ 2)
  EV = OV - UV

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


  pa = psych::fa.parallel(raw.res, fa = "pc")
  loadings = psych::pca(raw.res)$loadings

  list(PCA = out, loadings)
}