R/ppmc.R

Defines functions ppmc_item_pvalue pw_or ppmc_odds_ratio ppmc_conditional_probs ppmc_item_fit ppmc_rawscore_chisq ppmc_model_fit fit_ppmc

Documented in fit_ppmc

#' Posterior predictive model checks for assessing model fit
#'
#' For models estimated with `method = "mcmc"`, use the posterior distributions
#' to compute expected distributions for fit statistics and compare to values
#' in the observed data.
#'
#' @param model A [measrfit] object.
#' @param ndraws The number of posterior draws to base the checks on. Must be
#'   less than or equal to the total number of posterior draws retained in the
#'   estimated model. If `NULL` (the default) the total number from the
#'   estimated model is used.
#' @param probs The percentiles to be computed by the [stats::quantile()]
#'   function for summarizing the posterior distributions of the specified fit
#'   statistics.
#' @param return_draws Proportion of posterior draws for each specified fit
#'   statistic to be returned. This does not affect the calculation of the
#'   posterior predictive checks, but can be useful for visualizing the fit
#'   statistics. For example, if `ndraws = 500`, `return_draws = 0.2`, and
#'   `model_fit = "raw_score"`, then the raw score chi-square will be computed
#'   500 times (once for each draw) and 100 of those values (0.2 * 500) will be
#'   returned. If `0` (the default), only summaries of the posterior are
#'   returned (no individual samples).
#' @param model_fit The posterior predictive model checks to compute for an
#'   evaluation of model-level fit. If `NULL`, no model-level checks are
#'   computed. See details.
#' @param item_fit The posterior predictive model checks to compute for an
#'   evaluation of item-level fit. If `NULL`, no item-level checks are computed.
#'   Multiple checks can be provided in order to calculate more than one check
#'   simultaneously (e.g., `item_fit = c("conditional_prob", "odds_ratio")`).
#'   See details.
#' @param force If all requested \acronym{PPMC}s have already been added to the
#'   model object using [add_fit()], should they be recalculated. Default is
#'   `FALSE`.
#'
#' @details
#' Posterior predictive model checks (PPMCs) use the posterior distribution of
#' an estimated model to compute different statistics. This creates an expected
#' distribution of the given statistic, *if our estimated parameters are
#' correct*. We then compute the statistic in our observed data and compare the
#' observed value to the expected distribution. Observed values that fall
#' outside of the expected distributions indicate incompatibility between the
#' estimated model and the observed data.
#'
#' We currently support PPMCs at the model and item level. At the model level,
#' we calculate the expected raw score distribution (`model_fit = "raw_score"`),
#' as described by Thompson (2019) and Park et al. (2015).
#'
#' At the item level, we can calculate the conditional probability that a
#' respondent in each class provides a correct response (`item_fit =
#' "conditional_prob"`) as described by Thompson (2019) and Sinharay & Almond
#' (2007) or the overall proportion correct for an item (`item_fit = "pvalue"`),
#' as described by Thompson (2019). We can also calculate the odds ratio for
#' each pair of items (`item_fit = "odds_ratio"`) as described by Park et al.
#' (2015) and Sinharay et al. (2006).
#'
#' @return A list with two elements, "model_fit" and "item_fit". If either
#'   `model_fit = NULL` or `item_fit = NULL` in the function call, this will be
#'   a one-element list, with the null criteria excluded. Each list element, is
#'   itself a list with one element for each specified \acronym{PPMC} containing
#'   a [tibble][tibble::tibble-package]. For example if
#'   `item_fit = c("conditional_prob", "odds_ratio")`, the "item_fit" element
#'   will be a list of length two, where each element is a tibble containing the
#'   results of the \acronym{PPMC}. All tibbles follow the same general
#'   structure:
#'
#'   \itemize{
#'   \item `obs_{ppmc}`: The value of the relevant statistic in the observed
#'     data.
#'   \item `ppmc_mean`: The mean of the `ndraws` posterior samples calculated
#'     for the given statistic.
#'   \item Quantile columns: 1 column for each value of `probs`, providing the
#'     corresponding quantiles of the `ndraws` posterior samples calculated for
#'     the given statistic.
#'   \item `samples`: A list column, where each element contains a vector of
#'     length `(ndraws * return_draws)`, representing samples from the
#'     posterior distribution of the calculated statistic. This column is
#'     excluded if `return_draws = 0`.
#'   \item `ppp`: The posterior predictive p-value. This is the proportion of
#'     posterior samples for calculated statistic that are greater than the
#'     observed value. Values very close to 0 or 1 indicate incompatibility
#'     between the fitted model and the observed data.
#'   }
#'
#' @references Park, J. Y., Johnson, M. S., Lee, Y-S. (2015). Posterior
#'   predictive model checks for cognitive diagnostic models. *International
#'   Journal of Quantitative Research in Education, 2*(3-4), 244-264.
#'   \doi{10.1504/IJQRE.2015.071738}
#' @references Sinharay, S., & Almond, R. G. (2007). Assessing fit of cognitive
#'   diagnostic models. *Educational and Psychological Measurement, 67*(2),
#'   239-257. \doi{10.1177/0013164406292025}
#' @references Sinharay, S., Johnson, M. S., & Stern, H. S. (2006). Posterior
#'   predictive assessment of item response theory models. *Applied
#'   Psychological Measurement, 30*(4), 298-321.
#'   \doi{10.1177/0146621605285517}
#' @references Thompson, W. J. (2019). *Bayesian psychometrics for diagnostic
#'   assessments: A proof of concept* (Research Report No. 19-01). University
#'   of Kansas; Accessible Teaching, Learning, and Assessment Systems.
#'   \doi{10.35542/osf.io/jzqs8}
#'
#' @export
#' @examplesIf measr_examples()
#' mdm_dina <- measr_dcm(
#'   data = mdm_data, missing = NA, qmatrix = mdm_qmatrix,
#'   resp_id = "respondent", item_id = "item", type = "dina",
#'   method = "mcmc", seed = 63277, backend = "rstan",
#'   iter = 700, warmup = 500, chains = 2, refresh = 0
#' )
#'
#' fit_ppmc(mdm_dina, model_fit = "raw_score", item_fit = NULL)
fit_ppmc <- function(model, ndraws = NULL, probs = c(0.025, 0.975),
                     return_draws = 0,
                     model_fit = c("raw_score"),
                     item_fit = c("conditional_prob", "odds_ratio", "pvalue"),
                     force = FALSE) {
  model <- check_model(model, required_class = "measrdcm", name = "object")
  total_draws <- posterior::ndraws(posterior::as_draws(model))
  ndraws <- check_integer(ndraws, lb = 1, ub = total_draws,
                          name = "ndraws", allow_null = TRUE)
  probs <- check_double(probs, lb = 0, ub = 1, inclusive = TRUE, name = "probs")
  return_draws <- check_double(return_draws, lb = 0, ub = 1, inclusive = TRUE,
                               name = "return_draws")
  model_fit <- if (!is.null(model_fit)) {
    rlang::arg_match(model_fit, multiple = TRUE)
  }
  item_fit <- if (!is.null(item_fit)) {
    rlang::arg_match(item_fit, multiple = TRUE)
  }

  # check for existing fit analyses
  check_ppmc <- existing_ppmc_check(model = model, method = "ppmc",
                                    dots = list(ndraws = ndraws,
                                                probs = probs,
                                                return_draws = return_draws,
                                                model_fit = model_fit,
                                                item_fit = item_fit),
                                    overwrite = force)
  if (!check_ppmc$run && !force) {
    requested <- list(model_fit = model_fit,
                      item_fit = item_fit)
    requested[sapply(requested, is.null)] <- NULL

    res <- mapply(
      function(x, nm, model) {
        element <- model$fit$ppmc[[nm]][x]
        rlang::set_names(element, nm = x)
        return(element)
      },
      requested, names(requested),
      MoreArgs = list(model = model), SIMPLIFY = FALSE, USE.NAMES = TRUE
    )

    return(res)
  }


  if (length(model_fit) == 0 && length(item_fit) == 0) return(list())

  clean_qmatrix <- model$data$qmatrix %>%
    dplyr::select(-"item_id") %>%
    dplyr::rename_with(~glue::glue("att{1:(ncol(model$data$qmatrix) - 1)}"))
  stan_data <- create_stan_data(dat = model$data$data, qmat = clean_qmatrix,
                                type = model$type)

  stan_draws <- get_mcmc_draws(model, ndraws = ndraws)

  stan_pars <- create_stan_gqs_params(backend = model$backend,
                                      draws = stan_draws)
  stan_pars$data <- stan_data

  # compile model -----
  stan_mod <- create_stan_function(backend = model$backend,
                                   method = "gqs",
                                   code = gqs_script(full_data = TRUE),
                                   precompiled = stanmodels$gqs_ppmc,
                                   pars = stan_pars,
                                   silent = 2)
  out <- utils::capture.output( #nolint
    gqs_model <- do.call(stan_mod$func, stan_mod$pars)
  )
  all_draws <- posterior::as_draws_array(gqs_model)

  # organize ppmc data
  post_data <- posterior::subset_draws(all_draws, variable = "y_rep") %>%
    posterior::as_draws_df() %>%
    tibble::as_tibble() %>%
    tidyr::pivot_longer(cols = -c(".chain", ".iteration", ".draw")) %>%
    tidyr::separate_wider_regex(
      cols = "name",
      patterns = c("y_rep\\[", obs = "\\d+", "\\]")) %>%
    dplyr::mutate(obs = as.integer(.data$obs)) %>%
    dplyr::left_join(model$data$data %>%
                       dplyr::mutate(obs = seq_len(dplyr::n()),
                                     resp = as.integer(.data$resp_id),
                                     item = as.integer(.data$item_id)) %>%
                       dplyr::select("obs", "resp", "item"),
                     by = "obs", relationship = "many-to-one")

  # calculate raw score distribution
  model_level_fit <- if (!is.null(model_fit)) {
    ppmc_model_fit(model = model,
                   post_data = post_data,
                   probs = probs,
                   return_draws = return_draws,
                   type = check_ppmc$args$model_fit)
  } else {
    NULL
  }

  item_level_fit <- if (!is.null(item_fit)) {
    resp_prob <- extract_class_probs(model = gqs_model,
                                     attr = ncol(clean_qmatrix),
                                     method = model$method)
    pi_draws <- posterior::subset_draws(stan_draws, variable = "pi")

    ppmc_item_fit(model = model,
                  post_data = post_data,
                  attr = ncol(clean_qmatrix),
                  resp_prob = resp_prob,
                  pi_draws = pi_draws,
                  probs = probs,
                  return_draws = return_draws,
                  type = check_ppmc$args$item_fit)
  } else {
    NULL
  }

  ret_list <- if (is.null(model$fit$ppmc)) {
    list(model_fit = model_level_fit,
         item_fit = item_level_fit)
  } else {
    model$fit$ppmc$model_fit[names(model_level_fit)] <- NULL
    model$fit$ppmc$item_fit[names(item_level_fit)] <- NULL

    utils::modifyList(model$fit$ppmc,
                      list(model_fit = model_level_fit,
                           item_fit = item_level_fit))
  }

  ret_list <- list(model_fit = ret_list$model_fit[model_fit],
                   item_fit = ret_list$item_fit[item_fit])

  ret_list[sapply(ret_list, is.null)] <- NULL

  return(ret_list)
}


