From 68f507ecdc68214087f5b9ab3ef689e020205ce0 Mon Sep 17 00:00:00 2001 From: Jake Thompson Date: Tue, 25 Jun 2024 10:59:16 -0500 Subject: [PATCH 1/2] add ability to specify different thresholds for each attribute --- .lintr | 3 ++- R/model-evaluation.R | 5 +++-- R/reliability.R | 34 +++++++++++++++++++++++++++---- R/utils-reliability.R | 16 ++++++++++++--- man/model_evaluation.Rd | 2 +- man/reliability.Rd | 6 +++++- tests/testthat/test-reliability.R | 17 ++++++++++++++-- 7 files changed, 69 insertions(+), 14 deletions(-) diff --git a/.lintr b/.lintr index bc0a699..24eb3c0 100644 --- a/.lintr +++ b/.lintr @@ -1,5 +1,6 @@ linters: linters_with_defaults( - indentation_linter = NULL + indentation_linter = NULL, + return_linter = NULL ) exclusions: list( "R/stanmodels.R", diff --git a/R/model-evaluation.R b/R/model-evaluation.R index 32dbe8e..493a3a3 100644 --- a/R/model-evaluation.R +++ b/R/model-evaluation.R @@ -159,7 +159,8 @@ add_criterion <- function(x, criterion = c("loo", "waic"), overwrite = FALSE, #' @export #' @rdname model_evaluation -add_reliability <- function(x, overwrite = FALSE, save = TRUE) { +add_reliability <- function(x, threshold = 0.5, overwrite = FALSE, + save = TRUE) { model <- check_model(x, required_class = "measrfit", name = "x") overwrite <- check_logical(overwrite, name = "overwrite") save <- check_logical(save, name = "force_save") @@ -168,7 +169,7 @@ add_reliability <- function(x, overwrite = FALSE, save = TRUE) { run_reli <- length(model$reliability) == 0 || overwrite if (run_reli) { - model$reliability <- reliability(model) + model$reliability <- reliability(model, threshold = threshold, force = TRUE) } # re-save model object (if applicable) diff --git a/R/reliability.R b/R/reliability.R index fe44539..f370cfd 100644 --- a/R/reliability.R +++ b/R/reliability.R @@ -28,7 +28,11 @@ reliability <- function(model, ...) { #' #' @param threshold For `map_reliability`, the threshold applied to the #' attribute-level probabilities for determining the binary attribute -#' classifications. +#' classifications. Should be a numeric vector of length 1 (the same threshold +#' is applied to all attributes), or length equal to the number of attributes. +#' If a named vector is supplied, names should match the attribute names in the +#' Q-matrix used to estimate the model. If unnamed, thresholds should be in the +#' order the attributes were defined in the Q-matrix. #' #' @details The pattern-level reliability (`pattern_reliability`) statistics are #' described in Cui et al. (2012). Attribute-level classification reliability @@ -71,6 +75,31 @@ reliability <- function(model, ...) { #' #' reliability(rstn_mdm_lcdm) reliability.measrdcm <- function(model, ..., threshold = 0.5, force = FALSE) { + threshold <- check_double(threshold, lb = 0, ub = 1, inclusive = FALSE, + name = "threshold") + force <- check_logical(force, name = "force") + + att_names <- colnames(dplyr::select(model$data$qmatrix, -"item_id")) + if (length(threshold) == 1) { + threshold <- rep(threshold, times = length(att_names)) %>% + rlang::set_names(att_names) + } else if (length(threshold) == length(att_names)) { + if (is.null(names(threshold))) { + threshold <- rlang::set_names(threshold, att_names) + } else if (!all(names(threshold) %in% att_names)) { + bad_names <- setdiff(names(threshold), att_names) + rlang::abort( + message = glue::glue("Unknown attribute names provided: ", + "{paste(bad_names, collapse = ', ')}") + ) + } + } else { + rlang::abort( + message = glue::glue("`threshold` must be of length 1 or length ", + "{length(att_names)} (the number of attributes).") + ) + } + if ((!is.null(model$reliability) && length(model$reliability) > 0) && !force) { return(model$reliability) @@ -78,9 +107,6 @@ reliability.measrdcm <- function(model, ..., threshold = 0.5, force = FALSE) { # coerce model into a list of values required for reliability obj <- reli_list(model, threshold = threshold) - att_names <- model$data$qmatrix %>% - dplyr::select(-"item_id") %>% - colnames() tbl <- obj$acc p <- obj$prev diff --git a/R/utils-reliability.R b/R/utils-reliability.R index e52ba44..429be1b 100644 --- a/R/utils-reliability.R +++ b/R/utils-reliability.R @@ -159,9 +159,19 @@ reli_list <- function(model, threshold) { # map estimates binary_att <- attr_probs %>% - dplyr::mutate(dplyr::across(dplyr::everything(), - ~dplyr::case_when(.x >= threshold ~ 1L, - TRUE ~ 0L))) %>% + tibble::rowid_to_column(var = "resp_id") %>% + tidyr::pivot_longer(cols = -"resp_id", + names_to = "attribute", values_to = "probability") %>% + dplyr::left_join(tibble::enframe(threshold, name = "attribute", + value = "threshold"), + by = "attribute", + relationship = "many-to-one") %>% + dplyr::mutate(class = dplyr::case_when(.data$probability >= + .data$threshold ~ 1L, + .default = 0L)) %>% + dplyr::select("resp_id", "attribute", "class") %>% + tidyr::pivot_wider(names_from = "attribute", values_from = "class") %>% + dplyr::select(dplyr::all_of(names(threshold))) %>% as.matrix() %>% unname() %>% as.vector() diff --git a/man/model_evaluation.Rd b/man/model_evaluation.Rd index 79641e2..2ba86f7 100644 --- a/man/model_evaluation.Rd +++ b/man/model_evaluation.Rd @@ -17,7 +17,7 @@ add_criterion( r_eff = NA ) -add_reliability(x, overwrite = FALSE, save = TRUE) +add_reliability(x, threshold = 0.5, overwrite = FALSE, save = TRUE) add_fit( x, diff --git a/man/reliability.Rd b/man/reliability.Rd index b9c1407..530ab70 100644 --- a/man/reliability.Rd +++ b/man/reliability.Rd @@ -16,7 +16,11 @@ reliability(model, ...) \item{threshold}{For \code{map_reliability}, the threshold applied to the attribute-level probabilities for determining the binary attribute -classifications.} +classifications. Should be a numeric vector of length 1 (the same threshold +is applied to all attributes), or length equal to the number of attributes. +If a named vector is supplied, names should match the attribute names in the +Q-matrix used to estimate the model. If unnamed, thresholds should be in the +order the attributes were defined in the Q-matrix.} \item{force}{If reliability information has already been added to the model object with \code{\link[=add_reliability]{add_reliability()}}, should it be recalculated. Default is diff --git a/tests/testthat/test-reliability.R b/tests/testthat/test-reliability.R index d266a04..19be85a 100644 --- a/tests/testthat/test-reliability.R +++ b/tests/testthat/test-reliability.R @@ -1,5 +1,18 @@ test_that("dino reliability", { - dino_reli <- reliability(rstn_dino) + dino_reli <- reliability(rstn_dino, threshold = 0.5) + + # threshold errors + err <- rlang::catch_cnd(reliability(rstn_dino, + threshold = c("att1" = 0.5, + "asdf" = 0.8, + "test" = 0.7, + "att4" = 0.5, + "att5" = 0.3))) + expect_match(err$message, "Unknown attribute names") + err <- rlang::catch_cnd(reliability(rstn_dino, + threshold = c("att1" = 0.5, + "asdf" = 0.8))) + expect_match(err$message, "must be of length 1 or length 5") # list naming expect_equal(names(dino_reli), c("pattern_reliability", "map_reliability", @@ -35,7 +48,7 @@ test_that("reliability can be added to model object", { err <- rlang::catch_cnd(measr_extract(dina_mod, "probability_reliability")) expect_match(err$message, "Reliability information must be added to a model") - dina_mod <- add_reliability(dina_mod) + dina_mod <- add_reliability(dina_mod, threshold = rep(0.5, 5)) expect_equal(names(dina_mod$reliability), c("pattern_reliability", "map_reliability", "eap_reliability")) From 1c4678e51b0b913afd816d44d8749731178f69ae Mon Sep 17 00:00:00 2001 From: Jake Thompson Date: Tue, 25 Jun 2024 10:59:26 -0500 Subject: [PATCH 2/2] update docs --- man/loo.measrfit.Rd | 11 ++++++----- man/model_evaluation.Rd | 11 ++++++----- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/man/loo.measrfit.Rd b/man/loo.measrfit.Rd index f1c4204..c6f473c 100644 --- a/man/loo.measrfit.Rd +++ b/man/loo.measrfit.Rd @@ -14,13 +14,14 @@ \item{r_eff}{Vector of relative effective sample size estimates for the likelihood (\code{exp(log_lik)}) of each observation. This is related to the relative efficiency of estimating the normalizing term in -self-normalizing importance sampling when using posterior draws obtained +self-normalized importance sampling when using posterior draws obtained with MCMC. If MCMC draws are used and \code{r_eff} is not provided then the reported PSIS effective sample sizes and Monte Carlo error estimates -will be over-optimistic. If the posterior draws are independent then -\code{r_eff=1} and can be omitted. The warning message thrown when \code{r_eff} is -not specified can be disabled by setting \code{r_eff} to \code{NA}. See the -\code{\link[loo:relative_eff]{relative_eff()}} helper functions for computing \code{r_eff}.} +can be over-optimistic. If the posterior draws are (near) independent then +\code{r_eff=1} can be used. \code{r_eff} has to be a scalar (same value is used +for all observations) or a vector with length equal to the number of +observations. The default value is 1. See the \code{\link[loo:relative_eff]{relative_eff()}} helper +functions for help computing \code{r_eff}.} \item{force}{If the LOO criterion has already been added to the model object with \code{\link[=add_criterion]{add_criterion()}}, should it be recalculated. Default is \code{FALSE}.} diff --git a/man/model_evaluation.Rd b/man/model_evaluation.Rd index 2ba86f7..e361110 100644 --- a/man/model_evaluation.Rd +++ b/man/model_evaluation.Rd @@ -56,13 +56,14 @@ file will not be updated.} \item{r_eff}{Vector of relative effective sample size estimates for the likelihood (\code{exp(log_lik)}) of each observation. This is related to the relative efficiency of estimating the normalizing term in -self-normalizing importance sampling when using posterior draws obtained +self-normalized importance sampling when using posterior draws obtained with MCMC. If MCMC draws are used and \code{r_eff} is not provided then the reported PSIS effective sample sizes and Monte Carlo error estimates -will be over-optimistic. If the posterior draws are independent then -\code{r_eff=1} and can be omitted. The warning message thrown when \code{r_eff} is -not specified can be disabled by setting \code{r_eff} to \code{NA}. See the -\code{\link[loo:relative_eff]{relative_eff()}} helper functions for computing \code{r_eff}.} +can be over-optimistic. If the posterior draws are (near) independent then +\code{r_eff=1} can be used. \code{r_eff} has to be a scalar (same value is used +for all observations) or a vector with length equal to the number of +observations. The default value is 1. See the \code{\link[loo:relative_eff]{relative_eff()}} helper +functions for help computing \code{r_eff}.} \item{method}{A vector of model fit methods to evaluate and add to the model object.}