Nothing
######################################################################
# UMIT - Private University for Health Sciences,
# Medical Informatics and Technology
# Institute of Psychology
# Statistics and Psychometrics Working Group
#
# tcl_splitcr
#
# 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 22/04/2021
######################################################################
tcl_splitcr <- 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.
oldw <- getOption("warn") # addded AK 20-02-2022
options(warn = -1)
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 = " ")
}
}
#----------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 for the computation of GR,LR, and W due to inappropriate response patterns within subgroups:"),
# 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"
if(rlang::is_empty(del.pos)) del.pos<-NA # addded AK 20-02-022
result <- list(X=X.original, X.el=X.el, X.list=Xlist.n, model=object$model, del_pos=del.pos, call=call)
return(result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.