ppmc_model_fit <- function(model, post_data, probs, return_draws, type) {
  raw_score <- if ("raw_score" %in% type) {
    ppmc_rawscore_chisq(model = model, post_data = post_data,
                        probs = probs, return_draws = return_draws)
  } else {
    NULL
  }

  mod_res <- list(raw_score = raw_score)
  mod_res[sapply(mod_res, is.null)] <- NULL

  return(mod_res)
}

ppmc_rawscore_chisq <- function(model, post_data, probs, return_draws) {
  raw_score_post <- dtplyr::lazy_dt(post_data) %>%
    dplyr::summarize(raw_score = sum(.data$value), .by = c("resp", ".draw")) %>%
    dplyr::count(.data$.draw, .data$raw_score) %>%
    tibble::as_tibble() %>%
    tidyr::complete(.data$.draw, raw_score = 0:nrow(model$data$qmatrix),
                    fill = list(n = 0L))

  exp_raw_scores <- raw_score_post %>%
    dplyr::summarize(exp_resp = mean(.data$n), .by = "raw_score") %>%
    dplyr::mutate(exp_resp = sapply(.data$exp_resp,
                                    function(.x) max(.x, 0.0001)))

  chisq_ppmc <- raw_score_post %>%
    dplyr::left_join(exp_raw_scores, by = c("raw_score"),
                     relationship = "many-to-one") %>%
    dplyr::mutate(piece = ((.data$n - .data$exp_resp) ^ 2) / .data$exp_resp) %>%
    dplyr::summarize(chisq = sum(.data$piece), .by = ".draw")

  chisq_obs <- model$data$data %>%
    dplyr::summarize(raw_score = sum(.data$score), .by = "resp_id") %>%
    dplyr::count(.data$raw_score) %>%
    tidyr::complete(raw_score = 0:nrow(model$data$qmatrix),
                    fill = list(n = 0L)) %>%
    dplyr::left_join(exp_raw_scores, by = "raw_score",
                     relationship = "one-to-one") %>%
    dplyr::mutate(piece = ((.data$n - .data$exp_resp) ^ 2) / .data$exp_resp) %>%
    dplyr::summarize(chisq = sum(.data$piece)) %>%
    dplyr::pull("chisq")

  raw_score_res <- tibble::tibble(obs_chisq = chisq_obs,
                                  ppmc_mean = mean(chisq_ppmc$chisq),
                                  bounds = list(
                                    tibble::as_tibble_row(
                                      stats::quantile(chisq_ppmc$chisq,
                                                      probs = probs,
                                                      na.rm = TRUE)
                                    )
                                  ),
                                  ppp = mean(chisq_ppmc$chisq > chisq_obs)) %>%
    tidyr::unnest("bounds")

  if (return_draws > 0) {
    raw_score_res <- raw_score_res %>%
      dplyr::mutate(
        rawscore_samples = list(raw_score_post %>%
                                  tidyr::nest(raw_scores = -".draw") %>%
                                  dplyr::slice_sample(prop = return_draws) %>%
                                  dplyr::select(-".draw")),
        chisq_samples = list(chisq_ppmc %>%
                         dplyr::slice_sample(prop = return_draws) %>%
                         dplyr::pull("chisq")),
        .before = "ppp")
  }

  return(raw_score_res)
}

