R/infitcutoff_plot.R

In easyRasch2: Psychometric Analysis with Rasch Measurement Theory

#' Plot Distribution of Simulated Infit and Outfit MSQ Values
#'
#' Visualises the distribution of simulation-based conditional item fit values
#' from \code{\link{RMitemInfitCutoff}}, optionally overlaying observed item fit
#' from the original data.
#'
#' Uses `ggdist::stat_dotsinterval()` (when `data` is not supplied) or
#' `ggdist::stat_dots()` (when `data` is supplied) with
#' `point_interval = "median_hdci"`. The outer `.width` follows
#' `simfit$hdci_width`, so the shaded interval matches the one
#' [RMitemInfit()] tabulates.
#'
#' @param simfit The return value of \code{\link{RMitemInfitCutoff}} (a list with
#'   components `results`, `item_cutoffs`, `actual_iterations`, `sample_n`, and
#'   `item_names`).
#' @param data Optional. A data.frame or matrix of item responses for computing
#'   and overlaying observed conditional item fit values. Items must be scored
#'   starting at 0 (non-negative integers). When provided, the plot includes
#'   orange diamond markers for the observed infit/outfit MSQ alongside the
#'   simulated distribution, plus segment summaries from the cutoff intervals.
#' @param statistic Character string. Either `"infit"` (default) to show only the
#'   infit panel, `"outfit"` to show only the outfit panel, or `"both"` to show
#'   infit and outfit side by side (requires the `patchwork` package when
#'   `data` is supplied).
#'
#' @return A `ggplot` object (or a `patchwork` object when `statistic = "both"`
#'   and `data` is supplied).
#'
#' @details
#' When `data` is **not** supplied, the function plots the simulated MSQ
#' distributions as dot-interval plots using `ggdist::stat_dotsinterval()` with
#' median and Highest Density Continuous Interval (HDCI) summaries, faceted by
#' statistic (InfitMSQ / OutfitMSQ).
#'
#' When `data` **is** supplied, the function:
#' \enumerate{
#'   \item Fits a Rasch / Partial Credit model by CML via
#'     `psychotools::pcmodel()` (a dichotomous item is a 2-category PCM)
#'     and computes observed conditional infit and outfit MSQ via
#'     `iarm::out_infit()`.
#'   \item Overlays observed fit values as orange diamond markers on the
#'     simulated distributions.
#'   \item Shows per-item cutoff intervals (from `simfit$item_cutoffs`) as
#'     black line segments, with thicker segments for the 66% HDCI range and
#'     black dots for the median.
#' }
#'
#' The `ggplot2`, `ggdist`, and optionally `patchwork` packages must be
#' installed (they are in Suggests, not Imports).
#'
#' @seealso \code{\link{RMitemInfitCutoff}}, \code{\link{RMitemInfit}}
#'
#' @importFrom rlang .data
#' @export
#'
#' @examples
#' \donttest{
#' if (requireNamespace("iarm", quietly = TRUE) &&
#'     requireNamespace("ggdist", quietly = TRUE) &&
#'     requireNamespace("ggplot2", quietly = TRUE)) {
#'   set.seed(42)
#'   sim_data <- as.data.frame(
#'     matrix(sample(0:1, 200 * 10, replace = TRUE), nrow = 200, ncol = 10)
#'   )
#'   colnames(sim_data) <- paste0("Item", 1:10)
#'
#'   # Run simulation
#'   cutoff_res <- RMitemInfitCutoff(sim_data, iterations = 100,
#'                                   parallel = FALSE, seed = 42)
#'
#'   # Simulated distribution only (infit + outfit faceted)
#'   RMitemInfitPlot(cutoff_res)
#'
#'   # With observed fit overlaid (infit only, the default)
#'   RMitemInfitPlot(cutoff_res, data = sim_data)
#'
#'   # Both infit and outfit panels side by side
#'   if (requireNamespace("patchwork", quietly = TRUE)) {
#'     RMitemInfitPlot(cutoff_res, data = sim_data, statistic = "both")
#'   }
#' }
#' }
RMitemInfitPlot <- function(simfit, data, statistic = "infit") {
  # --- Check required packages ------------------------------------------------
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop(
      "Package 'ggplot2' is required for RMitemInfitPlot() but is not installed.\n",
      "Install it with: install.packages(\"ggplot2\")",
      call. = FALSE
    )
  }
  if (!requireNamespace("ggdist", quietly = TRUE)) {
    stop(
      "Package 'ggdist' is required for RMitemInfitPlot() but is not installed.\n",
      "Install it with: install.packages(\"ggdist\")",
      call. = FALSE
    )
  }

  statistic <- match.arg(statistic, c("infit", "outfit", "both"))

  # --- Validate simfit --------------------------------------------------------
  required_names <- c(
    "results",
    "item_cutoffs",
    "actual_iterations",
    "sample_n",
    "item_names"
  )
  missing_names <- setdiff(required_names, names(simfit))
  if (length(missing_names) > 0L) {
    stop(
      "`simfit` is missing required components: ",
      paste(missing_names, collapse = ", "),
      ".\nExpected the return value of RMitemInfitCutoff().",
      call. = FALSE
    )
  }

  results_df <- simfit$results
  item_cutoffs <- simfit$item_cutoffs
  actual_iterations <- simfit$actual_iterations
  sample_n <- simfit$sample_n
  item_names <- simfit$item_names

  # Standard sample-size clause for the (complete-case) simulation sample.
  # `sample_n_total` / `sample_has_na` are absent in cutoff objects made by
  # older versions, so fall back to the plain count.
  sample_clause <- .n_caption(
    sample_n,
    if (is.null(simfit$sample_n_total)) sample_n else simfit$sample_n_total,
    if (isTRUE(simfit$sample_has_na)) "complete cases" else character()
  )

  # Item factor levels (reversed for plotting top-to-bottom)
  item_levels <- rev(item_names)

  # Outer interval width. The plot used to draw a fixed 0.1st-to-99.9th
  # percentile whisker regardless of how the cutoffs were computed, so it
  # could show a wider interval than the one RMitemInfit() tabulates. It now
  # follows the cutoff object, which keeps the plot and the table describing
  # the same interval. Objects made by older versions carry no width, and
  # `cutoff_method = "quantile"` fixes the interval at the 2.5th/97.5th
  # percentiles, so both fall back to .95.
  outer_width <- if (
    identical(simfit$cutoff_method, "hdci") && !is.null(simfit$hdci_width)
  ) {
    simfit$hdci_width
  } else {
    0.95
  }
  outer_lo <- (1 - outer_width) / 2
  outer_hi <- 1 - outer_lo

  # --- Compute per-item summary intervals for segment overlays ----------------
  lo_hi <- do.call(
    rbind,
    lapply(item_names, function(item) {
      sub <- results_df[results_df$Item == item, ]
      data.frame(
        Item = item,
        min_infit_msq = stats::quantile(sub$InfitMSQ, outer_lo, na.rm = TRUE),
        max_infit_msq = stats::quantile(sub$InfitMSQ, outer_hi, na.rm = TRUE),
        p66lo_infit_msq = stats::quantile(sub$InfitMSQ, 0.167, na.rm = TRUE),
        p66hi_infit_msq = stats::quantile(sub$InfitMSQ, 0.833, na.rm = TRUE),
        median_infit = stats::median(sub$InfitMSQ, na.rm = TRUE),
        min_outfit_msq = stats::quantile(sub$OutfitMSQ, outer_lo, na.rm = TRUE),
        max_outfit_msq = stats::quantile(sub$OutfitMSQ, outer_hi, na.rm = TRUE),
        p66lo_outfit_msq = stats::quantile(sub$OutfitMSQ, 0.167, na.rm = TRUE),
        p66hi_outfit_msq = stats::quantile(sub$OutfitMSQ, 0.833, na.rm = TRUE),
        median_outfit = stats::median(sub$OutfitMSQ, na.rm = TRUE),
        stringsAsFactors = FALSE,
        row.names = NULL
      )
    })
  )
  rownames(lo_hi) <- NULL

  # --- Case 1: no observed data, show simulation distribution only ------------
  if (missing(data)) {
    # Pivot results to long format
    infit_long <- data.frame(
      Item = results_df$Item,
      statistic = "InfitMSQ",
      Value = results_df$InfitMSQ,
      stringsAsFactors = FALSE
    )
    outfit_long <- data.frame(
      Item = results_df$Item,
      statistic = "OutfitMSQ",
      Value = results_df$OutfitMSQ,
      stringsAsFactors = FALSE
    )
    results_long <- rbind(infit_long, outfit_long)
    results_long$Item <- factor(results_long$Item, levels = item_levels)

    p <- ggplot2::ggplot(
      results_long,
      ggplot2::aes(
        x = .data$Value,
        y = .data$Item
      )
    ) +
      ggdist::stat_dotsinterval(
        ggplot2::aes(slab_fill = ggplot2::after_stat(.data$level)),
        quantiles = actual_iterations,
        point_interval = "median_hdci",
        layout = "weave",
        slab_color = NA,
        .width = c(0.66, outer_width)
      ) +
      ggplot2::labs(
        x = "Conditional MSQ",
        y = "Item",
        caption = er2_caption(paste0(
          "Results from ",
          actual_iterations,
          " simulated datasets. ",
          sample_clause,
          " per dataset."
        ))
      ) +
      ggplot2::scale_color_manual(
        values = scales::brewer_pal()(3)[-1],
        aesthetics = "slab_fill",
        guide = "none"
      ) +
      ggplot2::facet_wrap(~statistic, ncol = 2) +
      ggplot2::scale_x_continuous(
        breaks = seq(0.5, 1.5, 0.1),
        minor_breaks = NULL
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme(panel.spacing = ggplot2::unit(0.7, "cm")) +
      er2_axis_margins() +
      er2_plot_caption()

    return(p)
  }

  # --- Case 2: observed data supplied -----------------------------------------
  if (!requireNamespace("iarm", quietly = TRUE)) {
    stop(
      "Package 'iarm' is required to compute observed item fit but is not installed.\n",
      "Install it with: install.packages(\"iarm\")",
      call. = FALSE
    )
  }
  validate_response_data(data)

  # Respondents with no responses at all break the CML fit below
  # (psychotools errors on all-NA rows); drop them, as in RMitemInfit().
  data <- .drop_empty_respondents(data)

  # rgl workaround
  old_rgl <- getOption("rgl.useNULL")
  options(rgl.useNULL = TRUE)
  on.exit(options(rgl.useNULL = old_rgl), add = TRUE)

  # CML fit via psychotools (a dichotomous item is a 2-category PCM);
  # iarm::out_infit() accepts pcmodel objects and matches the former
  # eRm::RM()/PCM() route to ~1e-6 (same pairing as RMitemInfitCutoff()).
  model_fit <- psychotools::pcmodel(data)
  cfit <- iarm::out_infit(model_fit)

  observed_df <- data.frame(
    Item = names(data),
    observed_infit = cfit$Infit,
    observed_outfit = cfit$Outfit,
    stringsAsFactors = FALSE
  )

  # --- Build infit data -------------------------------------------------------
  infit_sim <- data.frame(
    Item = results_df$Item,
    Value = results_df$InfitMSQ,
    stringsAsFactors = FALSE
  )
  infit_sim <- merge(
    infit_sim,
    observed_df[, c("Item", "observed_infit")],
    by = "Item",
    sort = FALSE
  )
  infit_sim$Item <- factor(infit_sim$Item, levels = item_levels)

  lo_hi$Item_f <- factor(lo_hi$Item, levels = item_levels)

  caption_text <- er2_caption(paste0(
    "Results from ",
    actual_iterations,
    " simulated datasets. ",
    sample_clause,
    " per dataset.\n",
    "Orange dots indicate observed conditional item fit. ",
    "Black dots indicate median fit from simulations."
  ))

  infit_p <- ggplot2::ggplot(
    infit_sim,
    ggplot2::aes(
      x = .data$Value,
      y = .data$Item
    )
  ) +
    ggdist::stat_dots(
      ggplot2::aes(slab_fill = ggplot2::after_stat(.data$level)),
      quantiles = actual_iterations,
      layout = "weave",
      slab_color = NA,
      .width = c(0.666, outer_width)
    ) +
    ggplot2::geom_segment(
      data = lo_hi,
      ggplot2::aes(
        x = .data$min_infit_msq,
        xend = .data$max_infit_msq,
        y = .data$Item_f,
        yend = .data$Item_f
      ),
      color = "black",
      linewidth = 0.7
    ) +
    ggplot2::geom_segment(
      data = lo_hi,
      ggplot2::aes(
        x = .data$p66lo_infit_msq,
        xend = .data$p66hi_infit_msq,
        y = .data$Item_f,
        yend = .data$Item_f
      ),
      color = "black",
      linewidth = 1.2
    ) +
    ggplot2::geom_point(
      data = lo_hi,
      ggplot2::aes(
        x = .data$median_infit,
        y = .data$Item_f
      ),
      size = 3.6
    ) +
    ggplot2::geom_point(
      ggplot2::aes(x = .data$observed_infit),
      color = "sienna2",
      shape = 18,
      position = ggplot2::position_nudge(y = -0.1),
      size = 4
    ) +
    ggplot2::labs(
      x = "Conditional Infit MSQ",
      y = "Item"
    ) +
    ggplot2::scale_color_manual(
      values = scales::brewer_pal()(3)[-1],
      aesthetics = "slab_fill",
      guide = "none"
    ) +
    ggplot2::scale_x_continuous(
      breaks = seq(0.5, 1.5, 0.1),
      minor_breaks = NULL
    ) +
    ggplot2::theme_minimal() +
    ggplot2::theme(panel.spacing = ggplot2::unit(0.7, "cm")) +
    er2_axis_margins() +
    er2_plot_caption()

  # --- Build outfit data ------------------------------------------------------
  outfit_sim <- data.frame(
    Item = results_df$Item,
    Value = results_df$OutfitMSQ,
    stringsAsFactors = FALSE
  )
  outfit_sim <- merge(
    outfit_sim,
    observed_df[, c("Item", "observed_outfit")],
    by = "Item",
    sort = FALSE
  )
  outfit_sim$Item <- factor(outfit_sim$Item, levels = item_levels)

  outfit_p <- ggplot2::ggplot(
    outfit_sim,
    ggplot2::aes(
      x = .data$Value,
      y = .data$Item
    )
  ) +
    ggdist::stat_dots(
      ggplot2::aes(slab_fill = ggplot2::after_stat(.data$level)),
      quantiles = actual_iterations,
      layout = "weave",
      slab_color = NA,
      .width = c(0.666, outer_width)
    ) +
    ggplot2::geom_segment(
      data = lo_hi,
      ggplot2::aes(
        x = .data$min_outfit_msq,
        xend = .data$max_outfit_msq,
        y = .data$Item_f,
        yend = .data$Item_f
      ),
      color = "black",
      linewidth = 0.7
    ) +
    ggplot2::geom_segment(
      data = lo_hi,
      ggplot2::aes(
        x = .data$p66lo_outfit_msq,
        xend = .data$p66hi_outfit_msq,
        y = .data$Item_f,
        yend = .data$Item_f
      ),
      color = "black",
      linewidth = 1.2
    ) +
    ggplot2::geom_point(
      data = lo_hi,
      ggplot2::aes(
        x = .data$median_outfit,
        y = .data$Item_f
      ),
      size = 3.6
    ) +
    ggplot2::geom_point(
      ggplot2::aes(x = .data$observed_outfit),
      color = "sienna2",
      shape = 18,
      position = ggplot2::position_nudge(y = -0.1),
      size = 4
    ) +
    ggplot2::labs(
      x = "Conditional Outfit MSQ",
      y = "Item",
      caption = caption_text
    ) +
    ggplot2::scale_color_manual(
      values = scales::brewer_pal()(3)[-1],
      aesthetics = "slab_fill",
      guide = "none"
    ) +
    ggplot2::scale_x_continuous(
      breaks = seq(0.5, 1.5, 0.1),
      minor_breaks = NULL
    ) +
    ggplot2::theme_minimal() +
    ggplot2::theme(panel.spacing = ggplot2::unit(0.7, "cm")) +
    er2_axis_margins() +
    er2_plot_caption()

  # --- Return the requested panel(s) ------------------------------------------
  if (statistic == "both") {
    if (!requireNamespace("patchwork", quietly = TRUE)) {
      stop(
        "Package 'patchwork' is required for statistic = \"both\" but is not installed.\n",
        "Install it with: install.packages(\"patchwork\")",
        call. = FALSE
      )
    }
    return(infit_p + outfit_p)
  } else if (statistic == "outfit") {
    return(outfit_p)
  } else {
    # default: infit only, add caption to infit panel
    infit_p <- infit_p +
      ggplot2::labs(caption = caption_text)
    return(infit_p)
  }
}

#' @rdname RMitemInfitPlot
#' @param ... Arguments passed on to \code{RMitemInfitPlot()}.
#' @details
#' `RMitemInfitCutoffPlot()` is a deprecated alias for `RMitemInfitPlot()`,
#' retained for backward compatibility with code written against easyRasch2
#' 0.8.0. It warns and forwards to `RMitemInfitPlot()`.
#' @export
RMitemInfitCutoffPlot <- function(...) {
  .Deprecated("RMitemInfitPlot")
  RMitemInfitPlot(...)
}