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_LRtest.Rm
#
# Part of R/tlc - Testing in Conditional Likelihood Context package
#
# This file contains a routine that is adopted from orignal eRm function
# "LRtest.Rm" by adding a normalization condition (sum0 = FALSE).
# 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) 2019, Last Modified 09/09/2019
######################################################################

tcl_LRtest.Rm <-
    function(object, splitcr = "median",
             se = TRUE, sum0 = FALSE) { # argument sum0=False added AK 2019-08-15
  # performs Andersen LR-test
  # object... object of class 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.
  # se...whether standard errors should be computed
  # "sum0" added in "objectg <-  RM ..." etc.

  call <- match.call()

  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 = " ")
    }

  }

  #----------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!")),
            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

  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 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


  result <- list(X=X.original, 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"

  return(result)

}