Merge pull request #54 from JeffreyCHoover/aic-bic

AIC and BIC

Author: wjakethompson

NEWS.md

@@ -4,6 +4,8 @@
 
 * Added new `measrfit()` function for creating measrfit objects from *Stan* models that were not originally created with measr.
 
+* Added `aic()` and `bic()` functions for calculating the Akaike and Bayesian information criteria, respectively, for models estimated with `method = "optim"`.
+
 # measr 1.0.0
 
 ## New documentation

R/aic-bic.R

@@ -0,0 +1,35 @@
+aic <- function(model) {
+  model <- check_model(model, required_class = "measrfit", name = "model")
+  log_lik <- model$model$value
+
+  num_params <- model$model$par %>%
+    tibble::as_tibble() %>%
+    dplyr::mutate(param = names(model$model$par)) %>%
+    dplyr::filter(!grepl("pi", .data$param),
+                  !grepl("log_Vc", .data$param)) %>%
+    nrow() - 1
+
+  aic <- (-2 * log_lik) + (2 * num_params)
+
+  return(aic)
+}
+
+bic <- function(model) {
+  model <- check_model(model, required_class = "measrfit", name = "model")
+  log_lik <- model$model$value
+
+  num_params <- model$model$par %>%
+    tibble::as_tibble() %>%
+    dplyr::mutate(param = names(model$model$par)) %>%
+    dplyr::filter(!grepl("pi", .data$param),
+                  !grepl("log_Vc", .data$param)) %>%
+    nrow() - 1
+
+  n <- model$data$data %>%
+    dplyr::distinct(.data$resp_id) %>%
+    nrow()
+
+  bic <- (-2 * log_lik) + (log(n) * num_params)
+
+  return(bic)
+}

R/model-evaluation.R

@@ -10,7 +10,8 @@
 #' @inheritParams dcm2::calc_m2
 #' @param x A [measrfit] object.
 #' @param criterion A vector of criteria to calculate and add to the model
-#'   object.
+#'   object. Must be one of `"loo"` or `"waic"` for models estimated with MCMC,
+#'   or one of `"aic"` or `"bic"` for model estimated with the optimizer.
 #' @param method A vector of model fit methods to evaluate and add to the model
 #'   object.
 #' @param probs The percentiles to be computed by the [stats::quantile()]
@@ -117,16 +118,23 @@ NULL
 
 #' @export
 #' @rdname model_evaluation