ppmc_item_fit <- function(model, post_data, attr, resp_prob, pi_draws, probs,
                          return_draws, type) {
  cond_prob <- if ("conditional_prob" %in% type) {
    ppmc_conditional_probs(model = model, attr = attr, resp_prob = resp_prob,
                           pi_draws = pi_draws, probs = probs,
                           return_draws = return_draws)
  } else {
    NULL
  }
  odds_ratio <- if ("odds_ratio" %in% type) {
    ppmc_odds_ratio(model = model, post_data = post_data, probs, return_draws)
  } else {
    NULL
  }
  pvalue <- if ("pvalue" %in% type) {
    ppmc_item_pvalue(model = model, post_data = post_data, probs, return_draws)
  } else {
    NULL
  }

  item_res <- list(conditional_prob = cond_prob,
                   odds_ratio = odds_ratio,
                   pvalue = pvalue)
  item_res[sapply(item_res, is.null)] <- NULL

  return(item_res)
}

ppmc_conditional_probs <- function(model, attr, resp_prob, pi_draws, probs,
                                   return_draws) {
  all_profiles <- profile_labels(attributes = attr)

  obs_class <- resp_prob %>%
    dplyr::mutate(dplyr::across(dplyr::where(posterior::is_rvar),
                                ~lapply(.x,
                                        function(x) {
                                          posterior::as_draws_df(x) %>%
                                            tibble::as_tibble()
                                        })
                                )) %>%
    tidyr::unnest(-"resp_id", names_sep = "_") %>%
    dplyr::select("resp_id",
                  dplyr::all_of(paste0(all_profiles$class[1], "_",
                                       c(".chain", ".iteration", ".draw"))),
                  dplyr::ends_with("_x")) %>%
    dplyr::rename_with(function(x) {
      x <- sub("_x", "", x)
      x <- sub("\\[[0-9,]*\\]_", "", x)
    }) %>%
    tidyr::pivot_longer(cols = -c("resp_id", ".chain", ".iteration", ".draw"),
                        names_to = "class_label",
                        values_to = "prob") %>%
    dtplyr::lazy_dt() %>%
    dplyr::mutate(max_class = .data$prob == max(.data$prob),
                  .by = c(".draw", "resp_id")) %>%
    dplyr::filter(.data$max_class) %>%
    dplyr::group_by(.data$.draw, .data$resp_id) %>%
    dplyr::slice_sample(n = 1) %>%
    dplyr::ungroup() %>%
    dplyr::left_join(all_profiles, by = c("class_label" = "class"),
                     relationship = "many-to-one") %>%
    dplyr::select(".draw", "resp_id", class = "class_id") %>%
    tibble::as_tibble()

  obs_cond_pval <- model$data$data %>%
    dplyr::mutate(resp_id = as.integer(.data$resp_id),
                  item_id = as.integer(.data$item_id)) %>%
    dplyr::left_join(obs_class, by = "resp_id",
                     relationship = "many-to-many") %>%
    dplyr::summarize(obs_cond_pval = mean(.data$score),
                     .by = c("item_id", "class")) %>%
    dplyr::arrange(.data$item_id, .data$class)

  cond_pval_res <- pi_draws %>%
    posterior::as_draws_df() %>%
    tibble::as_tibble() %>%
    tidyr::pivot_longer(-c(".chain", ".iteration", ".draw"),
                        names_to = c("i", "c"), values_to = "pi",
                        names_pattern = "pi\\[([0-9]*),([0-9]*)\\]",
                        names_transform = list(
                          i = ~as.integer(.x),
                          c = ~as.integer(.x)
                        )) %>%
    dplyr::select(-c(".chain", ".iteration", ".draw")) %>%
    tidyr::nest(cond_pval = "pi") %>%
    dplyr::arrange(.data$i, .data$c) %>%
    dplyr::left_join(obs_cond_pval, by = c("i" = "item_id",
                                           "c" = "class"),
                     relationship = "one-to-one") %>%
    dplyr::mutate(ppmc_mean = vapply(.data$cond_pval, function(.x) mean(.x$pi),
                                     double(1)),
                  bounds = lapply(.data$cond_pval,
                                  function(.x, probs) {
                                    tibble::as_tibble_row(
                                      stats::quantile(.x$pi,
                                                      probs = probs,
                                                      na.rm = TRUE)
                                    )
                                  },
                                  probs = probs),
                  ppp = mapply(function(exp, obs) {
                                 mean(exp$pi > obs)
                               },
                               .data$cond_pval, .data$obs_cond_pval)) %>%
    tidyr::unnest("bounds")

  if (return_draws > 0) {
    cond_pval_res <- cond_pval_res %>%
      dplyr::mutate(
        samples = lapply(.data$cond_pval,
                         function(.x) {
                           .x %>%
                             dplyr::slice_sample(prop = return_draws) %>%
                             dplyr::pull("pi")
                         }),
        .before = "ppp")
  }

  cond_pval_res <- cond_pval_res %>%
    dplyr::select(-"cond_pval") %>%
    dplyr::rename(item_id = "i", class_id = "c") %>%
    dplyr::left_join(all_profiles, by = "class_id",
                     relationship = "many-to-one") %>%
    dplyr::left_join(model$data$qmatrix %>%
                       dplyr::select(item = "item_id") %>%
                       dplyr::mutate(item_id = as.integer(.data$item)),
                     by = "item_id", relationship = "many-to-one") %>%
    dplyr::select(-"item_id", -"class_id") %>%
    dplyr::relocate("item", "class", .before = 1) %>%
    dplyr::rename(!!model$data$item_id := "item")

  return(cond_pval_res)
}

