R/reval_mst.R

In irtQ: Unidimensional Item Response Theory Modeling

#' Recursion-based MST evaluation method
#'
#' This function evaluates the measurement precision and bias in
#' Multistage-Adaptive Test (MST) panels using a recursion-based evaluation
#' method introduced by Lim et al. (2021). This function computes conditional
#' biases and standard errors of measurement (CSEMs) across a range of IRT
#' ability levels, facilitating efficient and accurate MST panel assessments
#' without extensive simulations.
#'
#' @param x A data frame containing the metadata for the item bank, which
#'   includes important information for each item such as the number of score
#'   categories and the IRT model applied. This metadata is essential for
#'   evaluating the MST panel, with items selected based on the specifications
#'   in the `module` argument. To construct this item metadata efficiently, the
#'   [irtQ::shape_df()] function is recommended. Further details on utilizing
#'   item bank metadata along with `module` for MST panel evaluation are
#'   provided below.
#' @param D A scaling constant used in IRT models to make the logistic function
#'   closely approximate the normal ogive function. A value of 1.7 is commonly
#'   used for this purpose. Default is 1.
#' @param route_map A binary square matrix that defines the MST structure,
#'   illustrating transitions between modules and stages. This concept and
#'   structure are inspired by the `transMatrix` argument in the `randomMST()`
#'   function from the \pkg{mstR} package (Magis et al., 2017), which provides a
#'   framework for representing MST pathways. For comprehensive understanding
#'   and examples of constructing `route_map`, refer to the \pkg{mstR} package
#'   (Magis et al., 2017) documentation. Also see below for details.
#' @param module A binary matrix that maps items from the item bank specified in
#'   `x` to modules within the MST framework. This parameter's structure is
#'   analogous to the `modules` argument in the `randomMST()`  function of the
#'   \pkg{mstR} package, enabling precise item-to-module assignments for MST
#'   configurations. For detailed instructions on creating and utilizing the
#'   `module` matrix effectively, consult the documentation of the \pkg{mstR}
#'   package (Magis et al., 2017). Also see below for details.
#' @param cut_score A list defining cut scores for routing test takers through
#'   MST stages. Each list element is a vector of cut scores for advancing
#'   participants to subsequent stage modules. In a 1-3-3 MST configuration, for
#'   example, `cut_score` might be defined as `cut_score = list(c(-0.5, 0.5),
#'   c(-0.6, 0.6))`, where `c(-0.5, 0.5)` are thresholds for routing from the
#'   first to the second stage, and `c(-0.6, 0.6)` for routing from the second
#'   to the third stage. This setup facilitates tailored test progression based
#'   on performance. Further examples and explanations are available below.
#' @param theta A vector of ability levels (theta) at which the MST panel's
#'   performance is assessed. This allows for the evaluation of measurement
#'   precision and bias across a continuum of ability levels. The default range
#'   is `theta = seq(-5, 5, 0.1)`.
#' @param intpol A logical value to enable linear interpolation in the inverse
#'   test characteristic curve (TCC) scoring, facilitating ability estimate
#'   approximation for observed sum scores not directly obtainable from the TCC,
#'   such as those below the sum of item guessing parameters. Default is TRUE,
#'   applying interpolation to bridge gaps in the TCC. Refer to
#'   [irtQ::est_score()] for more details and consult Lim et al. (2021) for
#'   insights into the interpolation technique within inverse TCC scoring.
#' @param range.tcc A vector to define the range of ability estimates for
#'   inverse TCC scoring, expressed as the two numeric values for lower and
#'   upper bounds. Default is to c(-7, 7).
#' @param tol A numeric value of the convergent tolerance for the  inverse TCC
#'   scoring. For the inverse TCC scoring, the bisection method is used for
#'   optimization. Default is 1e-4.
#'
#' @details The [irtQ::reval_mst()] function evaluates an MST panel by
#'   implementing a recursion-based method to assess measurement precision
#'   across IRT ability levels. This approach, detailed in Lim et al. (2021),
#'   enables the computation of conditional biases and CSEMs efficiently,
#'   bypassing the need for extensive simulations traditionally required for MST
#'   evaluation.
#'
#'   The `module` argument, used in conjunction with the item bank metadata `x`,
#'   systematically organizes items into modules for MST panel evaluation. Each
#'   row of `x` corresponds to an item, detailing its characteristics like score
#'   categories and IRT model. The `module` matrix, structured with the same
#'   number of rows as `x` and columns representing modules, indicates item
#'   assignments with 1s. This precise mapping enables the [irtQ::reval_mst()]
#'   function to evaluate the MST panel's performance by analyzing how items
#'   within each module contribute to measurement precision and bias, reflecting
#'   the tailored progression logic inherent in MST designs.
#'
#'   The `route_map` argument is essential for defining the MST's structure by
#'   indicating possible transitions between modules. Similar to the
#'   `transMatrix()` in the \pkg{mstR} package (Magis et al., 2017), `route_map`
#'   is a binary matrix that outlines which module transitions are possible
#'   within an MST design. Each "1" in the matrix represents a feasible
#'   transition from one module (row) to another (column), effectively mapping
#'   the flow of test takers through the MST based on their performance. For
#'   instance, a "1" at the intersection of row *i* and column *j* indicates the
#'   possibility for test takers to progress from the module corresponding to
#'   row *i* directly to the module denoted by column *j*. This structure allows
#'   [irtQ::reval_mst()] to simulate and evaluate the dynamic routing of test
#'   takers through various stages and modules of the MST panel.
#'
#'   To further detail the `cut_score` argument with an illustration: In a 1-3-3
#'   MST configuration, the list `cut_score = list(c(-0.5, 0.5), c(-0.6, 0.6))`
#'   operates as a decision guide at each stage. Initially, all test takers
#'   start in the first module. Upon completion, their scores determine their
#'   next stage module: scores below -0.5 route to the first module of the next
#'   stage, between -0.5 and 0.5 to the second, and above 0.5 to the third. This
#'   pattern allows for dynamic adaptation, tailoring the test path to
#'   individual performance levels.
#'
#' @return This function returns a list of seven internal objects. The four
#'   objects are:
#'
#' \item{panel.info}{A list of several sub-objects containing detailed information
#' about the MST panel configuration, including:
#'     \describe{
#'       \item{config}{A nested list indicating the arrangement of modules across
#'       stages, showing which modules are included in each stage. For example,
#'       the first stage includes module 1, the second stage includes
#'       modules 2 to 4, and so forth.}
#'       \item{pathway}{A matrix detailing all possible pathways through the MST panel.
#'       Each row represents a unique path a test taker might take, based on their
#'       performance and the cut scores defined.}
#'       \item{n.module}{A named vector indicating the number of modules available
#'       at each stage.}
#'       \item{n.stage}{A single numeric value representing the total number of
#'       stages in the MST panel.}
#'    }
#' }
#'
#' \item{item.by.mod}{A list where each entry represents a module
#' in the MST panel, detailing the item metadata within that module. Each module's
#' metadata includes item IDs, the number of categories, the IRT model
#' used (model), and the item parameters (e.g., par.1, par.2, par.3).}
#'
#' \item{item.by.path}{A list containing item metadata arranged
#' according to the paths through the MST structure. This detailed
#' breakdown allows for an analysis of item characteristics along specific
#' MST paths. Each list entry corresponds to a testing stage and path,
#' providing item metadata. This structure facilitates the examination of
#' how items function within the context of each unique path through the MST.}
#'
#' \item{eq.theta}{
#' Estimated ability levels (\eqn{\theta}) corresponding to the observed
#' scores, derived from the inverse TCC scoring method. This provides the
#' estimated \eqn{\theta} values for each potential pathway through the
#' MST stages. For each stage, \eqn{\theta} values are calculated for each
#' path, indicating the range of ability levels across the test takers.
#' For instance, in a three-stage MST, the `eq.theta` list may contain
#' \eqn{\theta} estimates for multiple paths within each stage, reflecting
#' the progression of ability estimates as participants move through the test.
#' The example below illustrates the structure of `eq.theta` output for a
#' 1-3-3 MST panel with varying paths:
#'      \describe{
#'        \item{stage.1}{`path.1` shows \eqn{\theta} estimates ranging from -7 to +7,
#'        demonstrating the initial spread of abilities.}
#'        \item{stage.2}{Multiple paths (`path.1`, `path.2`, ...) each with
#'        their own \eqn{\theta} estimates, indicating divergence in ability levels
#'        based on test performance.}
#'        \item{stage.3}{Further refinement of \eqn{\theta} estimates across paths,
#'        illustrating the final estimation of abilities after the last stage.}
#'   }
#' }
#'
#' \item{cdist.by.mod}{A list where each entry contains the conditional distributions
#' of the observed scores for each module given the ability levels.}
#'
#' \item{jdist.by.path}{Joint distributions of observed scores for different paths
#' at each stage in a MST panel. The example below outlines the organization of
#' `jdist.by.path` data in a hypothetical 1-3-3 MST panel:
#'      \describe{
#'        \item{stage.1}{Represents the distribution at the initial stage, indicating
#'        the broad spread of test-taker abilities at the outset.}
#'        \item{stage.2}{Represents the conditional joint distributions of the observed
#'        scores as test-takers move through different paths at the stage 2, based on
#'        their performance in earlier stages.}
#'        \item{stage.3}{Represents a further refinement of joint distribution of
#'        observed scores as test-takers move through different paths at the final
#'        stage 3, based on their performance in earlier stages.}
#'  }
#' }
#'
#' \item{eval.tb}{
#' A table summarizing the measurement precision of the MST panel. It contains the true
#' ability levels (`theta`) with the mean ability estimates (`mu`), variance
#' (`sigma2`), bias, and conditional standard error of measurement (CSEM) given
#' the true ability levels. This table highlights the MST panel's accuracy and
#' precision across different ability levels, providing insights into its effectiveness
#' in estimating test-taker abilities.}
#'
#' @author Hwanggyu Lim \email{hglim83@@gmail.com}
#'
#' @seealso [irtQ::shape_df()], [irtQ::est_score()]
#'
#' @references Lim, H., Davey, T., & Wells, C. S. (2021). A recursion-based
#'   analytical approach to evaluate the performance of MST. *Journal of
#'   Educational Measurement, 58*(2), 154-178. \doi{10.1111/jedm.12276}.
#'
#'   Magis, D., Yan, D., & von Davier, A. A. (2017). *Computerized adaptive and
#'   multistage testing with R: Using packages catR and mstR*. Springer.
#'   \doi{10.1007/978-3-319-69218-0}.
#'
#' @examples
#' \donttest{
#' ## ------------------------------------------------------------------------------
#' # Evaluation of a 1-3-3 MST panel using simMST data.
#' # This panel was assembled using module-level target TIFs and design
#' # constraints similar to those used in Lim et al.'s (2021) simulation study
#' # (it is not the identical dataset from that study).
#' # Details:
#' #    (a) Panel configuration: 1-3-3 MST panel
#' #    (b) Test length: 24 items (each module contains 8 items across all stages)
#' #    (c) IRT model: 3-parameter logistic model (3PLM)
#' ## ------------------------------------------------------------------------------
#' # Load the necessary library
#' library(dplyr)
#' library(tidyr)
#' library(ggplot2)
#'
#' # Import item bank metadata
#' x <- simMST$item_bank
#'
#' # Import module information
#' module <- simMST$module
#'
#' # Import routing map
#' route_map <- simMST$route_map
#'
#' # Import cut scores for routing to subsequent modules
#' cut_score <- simMST$cut_score
#'
#' # Import ability levels (theta) for evaluating measurement precision
#' theta <- simMST$theta
#'
#' # Evaluate MST panel using the reval_mst() function
#' eval <-
#'   reval_mst(x,
#'     D = 1.702, route_map = route_map, module = module,
#'     cut_score = cut_score, theta = theta, range.tcc = c(-7, 7)
#'   )
#'
#' # Review evaluation results
#' # The evaluation result table below includes conditional biases and
#' # standard errors of measurement (CSEMs) across ability levels
#' print(eval$eval.tb)
#'
#' # Generate plots for biases and CSEMs
#' p_eval <-
#'   eval$eval.tb %>%
#'   dplyr::select(theta, bias, csem) %>%
#'   tidyr::pivot_longer(
#'     cols = c(bias, csem),
#'     names_to = "criterion", values_to = "value"
#'   ) %>%
#'   ggplot2::ggplot(mapping = ggplot2::aes(x = theta, y = value)) +
#'   ggplot2::geom_point(mapping = ggplot2::aes(shape = criterion), size = 3) +
#'   ggplot2::geom_line(
#'     mapping = ggplot2::aes(
#'       color = criterion,
#'       linetype = criterion
#'     ),
#'     linewidth = 1.5
#'   ) +
#'   ggplot2::labs(x = expression(theta), y = NULL) +
#'   ggplot2::theme_classic() +
#'   ggplot2::theme_bw() +
#'   ggplot2::theme(legend.key.width = unit(1.5, "cm"))
#' print(p_eval)
#' }
#'
#' @export
#' @import dplyr
reval_mst <- function(x,
                      D = 1,
                      route_map,
                      module,
                      cut_score,
                      theta = seq(-5, 5, 1),
                      intpol = TRUE,
                      range.tcc = c(-7, 7),
                      tol = 1e-4) {

  ## -----------------------------------------------------
  # Extract all the panel information from the route map
  ## -----------------------------------------------------
  # Run the panel_info() function
  panel_data <- panel_info(route_map)

  # Pathways allowed
  pathway <- panel_data$pathway

  # Number of modules for each stage
  n.mod <- panel_data$n.module

  # Total number of modules across all stages
  tn.mod <- sum(n.mod)

  # Number of stages
  n.stg <- panel_data$n.stage

  ## -----------------------------------------------------
  # Create List of item metadata for modules and paths
  # across the stages
  ## -----------------------------------------------------
  # List containing item metadata of all modules
  meta_mod <-
    purrr::map(
      .x = 1:tn.mod,
      .f = ~ {
        x[module[, .x] == 1, ] %>%
          tibble::remove_rownames()
      }
    )
  names(meta_mod) <- paste0("m.", 1:tn.mod)

  # List containing item metadata for all possible (sub) pathways at each stage
  meta_path <-
    purrr::map(
      .x = 1:n.stg,
      .f = ~ {
        # Unique pathways upto the current stage
        path_uni <- unique(pathway[, 1:.x, drop = FALSE])

        # Item metadata for the unique pathways
        meta_uni <-
          purrr::map(
            .x = 1:nrow(path_uni),
            .f = ~ {
              meta_mod[path_uni[.x, ]] %>%
                dplyr::bind_rows()
            }
          )

        # Add the names and return the list
        names(meta_uni) <- paste0("path.", 1:length(meta_uni))
        return(meta_uni)
      }
    )
  names(meta_path) <- paste0("stage.", 1:n.stg)

  # Count the cumulative numbers of all possible (sub) pathways by each stage
  cum_path <- purrr::map_dbl(.x = meta_path,
                             .f = ~ {length(.x)})

  ## -----------------------------------------------------
  # Estimate the IRT theta scores using the Inverse-TCC method
  ## -----------------------------------------------------
  # Compute the IRT thetas corresponding to the all observed scores using
  # the Inverse-TCC scoring method for the all (sub) pathways
  eq_theta <-
    purrr::map(
      .x = meta_path,
      .f = ~ {
        # Implement the Inverse-TCC scoring for the current pathway
        tmp_meta <- .x
        tmp_eqtheta <-
          purrr::map(
            .x = tmp_meta,
            .f = ~ {
              inv_tcc_nr(
                x = .x, D = D,
                tol = tol, intpol = intpol,
                range.tcc = range.tcc
              )$est.theta
            }
          ) %>%
          do.call(what = "cbind")
      }
    )

  ## -----------------------------------------------------
  # Compute the conditional distributions of modules at each theta
  ## -----------------------------------------------------
  # List containing the conditional distributions of the modules given the ability levels across all modules.
  # Each element of the list includes the conditional distributions across all ability levels for each module.
  cdist_by_mod <-
    purrr::map(
      .x = meta_mod,
      .f = ~ {
        lwrc(x = .x, theta = theta, D = D)
      }
    )

  # Let's restructure the list above.
  # We need a list of conditional distributions of the modules across all ability levels.
  # Thus, each element of the list includes the conditional distributions for all
  # modules given each ability level.
  cdist_by_th <-
    purrr::map(
      .x = 1:length(theta),
      .f = ~ {
        num <- .x
        tmp_tb <-
          purrr::map(
            .x = cdist_by_mod,
            .f = ~ {
              .x[, num]
            }
          ) %>%
          bind.fill(type = "cbind", fill = 0L)
        colnames(tmp_tb) <- paste0("m.", 1:tn.mod)
        tmp_tb
      }
    )
  names(cdist_by_th) <- theta

  ## -----------------------------------------------------
  # Compute the joint conditional distributions at each theta
  ## -----------------------------------------------------
  # List of possible observed scores across all modules
  score_mod <-
    purrr::map(
      .x = meta_mod[pathway[1, ]],
      .f = ~ {
        0:sum(.x$cats - 1)
      }
    )
  names(score_mod) <- paste0("stage.", 1:n.stg)

  # Cumulative maximum sum scores across the stages
  score_maxcum <- cumsum(sapply(X = score_mod, FUN = "max"))

  # Empty list to include the cut scores used to route a test taker to next module across all stages
  cut4route <- vector("list", n.stg - 1)
  names(cut4route) <- paste0("stage.", 2:n.stg)

  # Empty list to include path labels across all stages
  path_label <- vector("list", n.stg - 1)
  names(path_label) <- paste0("stage.", 2:n.stg)

  # Empty list to include probabilities of the joint distributions across the stages
  # Each element of the list will contain the joint distributions across all ability levels
  joint_dist <- vector("list", n.stg)
  names(joint_dist) <- paste0("stage.", 1:n.stg)

  # For the routing module in the first stage,
  # assign the conditional distributions of the routing modules across the ability levels
  joint_dist[[1]] <-
    purrr::map(
      .x = cdist_by_th,
      .f = ~ {
        .x[, 1, drop = FALSE] %>%
          unname()
      }
    )
  names(joint_dist[[1]]) <- round(theta, 10)

  # Calculate the probabilities of the conditional joint distributions
  # given the ability levels across all stages
  for (s in 2:n.stg) {
    # The numbers of all possible (sub) pathways at the current stage
    n_path <- cum_path[s - 1]

    # Empty list to contain the cut scores used to route a test taker to next module
    cut4route[[s - 1]] <- vector("list", n_path)

    # A data.frame with the pathway by the current stage (first column)
    # and the modules (the current stage (second column)
    path_pre <-
      data.frame(pathway[, 1:s]) %>%
      tidyr::unite(col = "path", 1:(s - 1), sep = "_")

    # Unique pathways by the current stage
    uni_path_pre <- unique(path_pre$path)

    # Empty to contain the conditional joint distributions at the current stage
    jdist_temp <- vector("list", n_path)

    # Loop over the all possible (sub) pathways
    for (i in 1:n_path) {
      # Find cut scores for routing a next module
      # Firstly, find all unique next modules that will be assigned at the current pathway
      nmod_next <-
        path_pre %>%
        dplyr::filter(.data$path %in% uni_path_pre[i]) %>%
        dplyr::pull(2) %>%
        unique()

      # Secondly, find the indices which modules will be actually used among the all possible next modules
      idx_next <- match(nmod_next, unique(pathway[, s]))

      # Lastly, find the cut scores to be used to assign the next modules
      cut4route[[s - 1]][[i]] <- cut_score[[s - 1]][dplyr::lag(idx_next)[-1]]

      # Find the path labels to determine which module will be assigned to each
      # ability level based on the cut scores
      path_label[[s - 1]][[i]] <-
        give_path(
          score = unlist(eq_theta[[s - 1]][, i]),
          cut_sc = cut4route[[s - 1]][[i]]
        )$path

      # Cumulative maximum sum scores by the current stages
      max_score <- score_maxcum[s]

      # Number of all unique next modules (a.k.a. all unique possible routing numbers)
      n_route <- length(nmod_next)

      # All possible observed scores by the current stage
      score_cur <- 0:score_maxcum[s - 1]

      # Possible observed score for the next module
      score_next <- score_mod[[s]]

      # The conditional joint distributions given the ability levels by the current stage
      jdist_cur <-
        purrr::map(
          .x = joint_dist[[s - 1]],
          .f = ~ {
            .x[, i, drop = FALSE]
          }
        )

      # The conditional distribution of the next modules given the ability levels
      cdist_next <-
        purrr::map(
          .x = cdist_by_th,
          .f = ~ {
            .x[, nmod_next, drop = FALSE]
          }
        )

      # Column bind of the above two distributions given the ability levels
      cbind_dist <-
        purrr::map(
          .x = 1:length(theta),
          .f = ~ {
            bind.fill(
              List = list(jdist_cur[[.x]], cdist_next[[.x]]),
              type = "cbind",
              fill = 0L
            )
          }
        )

      # Compute the joint conditional distributions
      jdist_temp[[i]] <-
        jdist(
          cbind_dist = cbind_dist, path_lb = path_label[[s - 1]][[i]],
          n_route = n_route, max_score = max_score, score_cur = score_cur,
          score_next = score_next, theta = theta
        )
    }

    # Add to the list
    joint_dist[[s]] <-
      purrr::map(
        .x = 1:length(theta),
        .f = ~ {
          do.call("cbind", lapply(jdist_temp, "[[", .x))
        }
      )
    names(joint_dist[[s]]) <- round(theta, 10)
  }

  ## -----------------------------------------------------
  # Evaluate the measurement precision of the MST panel
  ## -----------------------------------------------------
  # Calculate the means and variances of the conditional joint distributions
  # across all ability levels
  cmoment <-
    lapply(
      X = joint_dist[[s]],
      FUN = cal_moment,
      node = eq_theta[[s]]
    ) %>%
    do.call(what = "rbind")

  # Create a table containing the moments and
  # conditional biases and standard errors of measurement (SEMs)
  eval_tb <-
    data.frame(theta = theta, cmoment) %>%
    dplyr::mutate(
      bias = .data$mu - .data$theta,
      csem = sqrt(.data$sigma2)
    ) %>%
    tibble::remove_rownames()


  ## -----------------------------------------------------
  # Return the results
  ## -----------------------------------------------------
  rst <- list(
    panel.info = panel_data,
    item.by.mod = meta_mod,
    item.by.path = meta_path,
    eq.theta = eq_theta,
    cdist.by.mod = cdist_by_th,
    jdist.by.path = joint_dist,
    eval.tb = eval_tb
  )
  rst
}