-add_criterion <- function(x, criterion = c("loo", "waic"), overwrite = FALSE,
-                          save = TRUE, ..., r_eff = NA) {
+add_criterion <- function(x, criterion = c("loo", "waic"),
+                          overwrite = FALSE, save = TRUE, ..., r_eff = NA) {
   model <- check_model(x, required_class = "measrfit", name = "x")
-  if (model$method != "mcmc") {
+  if (any(model$method != "mcmc" && any(criterion %in% c("loo", "waic")))) {
     rlang::abort("error_bad_method",
-                 message = glue::glue("Model criteria are only available for ",
-                                      "models estimated with ",
+                 message = glue::glue("LOO and WAIC model criteria are only ",
+                                      "available for models estimated with ",
                                       "`method = \"mcmc\"`."))
+  } else if (any(model$method != "optim" &&
+                 any(criterion %in% c("aic", "bic")))) {
+    rlang::abort("error_bad_method",
+                 message = glue::glue("AIC and BIC model criteria are only ",
+                                      "available for models estimated with ",
+                                      "`method = \"optim\"`."))
   }
-  criterion <- rlang::arg_match(criterion, values = c("loo", "waic"),
+  criterion <- rlang::arg_match(criterion,
+                                values = c("loo", "waic", "aic", "bic"),
                                 multiple = TRUE)
   overwrite <- check_logical(overwrite, name = "overwrite")
   save <- check_logical(save, name = "force_save")
@@ -140,7 +148,7 @@ add_criterion <- function(x, criterion = c("loo", "waic"), overwrite = FALSE,
   }
   all_criteria <- c(new_criteria, redo_criteria)
 
-  if (length(all_criteria) > 0) {
+  if (length(all_criteria) > 0 && (model$method == "mcmc")) {
     log_lik_array <- loglik_array(model)
   }
 
@@ -150,6 +158,12 @@ add_criterion <- function(x, criterion = c("loo", "waic"), overwrite = FALSE,
   if ("waic" %in% all_criteria) {
     model$criteria$waic <- waic(log_lik_array)
   }
+  if ("aic" %in% all_criteria) {
+    model$criteria$aic <- aic(model)
+  }
+  if ("bic" %in% all_criteria) {
+    model$criteria$bic <- bic(model)
+  }
 
   # re-save model object (if applicable)
   if (!is.null(model$file) && length(all_criteria) > 0 && save) {

man/model_evaluation.Rd

@@ -293,3 +293,14 @@ test_that("mcmc requirements error", {
   expect_s3_class(err, "error_bad_method")
   expect_match(err$message, "`method = \"mcmc\"`")
 })
+
+test_that("optim requirements error", {
+  skip_on_cran()
+
+  mcmc_mod <- cmds_ecpe_lcdm
+  mcmc_mod$method <- "mcmc"
+
+  err <- rlang::catch_cnd(add_criterion(mcmc_mod, criterion = c("aic", "bic")))
+  expect_s3_class(err, "error_bad_method")
+  expect_match(err$message, "`method = \"optim\"`")
+})

tests/testthat/test-ecpe.R

@@ -116,8 +116,7 @@ test_that("loo and waic can be added to model", {
   expect_equal(names(loo_model$criteria), "loo")
   expect_s3_class(loo_model$criteria$loo, "psis_loo")
 
-  lw_model <- add_criterion(loo_model, criterion = c("loo", "waic"),
-                            overwrite = TRUE)
+  lw_model <- add_criterion(loo_model, overwrite = TRUE)
   expect_equal(names(lw_model$criteria), c("loo", "waic"))
   expect_s3_class(lw_model$criteria$loo, "psis_loo")
   expect_s3_class(lw_model$criteria$waic, "waic")

tests/testthat/test-mcmc.R

@@ -0,0 +1,40 @@
+test_that("add criterion error messages work", {
+  err <- rlang::catch_cnd(add_criterion("test"))
+  expect_s3_class(err, "error_bad_argument")
+  expect_match(err$message, "must be an object with class measrfit")
+
+
+  err <- rlang::catch_cnd(add_criterion(rstn_dino, criterion = "waic"))
+  expect_s3_class(err, "error_bad_method")
+  expect_match(err$message, "LOO and WAIC model criteria are only available")
+
+  err <- rlang::catch_cnd(add_criterion(rstn_dino, criterion = "waic"))
+  expect_s3_class(err, "error_bad_method")
+  expect_match(err$message, "LOO and WAIC model criteria are only available")
+
+  test_dino <- rstn_dino
+  test_dino$method <- "mcmc"
+  err <- rlang::catch_cnd(add_criterion(test_dino, criterion = "aic"))
+  expect_s3_class(err, "error_bad_method")
+  expect_match(err$message, "AIC and BIC model criteria are only available")
+
+  err <- rlang::catch_cnd(add_criterion(test_dino, criterion = "bic"))
+  expect_s3_class(err, "error_bad_method")
+  expect_match(err$message, "AIC and BIC model criteria are only available")
+})
+
+test_that("AIC works", {
+  rstn_dino <- add_criterion(rstn_dino, criterion = "aic")
+
+  exp_aic <- 37151.96
+
+  expect_equal(rstn_dino$criteria$aic, exp_aic, tolerance = .0001)
+})
+
+test_that("BIC works", {
+  rstn_dino <- add_criterion(rstn_dino, criterion = "bic")
+
+  exp_bic <- 37647.64
+
+  expect_equal(rstn_dino$criteria$bic, exp_bic, tolerance = .0001)
+})

tests/testthat/test-model-evaluation.R

@@ -0,0 +1,20 @@
+test_that("aic works", {
+  num_params <- 101
+  log_lik <- -18474.98
+
+  exp_aic <- (-2 * log_lik) + (2 * num_params)
+  aic_val <- aic(rstn_dino)
+
+  expect_equal(exp_aic, aic_val)
+})
+
+test_that("bic works", {
+  num_params <- 101
+  n <- 1000
+  log_lik <- -18474.98
+
+  exp_bic <- (-2 * log_lik) + (log(n) * num_params)
+  bic_val <- bic(rstn_dino)
+
+  expect_equal(exp_bic, bic_val)
+})