ppmc_odds_ratio <- function(model, post_data, probs, return_draws) {
  obs_or <- model$data$data %>%
    dplyr::mutate(resp_id = as.integer(.data$resp_id),
                  item_id = as.integer(.data$item_id)) %>%
    tidyr::pivot_wider(names_from = "item_id", values_from = "score") %>%
    dplyr::select(-"resp_id") %>%
    pw_or() %>%
    dplyr::rename(obs_or = "or")

  or_res <- post_data %>%
    dplyr::select(-"obs") %>%
    tidyr::pivot_wider(names_from = "item",
                       values_from = "value") %>%
    dplyr::select(-"resp") %>%
    tidyr::nest(dat = !dplyr::starts_with(".")) %>%
    dplyr::mutate(dat = lapply(.data$dat, pw_or)) %>%
    tidyr::unnest("dat") %>%
    tidyr::nest(samples = -c("item_1", "item_2")) %>%
    dplyr::left_join(obs_or, by = c("item_1", "item_2"),
                     relationship = "one-to-one") %>%
    dplyr::mutate(ppmc_mean = vapply(.data$samples, function(.x) mean(.x$or),
                                     double(1)),
                  bounds = lapply(.data$samples,
                                  function(.x, probs) {
                                    tibble::as_tibble_row(
                                      stats::quantile(.x$or,
                                                      probs = probs,
                                                      na.rm = TRUE)
                                    )
                                  },
                                  probs = probs),
                  ppp = mapply(function(exp, obs) {
                    mean(exp$or > obs)
                  },
                  .data$samples, .data$obs_or)) %>%
    tidyr::unnest("bounds")

  if (return_draws > 0) {
    or_res <- or_res %>%
      dplyr::relocate("samples", .before = "ppp") %>%
      dplyr::mutate(
        samples = lapply(.data$samples,
                         function(x) {
                           x %>%
                             dplyr::slice_sample(prop = return_draws) %>%
                             dplyr::pull("or")
                         }))
  } else {
    or_res <- dplyr::select(or_res, -"samples")
  }

  or_res <- or_res %>%
    dplyr::left_join(model$data$qmatrix %>%
                       dplyr::select(item = "item_id") %>%
                       dplyr::mutate(item_id = as.integer(.data$item)),
                     by = c("item_1" = "item_id"),
                     relationship = "many-to-one") %>%
    dplyr::left_join(model$data$qmatrix %>%
                       dplyr::select(item = "item_id") %>%
                       dplyr::mutate(item_id = as.integer(.data$item)),
                     by = c("item_2" = "item_id"),
                     relationship = "many-to-one") %>%
    dplyr::select(-"item_1", -"item_2") %>%
    dplyr::rename(item_1 = "item.x", item_2 = "item.y") %>%
    dplyr::relocate("item_1", "item_2", .before = 1)

  return(or_res)
}