# This function computes the conditional joint distributions of the current pathways
# given the ability levels
jdist <- function(cbind_dist, path_lb, n_route, max_score, score_cur,
                  score_next, theta) {
  # Empty list across all ability levels
  cj_dist_all <- vector("list", length(theta))

  # Loop over the ability levels
  for (j in 1:length(theta)) {
    # Matrix to contain the conditional joint distribution at a given ability level
    cj_dist <- array(0, c(max_score + 1, n_route))
    rownames(cj_dist) <- 0:max_score

    # Loop over the observed scores of the (sub) pathway by the the current stage
    for (i in 1:length(score_cur)) {
      # Compute the conditional joint distribution
      cj_dist[i:(i + max(score_next)), path_lb[i]] <-
        cj_dist[i:(i + max(score_next)), path_lb[i]] +
        (cbind_dist[[j]][1:length(score_next), (path_lb[i] + 1)] * cbind_dist[[j]][i, 1])
    }

    # Add to the list
    cj_dist_all[[j]] <- cj_dist
  }

  # Return the results
  cj_dist_all
}

# This function determines the next path that each test taker with
# a given ability level needs to take
give_path <- function(score, cut_sc) {
  # Number of theta scores
  nscore <- length(score)

  # Number of categories (or labels) to which each theta will be given
  ncats <- length(cut_sc) + 1

  # A vector of labels to be given to each test taker
  labels <- 1:ncats

  # Create pathway vectors
  path <-
    cut(x = score, breaks = c(-Inf, cut_sc, Inf), labels = labels) %>%
    as.numeric()

  # Return the results
  rst <- list(path = path, score = score, cut.sc = cut_sc)
  rst
}