R/ggumProbability.R

In duckmayr/bggum: Bayesian Estimation of Generalized Graded Unfolding Model Parameters

#' GGUM Probability Function
#'
#' Calculate the probability of a response according to the GGUM
#'
#' The General Graded Unfolding Model (GGUM) is an item response model
#' designed to consider the possibility of disagreement for opposite reasons.
#' This function gives the probability of a respondent's response to a test
#' item given item and respondent parameters. The user can calculate the
#' probability of one particular response to an item, for any number of the
#' possible responses to the item, the probability of a vector of responses
#' (either responses by one person to multiple items, or by multiple people to
#' one item), or the probability of each response in a response matrix.
#'
#' The probability that respondent \eqn{i} chooses option \eqn{k} for item
#' \eqn{j} is given by
#' \deqn{\frac{\exp (\alpha_j [k (\theta_i - \delta_j) -
#' \sum_{m=0}^k \tau_{jm}]) + \exp (\alpha_j [(2K - k - 1)
#' (\theta_i - \delta_j) - \sum_{m=0}^k \tau_{jm}])}{%
#' \sum_{l=0}^{K-1} [\exp (\alpha_j [l (\theta_i - \delta_j) -
#' \sum_{m=0}^l \tau_{jm}]) + \exp (\alpha_j [(2K - l - 1)
#' (\theta_i - \delta_j) - \sum_{m=0}^l \tau_{jm}])]}}{%
#' (exp(\alpha_j [k(\theta_i-\delta_j) -
#' \sum_{m=0}^k \tau_{jm}]) + exp(\alpha_j [(2K - k - 1)
#' (\theta_i - \delta_j) - \sum_{m=0}^k \tau_{jm}])) /
#' (\sum_{l=0}^{K-1} [exp (\alpha_j [l (\theta_i - \delta_j) -
#' \sum_{m=0}^l \tau_{jm}]) + exp (\alpha_j [(2K - l - 1)
#' (\theta_i - \delta_j) - \sum_{m=0}^l \tau_{jm}])])},
#' where \eqn{\theta_i} is \eqn{i}'s latent trait parameter,
#' \eqn{\alpha_j} is the item's discrimination parameter,
#' \eqn{\delta_j} is the item's location parameter,
#' \eqn{\tau_{j0}, \ldots, \tau_{j(K-1)}} are the options' threshold
#' parameters, and \eqn{\tau_{j0}} is 0,
#' \eqn{K} is the number of options for item \eqn{j}, and
#' the options are indexed by \eqn{k = 0, \ldots, K-1}.
#'
#' @param response A numeric vector or matrix giving the response(s) for which
#'   probability should be calculated.
#' @param theta A numeric vector of latent trait score(s) for respondent(s)
#' @param alpha A numeric vector of discrimination parameter(s)
#' @param delta A numeric vector of location parameter(s)
#' @param tau A numeric vector (if responses to one item are given) or a list
#'   (if responses to multiple items are given); the tau parameters for each
#'   item is a numeric vector of length K (the number of possible responses)
#'   giving the options' threshold parameters; the first element of \code{tau}
#'   should be zero
#'
#' @return A matrix or vector of the same dimensions/length of \code{response}.
#'
#' @section Note:
#'   Please note that items' options should be zero-indexed.
#'
#' @references de la Torre, Jimmy, Stephen Stark, and Oleksandr S.
#'   Chernyshenko. 2006. \dQuote{Markov Chain Monte Carlo Estimation of Item
#'   Parameters for the Generalized Graded Unfolding Model.} \emph{Applied
#'   Psychological Measurement} 30(3): 216--232.
#' @references Roberts, James S., John R. Donoghue, and James E. Laughlin. 2000.
#'   \dQuote{A General Item Response Theory Model for Unfolding Unidimensional
#'   Polytomous Responses.} \emph{Applied Psychological Measurement}
#'   24(1): 3--32.
#'
#' @examples
#' ## What is the probability of a 1 response to a dichotomous item
#' ## with discrimination parameter 2, location parameter 0, and
#' ## option threshold vector (0, -1) for respondents at -1, 0, and 1
#' ## on the latent scale?
#' ggumProbability(response = rep(1, 3), theta = c(-1, 0, 1), alpha = 2,
#'                 delta = 0, tau = c(0, -1))
#' ## We can also use this function for getting the probability of all
#' ## observed responses given the data and item and person parameter estimtes.
#' ## Here's an example of that with some simulated data:
#' ## Simulate data with 10 items, each with four options, and 100 respondents
#' set.seed(123)
#' sim_data <- ggum_simulation(100, 10, 4)
#' head(ggumProbability(response = sim_data$response_matrix,
#'                      theta = sim_data$theta,
#'                      alpha = sim_data$alpha,
#'                      delta = sim_data$delta,
#'                      tau = sim_data$tau))
#'
#' @rdname ggumProbability
#' @export
ggumProbability <- function(response, theta, alpha, delta, tau) {
    if ( is.matrix(response) ) {
        n <- nrow(response)
        m <- ncol(response)
        if ( length(theta) != n | length(alpha) != m | length(delta) != m
             | length(tau) != m | !is.list(tau) ) {
            stop(paste("For a response matrix, provide parameter estimates",
                       "for all items and all respondents."), call. = FALSE)
        }
        result <- matrix(NA_real_, nrow = n, ncol = m)
        for ( j in 1:m ) {
            result[ , j] <- probCol(response[ , j], theta, alpha[j], delta[j],
                                    tau[[j]])
        }
    }
    else if ( length(theta) > 1 ) {
        if ( length(alpha) > 1 | length(delta) > 1 | is.list(tau) ) {
            stop(paste("For multiple items and respondents, provide a matrix",
                       "of responses."), call. = FALSE)
        }
        if ( length(response) == 1 ) {
            stop("Provide a response for each respondent.", call. = FALSE)
        }
        result <- probCol(response, theta, alpha, delta, tau)
    }
    else if ( length(alpha) > 1 ) {
        if ( length(response) == 1 ) {
            stop("Provide a response for each item.", call. = FALSE)
        }
        result <- probRow(response, theta, alpha, delta, tau)
    }
    else {
        if ( length(response) > 1 ) {
            response2 <- ifelse(is.na(response), 0, response)
            result <- sapply(response2, prob, theta, alpha, delta, tau)
        }
        else {
            result <- prob(ifelse(is.na(response), 0, response),
                           theta, alpha, delta, tau)
        }
    }
    result[which(is.na(response))] <- NA
    return(result)
}