pw_or <- function(dat) {
  dat <- as.matrix(dat)
  p <- ncol(dat)
  ind <- t(utils::combn(p, 2))
  nind <- nrow(ind)
  or <- numeric(nind)

  for (i in 1:nind) {
    xtab <- table(dat[, ind[i, 1]], dat[, ind[i, 2]])
    n00 <- xtab[1, 1]
    n01 <- xtab[1, 2]
    n10 <- xtab[2, 1]
    n11 <- xtab[2, 2]
    or[i] <- (n00 * n11) / (n01 * n10)
  }
  pwor <- tibble::as_tibble(ind, .name_repair = ~c("item_1", "item_2")) %>%
    dplyr::mutate(or = or)
  return(pwor)
}

ppmc_item_pvalue <- function(model, post_data, probs, return_draws) {
  obs_pvalue <- model$data$data %>%
    dplyr::mutate(item = as.numeric(.data$item_id)) %>%
    dplyr::summarize(obs_pvalue = mean(.data$score),
                     .by = c("item", "item_id"))

  pval_res <- post_data %>%
    dplyr::summarize(pvalue = mean(.data$value),
                     .by = c(".chain", ".iteration", ".draw", "item")) %>%
    tidyr::nest(samples = -"item") %>%
    dplyr::left_join(obs_pvalue, by = "item", relationship = "one-to-one") %>%
    dplyr::mutate(ppmc_mean = vapply(.data$samples,
                                     function(.x) mean(.x$pvalue),
                                     double(1)),
                  bounds = lapply(.data$samples,
                                  function(.x, probs) {
                                    tibble::as_tibble_row(
                                      stats::quantile(.x$pvalue,
                                                      probs = probs,
                                                      na.rm = TRUE)
                                    )
                                  },
                                  probs = probs),
                  ppp = mapply(function(exp, obs) {
                    mean(exp$pvalue > obs)
                  },
                  .data$samples, .data$obs_pvalue)) %>%
    tidyr::unnest("bounds")

  if (return_draws > 0) {
    pval_res <- pval_res %>%
      dplyr::relocate("samples", .before = "ppp") %>%
      dplyr::mutate(
        samples = lapply(.data$samples,
                         function(x) {
                           x %>%
                             dplyr::slice_sample(prop = return_draws) %>%
                             dplyr::pull("pvalue")
                         }))
  } else {
    pval_res <- dplyr::select(pval_res, -"samples")
  }

  pval_res %>%
    dplyr::select(-"item") %>%
    dplyr::rename(!!model$data$item_id := "item_id")
}
wjakethompson/measr documentation built on April 12, 2025, 9:46 p.m.