#' Item Response Theory Study
#'
#' Functions for estimating item response theory item and person parameters.
#'
#' @param x matrix or data.frame of scored item responses, one row per person,
#' one column per item.
#' @param model integer for the number of paramters to estimate per item. The
#' default of 1 specifies item intercepts or difficulties via the Rasch model,
#' 2 additionally specifies item slopes or discriminations, and 3
#' additionally specifies varying lower asymptotes.
#' @param complete logical with default \code{FALSE} indicating whether or not
#' \code{x} should be reduced to rows with complete data across all columns.
#' @param ip vector or matrix of item parameters, with a in the first column, b
#' in the second, and c in the third.
#' @param theta vector of thetas over which to estimate the corresponding function.
#' @param ... further arguments passed to or from other functions, including
#' \code{\link{glmer}}
#' @examples
#'
#' ritems <- c("r414q02", "r414q11", "r414q06", "r414q09",
#' "r452q03", "r452q04", "r452q06", "r452q07", "r458q01",
#' "r458q07", "r458q04")
#' rsitems <- paste0(ritems, "s")
#' pisagbr <- PISA09[PISA09$cnt == "GBR", rsitems]
#' irtgbr <- irtstudy(pisagbr)
#'
#' @export
irtstudy <- function(x, model = 1, complete = FALSE, ...) {
if (!all(unlist(x) %in% c(0, 1, NA)))
stop("'x' can only contain score values 0, 1, and NA.")
ni <- ncol(x)
if (complete)
x <- x[complete.cases(x), ]
else {
allna <- rowSums(is.na(x)) == ni
if (any(allna)) {
x <- x[!allna, ]
warning(sum(allna), " cases with missing data on all items removed.")
}
}
np <- nrow(x)
inames <- colnames(x)
if (is.null(inames))
inames <- paste0("item", 1:ni)
xl <- data.frame(score = c(unlist(x)), person = rep(1:np, ni),
item = rep(1:ni, each = np), row.names = NULL)
m <- lme4::glmer(score ~ -1 + (1 | item) + (1 | person),
data = xl, family = "binomial")
dat <- data.frame(x, theta = unlist(lme4::ranef(m)$person))
ip <- data.frame(a = 1, b = -unlist(lme4::ranef(m)$item), c = 0,
row.names = inames)
dat$se <- rtef(ip, dat$theta)$se
out <- list(data = dat, ip = ip, np = np, ni = ni,
vc = unlist(lme4::VarCorr(m)),
fit = c(AIC = AIC(m), BIC = BIC(m), logLik = logLik(m),
deviance = deviance(m), df.residual = df.residual(m)))
class(out) <- c("irtstudy", "list")
return(out)
}
#' @export
print.irtstudy <- function(x, ...) {
cat("\nItem Response Theory Study\n\n")
cat(x$np, "people,", x$ni, "items\n\n")
cat("Model fit with lme4::glmer\n")
print(x$fit)
cat("\nRandom effects\n")
print(data.frame(Std.Dev = sqrt(x$vc), Var = x$vc,
row.names = c("person", "item")))
}
# Item response function
#' @rdname irtstudy
#' @export
rirf <- function(ip, theta = seq(-4, 4, length = 100)){
if (NCOL(ip) == 1)
ip <- cbind(1, ip, 0)
else if (NCOL(ip) == 2)
ip <- cbind(ip, 0)
if (NCOL(ip) != 3)
stop("'ip' can only contain up to three parameters per item.")
ni <- NROW(ip)
out <- rbind(sapply(1:ni, function(i) ip[i, 3] + (1 - ip[i, 3]) /
(1 + exp(ip[i, 1] * (-theta + ip[i, 2])))))
colnames(out) <- rownames(ip)
out <- data.frame(theta = theta, out)
return(out)
}
# Item information function
#' @rdname irtstudy
#' @export
riif <- function(ip, theta = seq(-4, 4, length = 100)){
pr <- rirf(ip, theta)[, -1]
ia <- ip[, 1]
ic <- ip[, 3]
qr <- 1 - pr
info <- ia^2 * ((pr - ic)^2 / (1 - ic)^2) * (qr / pr)
out <- data.frame(theta = theta, info)
return(out)
}
# Item error function
#' @rdname irtstudy
#' @export
rief <- function(ip, theta = seq(-4, 4, length = 100)){
out <- data.frame(theta = theta, 1 / sqrt(riif(ip, theta)[, -1]))
colnames(out)[-1] <- rownames(ip)
return(out)
}
# Test response function
#' @rdname irtstudy
#' @export
rtrf <- function(ip, theta = seq(-4, 4, length = 100)){
out <- data.frame(theta = theta,
p = apply(rirf(ip, theta)[, -1], 1, sum))
return(out)
}
# Test information function
#' @rdname irtstudy
#' @export
rtif <- function(ip, theta = seq(-4, 4, length = 100)){
out <- data.frame(theta = theta,
i = apply(riif(ip, theta)[, -1], 1, sum))
return(out)
}
# Test error function
#' @rdname irtstudy
#' @export
rtef <- function(ip, theta = seq(-4, 4, length = 100)){
out <- data.frame(theta = theta,
se = 1 / sqrt(rtif(ip, theta)$i))
return(out)
}
# Simulate dichotomous responses
#' @rdname irtstudy
#' @export
irtsim <- function(ip, theta) {
p <- rirf(ip, theta)[, -1]
nr <- prod(dim(p))
out <- ifelse(p > runif(nr), 1, 0)
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.