######################################################################
# 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.