Fix and test brm_summary()

R/brm_summary.R

@@ -2,6 +2,8 @@
 #' @export
 #' @family results
 #' @description Summarize a basic MMRM model fit.
+#' @details Currently assumes the response variable is `CHG`
+#'   (change from baseline) and not `AVAL` (raw response).
 #' @return A `tibble` with summary statistics of the marginal posterior.
 #' @inheritParams brm_formula
 #' @param model Fitted `brms` model object from [brm_model()]
@@ -91,7 +93,7 @@ brm_summary <- function(
     data = data,
     emmeans_response = emmeans_response
   )
-  table_diff <- brm_summary_diff_change(
+  table_diff <- brm_summary_diff(
     data = data,
     emmeans_response = emmeans_response,
     group = group,
@@ -117,23 +119,25 @@ brm_summary_response <- function(data, emmeans_response) {
   out
 }
 
-brm_summary_diff_change <- function(
+brm_summary_diff <- function(
   data,
   emmeans_response,
   group,
   time,
   nuisance,
   control
 ) {
-  contrasts_diff <- list()
   reference <- tibble::as_tibble(emmeans_response)
+  contrasts_diff <- list()
   for (level_group in setdiff(sort(unique(reference[[group]])), control)) {
     for (level_time in sort(unique(reference[[time]]))) {
       contrast_treatment <- as.integer(
-        reference[[group]] == group & reference[[time]] == level_time
+        (reference[[group]] == level_group) &
+          (reference[[time]] == level_time)
       )
       contrast_control <- as.integer(
-        reference[[group]] == control & reference[[time]] == level_time
+        (reference[[group]] == control) & 
+          (reference[[time]] == level_time)
       )
       contrast <- contrast_treatment - contrast_control
       contrasts_diff[[length(contrasts_diff) + 1L]] <- contrast

tests/testthat/test-brm_summary.R

@@ -0,0 +1,62 @@
+test_that("brm_summary()", {
+  set.seed(0L)
+  sim <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L,
+    hyper_beta = 1,
+    hyper_sigma = 1,
+    hyper_correlation = 1
+  )
+  data <- sim$data
+  formula <- brm_formula(
+    response = "response",
+    group = "group",
+    time = "time",
+    patient = "patient",
+    effect_base = FALSE,
+    interaction_base = FALSE
+  )
+  tmp <- utils::capture.output(
+    suppressMessages(
+      suppressWarnings(
+        model <- brm_model(
+          data = data,
+          formula = formula,
+          chains = 1,
+          iter = 100,
+          refresh = 0
+        )
+      )
+    )
+  )
+  out <- brm_summary(
+    model = model,
+    group = "group",
+    time = "time",
+    patient = "patient",
+    control = 1
+  )
+  cols <- c(
+    "time",
+    "group",
+    "response_mean",
+    "response_lower",
+    "response_upper",
+    "diff_mean",
+    "diff_lower",
+    "diff_upper"
+  )
+  expect_equal(sort(colnames(out)), sort(cols))
+  expect_equal(nrow(out), 8L)
+  expect_equal(as.integer(out$time), rep(seq_len(4), each = 2))
+  expect_equal(as.integer(out$group), rep(seq_len(2), times = 4))
+  expect_true(all(out$response_mean > out$response_lower))
+  expect_true(all(out$response_mean > out$response_upper))
+  expect_true(all(out$diff_mean > out$diff_lower))
+  expect_true(all(out$diff_mean > out$diff_upper))
+  expect_equal(
+    out$diff_mean[seq_len(4)],
+    out$response_mean[out$group == 2] - out$response_mean[out$group == 1]
+  )
+})