R/tcl_splitcr_RStest.R

In tcl: Testing in Conditional Likelihood Context

######################################################################
# UMIT - Private University for Health Sciences,
#        Medical Informatics and Technology
#        Institute of Psychology
#        Statistics and Psychometrics Working Group
#
# tcl_splitcr_RStest
#
# Part of R/tlc - Testing in Conditional Likelihood Context package
#
# This file contains a routine that is adopted from orignal eRm function
# "LRtest.Rm"
# original URL: https://github.com/cran/eRm/blob/master/R/LRtest.Rm.R
#
# Licensed under the GNU General Public License Version 3 (June 2007)
# copyright (c) 2021, Last Modified 24/01/2022
######################################################################

tcl_splitcr_RStest <- function(X, splitcr = "median", model ="RM") {
  # splitcr... splitting criterion for LR-groups. "all.r" corresponds to a complete
  #            raw score split (r=1,...,k-1), "median" to a median raw score split,
  #            "mean" corresponds to the mean raw score split.
  #            optionally also a vector of length n for group split can be submitted.

  call <- match.call()

  colnames(X) <- paste("I",1:ncol(X),sep="")
  rownames(X) <- paste0("P", seq_len(nrow(X)))

  object <- list()
  object$X <- X
  object$model <- model

  spl.gr <- NULL

  X.original <- object$X
  if ((length(splitcr) > 1) & is.character(splitcr)) {
    # if splitcr is character vector, treated as factor
    splitcr <- as.factor(splitcr)
  }
  if (is.factor(splitcr)) {
    spl.nam <- deparse(substitute(splitcr))
    spl.lev <- levels(splitcr)
    spl.gr <- paste(spl.nam, spl.lev, sep = " ")
    splitcr <- unclass(splitcr)
  }

  numsplit <- is.numeric(splitcr)
  if (any(is.na(object$X))) {
    if (!numsplit && splitcr == "mean") {
      # mean split
      spl.gr <- c("Raw Scores < Mean", "Raw Scores >= Mean")
      X <- object$X
      # calculates index for NA groups from person.parameter.eRm
      dichX <- ifelse(is.na(X), 1, 0)
      strdata <- apply(dichX, 1, function(x) {
        paste(x, collapse = "")
      })
      gmemb <- as.vector(data.matrix(data.frame(strdata)))
      gindx <- unique(gmemb)
      rsum.all <- rowSums(X, na.rm = T)
      grmeans <- tapply(rsum.all, gmemb, mean)  #sorted
      ngr <- table(gmemb)  #sorted
      m.all <- rep(grmeans, ngr)  #sorted,expanded
      rsum.all <- rsum.all[order(gmemb)]
      spl <- ifelse(rsum.all < m.all, 1, 2)
      splitcr <- spl
      object$X <- X[order(gmemb), ]
    }
    if (!numsplit && splitcr == "median") {
      # median split
      spl.gr <- c("Raw Scores <= Median", "Raw Scores > Median")
      # cat('Warning message: Persons with median raw scores are
      # assigned to the lower raw score group!\n')
      X <- object$X
      # calculates index for NA groups from person.parameter.eRm
      dichX <- ifelse(is.na(X), 1, 0)
      strdata <- apply(dichX, 1, function(x) {
        paste(x, collapse = "")
      })
      gmemb <- as.vector(data.matrix(data.frame(strdata)))
      gindx <- unique(gmemb)
      rsum.all <- rowSums(X, na.rm = T)
      grmed <- tapply(rsum.all, gmemb, median)  #sorted
      ngr <- table(gmemb)  #sorted
      m.all <- rep(grmed, ngr)  #sorted,expanded
      rsum.all <- rsum.all[order(gmemb)]
      spl <- ifelse(rsum.all <= m.all, 1, 2)
      splitcr <- spl
      object$X <- X[order(gmemb), ]
    }
  }

  if (!is.numeric(splitcr)) {
    if (splitcr == "all.r")
    {
      # full raw score split ### begin MjM 2012-03-18
      rvind <- rowSums(object$X, na.rm = TRUE)  #person raw scoobject
      excl_0_k <- (rvind > 0) & (rvind < sum(apply(object$X,
                                                   2, max, na.rm = T)))
      Xlist <- by(object$X[excl_0_k, ], rvind[excl_0_k],
                  function(x) x)
      names(Xlist) <- as.list(paste("Raw Score =",
                                    sort(unique(rvind[excl_0_k]))))
      spl.gr <- unlist(names(Xlist))
    }  ### end MjM 2012-03-18

    if (splitcr == "median") {
      # median split
      spl.gr <- c("Raw Scores <= Median", "Raw Scores > Median")
      # removed rh 2010-12-17 cat('Warning message: Persons with
      # median raw scores are assigned to the lower raw score
      # group!\n')
      rv <- apply(object$X, 1, sum, na.rm = TRUE)
      rvsplit <- median(rv)
      rvind <- rep(0, length(rv))
      rvind[rv > rvsplit] <- 1  #group with highraw scoobject
      Xlist <- by(object$X, rvind, function(x) x)
      names(Xlist) <- list("low", "high")
    }

    if (splitcr == "mean") {
      # mean split
      spl.gr <- c("Raw Scores < Mean", "Raw Scores >= Mean")
      rv <- apply(object$X, 1, sum, na.rm = TRUE)
      rvsplit <- mean(rv)
      rvind <- rep(0, length(rv))
      rvind[rv > rvsplit] <- 1  #group with highraw scoobject
      Xlist <- by(object$X, rvind, function(x) x)
      names(Xlist) <- list("low", "high")
    }
  }

  if (is.numeric(splitcr)) {
    # manual raw score split
    spl.nam <- deparse(substitute(splitcr))
    if (length(splitcr) != dim(object$X)[1])
      stop("Mismatch between length of split vector and number of persons!")
    if (any(is.na(splitcr)))
      stop("Split vector should not contain NA's")

    rvind <- splitcr
    Xlist <- by(object$X, rvind, function(x) x)
    names(Xlist) <- as.list(sort(unique(splitcr)))
    if (is.null(spl.gr)) {
      spl.lev <- names(Xlist)
      spl.gr <- paste(spl.nam, spl.lev, sep = " ")
    }

  }

  ## code block commented out AK 24-01-2022
  #----------item to be deleted---------------
  # del.pos.l <- lapply(Xlist, function(x) {
  #   it.sub <- datcheck.LRtest(x, object$X, object$model)  #items to be removed within subgroup
  # })
  #
  # del.pos <- unique(unlist(del.pos.l))
  # if (length(del.pos) >= (ncol(object$X) - 1)) {
  #   stop("\nNo items with appropriate response patterns left to perform LR-test!\n")
  # }
  #
  # if (length(del.pos) > 0)
  # {
  #   ### begin MjM 2013-01-27
  #   warning(paste0("\n", prettyPaste("The following items were excluded due to inappropriate response patterns within subgroups:"),
  #                  "\n", paste(colnames(object$X)[del.pos], collapse = " "),
  #                  "\n\n", prettyPaste("Full and subgroup models are estimated without these items!")),
  #           call. = FALSE, immediate. = TRUE)
  # }  ### end MjM 2013-01-27
  #
  #
  # if (length(del.pos) > 0) {
  #   X.el <- object$X[, -(del.pos)]
  # } else {
  #   X.el <- object$X
  # }
  #
  # if (ifelse(length(splitcr) == 1, splitcr != "all.r", TRUE)) {
  #   ### begin MjM 2012-03-18 # for all cases except 'all.r'
  #   Xlist.n <- by(X.el, rvind, function(y) y)
  #   names(Xlist.n) <- names(Xlist)
  #   if (length(del.pos) > 0)
  #     Xlist.n <- c(Xlist.n, list(X.el))  # X.el added since we must refit whole group without del.pos items
  # } else {
  #   Xlist.n <- by(X.el[excl_0_k, ], rvind[excl_0_k], function(y) y)
  #   names(Xlist.n) <- names(Xlist)
  #   Xlist.n <- c(Xlist.n, list(X.el[excl_0_k, ]))  # X.el added since we must refit whole group without del.pos items
  # }  ### end MjM 2012-03-18

  ## code block commented out AK 24-04-2021
  # if (object$model == "RM") {
  #   likpar <- sapply(Xlist.n, function(x) {
  #     # matrix with loglik and npar for each subgroup
  #     # objectg <- RM(x,se=se) # original eRm code !!
  #     objectg <- eRm::RM(x, se = se, sum0 = sum0)  # argument "sum0" added AK 2019-08-15
  #     likg <- objectg$loglik
  #     nparg <- length(objectg$etapar)
  #     # betalab <- colnames(objectg$X)
  #     list(likg, nparg, objectg$betapar, objectg$etapar,
  #          objectg$se.beta, outobj = objectg)  # rh outobj added
  #     ### list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta)
  #     ### # rh outobj added
  #   })
  # }
  # if (object$model == "PCM") {
  #   likpar <- sapply(Xlist.n, function(x) {
  #     # matrix with loglik and npar for each subgroup
  #     # objectg <- PCM(x,se=se) # original eRm code !!
  #     objectg <- eRm::PCM(x, se = se, sum0 = sum0)   # argument "sum0" added added AK 2019-08-15
  #     likg <- objectg$loglik
  #     nparg <- length(objectg$etapar)
  #     list(likg, nparg, objectg$betapar, objectg$etapar,
  #          objectg$se.beta, outobj = objectg)  # rh outobj added
  #     ### list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta)
  #     ### # rh outobj added
  #   })
  # }
  # if (object$model == "RSM") {
  #   likpar <- sapply(Xlist.n, function(x) {
  #     # matrix with loglik and npar for each subgroup
  #     # objectg <- RSM(x,se=se) # original eRm code !!
  #     objectg <- eRm::RSM(x, se = se, sum0 = sum0)  # argument "sum0" added AK 2019-08-15
  #     likg <- objectg$loglik
  #     nparg <- length(objectg$etapar)
  #     list(likg, nparg, objectg$betapar, objectg$etapar,
  #          objectg$se.beta, outobj = objectg)  # rh outobj added
  #     ### list(likg,nparg,objectg$betapar,objectg$etapar,objectg$se.beta)
  #     ### # rh outobj added
  #   })
  # }
  #
  # ## extract fitted splitgroup models # rh 02-05-2010 begin MjM
  # ## 2012-03-18
  # if (ifelse(length(splitcr) == 1, splitcr != "all.r", TRUE)) {
  #   fitobj <- likpar[6, 1:length(unique(rvind))]
  # } else {
  #   fitobj <- likpar[6, 1:length(unique(rvind[excl_0_k]))]
  # }  ### end MjM 2012-03-18
  # likpar <- likpar[-6, ]
  #
  # if ((length(del.pos) > 0) | ifelse(length(splitcr) == 1,
  #                                    splitcr == "all.r", FALSE)) {
  #   # re-estimate full model ### MjM 2012-03-18
  #   pos <- length(Xlist.n)  #position of the full model
  #   loglik.all <- likpar[1, pos][[1]]  #loglik full model
  #   # etapar.all <- rep(0,likpar[2,pos]) #etapar full model
  #   # (filled with 0 for df computation)
  #   etapar.all <- rep(0, unlist(likpar[2, pos]))  #etapar full model (filled with 0 for df computation)
  #   likpar <- likpar[, -pos]
  #   Xlist.n <- Xlist.n[-pos]
  # } else {
  #   loglik.all <- object$loglik
  #   etapar.all <- object$etapar
  # }
  #
  # loglikg <- sum(unlist(likpar[1, ]))  #sum of likelihood value for subgroups
  # LR <- 2 * (abs(loglikg - loglik.all))  #LR value
  # df = sum(unlist(likpar[2, ])) - (length(etapar.all))  #final degrees of freedom
  # pvalue <- 1 - pchisq(LR, df)  #pvalue
  #
  # betalist <- likpar[3, ]  #organizing betalist

  # AK 24.04.2021 X.el added
  # result <- list(X=X.original, X.el=X.el, X.list=Xlist.n, model=object$model,LR=LR,
  #                df=df, pvalue=pvalue, likgroup=unlist(likpar[1,],use.names=FALSE),
  #                betalist=betalist, etalist=likpar[4,],selist=likpar[5,], spl.gr=spl.gr, call=call, fitobj=fitobj)  ## rh fitobj added
  # class(result) <- "LR"
  # result <- list(X=X.original, X.el=X.el, X.list=Xlist.n, model=object$model, call=call) # AK 24.01.2022
  result <- list(X=X.original, X_list=Xlist, model=object$model,call=call)

  return(result)

}