Nothing
#' Get sufficient statistics
#'
#' Find sufficient statistics from a data set
#'
#' @param dat Data set with items as columns and examinees as rows. Missing
#' responses should be coded as NA.
#'
#' @return Vector of sufficient statistics
#'
#' @examples
#'
#' set.seed(2524)
#'
#' diff_1 <- rnorm(10)
#' diff_2 <- rnorm(10)
#'
#' N <- 500
#'
#' th <- MASS::mvrnorm(N, mu = c(0, -1),
#' Sigma = matrix(c(1, .5 * 2, .5 * 2, 4), nrow = 2))
#'
#' probs_1 <- 1 / (1 + exp(-outer(th[, 1], diff_1, "-")))
#' probs_2 <- 1 / (1 + exp(-outer(th[, 2], diff_2, "-")))
#'
#' probs <- cbind(probs_1, probs_2)
#'
#' dat <- apply(probs, 2, function(p) as.numeric(p > runif(N)))
#'
#' get_ss(dat)
#'
#' @export
get_ss <- function(dat){
# number of items
nitems <- ncol(dat)
# number of categories for each item
ncats <- apply(dat, 2, function(x) length(unique(x[!is.na(x)])))
# number of observations for each item
nobs <- apply(dat, 2, function(x) sum(!is.na(x)))
# find the proportion of subjects who have reached at least each level
ss <- as.numeric(unlist(sapply(1:nitems, function(i){
(nobs[i] - cumsum(summary(factor(dat[, i])))[1:(ncats[i] - 1)]) / nobs[i]
})))
ss
}
#' Get Thurstone thresholds
#'
#' Find Thurstone thresholds from a fitted model.
#'
#' @param pars Vector of estimated item parameters
#' @param itemtype Item type: "1PL", "PCM", or "PCM2".
#' @param item_ind Vector with one element for each parameter indicating which
#' item each parameter is associated with.
#' @param alpha Vector of item steepnesses, with one element for each item.
#' Recycled if of length 1.
#'
#' @return Vector of Thurstone thresholds
#'
#' @examples
#'
#' if(require(TAM)){
#'
#' set.seed(2524)
#'
#' diff <- rnorm(10)
#'
#' N <- 500
#'
#' th <- rnorm(N)
#'
#' probs <- 1 / (1 + exp(-outer(th, diff, "-")))
#'
#' dat <- apply(probs, 2, function(p) as.numeric(p > runif(N)))
#'
#' # fit the model
#'
#' mod <- TAM::tam.mml(resp = dat, irtmodel = "1PL")
#'
#' get_thresh(mod$xsi$xsi, itemtype = "1PL", item_ind = 1:10)
#'
#' }
#'
#' @export
get_thresh <- function(pars, itemtype, item_ind, alpha = 1){
ncats <- tapply(item_ind, item_ind, function(x) length(x) + 1)
if (length(alpha) == 1) alpha <- rep(alpha, max(item_ind))
if (itemtype == "PCM2"){
## transform to Masters' notation
out <- lapply(1:max(item_ind), function(i){
delta <- pars[item_ind == i][1]
tau <- pars[item_ind == i][-1]
tau <- c(tau, -sum(tau))
cbind(delta + tau, i)
})
out <- do.call(rbind, out)
pars <- out[, 1]
item_ind <- out[, 2]
}
# category response function - uses Masters' notation (not CQ notation)
prob <- function(theta, xsi, alpha){
xsi <- c(-sum(xsi), xsi)
nums <- outer(alpha * theta, xsi, "-")
nums <- t(exp(apply(nums, 1, cumsum)))
den <- rowSums(nums)
nums / den
}
f <- function(theta, pars, k, alpha){
sum(prob(theta = theta, xsi = pars, alpha = alpha)[1:k]) - .5
}
item <- rep(1:length(ncats), ncats - 1)
out <- lapply(1: length(ncats), function(i){
pars_i <- pars[item == i]
k <- length(pars_i)
sapply(1:k, function(x){
stats::uniroot(f = f, interval = c(-6, 6), pars = pars_i, k = x,
alpha = alpha[i], extendInt = "yes")$root
})
})
unlist(out)
}
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.