From 80940465b1fa0b0a151a895b514077cf91adb4f2 Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Fri, 2 Jun 2023 13:42:07 -0400
Subject: [PATCH 01/17] Subtract draws

---
 R/brm_summary.R    | 184 ++++++++++++++++++++++++++++-----------------
 R/utils_assert.R   |   2 +-
 man/brm_summary.Rd |  23 +++++-
 3 files changed, 137 insertions(+), 72 deletions(-)

diff --git a/R/brm_summary.R b/R/brm_summary.R
index 2d1bb4d5..e7722fb8 100644
--- a/R/brm_summary.R
+++ b/R/brm_summary.R
@@ -7,12 +7,21 @@
 #' @return A `tibble` with summary statistics of the marginal posterior.
 #' @inheritParams brm_formula
 #' @param model Fitted `brms` model object from [brm_model()].
-#' @param control Character of length 1, name of the control arm
-#'   in the `group` column in the data.
+#' @param control Element of the `group` column in the data which indicates
+#'   the control group for the purposes of calculating treatment differences.
+#' @param baseline Element of the `time` column in the data
+#'   which indicates the baseline time for the purposes of calculating
+#'   change from baseline. Ignored if `response_type = "change"`.
+#' @param response_type Character of length 1, `"response"` if the
+#'   response variable is the raw outcome variable (such as AVAL)
+#'   or `"change"` if the response variable is change from baseline
+#'   (e.g. CHG).
 #' @examples
 #' set.seed(0L)
 #' sim <- brm_simulate()
 #' data <- sim$data
+#' data$group <- paste("treatment", data$group)
+#' data$time <- paste("visit", data$time)
 #' formula <- brm_formula(
 #'   response = "response",
 #'   group = "group",
@@ -39,7 +48,9 @@
 #'   group = "group",
 #'   time = "time",
 #'   patient = "patient",
-#'   control = 1
+#'   control = "treatment 1",
+#'   baseline = "visit 1",
+#'   response_type = "response"
 #' )
 brm_summary <- function(
   model,
@@ -48,12 +59,22 @@ brm_summary <- function(
   time = "AVISIT",
   patient = "USUBJID",
   covariates = character(0),
-  control = "Placebo"
+  control = "Placebo",
+  baseline = "Baseline",
+  response_type = "change"
 ) {
   assert_chr(base, "base arg must be a nonempty character string")
   assert_chr(group, "group arg must be a nonempty character string")
   assert_chr(time, "time arg must be a nonempty character string")
   assert_chr(patient, "patient arg must be a nonempty character string")
+  assert_chr(
+    response_type,
+    "response_type arg must be a nonempty character string"
+  )
+  assert(
+    response_type %in% c("response", "change"),
+    message = "response_type must be either \"response\" or \"change\""
+  )
   assert_chr_vec(covariates, "covariates arg must be a character vector")
   assert(
     control,
@@ -62,6 +83,13 @@ brm_summary <- function(
     !anyNA(.),
     message = "control arg must be a length-1 non-missing atomic value"
   )
+  assert(
+    baseline,
+    is.atomic(.),
+    length(.) == 1L,
+    !anyNA(.),
+    message = "baseline arg must be a length-1 non-missing atomic value"
+  )
   assert(is.data.frame(model$data))
   data <- model$data
   assert(
@@ -85,82 +113,104 @@ brm_summary <- function(
     message = "control arg must be in data[[group]]"
   )
   nuisance <- c(base, patient, covariates)
-  emmeans_response <- emmeans::emmeans(
+  emmeans <- emmeans::emmeans(
     object = model,
-    specs = as.formula(sprintf("~%s:%s", time, group)),
+    specs = as.formula(sprintf("~%s:%s", group, time)),
     weights = "proportional",
     nuisance = nuisance
   )
-  table_response <- brm_summary_response(
-    data = data,
-    emmeans_response = emmeans_response
+  draws_response <- posterior::as_draws_df(as.mcmc(emmeans))
+  .chain <- draws_response[[".chain"]]
+  .iteration <- draws_response[[".iteration"]]
+  .draw <- draws_response[[".draw"]]
+  draws_response[[".chain"]] <- NULL
+  draws_response[[".iteration"]] <- NULL
+  draws_response[[".draw"]] <- NULL
+  colnames(draws_response) <- gsub(
+    pattern = sprintf("^%s ", group),
+    replacement = "",
+    x = colnames(draws_response)
   )
-  table_diff <- brm_summary_diff(
-    data = data,
-    emmeans_response = emmeans_response,
-    group = group,
-    time = time,
-    nuisance = nuisance,
-    control = control
+  colnames(draws_response) <- gsub(
+    pattern = sprintf(", %s ", time),
+    replacement = ", ",
+    x = colnames(draws_response)
   )
-  dplyr::left_join(
-    x = table_response,
-    y = table_diff,
-    by = c(group, time)
+  groups <- unique(gsub(",.*$", "", colnames(draws_response)))
+  times <- unique(gsub("^.*, ", "", colnames(draws_response)))
+  draws_response[[".chain"]] <- .chain
+  draws_response[[".iteration"]] <- .iteration
+  draws_response[[".draw"]] <- .draw
+  control <- as.character(control)
+  time <- as.character(time)
+  assert(
+    control %in% groups,
+    message = sprintf(
+      "control argument \"%s\" is not in one of the treatment groups: %s",
+      control,
+      paste(groups, collapse = ", ")
+    )
   )
+  if (response_type == "response") {
+    assert(
+      baseline %in% times,
+      message = sprintf(
+        "baseline argument \"%s\" is not in one of the time points: %s",
+        baseline,
+        paste(times, collapse = ", ")
+      )
+    )
+  }
+  if (response_type == "response") {
+    draws_change <- subtract_baseline(
+      draws = draws_response,
+      groups = groups,
+      times = times,
+      baseline = baseline
+    )
+    draws_diff <- subtract_control(
+      draws = draws_change,
+      groups = groups,
+      times = setdiff(times, baseline),
+      control = control
+    )
+  } else {
+    draws_diff <- subtract_control(
+      draws = draws_response,
+      groups = groups,
+      times = times,
+      control = control
+    )
+  }
+  
+  browser()
+  
 }
 
-brm_summary_response <- function(data, emmeans_response) {
-  out <- tibble::as_tibble(emmeans_response)
-  out$response_mean <- out$emmean
-  out$response_lower <- out$lower.HPD
-  out$response_upper <- out$upper.HPD
-  out$emmean <- NULL
-  out$lower.HPD <- NULL
-  out$upper.HPD <- NULL
+subtract_baseline <- function(draws, groups, times, baseline) {
+  out <- draws[, c(".chain", ".iteration", ".draw")]
+  for (group in groups) {
+    for (time in setdiff(times, baseline)) {
+      name1 <- marginal_name(group, baseline)
+      name2 <- marginal_name(group, time)
+      out[[name2]] <- draws[[name2]] - draws[[name1]]
+    }
+  }
   out
 }
 
-brm_summary_diff <- function(
-  data,
-  emmeans_response,
-  group,
-  time,
-  nuisance,
-  control
-) {
-  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]] == level_group) &
-          (reference[[time]] == level_time)
-      )
-      contrast_control <- as.integer(
-        (reference[[group]] == control) &
-          (reference[[time]] == level_time)
-      )
-      contrast <- contrast_treatment - contrast_control
-      contrasts_diff[[length(contrasts_diff) + 1L]] <- contrast
+subtract_control <- function(draws, groups, times, control) {
+  out <- draws[, c(".chain", ".iteration", ".draw")]
+  for (group in setdiff(groups, control)) {
+    for (time in times) {
+      name1 <- marginal_name(control, time)
+      name2 <- marginal_name(group, time)
+      out[[name2]] <- draws[[name2]] - draws[[name1]]
     }
   }
-  emmeans_diff <- emmeans::contrast(
-    emmeans_response,
-    method = contrasts_diff,
-    adjust = "none",
-    nuisance = nuisance
-  )
-  subset_reference <- reference[reference[[group]] != control, ]
-  out <- tibble::as_tibble(emmeans_diff)
-  out[[group]] <- subset_reference[[group]]
-  out[[time]] <- subset_reference[[time]]
-  out$contrast <- NULL
-  out$diff_mean <- out$estimate
-  out$diff_lower <- out$lower.HPD
-  out$diff_upper <- out$upper.HPD
-  out$estimate <- NULL
-  out$lower.HPD <- NULL
-  out$upper.HPD <- NULL
   out
 }
+
+marginal_name <- function(group, time) {
+  sprintf("%s, %s", group , time)
+}
diff --git a/R/utils_assert.R b/R/utils_assert.R
index 7b2de1a9..1bb5efa6 100644
--- a/R/utils_assert.R
+++ b/R/utils_assert.R
@@ -33,7 +33,7 @@ assert_chr_vec <- function(value, message = NULL) {
 
 assert_chr <- function(value, message = NULL) {
   assert_chr_vec(value, message = message)
-  assert(value, length(.) == 1L)
+  assert(value, length(.) == 1L, message = message)
 }
 
 assert_lgl <- function(value, message = NULL) {
diff --git a/man/brm_summary.Rd b/man/brm_summary.Rd
index fa464928..3f68d289 100644
--- a/man/brm_summary.Rd
+++ b/man/brm_summary.Rd
@@ -11,7 +11,9 @@ brm_summary(
   time = "AVISIT",
   patient = "USUBJID",
   covariates = character(0),
-  control = "Placebo"
+  control = "Placebo",
+  baseline = "Baseline",
+  response_type = "change"
 )
 }
 \arguments{
@@ -32,8 +34,17 @@ in the data.}
 \item{covariates}{Character vector of names of other covariates
 in the data.}
 
-\item{control}{Character of length 1, name of the control arm
-in the \code{group} column in the data.}
+\item{control}{Element of the \code{group} column in the data which indicates
+the control group for the purposes of calculating treatment differences.}
+
+\item{baseline}{Element of the \code{time} column in the data
+which indicates the baseline time for the purposes of calculating
+change from baseline. Ignored if \code{response_type = "change"}.}
+
+\item{response_type}{Character of length 1, \code{"response"} if the
+response variable is the raw outcome variable (such as AVAL)
+or \code{"change"} if the response variable is change from baseline
+(e.g. CHG).}
 }
 \value{
 A \code{tibble} with summary statistics of the marginal posterior.
@@ -49,6 +60,8 @@ Currently assumes the response variable is \code{CHG}
 set.seed(0L)
 sim <- brm_simulate()
 data <- sim$data
+data$group <- paste("treatment", data$group)
+data$time <- paste("visit", data$time)
 formula <- brm_formula(
   response = "response",
   group = "group",
@@ -75,7 +88,9 @@ brm_summary(
   group = "group",
   time = "time",
   patient = "patient",
-  control = 1
+  control = "treatment 1",
+  baseline = "visit 1",
+  response_type = "response"
 )
 }
 \concept{results}

From 00f3f13d0584844225a7cd6f886d367e4d2fde41 Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Fri, 2 Jun 2023 14:33:03 -0400
Subject: [PATCH 02/17] Sketch neat summaries

---
 R/brm_summary.R    | 66 +++++++++++++++++++++++++++++++++++++++++++---
 man/brm_summary.Rd |  6 ++++-
 2 files changed, 67 insertions(+), 5 deletions(-)

diff --git a/R/brm_summary.R b/R/brm_summary.R
index e7722fb8..91e286b7 100644
--- a/R/brm_summary.R
+++ b/R/brm_summary.R
@@ -16,6 +16,8 @@
 #'   response variable is the raw outcome variable (such as AVAL)
 #'   or `"change"` if the response variable is change from baseline
 #'   (e.g. CHG).
+#' @param level Numeric of length 1 between 0 and 1, credible level
+#'   for the credible intervals in the output.
 #' @examples
 #' set.seed(0L)
 #' sim <- brm_simulate()
@@ -61,7 +63,8 @@ brm_summary <- function(
   covariates = character(0),
   control = "Placebo",
   baseline = "Baseline",
-  response_type = "change"
+  response_type = "change",
+  level = 0.95
 ) {
   assert_chr(base, "base arg must be a nonempty character string")
   assert_chr(group, "group arg must be a nonempty character string")
@@ -90,6 +93,8 @@ brm_summary <- function(
     !anyNA(.),
     message = "baseline arg must be a length-1 non-missing atomic value"
   )
+  assert_num(level, "level arg must be a length-1 numeric between 0 and 1")
+  assert(level, . >= 0, . <= 1, message = "level arg must be between 0 and 1")
   assert(is.data.frame(model$data))
   data <- model$data
   assert(
@@ -182,9 +187,19 @@ brm_summary <- function(
       control = control
     )
   }
-  
-  browser()
-  
+  table_response <- summarize_marginals(draws_response, level)
+  table_change <- if_any(
+    response_type == "response",
+    summarize_marginals(draws_change, level),
+    NULL
+  )
+  table_diff <- summarize_marginals(draws_diff, level)
+  dplyr::bind_rows(
+    response = table_response,
+    change = table_change,
+    difference = table_diff,
+    .id = "marginal"
+  )
 }
 
 subtract_baseline <- function(draws, groups, times, baseline) {
@@ -214,3 +229,46 @@ subtract_control <- function(draws, groups, times, control) {
 marginal_name <- function(group, time) {
   sprintf("%s, %s", group , time)
 }
+
+summarize_marginals <- function(draws, level) {
+  level_lower <- (1 - level) / 2
+  level_upper <- 1 - level_lower
+  draws[[".chain"]] <- NULL
+  draws[[".iteration"]] <- NULL
+  draws[[".draw"]] <- NULL
+  value <- tibble::tibble(
+    group = gsub(",.*$", "", colnames(draws)),
+    time = gsub("^.*, ", "", colnames(draws)),
+    mean = purrr::map_dbl(draws, mean),
+    median = purrr::map_dbl(draws, median),
+    sd = purrr::map_dbl(draws, sd),
+    lower = purrr::map_dbl(draws, ~quantile(.x, level_lower)),
+    upper = purrr::map_dbl(draws, ~quantile(.x, level_upper))
+  )
+  mcse <- tibble::tibble(
+    group = gsub(",.*$", "", colnames(draws)),
+    time = gsub("^.*, ", "", colnames(draws)),
+    mean = purrr::map_dbl(draws, posterior::mcse_mean),
+    median = purrr::map_dbl(draws, posterior::mcse_median),
+    sd = purrr::map_dbl(draws, posterior::mcse_sd),
+    lower = purrr::map_dbl(draws, ~posterior::mcse_quantile(.x, level_lower)),
+    upper = purrr::map_dbl(draws, ~posterior::mcse_quantile(.x, level_upper))
+  )
+  value <- tidyr::pivot_longer(
+    data = value,
+    cols = -any_of(c("group", "time")),
+    names_to = "statistic",
+    values_to = "value"
+  )
+  mcse <- tidyr::pivot_longer(
+    data = mcse,
+    cols = -any_of(c("group", "time")),
+    names_to = "statistic",
+    values_to = "mcse"
+  )
+  dplyr::left_join(
+    x = value,
+    y = mcse,
+    by = c("group", "time", "statistic")
+  )
+}
diff --git a/man/brm_summary.Rd b/man/brm_summary.Rd
index 3f68d289..56bf8307 100644
--- a/man/brm_summary.Rd
+++ b/man/brm_summary.Rd
@@ -13,7 +13,8 @@ brm_summary(
   covariates = character(0),
   control = "Placebo",
   baseline = "Baseline",
-  response_type = "change"
+  response_type = "change",
+  level = 0.95
 )
 }
 \arguments{
@@ -45,6 +46,9 @@ change from baseline. Ignored if \code{response_type = "change"}.}
 response variable is the raw outcome variable (such as AVAL)
 or \code{"change"} if the response variable is change from baseline
 (e.g. CHG).}
+
+\item{level}{Numeric of length 1 between 0 and 1, credible level
+for the credible intervals in the output.}
 }
 \value{
 A \code{tibble} with summary statistics of the marginal posterior.

From 97cf007c520cb3d2b32109eb20745cc57416560a Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Fri, 2 Jun 2023 15:39:32 -0400
Subject: [PATCH 03/17] Implement #8

---
 R/brm_summary.R    | 135 ++++++++++++++++++++++++++++++++++-----------
 man/brm_summary.Rd |  39 +++++++++----
 2 files changed, 131 insertions(+), 43 deletions(-)

diff --git a/R/brm_summary.R b/R/brm_summary.R
index 91e286b7..73a25ea7 100644
--- a/R/brm_summary.R
+++ b/R/brm_summary.R
@@ -4,20 +4,33 @@
 #' @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.
+#' @return A list of two data frames, one with summary statistics on the
+#'   marginal posterior and another with posterior probabilities
+#'   on the treatment effects.
 #' @inheritParams brm_formula
 #' @param model Fitted `brms` model object from [brm_model()].
+#' @param outcome Character of length 1, `"response"` if the
+#'   response variable is the raw outcome variable (such as AVAL)
+#'   or `"change"` if the response variable is change from baseline
+#'   (e.g. CHG).
 #' @param control Element of the `group` column in the data which indicates
 #'   the control group for the purposes of calculating treatment differences.
 #' @param baseline Element of the `time` column in the data
 #'   which indicates the baseline time for the purposes of calculating
-#'   change from baseline. Ignored if `response_type = "change"`.
-#' @param response_type Character of length 1, `"response"` if the
-#'   response variable is the raw outcome variable (such as AVAL)
-#'   or `"change"` if the response variable is change from baseline
-#'   (e.g. CHG).
+#'   change from baseline.
 #' @param level Numeric of length 1 between 0 and 1, credible level
-#'   for the credible intervals in the output.
+#'   for the credible intervals. Only relevant when `return = "marginals"`.
+#' @param direction Character vector of the same length as `threshold`.
+#'   `"greater"` to compute the marginal posterior probability that the
+#'   treatment effect is greater than the threshold,
+#'   `"less"` to compute the marginal posterior probability that the
+#'   treatment effect is less than the threshold.
+#'   Each element `direction[i]` corresponds to `threshold[i]`
+#'   for all `i` from 1 to `length(direction)`.
+#' @param threshold Numeric vector of the same length as `direction`,
+#'   treatment effect threshold for computing posterior probabilities.
+#'   Each element `direction[i]` corresponds to `threshold[i]` for
+#'   all `i` from 1 to `length(direction)`.
 #' @examples
 #' set.seed(0L)
 #' sim <- brm_simulate()
@@ -52,7 +65,7 @@
 #'   patient = "patient",
 #'   control = "treatment 1",
 #'   baseline = "visit 1",
-#'   response_type = "response"
+#'   outcome = "response"
 #' )
 brm_summary <- function(
   model,
@@ -61,22 +74,24 @@ brm_summary <- function(
   time = "AVISIT",
   patient = "USUBJID",
   covariates = character(0),
+  outcome = "change",
   control = "Placebo",
   baseline = "Baseline",
-  response_type = "change",
-  level = 0.95
+  level = 0.95,
+  direction = "greater",
+  threshold = 0
 ) {
   assert_chr(base, "base arg must be a nonempty character string")
   assert_chr(group, "group arg must be a nonempty character string")
   assert_chr(time, "time arg must be a nonempty character string")
   assert_chr(patient, "patient arg must be a nonempty character string")
   assert_chr(
-    response_type,
-    "response_type arg must be a nonempty character string"
+    outcome,
+    "outcome arg must be a nonempty character string"
   )
   assert(
-    response_type %in% c("response", "change"),
-    message = "response_type must be either \"response\" or \"change\""
+    outcome %in% c("response", "change"),
+    message = "outcome must be either \"response\" or \"change\""
   )
   assert_chr_vec(covariates, "covariates arg must be a character vector")
   assert(
@@ -95,6 +110,24 @@ brm_summary <- function(
   )
   assert_num(level, "level arg must be a length-1 numeric between 0 and 1")
   assert(level, . >= 0, . <= 1, message = "level arg must be between 0 and 1")
+  assert(
+    direction,
+    is.character(.),
+    !anyNA(.),
+    nzchar(.),
+    . %in% c("greater", "less"),
+    message = "elements of the direction arg must be \"greater\" or \"less\""
+  )
+  assert(
+    threshold,
+    is.numeric(.),
+    is.finite(.),
+    message = "threshold arg must be a numeric vector"
+  )
+  assert(
+    length(direction) == length(threshold),
+    message = "direction and threshold must have the same length"
+  )
   assert(is.data.frame(model$data))
   data <- model$data
   assert(
@@ -141,8 +174,8 @@ brm_summary <- function(
     replacement = ", ",
     x = colnames(draws_response)
   )
-  groups <- unique(gsub(",.*$", "", colnames(draws_response)))
-  times <- unique(gsub("^.*, ", "", colnames(draws_response)))
+  groups <- unique(group_names(draws_response))
+  times <- unique(time_names(draws_response))
   draws_response[[".chain"]] <- .chain
   draws_response[[".iteration"]] <- .iteration
   draws_response[[".draw"]] <- .draw
@@ -156,7 +189,7 @@ brm_summary <- function(
       paste(groups, collapse = ", ")
     )
   )
-  if (response_type == "response") {
+  if (outcome == "response") {
     assert(
       baseline %in% times,
       message = sprintf(
@@ -166,21 +199,21 @@ brm_summary <- function(
       )
     )
   }
-  if (response_type == "response") {
+  if (outcome == "response") {
     draws_change <- subtract_baseline(
       draws = draws_response,
       groups = groups,
       times = times,
       baseline = baseline
     )
-    draws_diff <- subtract_control(
+    draws_difference <- subtract_control(
       draws = draws_change,
       groups = groups,
       times = setdiff(times, baseline),
       control = control
     )
   } else {
-    draws_diff <- subtract_control(
+    draws_difference <- subtract_control(
       draws = draws_response,
       groups = groups,
       times = times,
@@ -189,17 +222,23 @@ brm_summary <- function(
   }
   table_response <- summarize_marginals(draws_response, level)
   table_change <- if_any(
-    response_type == "response",
+    outcome == "response",
     summarize_marginals(draws_change, level),
     NULL
   )
-  table_diff <- summarize_marginals(draws_diff, level)
-  dplyr::bind_rows(
+  table_difference <- summarize_marginals(draws_difference, level)
+  marginals <- dplyr::bind_rows(
     response = table_response,
     change = table_change,
-    difference = table_diff,
+    difference = table_difference,
     .id = "marginal"
   )
+  probabilities <- summarize_probabilities(
+    draws = draws_difference,
+    direction = direction,
+    threshold = threshold
+  )
+  list(marginals = marginals, probabilities = probabilities)
 }
 
 subtract_baseline <- function(draws, groups, times, baseline) {
@@ -226,10 +265,6 @@ subtract_control <- function(draws, groups, times, control) {
   out
 }
 
-marginal_name <- function(group, time) {
-  sprintf("%s, %s", group , time)
-}
-
 summarize_marginals <- function(draws, level) {
   level_lower <- (1 - level) / 2
   level_upper <- 1 - level_lower
@@ -237,8 +272,8 @@ summarize_marginals <- function(draws, level) {
   draws[[".iteration"]] <- NULL
   draws[[".draw"]] <- NULL
   value <- tibble::tibble(
-    group = gsub(",.*$", "", colnames(draws)),
-    time = gsub("^.*, ", "", colnames(draws)),
+    group = group_names(draws),
+    time = time_names(draws),
     mean = purrr::map_dbl(draws, mean),
     median = purrr::map_dbl(draws, median),
     sd = purrr::map_dbl(draws, sd),
@@ -246,8 +281,8 @@ summarize_marginals <- function(draws, level) {
     upper = purrr::map_dbl(draws, ~quantile(.x, level_upper))
   )
   mcse <- tibble::tibble(
-    group = gsub(",.*$", "", colnames(draws)),
-    time = gsub("^.*, ", "", colnames(draws)),
+    group = group_names(draws),
+    time = time_names(draws),
     mean = purrr::map_dbl(draws, posterior::mcse_mean),
     median = purrr::map_dbl(draws, posterior::mcse_median),
     sd = purrr::map_dbl(draws, posterior::mcse_sd),
@@ -272,3 +307,39 @@ summarize_marginals <- function(draws, level) {
     by = c("group", "time", "statistic")
   )
 }
+
+summarize_probabilities <- function(draws, direction, threshold) {
+  draws[[".chain"]] <- NULL
+  draws[[".iteration"]] <- NULL
+  draws[[".draw"]] <- NULL
+  tibble::tibble(
+    group = group_names(draws),
+    time = time_names(draws),
+    direction = direction,
+    theshold = threshold,
+    value = purrr::map_dbl(
+      draws,
+      ~marginal_probability(.x, direction, threshold)
+    )
+  )
+}
+
+marginal_probability <- function(difference, direction, threshold) {
+  if_any(
+    direction == "greater",
+    mean(difference > threshold),
+    mean(difference < threshold)
+  )
+}
+
+marginal_name <- function(group, time) {
+  sprintf("%s, %s", group , time)
+}
+
+group_names <- function(draws) {
+  gsub(",.*$", "", colnames(draws))
+}
+
+time_names <- function(draws) {
+  gsub("^.*, ", "", colnames(draws))
+}
diff --git a/man/brm_summary.Rd b/man/brm_summary.Rd
index 56bf8307..76b50246 100644
--- a/man/brm_summary.Rd
+++ b/man/brm_summary.Rd
@@ -11,10 +11,12 @@ brm_summary(
   time = "AVISIT",
   patient = "USUBJID",
   covariates = character(0),
+  outcome = "change",
   control = "Placebo",
   baseline = "Baseline",
-  response_type = "change",
-  level = 0.95
+  level = 0.95,
+  direction = "greater",
+  threshold = 0
 )
 }
 \arguments{
@@ -35,23 +37,38 @@ in the data.}
 \item{covariates}{Character vector of names of other covariates
 in the data.}
 
+\item{outcome}{Character of length 1, \code{"response"} if the
+response variable is the raw outcome variable (such as AVAL)
+or \code{"change"} if the response variable is change from baseline
+(e.g. CHG).}
+
 \item{control}{Element of the \code{group} column in the data which indicates
 the control group for the purposes of calculating treatment differences.}
 
 \item{baseline}{Element of the \code{time} column in the data
 which indicates the baseline time for the purposes of calculating
-change from baseline. Ignored if \code{response_type = "change"}.}
-
-\item{response_type}{Character of length 1, \code{"response"} if the
-response variable is the raw outcome variable (such as AVAL)
-or \code{"change"} if the response variable is change from baseline
-(e.g. CHG).}
+change from baseline.}
 
 \item{level}{Numeric of length 1 between 0 and 1, credible level
-for the credible intervals in the output.}
+for the credible intervals. Only relevant when \code{return = "marginals"}.}
+
+\item{direction}{Character vector of the same length as \code{threshold}.
+\code{"greater"} to compute the marginal posterior probability that the
+treatment effect is greater than the threshold,
+\code{"less"} to compute the marginal posterior probability that the
+treatment effect is less than the threshold.
+Each element \code{direction[i]} corresponds to \code{threshold[i]}
+for all \code{i} from 1 to \code{length(direction)}.}
+
+\item{threshold}{Numeric vector of the same length as \code{direction},
+treatment effect threshold for computing posterior probabilities.
+Each element \code{direction[i]} corresponds to \code{threshold[i]} for
+all \code{i} from 1 to \code{length(direction)}.}
 }
 \value{
-A \code{tibble} with summary statistics of the marginal posterior.
+A list of two data frames, one with summary statistics on the
+marginal posterior and another with posterior probabilities
+on the treatment effects.
 }
 \description{
 Summarize a basic MMRM model fit.
@@ -94,7 +111,7 @@ brm_summary(
   patient = "patient",
   control = "treatment 1",
   baseline = "visit 1",
-  response_type = "response"
+  outcome = "response"
 )
 }
 \concept{results}

From 3c2573e292dd875cfb00879d0de926c0af3930ad Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 10:46:57 -0400
Subject: [PATCH 04/17] Add brm_marginals()

---
 NAMESPACE            |   1 +
 R/brm_marginals.R    | 224 +++++++++++++++++++++++++++++++++++++++++
 R/brm_summary.R      | 230 +++++++------------------------------------
 _pkgdown.yml         |   2 +-
 man/brm_marginals.Rd | 114 +++++++++++++++++++++
 man/brm_summary.Rd   |  79 ++++++---------
 6 files changed, 407 insertions(+), 243 deletions(-)
 create mode 100644 R/brm_marginals.R
 create mode 100644 man/brm_marginals.Rd

diff --git a/NAMESPACE b/NAMESPACE
index c166e393..960f2ed3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(brm_formula)
+export(brm_marginals)
 export(brm_model)
 export(brm_simulate)
 export(brm_summary)
diff --git a/R/brm_marginals.R b/R/brm_marginals.R
new file mode 100644
index 00000000..df1524ac
--- /dev/null
+++ b/R/brm_marginals.R
@@ -0,0 +1,224 @@
+#' @title MMRM marginal posterior samples.
+#' @export
+#' @family results
+#' @description Get marginal posteior samples from a fitted MMRM.
+#' @details Currently assumes the response variable is `CHG`
+#'   (change from baseline) and not `AVAL` (raw response).
+#' @return A named list of tibbles of MCMC samples of the marginal posterior
+#'   distribution of each treatment group and time point:
+#'   * `response`: on the scale of the response variable.
+#'   * `change`: change from baseline, where the `baseline` argument determines
+#'     the time point at baseline. Only returned if the `outcome` argument is
+#'     `"response"`. (If `outcome` is `"change"`, then `response` already
+#'     represents change from baseline.)
+#'   * `difference`: treatment effect of change from baseline, where the
+#'     `control` argument identifies the placebo or active control group.
+#'   In each tibble, there is 1 row per posterior sample and one column for
+#'   each type of marginal distribution (i.e. each combination of treatment
+#'   group and discrete time point.
+#'   Treatment and time are comma-delimited in the column names.
+#' @inheritParams brm_formula
+#' @param model Fitted `brms` model object from [brm_model()].
+#' @param outcome Character of length 1, `"response"` if the
+#'   response variable is the raw outcome variable (such as AVAL)
+#'   or `"change"` if the response variable is change from baseline
+#'   (e.g. CHG).
+#' @param control Element of the `group` column in the data which indicates
+#'   the control group for the purposes of calculating treatment differences.
+#' @param baseline Element of the `time` column in the data
+#'   which indicates the baseline time for the purposes of calculating
+#'   change from baseline.
+#' @examples
+#' set.seed(0L)
+#' sim <- brm_simulate()
+#' data <- sim$data
+#' data$group <- paste("treatment", data$group)
+#' data$time <- paste("visit", data$time)
+#' 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
+#'       )
+#'     )
+#'   )
+#' )
+#' brm_marginals(
+#'   model = model,
+#'   group = "group",
+#'   time = "time",
+#'   patient = "patient",
+#'   control = "treatment 1",
+#'   baseline = "visit 1",
+#'   outcome = "response"
+#' )
+brm_marginals <- function(
+  model,
+  base = "BASE",
+  group = "TRT01P",
+  time = "AVISIT",
+  patient = "USUBJID",
+  covariates = character(0),
+  outcome = "change",
+  control = "Placebo",
+  baseline = "Baseline"
+) {
+  assert_chr(base, "base arg must be a nonempty character string")
+  assert_chr(group, "group arg must be a nonempty character string")
+  assert_chr(time, "time arg must be a nonempty character string")
+  assert_chr(patient, "patient arg must be a nonempty character string")
+  assert_chr(
+    outcome,
+    "outcome arg must be a nonempty character string"
+  )
+  assert(
+    outcome %in% c("response", "change"),
+    message = "outcome must be either \"response\" or \"change\""
+  )
+  assert_chr_vec(covariates, "covariates arg must be a character vector")
+  assert(
+    control,
+    is.atomic(.),
+    length(.) == 1L,
+    !anyNA(.),
+    message = "control arg must be a length-1 non-missing atomic value"
+  )
+  assert(
+    baseline,
+    is.atomic(.),
+    length(.) == 1L,
+    !anyNA(.),
+    message = "baseline arg must be a length-1 non-missing atomic value"
+  )
+  assert(is.data.frame(model$data))
+  data <- model$data
+  assert(
+    group %in% colnames(data),
+    message = "group arg must be a data column name"
+  )
+  assert(
+    time %in% colnames(data),
+    message = "time arg must be a data column name"
+  )
+  assert(
+    patient %in% colnames(data),
+    message = "patient arg must be a data column name"
+  )
+  assert(
+    covariates %in% colnames(data),
+    message = "all covariates must be data column names"
+  )
+  assert(
+    control %in% data[[group]],
+    message = "control arg must be in data[[group]]"
+  )
+  nuisance <- c(base, patient, covariates)
+  emmeans <- emmeans::emmeans(
+    object = model,
+    specs = as.formula(sprintf("~%s:%s", group, time)),
+    weights = "proportional",
+    nuisance = nuisance
+  )
+  mcmc <- coda::as.mcmc(emmeans, fixed = TRUE, names = FALSE)
+  draws_response <- posterior::as_draws_df(mcmc)
+  groups <- unique(names_group(draws_response))
+  times <- unique(names_time(draws_response))
+  control <- as.character(control)
+  time <- as.character(time)
+  assert(
+    control %in% groups,
+    message = sprintf(
+      "control argument \"%s\" is not in one of the treatment groups: %s",
+      control,
+      paste(groups, collapse = ", ")
+    )
+  )
+  if (identical(outcome, "response")) {
+    assert(
+      baseline %in% times,
+      message = sprintf(
+        "baseline argument \"%s\" is not in one of the time points: %s",
+        baseline,
+        paste(times, collapse = ", ")
+      )
+    )
+    draws_change <- subtract_baseline(
+      draws = draws_response,
+      groups = groups,
+      times = times,
+      baseline = baseline
+    )
+    draws_difference <- subtract_control(
+      draws = draws_change,
+      groups = groups,
+      times = setdiff(times, baseline),
+      control = control
+    )
+  } else {
+    draws_difference <- subtract_control(
+      draws = draws_response,
+      groups = groups,
+      times = times,
+      control = control
+    )
+  }
+  out <- list()
+  out$response <- draws_response
+  if (identical(outcome, "response")) {
+    out$change <- draws_change
+  }
+  out$difference <- draws_difference
+  out
+}
+
+subtract_baseline <- function(draws, groups, times, baseline) {
+  out <- draws[, names_mcmc]
+  for (group in groups) {
+    for (time in setdiff(times, baseline)) {
+      name1 <- name_marginal(group, baseline)
+      name2 <- name_marginal(group, time)
+      out[[name2]] <- draws[[name2]] - draws[[name1]]
+    }
+  }
+  out
+}
+
+subtract_control <- function(draws, groups, times, control) {
+  out <- draws[, names_mcmc]
+  for (group in setdiff(groups, control)) {
+    for (time in times) {
+      name1 <- name_marginal(control, time)
+      name2 <- name_marginal(group, time)
+      out[[name2]] <- draws[[name2]] - draws[[name1]]
+    }
+  }
+  out
+}
+
+name_marginal <- function(group, time) {
+  sprintf("%s, %s", group , time)
+}
+
+names_group <- function(draws) {
+  names <- setdiff(colnames(draws), names_mcmc)
+  gsub(",.*$", "", names)
+}
+
+names_time <- function(draws) {
+  names <- setdiff(colnames(draws), names_mcmc)
+  gsub("^.*, ", "", names)
+}
+
+names_mcmc <- c(".chain", ".draw", ".iteration")
diff --git a/R/brm_summary.R b/R/brm_summary.R
index 73a25ea7..9f3101a5 100644
--- a/R/brm_summary.R
+++ b/R/brm_summary.R
@@ -2,22 +2,26 @@
 #' @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 list of two data frames, one with summary statistics on the
-#'   marginal posterior and another with posterior probabilities
-#'   on the treatment effects.
-#' @inheritParams brm_formula
-#' @param model Fitted `brms` model object from [brm_model()].
-#' @param outcome Character of length 1, `"response"` if the
-#'   response variable is the raw outcome variable (such as AVAL)
-#'   or `"change"` if the response variable is change from baseline
-#'   (e.g. CHG).
-#' @param control Element of the `group` column in the data which indicates
-#'   the control group for the purposes of calculating treatment differences.
-#' @param baseline Element of the `time` column in the data
-#'   which indicates the baseline time for the purposes of calculating
-#'   change from baseline.
+#' @return A named list of two data frames with one row per summary
+#'   statistic:
+#'   * `means`: summary statistics of the marginal posterior mean of
+#'     each treatment group and time point. Columns include `marginal`
+#'     for the type of marginal distribution (see the "Value" section of
+#'     the [brm_marginals()] help file for details), `group` for treatment
+#'     group, `time` for discrete time point, `statistic` for the type of
+#'     summary statistic, `value` for the numeric value of the estimate,
+#'     and `mcse` for the Monte Carlo standard error of the estimate.
+#'   * `probabilities`: marginal posterior probabilities of the form
+#'     `Prob(treatment effect > threshold | data)` and/or
+#'     `Prob(treatment effect < threshold | data)`. Columns include
+#'     `group` for the treatment group, `time` for the discrete time point,
+#'     `direction` to indicate the direction of the comparison
+#'     (`"greater"` for `>`, `"less"` for `<`), `threshold` for the
+#'     treatment effect threshold in the probability statement,
+#'     and `value` for the numberic value of the estimate of the
+#'     probability.
+#' @param marginals Posterior draws of the marginal posterior
+#'   obtained from [brm_marginals()].
 #' @param level Numeric of length 1 between 0 and 1, credible level
 #'   for the credible intervals. Only relevant when `return = "marginals"`.
 #' @param direction Character vector of the same length as `threshold`.
@@ -58,7 +62,7 @@
 #'     )
 #'   )
 #' )
-#' brm_summary(
+#' marginals <- brm_marginals(
 #'   model = model,
 #'   group = "group",
 #'   time = "time",
@@ -67,46 +71,16 @@
 #'   baseline = "visit 1",
 #'   outcome = "response"
 #' )
+#' brm_summary(marginals)
 brm_summary <- function(
-  model,
-  base = "BASE",
-  group = "TRT01P",
-  time = "AVISIT",
-  patient = "USUBJID",
-  covariates = character(0),
-  outcome = "change",
-  control = "Placebo",
-  baseline = "Baseline",
+  marginals,
   level = 0.95,
   direction = "greater",
   threshold = 0
 ) {
-  assert_chr(base, "base arg must be a nonempty character string")
-  assert_chr(group, "group arg must be a nonempty character string")
-  assert_chr(time, "time arg must be a nonempty character string")
-  assert_chr(patient, "patient arg must be a nonempty character string")
-  assert_chr(
-    outcome,
-    "outcome arg must be a nonempty character string"
-  )
-  assert(
-    outcome %in% c("response", "change"),
-    message = "outcome must be either \"response\" or \"change\""
-  )
-  assert_chr_vec(covariates, "covariates arg must be a character vector")
-  assert(
-    control,
-    is.atomic(.),
-    length(.) == 1L,
-    !anyNA(.),
-    message = "control arg must be a length-1 non-missing atomic value"
-  )
   assert(
-    baseline,
-    is.atomic(.),
-    length(.) == 1L,
-    !anyNA(.),
-    message = "baseline arg must be a length-1 non-missing atomic value"
+    is.list(marginals),
+    message = "marginals arg must be a named list from brm_marginals()"
   )
   assert_num(level, "level arg must be a length-1 numeric between 0 and 1")
   assert(level, . >= 0, . <= 1, message = "level arg must be between 0 and 1")
@@ -128,141 +102,25 @@ brm_summary <- function(
     length(direction) == length(threshold),
     message = "direction and threshold must have the same length"
   )
-  assert(is.data.frame(model$data))
-  data <- model$data
-  assert(
-    group %in% colnames(data),
-    message = "group arg must be a data column name"
-  )
-  assert(
-    time %in% colnames(data),
-    message = "time arg must be a data column name"
-  )
-  assert(
-    patient %in% colnames(data),
-    message = "patient arg must be a data column name"
-  )
-  assert(
-    covariates %in% colnames(data),
-    message = "all covariates must be data column names"
-  )
-  assert(
-    control %in% data[[group]],
-    message = "control arg must be in data[[group]]"
-  )
-  nuisance <- c(base, patient, covariates)
-  emmeans <- emmeans::emmeans(
-    object = model,
-    specs = as.formula(sprintf("~%s:%s", group, time)),
-    weights = "proportional",
-    nuisance = nuisance
-  )
-  draws_response <- posterior::as_draws_df(as.mcmc(emmeans))
-  .chain <- draws_response[[".chain"]]
-  .iteration <- draws_response[[".iteration"]]
-  .draw <- draws_response[[".draw"]]
-  draws_response[[".chain"]] <- NULL
-  draws_response[[".iteration"]] <- NULL
-  draws_response[[".draw"]] <- NULL
-  colnames(draws_response) <- gsub(
-    pattern = sprintf("^%s ", group),
-    replacement = "",
-    x = colnames(draws_response)
-  )
-  colnames(draws_response) <- gsub(
-    pattern = sprintf(", %s ", time),
-    replacement = ", ",
-    x = colnames(draws_response)
-  )
-  groups <- unique(group_names(draws_response))
-  times <- unique(time_names(draws_response))
-  draws_response[[".chain"]] <- .chain
-  draws_response[[".iteration"]] <- .iteration
-  draws_response[[".draw"]] <- .draw
-  control <- as.character(control)
-  time <- as.character(time)
-  assert(
-    control %in% groups,
-    message = sprintf(
-      "control argument \"%s\" is not in one of the treatment groups: %s",
-      control,
-      paste(groups, collapse = ", ")
-    )
-  )
-  if (outcome == "response") {
-    assert(
-      baseline %in% times,
-      message = sprintf(
-        "baseline argument \"%s\" is not in one of the time points: %s",
-        baseline,
-        paste(times, collapse = ", ")
-      )
-    )
-  }
-  if (outcome == "response") {
-    draws_change <- subtract_baseline(
-      draws = draws_response,
-      groups = groups,
-      times = times,
-      baseline = baseline
-    )
-    draws_difference <- subtract_control(
-      draws = draws_change,
-      groups = groups,
-      times = setdiff(times, baseline),
-      control = control
-    )
-  } else {
-    draws_difference <- subtract_control(
-      draws = draws_response,
-      groups = groups,
-      times = times,
-      control = control
-    )
-  }
-  table_response <- summarize_marginals(draws_response, level)
+  table_response <- summarize_marginals(marginals$response, level)
   table_change <- if_any(
-    outcome == "response",
-    summarize_marginals(draws_change, level),
+    "change" %in% names(marginals),
+    summarize_marginals(marginals$change, level),
     NULL
   )
-  table_difference <- summarize_marginals(draws_difference, level)
-  marginals <- dplyr::bind_rows(
+  table_difference <- summarize_marginals(marginals$difference, level)
+  means <- dplyr::bind_rows(
     response = table_response,
     change = table_change,
     difference = table_difference,
     .id = "marginal"
   )
   probabilities <- summarize_probabilities(
-    draws = draws_difference,
+    draws = marginals$difference,
     direction = direction,
     threshold = threshold
   )
-  list(marginals = marginals, probabilities = probabilities)
-}
-
-subtract_baseline <- function(draws, groups, times, baseline) {
-  out <- draws[, c(".chain", ".iteration", ".draw")]
-  for (group in groups) {
-    for (time in setdiff(times, baseline)) {
-      name1 <- marginal_name(group, baseline)
-      name2 <- marginal_name(group, time)
-      out[[name2]] <- draws[[name2]] - draws[[name1]]
-    }
-  }
-  out
-}
-
-subtract_control <- function(draws, groups, times, control) {
-  out <- draws[, c(".chain", ".iteration", ".draw")]
-  for (group in setdiff(groups, control)) {
-    for (time in times) {
-      name1 <- marginal_name(control, time)
-      name2 <- marginal_name(group, time)
-      out[[name2]] <- draws[[name2]] - draws[[name1]]
-    }
-  }
-  out
+  list(means = means, probabilities = probabilities)
 }
 
 summarize_marginals <- function(draws, level) {
@@ -272,8 +130,8 @@ summarize_marginals <- function(draws, level) {
   draws[[".iteration"]] <- NULL
   draws[[".draw"]] <- NULL
   value <- tibble::tibble(
-    group = group_names(draws),
-    time = time_names(draws),
+    group = names_group(draws),
+    time = names_time(draws),
     mean = purrr::map_dbl(draws, mean),
     median = purrr::map_dbl(draws, median),
     sd = purrr::map_dbl(draws, sd),
@@ -281,8 +139,8 @@ summarize_marginals <- function(draws, level) {
     upper = purrr::map_dbl(draws, ~quantile(.x, level_upper))
   )
   mcse <- tibble::tibble(
-    group = group_names(draws),
-    time = time_names(draws),
+    group = names_group(draws),
+    time = names_time(draws),
     mean = purrr::map_dbl(draws, posterior::mcse_mean),
     median = purrr::map_dbl(draws, posterior::mcse_median),
     sd = purrr::map_dbl(draws, posterior::mcse_sd),
@@ -313,8 +171,8 @@ summarize_probabilities <- function(draws, direction, threshold) {
   draws[[".iteration"]] <- NULL
   draws[[".draw"]] <- NULL
   tibble::tibble(
-    group = group_names(draws),
-    time = time_names(draws),
+    group = names_group(draws),
+    time = names_time(draws),
     direction = direction,
     theshold = threshold,
     value = purrr::map_dbl(
@@ -331,15 +189,3 @@ marginal_probability <- function(difference, direction, threshold) {
     mean(difference < threshold)
   )
 }
-
-marginal_name <- function(group, time) {
-  sprintf("%s, %s", group , time)
-}
-
-group_names <- function(draws) {
-  gsub(",.*$", "", colnames(draws))
-}
-
-time_names <- function(draws) {
-  gsub("^.*, ", "", colnames(draws))
-}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 98d2c2af..7c452d8c 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -17,5 +17,5 @@ reference:
 - title: Results
   contents:
   - '`brm_convergence`'
-  - '`brm_plot`'
+  - '`brm_marginals`'
   - '`brm_summary`'
diff --git a/man/brm_marginals.Rd b/man/brm_marginals.Rd
new file mode 100644
index 00000000..8f4d77e7
--- /dev/null
+++ b/man/brm_marginals.Rd
@@ -0,0 +1,114 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/brm_marginals.R
+\name{brm_marginals}
+\alias{brm_marginals}
+\title{MMRM marginal posterior samples.}
+\usage{
+brm_marginals(
+  model,
+  base = "BASE",
+  group = "TRT01P",
+  time = "AVISIT",
+  patient = "USUBJID",
+  covariates = character(0),
+  outcome = "change",
+  control = "Placebo",
+  baseline = "Baseline"
+)
+}
+\arguments{
+\item{model}{Fitted \code{brms} model object from \code{\link[=brm_model]{brm_model()}}.}
+
+\item{base}{Character of length 1, name of the baseline response variable
+in the data.}
+
+\item{group}{Character of length 1, name of the treatment group
+variable in the data.}
+
+\item{time}{Character of length 1, name of the discrete time variable
+in the data.}
+
+\item{patient}{Character of length 1, name of the patient ID variable
+in the data.}
+
+\item{covariates}{Character vector of names of other covariates
+in the data.}
+
+\item{outcome}{Character of length 1, \code{"response"} if the
+response variable is the raw outcome variable (such as AVAL)
+or \code{"change"} if the response variable is change from baseline
+(e.g. CHG).}
+
+\item{control}{Element of the \code{group} column in the data which indicates
+the control group for the purposes of calculating treatment differences.}
+
+\item{baseline}{Element of the \code{time} column in the data
+which indicates the baseline time for the purposes of calculating
+change from baseline.}
+}
+\value{
+A named list of tibbles of MCMC samples of the marginal posterior
+distribution of each treatment group and time point:
+\itemize{
+\item \code{response}: on the scale of the response variable.
+\item \code{change}: change from baseline, where the \code{baseline} argument determines
+the time point at baseline. Only returned if the \code{outcome} argument is
+\code{"response"}. (If \code{outcome} is \code{"change"}, then \code{response} already
+represents change from baseline.)
+\item \code{difference}: treatment effect of change from baseline, where the
+\code{control} argument identifies the placebo or active control group.
+In each tibble, there is 1 row per posterior sample and one column for
+each type of marginal distribution (i.e. each combination of treatment
+group and discrete time point.
+Treatment and time are comma-delimited in the column names.
+}
+}
+\description{
+Get marginal posteior samples from a fitted MMRM.
+}
+\details{
+Currently assumes the response variable is \code{CHG}
+(change from baseline) and not \code{AVAL} (raw response).
+}
+\examples{
+set.seed(0L)
+sim <- brm_simulate()
+data <- sim$data
+data$group <- paste("treatment", data$group)
+data$time <- paste("visit", data$time)
+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
+      )
+    )
+  )
+)
+brm_marginals(
+  model = model,
+  group = "group",
+  time = "time",
+  patient = "patient",
+  control = "treatment 1",
+  baseline = "visit 1",
+  outcome = "response"
+)
+}
+\seealso{
+Other results: 
+\code{\link{brm_summary}()}
+}
+\concept{results}
diff --git a/man/brm_summary.Rd b/man/brm_summary.Rd
index 76b50246..a3ea1735 100644
--- a/man/brm_summary.Rd
+++ b/man/brm_summary.Rd
@@ -4,50 +4,11 @@
 \alias{brm_summary}
 \title{Summarize an MMRM.}
 \usage{
-brm_summary(
-  model,
-  base = "BASE",
-  group = "TRT01P",
-  time = "AVISIT",
-  patient = "USUBJID",
-  covariates = character(0),
-  outcome = "change",
-  control = "Placebo",
-  baseline = "Baseline",
-  level = 0.95,
-  direction = "greater",
-  threshold = 0
-)
+brm_summary(marginals, level = 0.95, direction = "greater", threshold = 0)
 }
 \arguments{
-\item{model}{Fitted \code{brms} model object from \code{\link[=brm_model]{brm_model()}}.}
-
-\item{base}{Character of length 1, name of the baseline response variable
-in the data.}
-
-\item{group}{Character of length 1, name of the treatment group
-variable in the data.}
-
-\item{time}{Character of length 1, name of the discrete time variable
-in the data.}
-
-\item{patient}{Character of length 1, name of the patient ID variable
-in the data.}
-
-\item{covariates}{Character vector of names of other covariates
-in the data.}
-
-\item{outcome}{Character of length 1, \code{"response"} if the
-response variable is the raw outcome variable (such as AVAL)
-or \code{"change"} if the response variable is change from baseline
-(e.g. CHG).}
-
-\item{control}{Element of the \code{group} column in the data which indicates
-the control group for the purposes of calculating treatment differences.}
-
-\item{baseline}{Element of the \code{time} column in the data
-which indicates the baseline time for the purposes of calculating
-change from baseline.}
+\item{marginals}{Posterior draws of the marginal posterior
+obtained from \code{\link[=brm_marginals]{brm_marginals()}}.}
 
 \item{level}{Numeric of length 1 between 0 and 1, credible level
 for the credible intervals. Only relevant when \code{return = "marginals"}.}
@@ -66,17 +27,30 @@ Each element \code{direction[i]} corresponds to \code{threshold[i]} for
 all \code{i} from 1 to \code{length(direction)}.}
 }
 \value{
-A list of two data frames, one with summary statistics on the
-marginal posterior and another with posterior probabilities
-on the treatment effects.
+A named list of two data frames with one row per summary
+statistic:
+\itemize{
+\item \code{means}: summary statistics of the marginal posterior mean of
+each treatment group and time point. Columns include \code{marginal}
+for the type of marginal distribution (see the "Value" section of
+the \code{\link[=brm_marginals]{brm_marginals()}} help file for details), \code{group} for treatment
+group, \code{time} for discrete time point, \code{statistic} for the type of
+summary statistic, \code{value} for the numeric value of the estimate,
+and \code{mcse} for the Monte Carlo standard error of the estimate.
+\item \code{probabilities}: marginal posterior probabilities of the form
+\verb{Prob(treatment effect > threshold | data)} and/or
+\verb{Prob(treatment effect < threshold | data)}. Columns include
+\code{group} for the treatment group, \code{time} for the discrete time point,
+\code{direction} to indicate the direction of the comparison
+(\code{"greater"} for \code{>}, \code{"less"} for \code{<}), \code{threshold} for the
+treatment effect threshold in the probability statement,
+and \code{value} for the numberic value of the estimate of the
+probability.
+}
 }
 \description{
 Summarize a basic MMRM model fit.
 }
-\details{
-Currently assumes the response variable is \code{CHG}
-(change from baseline) and not \code{AVAL} (raw response).
-}
 \examples{
 set.seed(0L)
 sim <- brm_simulate()
@@ -104,7 +78,7 @@ tmp <- utils::capture.output(
     )
   )
 )
-brm_summary(
+marginals <- brm_marginals(
   model = model,
   group = "group",
   time = "time",
@@ -113,5 +87,10 @@ brm_summary(
   baseline = "visit 1",
   outcome = "response"
 )
+brm_summary(marginals)
+}
+\seealso{
+Other results: 
+\code{\link{brm_marginals}()}
 }
 \concept{results}

From 5e7864d43d8a5bfa9682f43a1509683bfb831341 Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 10:49:42 -0400
Subject: [PATCH 05/17] lint

---
 R/brm_marginals.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/brm_marginals.R b/R/brm_marginals.R
index df1524ac..0a028521 100644
--- a/R/brm_marginals.R
+++ b/R/brm_marginals.R
@@ -208,7 +208,7 @@ subtract_control <- function(draws, groups, times, control) {
 }
 
 name_marginal <- function(group, time) {
-  sprintf("%s, %s", group , time)
+  sprintf("%s, %s", group, time)
 }
 
 names_group <- function(draws) {

From dbc3f555d42079aa4fec030558fd4bf7eecdba06 Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 11:12:26 -0400
Subject: [PATCH 06/17] conciseness

---
 R/brm_summary.R | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/R/brm_summary.R b/R/brm_summary.R
index 9f3101a5..69e9a167 100644
--- a/R/brm_summary.R
+++ b/R/brm_summary.R
@@ -126,9 +126,7 @@ brm_summary <- function(
 summarize_marginals <- function(draws, level) {
   level_lower <- (1 - level) / 2
   level_upper <- 1 - level_lower
-  draws[[".chain"]] <- NULL
-  draws[[".iteration"]] <- NULL
-  draws[[".draw"]] <- NULL
+  draws[names_mcmc] <- NULL
   value <- tibble::tibble(
     group = names_group(draws),
     time = names_time(draws),
@@ -167,9 +165,7 @@ summarize_marginals <- function(draws, level) {
 }
 
 summarize_probabilities <- function(draws, direction, threshold) {
-  draws[[".chain"]] <- NULL
-  draws[[".iteration"]] <- NULL
-  draws[[".draw"]] <- NULL
+  draws[names_mcmc] <- NULL
   tibble::tibble(
     group = names_group(draws),
     time = names_time(draws),

From 0992d7875bddfa02ebdea76e3f29a7528909ccea Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 11:24:58 -0400
Subject: [PATCH 07/17] Add a test

---
 tests/testthat/test-brm_marginals.R | 90 +++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)
 create mode 100644 tests/testthat/test-brm_marginals.R

diff --git a/tests/testthat/test-brm_marginals.R b/tests/testthat/test-brm_marginals.R
new file mode 100644
index 00000000..914620da
--- /dev/null
+++ b/tests/testthat/test-brm_marginals.R
@@ -0,0 +1,90 @@
+test_that("brm_marginals() on response", {
+  set.seed(0L)
+  sim <- brm_simulate()
+  data <- sim$data
+  data$group <- paste("treatment", data$group)
+  data$time <- paste("visit", data$time)
+  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
+        )
+      )
+    )
+  )
+  suppressWarnings(
+    out <- brm_marginals(
+      model = model,
+      group = "group",
+      time = "time",
+      patient = "patient",
+      control = "treatment 1",
+      baseline = "visit 1",
+      outcome = "response"
+    )
+  )
+  fields <- c("response", "change", "difference")
+  columns_df <- expand.grid(
+    group = sort(unique(data$group)),
+    time = sort(unique(data$time)),
+    stringsAsFactors = FALSE
+  )
+  columns <- paste(columns_df$group, columns_df$time, sep = ", ")
+  expect_equal(sort(names(out)), sort(fields))
+  for (field in fields) {
+    x <- out[[field]]
+    expect_true(tibble::is_tibble(x))
+    expect_true(all(colnames(x) %in% c(columns, names_mcmc)))
+    expect_false(any(unlist(lapply(x, anyNA))))
+    expect_equal(nrow(x), 50)
+  }
+  expect_equal(
+    sort(colnames(out$response)),
+    sort(c(columns, names_mcmc))
+  )
+  columns_df <- columns_df[columns_df$time != "visit 1", ]
+  columns <- paste(columns_df$group, columns_df$time, sep = ", ")
+  expect_equal(
+    sort(colnames(out$change)),
+    sort(c(columns, names_mcmc))
+  )
+  columns_df <- columns_df[columns_df$group != "treatment 1", ]
+  columns <- paste(columns_df$group, columns_df$time, sep = ", ")
+  expect_equal(
+    sort(colnames(out$difference)),
+    sort(c(columns, names_mcmc))
+  )
+  for (group in setdiff(unique(data$group), "treatment 1")) {
+    for (time in setdiff(unique(data$time), "visit 1")) {
+      name1 <- paste("treatment 1", time, sep = ", ")
+      name2 <- paste(group, time, sep = ", ")
+      expect_equal(
+        out$difference[[name2]],
+        out$change[[name2]] - out$change[[name1]]
+      )
+    }
+  }
+  for (group in unique(data$group)) {
+    for (time in setdiff(unique(data$time), "visit 1")) {
+      name1 <- paste(group, "visit 1", sep = ", ")
+      name2 <- paste(group, time, sep = ", ")
+      expect_equal(
+        out$change[[name2]],
+        out$response[[name2]] - out$response[[name1]]
+      )
+    }
+  }
+})

From 6108d34c87db7475f518ee6a45a8b454593bc6dc Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 11:26:38 -0400
Subject: [PATCH 08/17] finish testing brm_marginals()

---
 tests/testthat/test-brm_marginals.R | 75 +++++++++++++++++++++++++++++
 1 file changed, 75 insertions(+)

diff --git a/tests/testthat/test-brm_marginals.R b/tests/testthat/test-brm_marginals.R
index 914620da..d601f719 100644
--- a/tests/testthat/test-brm_marginals.R
+++ b/tests/testthat/test-brm_marginals.R
@@ -88,3 +88,78 @@ test_that("brm_marginals() on response", {
     }
   }
 })
+
+test_that("brm_marginals() on change", {
+  set.seed(0L)
+  sim <- brm_simulate()
+  data <- sim$data
+  data$group <- paste("treatment", data$group)
+  data$time <- paste("visit", data$time)
+  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
+        )
+      )
+    )
+  )
+  suppressWarnings(
+    out <- brm_marginals(
+      model = model,
+      group = "group",
+      time = "time",
+      patient = "patient",
+      control = "treatment 1",
+      baseline = "visit 1",
+      outcome = "change"
+    )
+  )
+  fields <- c("response", "difference")
+  columns_df <- expand.grid(
+    group = sort(unique(data$group)),
+    time = sort(unique(data$time)),
+    stringsAsFactors = FALSE
+  )
+  columns <- paste(columns_df$group, columns_df$time, sep = ", ")
+  expect_equal(sort(names(out)), sort(fields))
+  for (field in fields) {
+    x <- out[[field]]
+    expect_true(tibble::is_tibble(x))
+    expect_true(all(colnames(x) %in% c(columns, names_mcmc)))
+    expect_false(any(unlist(lapply(x, anyNA))))
+    expect_equal(nrow(x), 50)
+  }
+  expect_equal(
+    sort(colnames(out$response)),
+    sort(c(columns, names_mcmc))
+  )
+  columns_df <- columns_df[columns_df$group != "treatment 1", ]
+  columns <- paste(columns_df$group, columns_df$time, sep = ", ")
+  expect_equal(
+    sort(colnames(out$difference)),
+    sort(c(columns, names_mcmc))
+  )
+  for (group in setdiff(unique(data$group), "treatment 1")) {
+    for (time in unique(data$time)) {
+      name1 <- paste("treatment 1", time, sep = ", ")
+      name2 <- paste(group, time, sep = ", ")
+      expect_equal(
+        out$difference[[name2]],
+        out$response[[name2]] - out$response[[name1]]
+      )
+    }
+  }
+})

From f1f983bfecbcbec753c8f00ff0436d5a0f4a22bf Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 11:50:24 -0400
Subject: [PATCH 09/17] Add a test

---
 R/brm_summary.R                     |   8 +-
 R/utils_data.R                      |   6 ++
 tests/testthat/test-brm_marginals.R |  36 +++----
 tests/testthat/test-brm_summary.R   | 149 ++++++++++++++++++++++------
 4 files changed, 143 insertions(+), 56 deletions(-)
 create mode 100644 R/utils_data.R

diff --git a/R/brm_summary.R b/R/brm_summary.R
index 69e9a167..18c5d858 100644
--- a/R/brm_summary.R
+++ b/R/brm_summary.R
@@ -157,25 +157,27 @@ summarize_marginals <- function(draws, level) {
     names_to = "statistic",
     values_to = "mcse"
   )
-  dplyr::left_join(
+  out <- dplyr::left_join(
     x = value,
     y = mcse,
     by = c("group", "time", "statistic")
   )
+  unname_df(out)
 }
 
 summarize_probabilities <- function(draws, direction, threshold) {
   draws[names_mcmc] <- NULL
-  tibble::tibble(
+  out <- tibble::tibble(
     group = names_group(draws),
     time = names_time(draws),
     direction = direction,
-    theshold = threshold,
+    threshold = threshold,
     value = purrr::map_dbl(
       draws,
       ~marginal_probability(.x, direction, threshold)
     )
   )
+  unname_df(out)
 }
 
 marginal_probability <- function(difference, direction, threshold) {
diff --git a/R/utils_data.R b/R/utils_data.R
new file mode 100644
index 00000000..c1e387a9
--- /dev/null
+++ b/R/utils_data.R
@@ -0,0 +1,6 @@
+unname_df <- function(x) {
+  for (i in seq_along(x)) {
+    x[[i]] <- unname(x[[i]])
+  }
+  x
+}
diff --git a/tests/testthat/test-brm_marginals.R b/tests/testthat/test-brm_marginals.R
index d601f719..462c158c 100644
--- a/tests/testthat/test-brm_marginals.R
+++ b/tests/testthat/test-brm_marginals.R
@@ -25,16 +25,14 @@ test_that("brm_marginals() on response", {
       )
     )
   )
-  suppressWarnings(
-    out <- brm_marginals(
-      model = model,
-      group = "group",
-      time = "time",
-      patient = "patient",
-      control = "treatment 1",
-      baseline = "visit 1",
-      outcome = "response"
-    )
+  out <- brm_marginals(
+    model = model,
+    group = "group",
+    time = "time",
+    patient = "patient",
+    control = "treatment 1",
+    baseline = "visit 1",
+    outcome = "response"
   )
   fields <- c("response", "change", "difference")
   columns_df <- expand.grid(
@@ -116,16 +114,14 @@ test_that("brm_marginals() on change", {
       )
     )
   )
-  suppressWarnings(
-    out <- brm_marginals(
-      model = model,
-      group = "group",
-      time = "time",
-      patient = "patient",
-      control = "treatment 1",
-      baseline = "visit 1",
-      outcome = "change"
-    )
+  out <- brm_marginals(
+    model = model,
+    group = "group",
+    time = "time",
+    patient = "patient",
+    control = "treatment 1",
+    baseline = "visit 1",
+    outcome = "change"
   )
   fields <- c("response", "difference")
   columns_df <- expand.grid(
diff --git a/tests/testthat/test-brm_summary.R b/tests/testthat/test-brm_summary.R
index ce7aa04c..8b07ae6d 100644
--- a/tests/testthat/test-brm_summary.R
+++ b/tests/testthat/test-brm_summary.R
@@ -1,14 +1,9 @@
-test_that("brm_summary()", {
+test_that("brm_summary() on response", {
   set.seed(0L)
-  sim <- brm_simulate(
-    n_group = 2L,
-    n_patient = 100L,
-    n_time = 4L,
-    hyper_beta = 1,
-    hyper_sigma = 1,
-    hyper_correlation = 1
-  )
+  sim <- brm_simulate()
   data <- sim$data
+  data$group <- paste("treatment", data$group)
+  data$time <- paste("visit", data$time)
   formula <- brm_formula(
     response = "response",
     group = "group",
@@ -30,34 +25,122 @@ test_that("brm_summary()", {
       )
     )
   )
-  out <- brm_summary(
+  marginals <- brm_marginals(
     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), times = 2))
-  expect_equal(as.integer(out$group), rep(seq_len(2), each = 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, na.rm = TRUE))
-  expect_true(all(out$diff_mean < out$diff_upper, na.rm = TRUE))
+    control = "treatment 1",
+    baseline = "visit 1",
+    outcome = "response"
+  )
+  suppressWarnings(
+    out <- brm_summary(
+      marginals,
+      level = 0.95,
+      threshold = 0,
+      direction = "greater"
+    )
+  )
+  expect_equal(sort(names(out)), sort(c("means", "probabilities")))
+  for (field in names(out)) {
+    x <- out[[field]]
+    expect_true(tibble::is_tibble(x))
+    expect_false(any(unlist(lapply(x, anyNA))))
+  }
+  x <- out$means
+  expect_equal(
+    sort(colnames(x)),
+    sort(c("marginal", "group", "time", "statistic", "value", "mcse"))
+  )
+  for (marginal in c("response", "change", "difference")) {
+    groups <- unique(data$group)
+    times <- unique(data$time)
+    if (marginal %in% c("change", "difference")) {
+      times <- setdiff(times, "visit 1")
+    }
+    if (identical(marginal, "difference")) {
+      groups <- setdiff(groups, "difference")
+    }
+    for (group in groups) {
+      for (time in times) {
+        name <- paste(group, time, sep = ", ")
+        index <- x$marginal == "response" & x$group == group & x$time == time
+        subset <- x[index, ]
+        expect_equal(
+          unname(subset$value[subset$statistic == "mean"]),
+          mean(marginals$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "median"]),
+          median(marginals$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "sd"]),
+          sd(marginals$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "lower"]),
+          unname(quantile(marginals$response[[name]], probs = 0.025))
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "upper"]),
+          unname(quantile(marginals$response[[name]], probs = 0.975))
+        )
+        suppressWarnings({
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "mean"]),
+            posterior::mcse_mean(marginals$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "median"]),
+            posterior::mcse_median(marginals$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "sd"]),
+            posterior::mcse_sd(marginals$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "lower"]),
+            unname(
+              posterior::mcse_quantile(
+                marginals$response[[name]],
+                probs = 0.025
+              )
+            )
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "upper"]),
+            unname(
+              posterior::mcse_quantile(
+                marginals$response[[name]],
+                probs = 0.975
+              )
+            )
+          )
+        })
+      }
+    }
+  }
+  x <- out$probabilities
+  expect_equal(
+    sort(colnames(x)),
+    sort(c("group", "time", "direction", "threshold", "value"))
+  )
+  expect_equal(x$group, rep("treatment 2", 3))
+  expect_equal(x$time, paste("visit", seq(2, 4)))
+  expect_equal(x$direction, rep("greater", 3))
+  expect_equal(x$threshold, rep(0, 3))
+  expect_equal(
+    x$value[1L],
+    mean(marginals$difference[["treatment 2, visit 2"]] > 0)
+  )
+  expect_equal(
+    x$value[2L],
+    mean(marginals$difference[["treatment 2, visit 3"]] > 0)
+  )
   expect_equal(
-    out$diff_mean[seq(5, 8)],
-    out$response_mean[out$group == 2] - out$response_mean[out$group == 1],
-    tolerance = 0.01
+    x$value[3L],
+    mean(marginals$difference[["treatment 2, visit 4"]] > 0)
   )
 })

From 9a31f2a2567b2e2b3f6880260ae3a66e855545ff Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 11:58:00 -0400
Subject: [PATCH 10/17] work on testing

---
 tests/testthat/test-utils_data.R | 5 +++++
 1 file changed, 5 insertions(+)
 create mode 100644 tests/testthat/test-utils_data.R

diff --git a/tests/testthat/test-utils_data.R b/tests/testthat/test-utils_data.R
new file mode 100644
index 00000000..75aa4185
--- /dev/null
+++ b/tests/testthat/test-utils_data.R
@@ -0,0 +1,5 @@
+test_that("unname_df()", {
+  x <- unname_df(tibble::tibble(x = c(a = 1, b = 2), y = c(c = 3, d = 4)))
+  expect_null(names(x$x))
+  expect_null(names(x$y))
+})

From 7ffe2ab4f6d5feabfd435530aa0c854ac6db44ae Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 12:45:06 -0400
Subject: [PATCH 11/17] function names and tests

---
 NAMESPACE                                     |   5 +-
 R/{brm_marginals.R => brm_marginal_draws.R}   |  16 +-
 R/brm_marginal_probabilities.R                | 119 +++++++++
 R/brm_marginal_summaries.R                    | 125 +++++++++
 R/brm_summary.R                               | 189 --------------
 ...brm_marginals.Rd => brm_marginal_draws.Rd} |  27 +-
 man/brm_marginal_probabilities.Rd             |  84 +++++++
 man/brm_marginal_summaries.Rd                 |  80 ++++++
 man/brm_summary.Rd                            |  96 -------
 ..._marginals.R => test-brm_marginal_draws.R} |   8 +-
 ...ry.R => test-brm_marginal_probabilities.R} |   2 +-
 tests/testthat/test-brm_marginal_summaries.R  | 238 ++++++++++++++++++
 12 files changed, 673 insertions(+), 316 deletions(-)
 rename R/{brm_marginals.R => brm_marginal_draws.R} (94%)
 create mode 100644 R/brm_marginal_probabilities.R
 create mode 100644 R/brm_marginal_summaries.R
 delete mode 100644 R/brm_summary.R
 rename man/{brm_marginals.Rd => brm_marginal_draws.Rd} (85%)
 create mode 100644 man/brm_marginal_probabilities.Rd
 create mode 100644 man/brm_marginal_summaries.Rd
 delete mode 100644 man/brm_summary.Rd
 rename tests/testthat/{test-brm_marginals.R => test-brm_marginal_draws.R} (96%)
 rename tests/testthat/{test-brm_summary.R => test-brm_marginal_probabilities.R} (99%)
 create mode 100644 tests/testthat/test-brm_marginal_summaries.R

diff --git a/NAMESPACE b/NAMESPACE
index 960f2ed3..bd32fc12 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,10 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
 export(brm_formula)
-export(brm_marginals)
+export(brm_marginal_draws)
+export(brm_marginal_probabilities)
+export(brm_marginal_summaries)
 export(brm_model)
 export(brm_simulate)
-export(brm_summary)
 importFrom(MASS,mvrnorm)
 importFrom(brms,brm)
 importFrom(brms,brmsformula)
diff --git a/R/brm_marginals.R b/R/brm_marginal_draws.R
similarity index 94%
rename from R/brm_marginals.R
rename to R/brm_marginal_draws.R
index 0a028521..d731636f 100644
--- a/R/brm_marginals.R
+++ b/R/brm_marginal_draws.R
@@ -1,10 +1,8 @@
-#' @title MMRM marginal posterior samples.
+#' @title MCMC draws from the marginal posterior of an MMRM
 #' @export
-#' @family results
-#' @description Get marginal posteior samples from a fitted MMRM.
-#' @details Currently assumes the response variable is `CHG`
-#'   (change from baseline) and not `AVAL` (raw response).
-#' @return A named list of tibbles of MCMC samples of the marginal posterior
+#' @family marginals
+#' @description Get marginal posterior draws from a fitted MMRM.
+#' @return A named list of tibbles of MCMC draws of the marginal posterior
 #'   distribution of each treatment group and time point:
 #'   * `response`: on the scale of the response variable.
 #'   * `change`: change from baseline, where the `baseline` argument determines
@@ -55,7 +53,7 @@
 #'     )
 #'   )
 #' )
-#' brm_marginals(
+#' brm_marginal_draws(
 #'   model = model,
 #'   group = "group",
 #'   time = "time",
@@ -64,7 +62,7 @@
 #'   baseline = "visit 1",
 #'   outcome = "response"
 #' )
-brm_marginals <- function(
+brm_marginal_draws <- function(
   model,
   base = "BASE",
   group = "TRT01P",
@@ -208,7 +206,7 @@ subtract_control <- function(draws, groups, times, control) {
 }
 
 name_marginal <- function(group, time) {
-  sprintf("%s, %s", group, time)
+  sprintf("%s, %s", group , time)
 }
 
 names_group <- function(draws) {
diff --git a/R/brm_marginal_probabilities.R b/R/brm_marginal_probabilities.R
new file mode 100644
index 00000000..053c8c43
--- /dev/null
+++ b/R/brm_marginal_probabilities.R
@@ -0,0 +1,119 @@
+#' @title Marginal probabilities on the treatment effect for an MMRM.
+#' @export
+#' @family marginals
+#' @description Marginal probabilities on the treatment effect for an MMRM.
+#' @return A data frame of probabilities of the form
+#'   `Prob(treatment effect > threshold | data)` and/or
+#'   `Prob(treatment effect < threshold | data)`. It has one row per
+#'   probability and the following columns:
+#'     * `group`: treatment group.
+#'     * `time`: discrete time point,
+#'     * `direction`: direction of the comparison in the marginal probability:
+#'       `"greater"` for `>`, `"less"` for `<`
+#'     * `threshold`: treatment effect threshold in the probability statement.
+#'     * `value`: numeric value of the estimate of the probability.
+#' @inheritParams brm_marginal_summaries
+#' @param direction Character vector of the same length as `threshold`.
+#'   `"greater"` to compute the marginal posterior probability that the
+#'   treatment effect is greater than the threshold,
+#'   `"less"` to compute the marginal posterior probability that the
+#'   treatment effect is less than the threshold.
+#'   Each element `direction[i]` corresponds to `threshold[i]`
+#'   for all `i` from 1 to `length(direction)`.
+#' @param threshold Numeric vector of the same length as `direction`,
+#'   treatment effect threshold for computing posterior probabilities.
+#'   Each element `direction[i]` corresponds to `threshold[i]` for
+#'   all `i` from 1 to `length(direction)`.
+#' @examples
+#' set.seed(0L)
+#' sim <- brm_simulate()
+#' data <- sim$data
+#' data$group <- paste("treatment", data$group)
+#' data$time <- paste("visit", data$time)
+#' 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
+#'       )
+#'     )
+#'   )
+#' )
+#' draws <- brm_marginal_draws(
+#'   model = model,
+#'   group = "group",
+#'   time = "time",
+#'   patient = "patient",
+#'   control = "treatment 1",
+#'   baseline = "visit 1",
+#'   outcome = "response"
+#' )
+#' brm_marginal_probabilities(draws)
+brm_marginal_probabilities <- function(
+  draws,
+  direction = "greater",
+  threshold = 0
+) {
+  assert(
+    is.list(draws),
+    message = "draws arg must be a named list from brm_marginal_draws()"
+  )
+  assert(
+    direction,
+    is.character(.),
+    !anyNA(.),
+    nzchar(.),
+    . %in% c("greater", "less"),
+    message = "elements of the direction arg must be \"greater\" or \"less\""
+  )
+  assert(
+    threshold,
+    is.numeric(.),
+    is.finite(.),
+    message = "threshold arg must be a numeric vector"
+  )
+  assert(
+    length(direction) == length(threshold),
+    message = "direction and threshold must have the same length"
+  )
+  summarize_probabilities(
+    draws = draws$difference,
+    direction = direction,
+    threshold = threshold
+  )
+}
+
+summarize_probabilities <- function(draws, direction, threshold) {
+  draws[names_mcmc] <- NULL
+  out <- tibble::tibble(
+    group = names_group(draws),
+    time = names_time(draws),
+    direction = direction,
+    threshold = threshold,
+    value = purrr::map_dbl(
+      draws,
+      ~marginal_probability(.x, direction, threshold)
+    )
+  )
+  unname_df(out)
+}
+
+marginal_probability <- function(difference, direction, threshold) {
+  if_any(
+    direction == "greater",
+    mean(difference > threshold),
+    mean(difference < threshold)
+  )
+}
diff --git a/R/brm_marginal_summaries.R b/R/brm_marginal_summaries.R
new file mode 100644
index 00000000..c5ef9d39
--- /dev/null
+++ b/R/brm_marginal_summaries.R
@@ -0,0 +1,125 @@
+#' @title Summary statistics of the marginal posterior of an MMRM.
+#' @export
+#' @family marginals
+#' @description Summary statistics of the marginal posterior of an MMRM.
+#' @return A data frame with one row per summary statistic and the following
+#'   columns:
+#'    * `marginal`: type of marginal distribution. If `outcome` was `"response"`
+#'      in [brm_marginal_draws()], then possible values include
+#'      `"response"` for the response on the raw scale, `"change"` for
+#'      change from baseline, and `"difference"` for treatment difference
+#'      in terms of change from baseline. If `outcome` was `"change"`,
+#'      then possible values include `"response"` for the respons one the
+#'      change from baseline scale and `"difference"` for treatment difference.
+#'    * `group`: treatment group.
+#'    * `time`: discrete time point.
+#'    * `statistic`: type of summary statistic.
+#'    * `value`: numeric value of the estimate.
+#'    * `mcse`: Monte Carlo standard error of the estimate.
+#' @param draws Posterior draws of the marginal posterior
+#'   obtained from [brm_marginal_draws()].
+#' @param level Numeric of length 1 between 0 and 1, credible level
+#'   for the credible intervals.
+#' @examples
+#' set.seed(0L)
+#' sim <- brm_simulate()
+#' data <- sim$data
+#' data$group <- paste("treatment", data$group)
+#' data$time <- paste("visit", data$time)
+#' 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
+#'       )
+#'     )
+#'   )
+#' )
+#' draws <- brm_marginal_draws(
+#'   model = model,
+#'   group = "group",
+#'   time = "time",
+#'   patient = "patient",
+#'   control = "treatment 1",
+#'   baseline = "visit 1",
+#'   outcome = "response"
+#' )
+#' brm_marginal_summaries(draws)
+brm_marginal_summaries <- function(
+  draws,
+  level = 0.95
+) {
+  assert(
+    is.list(draws),
+    message = "marginals arg must be a named list from brm_marginal_draws()"
+  )
+  assert_num(level, "level arg must be a length-1 numeric between 0 and 1")
+  assert(level, . >= 0, . <= 1, message = "level arg must be between 0 and 1")
+  table_response <- summarize_marginals(draws$response, level)
+  table_change <- if_any(
+    "change" %in% names(marginals),
+    summarize_marginals(draws$change, level),
+    NULL
+  )
+  table_difference <- summarize_marginals(draws$difference, level)
+  dplyr::bind_rows(
+    response = table_response,
+    change = table_change,
+    difference = table_difference,
+    .id = "marginal"
+  )
+}
+
+summarize_marginals <- function(draws, level) {
+  level_lower <- (1 - level) / 2
+  level_upper <- 1 - level_lower
+  draws[names_mcmc] <- NULL
+  value <- tibble::tibble(
+    group = names_group(draws),
+    time = names_time(draws),
+    mean = purrr::map_dbl(draws, mean),
+    median = purrr::map_dbl(draws, median),
+    sd = purrr::map_dbl(draws, sd),
+    lower = purrr::map_dbl(draws, ~quantile(.x, level_lower)),
+    upper = purrr::map_dbl(draws, ~quantile(.x, level_upper))
+  )
+  mcse <- tibble::tibble(
+    group = names_group(draws),
+    time = names_time(draws),
+    mean = purrr::map_dbl(draws, posterior::mcse_mean),
+    median = purrr::map_dbl(draws, posterior::mcse_median),
+    sd = purrr::map_dbl(draws, posterior::mcse_sd),
+    lower = purrr::map_dbl(draws, ~posterior::mcse_quantile(.x, level_lower)),
+    upper = purrr::map_dbl(draws, ~posterior::mcse_quantile(.x, level_upper))
+  )
+  value <- tidyr::pivot_longer(
+    data = value,
+    cols = -any_of(c("group", "time")),
+    names_to = "statistic",
+    values_to = "value"
+  )
+  mcse <- tidyr::pivot_longer(
+    data = mcse,
+    cols = -any_of(c("group", "time")),
+    names_to = "statistic",
+    values_to = "mcse"
+  )
+  out <- dplyr::left_join(
+    x = value,
+    y = mcse,
+    by = c("group", "time", "statistic")
+  )
+  unname_df(out)
+}
diff --git a/R/brm_summary.R b/R/brm_summary.R
deleted file mode 100644
index 18c5d858..00000000
--- a/R/brm_summary.R
+++ /dev/null
@@ -1,189 +0,0 @@
-#' @title Summarize an MMRM.
-#' @export
-#' @family results
-#' @description Summarize a basic MMRM model fit.
-#' @return A named list of two data frames with one row per summary
-#'   statistic:
-#'   * `means`: summary statistics of the marginal posterior mean of
-#'     each treatment group and time point. Columns include `marginal`
-#'     for the type of marginal distribution (see the "Value" section of
-#'     the [brm_marginals()] help file for details), `group` for treatment
-#'     group, `time` for discrete time point, `statistic` for the type of
-#'     summary statistic, `value` for the numeric value of the estimate,
-#'     and `mcse` for the Monte Carlo standard error of the estimate.
-#'   * `probabilities`: marginal posterior probabilities of the form
-#'     `Prob(treatment effect > threshold | data)` and/or
-#'     `Prob(treatment effect < threshold | data)`. Columns include
-#'     `group` for the treatment group, `time` for the discrete time point,
-#'     `direction` to indicate the direction of the comparison
-#'     (`"greater"` for `>`, `"less"` for `<`), `threshold` for the
-#'     treatment effect threshold in the probability statement,
-#'     and `value` for the numberic value of the estimate of the
-#'     probability.
-#' @param marginals Posterior draws of the marginal posterior
-#'   obtained from [brm_marginals()].
-#' @param level Numeric of length 1 between 0 and 1, credible level
-#'   for the credible intervals. Only relevant when `return = "marginals"`.
-#' @param direction Character vector of the same length as `threshold`.
-#'   `"greater"` to compute the marginal posterior probability that the
-#'   treatment effect is greater than the threshold,
-#'   `"less"` to compute the marginal posterior probability that the
-#'   treatment effect is less than the threshold.
-#'   Each element `direction[i]` corresponds to `threshold[i]`
-#'   for all `i` from 1 to `length(direction)`.
-#' @param threshold Numeric vector of the same length as `direction`,
-#'   treatment effect threshold for computing posterior probabilities.
-#'   Each element `direction[i]` corresponds to `threshold[i]` for
-#'   all `i` from 1 to `length(direction)`.
-#' @examples
-#' set.seed(0L)
-#' sim <- brm_simulate()
-#' data <- sim$data
-#' data$group <- paste("treatment", data$group)
-#' data$time <- paste("visit", data$time)
-#' 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
-#'       )
-#'     )
-#'   )
-#' )
-#' marginals <- brm_marginals(
-#'   model = model,
-#'   group = "group",
-#'   time = "time",
-#'   patient = "patient",
-#'   control = "treatment 1",
-#'   baseline = "visit 1",
-#'   outcome = "response"
-#' )
-#' brm_summary(marginals)
-brm_summary <- function(
-  marginals,
-  level = 0.95,
-  direction = "greater",
-  threshold = 0
-) {
-  assert(
-    is.list(marginals),
-    message = "marginals arg must be a named list from brm_marginals()"
-  )
-  assert_num(level, "level arg must be a length-1 numeric between 0 and 1")
-  assert(level, . >= 0, . <= 1, message = "level arg must be between 0 and 1")
-  assert(
-    direction,
-    is.character(.),
-    !anyNA(.),
-    nzchar(.),
-    . %in% c("greater", "less"),
-    message = "elements of the direction arg must be \"greater\" or \"less\""
-  )
-  assert(
-    threshold,
-    is.numeric(.),
-    is.finite(.),
-    message = "threshold arg must be a numeric vector"
-  )
-  assert(
-    length(direction) == length(threshold),
-    message = "direction and threshold must have the same length"
-  )
-  table_response <- summarize_marginals(marginals$response, level)
-  table_change <- if_any(
-    "change" %in% names(marginals),
-    summarize_marginals(marginals$change, level),
-    NULL
-  )
-  table_difference <- summarize_marginals(marginals$difference, level)
-  means <- dplyr::bind_rows(
-    response = table_response,
-    change = table_change,
-    difference = table_difference,
-    .id = "marginal"
-  )
-  probabilities <- summarize_probabilities(
-    draws = marginals$difference,
-    direction = direction,
-    threshold = threshold
-  )
-  list(means = means, probabilities = probabilities)
-}
-
-summarize_marginals <- function(draws, level) {
-  level_lower <- (1 - level) / 2
-  level_upper <- 1 - level_lower
-  draws[names_mcmc] <- NULL
-  value <- tibble::tibble(
-    group = names_group(draws),
-    time = names_time(draws),
-    mean = purrr::map_dbl(draws, mean),
-    median = purrr::map_dbl(draws, median),
-    sd = purrr::map_dbl(draws, sd),
-    lower = purrr::map_dbl(draws, ~quantile(.x, level_lower)),
-    upper = purrr::map_dbl(draws, ~quantile(.x, level_upper))
-  )
-  mcse <- tibble::tibble(
-    group = names_group(draws),
-    time = names_time(draws),
-    mean = purrr::map_dbl(draws, posterior::mcse_mean),
-    median = purrr::map_dbl(draws, posterior::mcse_median),
-    sd = purrr::map_dbl(draws, posterior::mcse_sd),
-    lower = purrr::map_dbl(draws, ~posterior::mcse_quantile(.x, level_lower)),
-    upper = purrr::map_dbl(draws, ~posterior::mcse_quantile(.x, level_upper))
-  )
-  value <- tidyr::pivot_longer(
-    data = value,
-    cols = -any_of(c("group", "time")),
-    names_to = "statistic",
-    values_to = "value"
-  )
-  mcse <- tidyr::pivot_longer(
-    data = mcse,
-    cols = -any_of(c("group", "time")),
-    names_to = "statistic",
-    values_to = "mcse"
-  )
-  out <- dplyr::left_join(
-    x = value,
-    y = mcse,
-    by = c("group", "time", "statistic")
-  )
-  unname_df(out)
-}
-
-summarize_probabilities <- function(draws, direction, threshold) {
-  draws[names_mcmc] <- NULL
-  out <- tibble::tibble(
-    group = names_group(draws),
-    time = names_time(draws),
-    direction = direction,
-    threshold = threshold,
-    value = purrr::map_dbl(
-      draws,
-      ~marginal_probability(.x, direction, threshold)
-    )
-  )
-  unname_df(out)
-}
-
-marginal_probability <- function(difference, direction, threshold) {
-  if_any(
-    direction == "greater",
-    mean(difference > threshold),
-    mean(difference < threshold)
-  )
-}
diff --git a/man/brm_marginals.Rd b/man/brm_marginal_draws.Rd
similarity index 85%
rename from man/brm_marginals.Rd
rename to man/brm_marginal_draws.Rd
index 8f4d77e7..28e79683 100644
--- a/man/brm_marginals.Rd
+++ b/man/brm_marginal_draws.Rd
@@ -1,10 +1,10 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/brm_marginals.R
-\name{brm_marginals}
-\alias{brm_marginals}
-\title{MMRM marginal posterior samples.}
+% Please edit documentation in R/brm_marginal_draws.R
+\name{brm_marginal_draws}
+\alias{brm_marginal_draws}
+\title{MCMC draws from the marginal posterior of an MMRM}
 \usage{
-brm_marginals(
+brm_marginal_draws(
   model,
   base = "BASE",
   group = "TRT01P",
@@ -47,7 +47,7 @@ which indicates the baseline time for the purposes of calculating
 change from baseline.}
 }
 \value{
-A named list of tibbles of MCMC samples of the marginal posterior
+A named list of tibbles of MCMC draws of the marginal posterior
 distribution of each treatment group and time point:
 \itemize{
 \item \code{response}: on the scale of the response variable.
@@ -64,11 +64,7 @@ Treatment and time are comma-delimited in the column names.
 }
 }
 \description{
-Get marginal posteior samples from a fitted MMRM.
-}
-\details{
-Currently assumes the response variable is \code{CHG}
-(change from baseline) and not \code{AVAL} (raw response).
+Get marginal posterior draws from a fitted MMRM.
 }
 \examples{
 set.seed(0L)
@@ -97,7 +93,7 @@ tmp <- utils::capture.output(
     )
   )
 )
-brm_marginals(
+brm_marginal_draws(
   model = model,
   group = "group",
   time = "time",
@@ -108,7 +104,8 @@ brm_marginals(
 )
 }
 \seealso{
-Other results: 
-\code{\link{brm_summary}()}
+Other marginals: 
+\code{\link{brm_marginal_probabilities}()},
+\code{\link{brm_marginal_summaries}()}
 }
-\concept{results}
+\concept{marginals}
diff --git a/man/brm_marginal_probabilities.Rd b/man/brm_marginal_probabilities.Rd
new file mode 100644
index 00000000..e6ce0831
--- /dev/null
+++ b/man/brm_marginal_probabilities.Rd
@@ -0,0 +1,84 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/brm_marginal_probabilities.R
+\name{brm_marginal_probabilities}
+\alias{brm_marginal_probabilities}
+\title{Marginal probabilities on the treatment effect for an MMRM.}
+\usage{
+brm_marginal_probabilities(draws, direction = "greater", threshold = 0)
+}
+\arguments{
+\item{draws}{Posterior draws of the marginal posterior
+obtained from \code{\link[=brm_marginal_draws]{brm_marginal_draws()}}.}
+
+\item{direction}{Character vector of the same length as \code{threshold}.
+\code{"greater"} to compute the marginal posterior probability that the
+treatment effect is greater than the threshold,
+\code{"less"} to compute the marginal posterior probability that the
+treatment effect is less than the threshold.
+Each element \code{direction[i]} corresponds to \code{threshold[i]}
+for all \code{i} from 1 to \code{length(direction)}.}
+
+\item{threshold}{Numeric vector of the same length as \code{direction},
+treatment effect threshold for computing posterior probabilities.
+Each element \code{direction[i]} corresponds to \code{threshold[i]} for
+all \code{i} from 1 to \code{length(direction)}.}
+}
+\value{
+A data frame of probabilities of the form
+\verb{Prob(treatment effect > threshold | data)} and/or
+\verb{Prob(treatment effect < threshold | data)}. It has one row per
+probability and the following columns:
+* \code{group}: treatment group.
+* \code{time}: discrete time point,
+* \code{direction}: direction of the comparison in the marginal probability:
+\code{"greater"} for \code{>}, \code{"less"} for \code{<}
+* \code{threshold}: treatment effect threshold in the probability statement.
+* \code{value}: numeric value of the estimate of the probability.
+}
+\description{
+Marginal probabilities on the treatment effect for an MMRM.
+}
+\examples{
+set.seed(0L)
+sim <- brm_simulate()
+data <- sim$data
+data$group <- paste("treatment", data$group)
+data$time <- paste("visit", data$time)
+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
+      )
+    )
+  )
+)
+draws <- brm_marginal_draws(
+  model = model,
+  group = "group",
+  time = "time",
+  patient = "patient",
+  control = "treatment 1",
+  baseline = "visit 1",
+  outcome = "response"
+)
+brm_marginal_probabilities(draws)
+}
+\seealso{
+Other marginals: 
+\code{\link{brm_marginal_draws}()},
+\code{\link{brm_marginal_summaries}()}
+}
+\concept{marginals}
diff --git a/man/brm_marginal_summaries.Rd b/man/brm_marginal_summaries.Rd
new file mode 100644
index 00000000..ace6cd94
--- /dev/null
+++ b/man/brm_marginal_summaries.Rd
@@ -0,0 +1,80 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/brm_marginal_summaries.R
+\name{brm_marginal_summaries}
+\alias{brm_marginal_summaries}
+\title{Summary statistics of the marginal posterior of an MMRM.}
+\usage{
+brm_marginal_summaries(draws, level = 0.95)
+}
+\arguments{
+\item{draws}{Posterior draws of the marginal posterior
+obtained from \code{\link[=brm_marginal_draws]{brm_marginal_draws()}}.}
+
+\item{level}{Numeric of length 1 between 0 and 1, credible level
+for the credible intervals.}
+}
+\value{
+A data frame with one row per summary statistic and the following
+columns:
+\itemize{
+\item \code{marginal}: type of marginal distribution. If \code{outcome} was \code{"response"}
+in \code{\link[=brm_marginal_draws]{brm_marginal_draws()}}, then possible values include
+\code{"response"} for the response on the raw scale, \code{"change"} for
+change from baseline, and \code{"difference"} for treatment difference
+in terms of change from baseline. If \code{outcome} was \code{"change"},
+then possible values include \code{"response"} for the respons one the
+change from baseline scale and \code{"difference"} for treatment difference.
+\item \code{group}: treatment group.
+\item \code{time}: discrete time point.
+\item \code{statistic}: type of summary statistic.
+\item \code{value}: numeric value of the estimate.
+\item \code{mcse}: Monte Carlo standard error of the estimate.
+}
+}
+\description{
+Summary statistics of the marginal posterior of an MMRM.
+}
+\examples{
+set.seed(0L)
+sim <- brm_simulate()
+data <- sim$data
+data$group <- paste("treatment", data$group)
+data$time <- paste("visit", data$time)
+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
+      )
+    )
+  )
+)
+draws <- brm_marginal_draws(
+  model = model,
+  group = "group",
+  time = "time",
+  patient = "patient",
+  control = "treatment 1",
+  baseline = "visit 1",
+  outcome = "response"
+)
+brm_marginal_summaries(draws)
+}
+\seealso{
+Other marginals: 
+\code{\link{brm_marginal_draws}()},
+\code{\link{brm_marginal_probabilities}()}
+}
+\concept{marginals}
diff --git a/man/brm_summary.Rd b/man/brm_summary.Rd
deleted file mode 100644
index a3ea1735..00000000
--- a/man/brm_summary.Rd
+++ /dev/null
@@ -1,96 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/brm_summary.R
-\name{brm_summary}
-\alias{brm_summary}
-\title{Summarize an MMRM.}
-\usage{
-brm_summary(marginals, level = 0.95, direction = "greater", threshold = 0)
-}
-\arguments{
-\item{marginals}{Posterior draws of the marginal posterior
-obtained from \code{\link[=brm_marginals]{brm_marginals()}}.}
-
-\item{level}{Numeric of length 1 between 0 and 1, credible level
-for the credible intervals. Only relevant when \code{return = "marginals"}.}
-
-\item{direction}{Character vector of the same length as \code{threshold}.
-\code{"greater"} to compute the marginal posterior probability that the
-treatment effect is greater than the threshold,
-\code{"less"} to compute the marginal posterior probability that the
-treatment effect is less than the threshold.
-Each element \code{direction[i]} corresponds to \code{threshold[i]}
-for all \code{i} from 1 to \code{length(direction)}.}
-
-\item{threshold}{Numeric vector of the same length as \code{direction},
-treatment effect threshold for computing posterior probabilities.
-Each element \code{direction[i]} corresponds to \code{threshold[i]} for
-all \code{i} from 1 to \code{length(direction)}.}
-}
-\value{
-A named list of two data frames with one row per summary
-statistic:
-\itemize{
-\item \code{means}: summary statistics of the marginal posterior mean of
-each treatment group and time point. Columns include \code{marginal}
-for the type of marginal distribution (see the "Value" section of
-the \code{\link[=brm_marginals]{brm_marginals()}} help file for details), \code{group} for treatment
-group, \code{time} for discrete time point, \code{statistic} for the type of
-summary statistic, \code{value} for the numeric value of the estimate,
-and \code{mcse} for the Monte Carlo standard error of the estimate.
-\item \code{probabilities}: marginal posterior probabilities of the form
-\verb{Prob(treatment effect > threshold | data)} and/or
-\verb{Prob(treatment effect < threshold | data)}. Columns include
-\code{group} for the treatment group, \code{time} for the discrete time point,
-\code{direction} to indicate the direction of the comparison
-(\code{"greater"} for \code{>}, \code{"less"} for \code{<}), \code{threshold} for the
-treatment effect threshold in the probability statement,
-and \code{value} for the numberic value of the estimate of the
-probability.
-}
-}
-\description{
-Summarize a basic MMRM model fit.
-}
-\examples{
-set.seed(0L)
-sim <- brm_simulate()
-data <- sim$data
-data$group <- paste("treatment", data$group)
-data$time <- paste("visit", data$time)
-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
-      )
-    )
-  )
-)
-marginals <- brm_marginals(
-  model = model,
-  group = "group",
-  time = "time",
-  patient = "patient",
-  control = "treatment 1",
-  baseline = "visit 1",
-  outcome = "response"
-)
-brm_summary(marginals)
-}
-\seealso{
-Other results: 
-\code{\link{brm_marginals}()}
-}
-\concept{results}
diff --git a/tests/testthat/test-brm_marginals.R b/tests/testthat/test-brm_marginal_draws.R
similarity index 96%
rename from tests/testthat/test-brm_marginals.R
rename to tests/testthat/test-brm_marginal_draws.R
index 462c158c..92927842 100644
--- a/tests/testthat/test-brm_marginals.R
+++ b/tests/testthat/test-brm_marginal_draws.R
@@ -1,4 +1,4 @@
-test_that("brm_marginals() on response", {
+test_that("brm_marginal_draws() on response", {
   set.seed(0L)
   sim <- brm_simulate()
   data <- sim$data
@@ -25,7 +25,7 @@ test_that("brm_marginals() on response", {
       )
     )
   )
-  out <- brm_marginals(
+  out <- brm_marginal_draws(
     model = model,
     group = "group",
     time = "time",
@@ -87,7 +87,7 @@ test_that("brm_marginals() on response", {
   }
 })
 
-test_that("brm_marginals() on change", {
+test_that("brm_marginal_draws() on change", {
   set.seed(0L)
   sim <- brm_simulate()
   data <- sim$data
@@ -114,7 +114,7 @@ test_that("brm_marginals() on change", {
       )
     )
   )
-  out <- brm_marginals(
+  out <- brm_marginal_draws(
     model = model,
     group = "group",
     time = "time",
diff --git a/tests/testthat/test-brm_summary.R b/tests/testthat/test-brm_marginal_probabilities.R
similarity index 99%
rename from tests/testthat/test-brm_summary.R
rename to tests/testthat/test-brm_marginal_probabilities.R
index 8b07ae6d..9b2a08ff 100644
--- a/tests/testthat/test-brm_summary.R
+++ b/tests/testthat/test-brm_marginal_probabilities.R
@@ -25,7 +25,7 @@ test_that("brm_summary() on response", {
       )
     )
   )
-  marginals <- brm_marginals(
+  marginals <- brm_marginal_draws(
     model = model,
     group = "group",
     time = "time",
diff --git a/tests/testthat/test-brm_marginal_summaries.R b/tests/testthat/test-brm_marginal_summaries.R
new file mode 100644
index 00000000..e52995b7
--- /dev/null
+++ b/tests/testthat/test-brm_marginal_summaries.R
@@ -0,0 +1,238 @@
+test_that("brm_marginal_summaries() on response", {
+  set.seed(0L)
+  sim <- brm_simulate()
+  data <- sim$data
+  data$group <- paste("treatment", data$group)
+  data$time <- paste("visit", data$time)
+  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
+        )
+      )
+    )
+  )
+  draws <- brm_marginal_draws(
+    model = model,
+    group = "group",
+    time = "time",
+    patient = "patient",
+    control = "treatment 1",
+    baseline = "visit 1",
+    outcome = "response"
+  )
+  suppressWarnings(
+    x <- brm_marginal_summaries(
+      draws,
+      level = 0.95
+    )
+  )
+  expect_equal(
+    sort(colnames(x)),
+    sort(c("marginal", "group", "time", "statistic", "value", "mcse"))
+  )
+  expect_equal(
+    sort(unique(x$marginal)),
+    sort(c("response", "change", "difference"))
+  )
+  for (marginal in c("response", "change", "difference")) {
+    groups <- unique(data$group)
+    times <- unique(data$time)
+    if (marginal %in% c("change", "difference")) {
+      times <- setdiff(times, "visit 1")
+    }
+    if (identical(marginal, "difference")) {
+      groups <- setdiff(groups, "difference")
+    }
+    for (group in groups) {
+      for (time in times) {
+        name <- paste(group, time, sep = ", ")
+        index <- x$marginal == "response" & x$group == group & x$time == time
+        subset <- x[index, ]
+        expect_equal(
+          unname(subset$value[subset$statistic == "mean"]),
+          mean(draws$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "median"]),
+          median(draws$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "sd"]),
+          sd(draws$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "lower"]),
+          unname(quantile(draws$response[[name]], probs = 0.025))
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "upper"]),
+          unname(quantile(draws$response[[name]], probs = 0.975))
+        )
+        suppressWarnings({
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "mean"]),
+            posterior::mcse_mean(draws$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "median"]),
+            posterior::mcse_median(draws$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "sd"]),
+            posterior::mcse_sd(draws$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "lower"]),
+            unname(
+              posterior::mcse_quantile(
+                draws$response[[name]],
+                probs = 0.025
+              )
+            )
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "upper"]),
+            unname(
+              posterior::mcse_quantile(
+                draws$response[[name]],
+                probs = 0.975
+              )
+            )
+          )
+        })
+      }
+    }
+  }
+})
+
+test_that("brm_marginal_summaries() on change", {
+  set.seed(0L)
+  sim <- brm_simulate()
+  data <- sim$data
+  data$group <- paste("treatment", data$group)
+  data$time <- paste("visit", data$time)
+  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
+        )
+      )
+    )
+  )
+  draws <- brm_marginal_draws(
+    model = model,
+    group = "group",
+    time = "time",
+    patient = "patient",
+    control = "treatment 1",
+    baseline = "visit 1",
+    outcome = "change"
+  )
+  suppressWarnings(
+    x <- brm_marginal_summaries(
+      draws,
+      level = 0.95
+    )
+  )
+  expect_equal(
+    sort(colnames(x)),
+    sort(c("marginal", "group", "time", "statistic", "value", "mcse"))
+  )
+  expect_equal(
+    sort(unique(x$marginal)),
+    sort(c("response", "difference"))
+  )
+  for (marginal in c("response", "difference")) {
+    groups <- unique(data$group)
+    times <- unique(data$time)
+    if (identical(marginal, "difference")) {
+      groups <- setdiff(groups, "difference")
+    }
+    for (group in groups) {
+      for (time in times) {
+        name <- paste(group, time, sep = ", ")
+        index <- x$marginal == "response" & x$group == group & x$time == time
+        subset <- x[index, ]
+        expect_equal(
+          unname(subset$value[subset$statistic == "mean"]),
+          mean(draws$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "median"]),
+          median(draws$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "sd"]),
+          sd(draws$response[[name]])
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "lower"]),
+          unname(quantile(draws$response[[name]], probs = 0.025))
+        )
+        expect_equal(
+          unname(subset$value[subset$statistic == "upper"]),
+          unname(quantile(draws$response[[name]], probs = 0.975))
+        )
+        suppressWarnings({
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "mean"]),
+            posterior::mcse_mean(draws$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "median"]),
+            posterior::mcse_median(draws$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "sd"]),
+            posterior::mcse_sd(draws$response[[name]])
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "lower"]),
+            unname(
+              posterior::mcse_quantile(
+                draws$response[[name]],
+                probs = 0.025
+              )
+            )
+          )
+          expect_equal(
+            unname(subset$mcse[subset$statistic == "upper"]),
+            unname(
+              posterior::mcse_quantile(
+                draws$response[[name]],
+                probs = 0.975
+              )
+            )
+          )
+        })
+      }
+    }
+  }
+})

From 8b9afda1075d661a34856ff06064faaa491c70a1 Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 13:16:05 -0400
Subject: [PATCH 12/17] Test brm_marginal_probabilities()

---
 R/brm_marginal_probabilities.R                |  35 ++--
 R/brm_marginal_summaries.R                    |  16 +-
 man/brm_marginal_probabilities.Rd             |   2 +-
 man/brm_marginal_summaries.Rd                 |   4 +-
 .../test-brm_marginal_probabilities.R         | 186 +++++++++---------
 tests/testthat/test-brm_marginal_summaries.R  |  10 +-
 6 files changed, 137 insertions(+), 116 deletions(-)

diff --git a/R/brm_marginal_probabilities.R b/R/brm_marginal_probabilities.R
index 053c8c43..70611cbb 100644
--- a/R/brm_marginal_probabilities.R
+++ b/R/brm_marginal_probabilities.R
@@ -60,7 +60,7 @@
 #'   baseline = "visit 1",
 #'   outcome = "response"
 #' )
-#' brm_marginal_probabilities(draws)
+#' brm_marginal_probabilities(draws, direction = "greater", threshold = 0)
 brm_marginal_probabilities <- function(
   draws,
   direction = "greater",
@@ -88,26 +88,39 @@ brm_marginal_probabilities <- function(
     length(direction) == length(threshold),
     message = "direction and threshold must have the same length"
   )
-  summarize_probabilities(
-    draws = draws$difference,
-    direction = direction,
-    threshold = threshold
+  draws <- tibble::as_tibble(draws$difference)
+  for (name in names_mcmc) {
+    draws[[name]] <- NULL
+  }
+  out <- purrr::map2_df(
+    .x = direction,
+    .y = threshold,
+    .f = ~summarize_probabilities(
+      draws = draws,
+      direction = .x,
+      threshold = .y
+    )
   )
+  columns <- c("direction", "threshold", "group", "time", "value")
+  out <- out[, columns]
+  args <- lapply(columns, as.symbol)
+  args$.data <- out
+  do.call(what = dplyr::arrange, args = args)
 }
 
 summarize_probabilities <- function(draws, direction, threshold) {
-  draws[names_mcmc] <- NULL
+  values <- purrr::map_dbl(
+    draws,
+    ~marginal_probability(.x, direction, threshold)
+  )
   out <- tibble::tibble(
     group = names_group(draws),
     time = names_time(draws),
     direction = direction,
     threshold = threshold,
-    value = purrr::map_dbl(
-      draws,
-      ~marginal_probability(.x, direction, threshold)
-    )
+    value = values
   )
-  unname_df(out)
+  out <- unname_df(out)
 }
 
 marginal_probability <- function(difference, direction, threshold) {
diff --git a/R/brm_marginal_summaries.R b/R/brm_marginal_summaries.R
index c5ef9d39..201fe825 100644
--- a/R/brm_marginal_summaries.R
+++ b/R/brm_marginal_summaries.R
@@ -4,6 +4,8 @@
 #' @description Summary statistics of the marginal posterior of an MMRM.
 #' @return A data frame with one row per summary statistic and the following
 #'   columns:
+#'    * `group`: treatment group.
+#'    * `time`: discrete time point.
 #'    * `marginal`: type of marginal distribution. If `outcome` was `"response"`
 #'      in [brm_marginal_draws()], then possible values include
 #'      `"response"` for the response on the raw scale, `"change"` for
@@ -11,8 +13,6 @@
 #'      in terms of change from baseline. If `outcome` was `"change"`,
 #'      then possible values include `"response"` for the respons one the
 #'      change from baseline scale and `"difference"` for treatment difference.
-#'    * `group`: treatment group.
-#'    * `time`: discrete time point.
 #'    * `statistic`: type of summary statistic.
 #'    * `value`: numeric value of the estimate.
 #'    * `mcse`: Monte Carlo standard error of the estimate.
@@ -74,18 +74,26 @@ brm_marginal_summaries <- function(
     NULL
   )
   table_difference <- summarize_marginals(draws$difference, level)
-  dplyr::bind_rows(
+  out <- dplyr::bind_rows(
     response = table_response,
     change = table_change,
     difference = table_difference,
     .id = "marginal"
   )
+  columns <- c("marginal", "statistic", "group", "time", "value", "mcse")
+  out <- out[, columns]
+  args <- lapply(columns, as.symbol)
+  args$.data <- out
+  do.call(what = dplyr::arrange, args = args)
 }
 
 summarize_marginals <- function(draws, level) {
+  draws <- tibble::as_tibble(draws)
+  for (name in names_mcmc) {
+    draws[[name]] <- NULL
+  }
   level_lower <- (1 - level) / 2
   level_upper <- 1 - level_lower
-  draws[names_mcmc] <- NULL
   value <- tibble::tibble(
     group = names_group(draws),
     time = names_time(draws),
diff --git a/man/brm_marginal_probabilities.Rd b/man/brm_marginal_probabilities.Rd
index e6ce0831..00937662 100644
--- a/man/brm_marginal_probabilities.Rd
+++ b/man/brm_marginal_probabilities.Rd
@@ -74,7 +74,7 @@ draws <- brm_marginal_draws(
   baseline = "visit 1",
   outcome = "response"
 )
-brm_marginal_probabilities(draws)
+brm_marginal_probabilities(draws, direction = "greater", threshold = 0)
 }
 \seealso{
 Other marginals: 
diff --git a/man/brm_marginal_summaries.Rd b/man/brm_marginal_summaries.Rd
index ace6cd94..ee9fa81d 100644
--- a/man/brm_marginal_summaries.Rd
+++ b/man/brm_marginal_summaries.Rd
@@ -17,6 +17,8 @@ for the credible intervals.}
 A data frame with one row per summary statistic and the following
 columns:
 \itemize{
+\item \code{group}: treatment group.
+\item \code{time}: discrete time point.
 \item \code{marginal}: type of marginal distribution. If \code{outcome} was \code{"response"}
 in \code{\link[=brm_marginal_draws]{brm_marginal_draws()}}, then possible values include
 \code{"response"} for the response on the raw scale, \code{"change"} for
@@ -24,8 +26,6 @@ change from baseline, and \code{"difference"} for treatment difference
 in terms of change from baseline. If \code{outcome} was \code{"change"},
 then possible values include \code{"response"} for the respons one the
 change from baseline scale and \code{"difference"} for treatment difference.
-\item \code{group}: treatment group.
-\item \code{time}: discrete time point.
 \item \code{statistic}: type of summary statistic.
 \item \code{value}: numeric value of the estimate.
 \item \code{mcse}: Monte Carlo standard error of the estimate.
diff --git a/tests/testthat/test-brm_marginal_probabilities.R b/tests/testthat/test-brm_marginal_probabilities.R
index 9b2a08ff..80ba8584 100644
--- a/tests/testthat/test-brm_marginal_probabilities.R
+++ b/tests/testthat/test-brm_marginal_probabilities.R
@@ -1,4 +1,4 @@
-test_that("brm_summary() on response", {
+test_that("brm_marginal_probabilities() on response", {
   set.seed(0L)
   sim <- brm_simulate()
   data <- sim$data
@@ -25,7 +25,7 @@ test_that("brm_summary() on response", {
       )
     )
   )
-  marginals <- brm_marginal_draws(
+  draws <- brm_marginal_draws(
     model = model,
     group = "group",
     time = "time",
@@ -34,95 +34,11 @@ test_that("brm_summary() on response", {
     baseline = "visit 1",
     outcome = "response"
   )
-  suppressWarnings(
-    out <- brm_summary(
-      marginals,
-      level = 0.95,
-      threshold = 0,
-      direction = "greater"
-    )
+  x <- brm_marginal_probabilities(
+    draws,
+    threshold = 0,
+    direction = "greater"
   )
-  expect_equal(sort(names(out)), sort(c("means", "probabilities")))
-  for (field in names(out)) {
-    x <- out[[field]]
-    expect_true(tibble::is_tibble(x))
-    expect_false(any(unlist(lapply(x, anyNA))))
-  }
-  x <- out$means
-  expect_equal(
-    sort(colnames(x)),
-    sort(c("marginal", "group", "time", "statistic", "value", "mcse"))
-  )
-  for (marginal in c("response", "change", "difference")) {
-    groups <- unique(data$group)
-    times <- unique(data$time)
-    if (marginal %in% c("change", "difference")) {
-      times <- setdiff(times, "visit 1")
-    }
-    if (identical(marginal, "difference")) {
-      groups <- setdiff(groups, "difference")
-    }
-    for (group in groups) {
-      for (time in times) {
-        name <- paste(group, time, sep = ", ")
-        index <- x$marginal == "response" & x$group == group & x$time == time
-        subset <- x[index, ]
-        expect_equal(
-          unname(subset$value[subset$statistic == "mean"]),
-          mean(marginals$response[[name]])
-        )
-        expect_equal(
-          unname(subset$value[subset$statistic == "median"]),
-          median(marginals$response[[name]])
-        )
-        expect_equal(
-          unname(subset$value[subset$statistic == "sd"]),
-          sd(marginals$response[[name]])
-        )
-        expect_equal(
-          unname(subset$value[subset$statistic == "lower"]),
-          unname(quantile(marginals$response[[name]], probs = 0.025))
-        )
-        expect_equal(
-          unname(subset$value[subset$statistic == "upper"]),
-          unname(quantile(marginals$response[[name]], probs = 0.975))
-        )
-        suppressWarnings({
-          expect_equal(
-            unname(subset$mcse[subset$statistic == "mean"]),
-            posterior::mcse_mean(marginals$response[[name]])
-          )
-          expect_equal(
-            unname(subset$mcse[subset$statistic == "median"]),
-            posterior::mcse_median(marginals$response[[name]])
-          )
-          expect_equal(
-            unname(subset$mcse[subset$statistic == "sd"]),
-            posterior::mcse_sd(marginals$response[[name]])
-          )
-          expect_equal(
-            unname(subset$mcse[subset$statistic == "lower"]),
-            unname(
-              posterior::mcse_quantile(
-                marginals$response[[name]],
-                probs = 0.025
-              )
-            )
-          )
-          expect_equal(
-            unname(subset$mcse[subset$statistic == "upper"]),
-            unname(
-              posterior::mcse_quantile(
-                marginals$response[[name]],
-                probs = 0.975
-              )
-            )
-          )
-        })
-      }
-    }
-  }
-  x <- out$probabilities
   expect_equal(
     sort(colnames(x)),
     sort(c("group", "time", "direction", "threshold", "value"))
@@ -133,14 +49,98 @@ test_that("brm_summary() on response", {
   expect_equal(x$threshold, rep(0, 3))
   expect_equal(
     x$value[1L],
-    mean(marginals$difference[["treatment 2, visit 2"]] > 0)
+    mean(draws$difference[["treatment 2, visit 2"]] > 0)
+  )
+  expect_equal(
+    x$value[2L],
+    mean(draws$difference[["treatment 2, visit 3"]] > 0)
+  )
+  expect_equal(
+    x$value[3L],
+    mean(draws$difference[["treatment 2, visit 4"]] > 0)
+  )
+})
+
+test_that("brm_marginal_probabilities() on change and multiple probs", {
+  set.seed(0L)
+  sim <- brm_simulate()
+  data <- sim$data
+  data$group <- paste("treatment", data$group)
+  data$time <- paste("visit", data$time)
+  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
+        )
+      )
+    )
+  )
+  draws <- brm_marginal_draws(
+    model = model,
+    group = "group",
+    time = "time",
+    patient = "patient",
+    control = "treatment 1",
+    baseline = "visit 1",
+    outcome = "change"
+  )
+  x <- brm_marginal_probabilities(
+    draws,
+    direction = c("less", "greater"),
+    threshold = c(-1.55, -1.7)
+  )
+  expect_equal(
+    sort(colnames(x)),
+    sort(c("group", "time", "direction", "threshold", "value"))
+  )
+  expect_equal(x$group, rep("treatment 2", 8))
+  expect_equal(x$time, rep(paste("visit", seq(1, 4)), times = 2))
+  expect_equal(x$direction, rep(c("greater", "less"), each = 4))
+  expect_equal(x$threshold, c(rep(-1.7, 4), rep(-1.55, 4)))
+  expect_equal(
+    x$value[1L],
+    mean(draws$difference[["treatment 2, visit 1"]] > -1.7)
   )
   expect_equal(
     x$value[2L],
-    mean(marginals$difference[["treatment 2, visit 3"]] > 0)
+    mean(draws$difference[["treatment 2, visit 2"]] > -1.7)
   )
   expect_equal(
     x$value[3L],
-    mean(marginals$difference[["treatment 2, visit 4"]] > 0)
+    mean(draws$difference[["treatment 2, visit 3"]] > -1.7)
+  )
+  expect_equal(
+    x$value[4L],
+    mean(draws$difference[["treatment 2, visit 4"]] > -1.7)
+  )
+  
+  expect_equal(
+    x$value[5L],
+    mean(draws$difference[["treatment 2, visit 1"]] < -1.55)
+  )
+  expect_equal(
+    x$value[6L],
+    mean(draws$difference[["treatment 2, visit 2"]] < -1.55)
+  )
+  expect_equal(
+    x$value[7L],
+    mean(draws$difference[["treatment 2, visit 3"]] < -1.55)
+  )
+  expect_equal(
+    x$value[8L],
+    mean(draws$difference[["treatment 2, visit 4"]] < -1.55)
   )
 })
diff --git a/tests/testthat/test-brm_marginal_summaries.R b/tests/testthat/test-brm_marginal_summaries.R
index e52995b7..d31f243d 100644
--- a/tests/testthat/test-brm_marginal_summaries.R
+++ b/tests/testthat/test-brm_marginal_summaries.R
@@ -158,7 +158,7 @@ test_that("brm_marginal_summaries() on change", {
   suppressWarnings(
     x <- brm_marginal_summaries(
       draws,
-      level = 0.95
+      level = 0.9
     )
   )
   expect_equal(
@@ -194,11 +194,11 @@ test_that("brm_marginal_summaries() on change", {
         )
         expect_equal(
           unname(subset$value[subset$statistic == "lower"]),
-          unname(quantile(draws$response[[name]], probs = 0.025))
+          unname(quantile(draws$response[[name]], probs = 0.05))
         )
         expect_equal(
           unname(subset$value[subset$statistic == "upper"]),
-          unname(quantile(draws$response[[name]], probs = 0.975))
+          unname(quantile(draws$response[[name]], probs = 0.95))
         )
         suppressWarnings({
           expect_equal(
@@ -218,7 +218,7 @@ test_that("brm_marginal_summaries() on change", {
             unname(
               posterior::mcse_quantile(
                 draws$response[[name]],
-                probs = 0.025
+                probs = 0.05
               )
             )
           )
@@ -227,7 +227,7 @@ test_that("brm_marginal_summaries() on change", {
             unname(
               posterior::mcse_quantile(
                 draws$response[[name]],
-                probs = 0.975
+                probs = 0.95
               )
             )
           )

From a3c0f6628dd5fd0b84de5d80e5ac75396af3fa5e Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 13:16:42 -0400
Subject: [PATCH 13/17] pkgdown

---
 _pkgdown.yml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/_pkgdown.yml b/_pkgdown.yml
index 7c452d8c..7050abfc 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -14,8 +14,8 @@ reference:
   contents:
   - '`brm_formula`'
   - '`brm_model`'
-- title: Results
+- title: Marginals
   contents:
-  - '`brm_convergence`'
-  - '`brm_marginals`'
-  - '`brm_summary`'
+  - '`brm_marginal_draws`'
+  - '`brm_marginal_summaries`'
+  - '`brm_marginal_probabilities`'

From f677ff5a83a80abd084ea1f1745ce0912a6ea46e Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 13:22:07 -0400
Subject: [PATCH 14/17] lints and spelling

---
 R/brm_marginal_draws.R                        |  2 +-
 R/brm_marginal_probabilities.R                |  2 +-
 R/brm_marginal_summaries.R                    | 20 ++++++-------------
 inst/WORDLIST                                 |  2 ++
 man/brm_marginal_probabilities.Rd             |  2 +-
 man/brm_marginal_summaries.Rd                 |  8 ++++----
 .../test-brm_marginal_probabilities.R         |  1 -
 7 files changed, 15 insertions(+), 22 deletions(-)

diff --git a/R/brm_marginal_draws.R b/R/brm_marginal_draws.R
index d731636f..52d4091c 100644
--- a/R/brm_marginal_draws.R
+++ b/R/brm_marginal_draws.R
@@ -206,7 +206,7 @@ subtract_control <- function(draws, groups, times, control) {
 }
 
 name_marginal <- function(group, time) {
-  sprintf("%s, %s", group , time)
+  sprintf("%s, %s", group, time)
 }
 
 names_group <- function(draws) {
diff --git a/R/brm_marginal_probabilities.R b/R/brm_marginal_probabilities.R
index 70611cbb..bf565b1c 100644
--- a/R/brm_marginal_probabilities.R
+++ b/R/brm_marginal_probabilities.R
@@ -2,7 +2,7 @@
 #' @export
 #' @family marginals
 #' @description Marginal probabilities on the treatment effect for an MMRM.
-#' @return A data frame of probabilities of the form
+#' @return A tibble of probabilities of the form
 #'   `Prob(treatment effect > threshold | data)` and/or
 #'   `Prob(treatment effect < threshold | data)`. It has one row per
 #'   probability and the following columns:
diff --git a/R/brm_marginal_summaries.R b/R/brm_marginal_summaries.R
index 201fe825..20eb2a12 100644
--- a/R/brm_marginal_summaries.R
+++ b/R/brm_marginal_summaries.R
@@ -2,17 +2,17 @@
 #' @export
 #' @family marginals
 #' @description Summary statistics of the marginal posterior of an MMRM.
-#' @return A data frame with one row per summary statistic and the following
+#' @return A tibble with one row per summary statistic and the following
 #'   columns:
-#'    * `group`: treatment group.
-#'    * `time`: discrete time point.
 #'    * `marginal`: type of marginal distribution. If `outcome` was `"response"`
 #'      in [brm_marginal_draws()], then possible values include
 #'      `"response"` for the response on the raw scale, `"change"` for
 #'      change from baseline, and `"difference"` for treatment difference
 #'      in terms of change from baseline. If `outcome` was `"change"`,
-#'      then possible values include `"response"` for the respons one the
+#'      then possible values include `"response"` for the response one the
 #'      change from baseline scale and `"difference"` for treatment difference.
+#'    * `group`: treatment group.
+#'    * `time`: discrete time point.
 #'    * `statistic`: type of summary statistic.
 #'    * `value`: numeric value of the estimate.
 #'    * `mcse`: Monte Carlo standard error of the estimate.
@@ -74,26 +74,18 @@ brm_marginal_summaries <- function(
     NULL
   )
   table_difference <- summarize_marginals(draws$difference, level)
-  out <- dplyr::bind_rows(
+  dplyr::bind_rows(
     response = table_response,
     change = table_change,
     difference = table_difference,
     .id = "marginal"
   )
-  columns <- c("marginal", "statistic", "group", "time", "value", "mcse")
-  out <- out[, columns]
-  args <- lapply(columns, as.symbol)
-  args$.data <- out
-  do.call(what = dplyr::arrange, args = args)
 }
 
 summarize_marginals <- function(draws, level) {
-  draws <- tibble::as_tibble(draws)
-  for (name in names_mcmc) {
-    draws[[name]] <- NULL
-  }
   level_lower <- (1 - level) / 2
   level_upper <- 1 - level_lower
+  draws[names_mcmc] <- NULL
   value <- tibble::tibble(
     group = names_group(draws),
     time = names_time(draws),
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 50da38b6..1f44a981 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -22,3 +22,5 @@ Schnell
 TBD
 Ther
 LKJ
+AVAL
+tibble
diff --git a/man/brm_marginal_probabilities.Rd b/man/brm_marginal_probabilities.Rd
index 00937662..b86de907 100644
--- a/man/brm_marginal_probabilities.Rd
+++ b/man/brm_marginal_probabilities.Rd
@@ -24,7 +24,7 @@ Each element \code{direction[i]} corresponds to \code{threshold[i]} for
 all \code{i} from 1 to \code{length(direction)}.}
 }
 \value{
-A data frame of probabilities of the form
+A tibble of probabilities of the form
 \verb{Prob(treatment effect > threshold | data)} and/or
 \verb{Prob(treatment effect < threshold | data)}. It has one row per
 probability and the following columns:
diff --git a/man/brm_marginal_summaries.Rd b/man/brm_marginal_summaries.Rd
index ee9fa81d..ac70a368 100644
--- a/man/brm_marginal_summaries.Rd
+++ b/man/brm_marginal_summaries.Rd
@@ -14,18 +14,18 @@ obtained from \code{\link[=brm_marginal_draws]{brm_marginal_draws()}}.}
 for the credible intervals.}
 }
 \value{
-A data frame with one row per summary statistic and the following
+A tibble with one row per summary statistic and the following
 columns:
 \itemize{
-\item \code{group}: treatment group.
-\item \code{time}: discrete time point.
 \item \code{marginal}: type of marginal distribution. If \code{outcome} was \code{"response"}
 in \code{\link[=brm_marginal_draws]{brm_marginal_draws()}}, then possible values include
 \code{"response"} for the response on the raw scale, \code{"change"} for
 change from baseline, and \code{"difference"} for treatment difference
 in terms of change from baseline. If \code{outcome} was \code{"change"},
-then possible values include \code{"response"} for the respons one the
+then possible values include \code{"response"} for the response one the
 change from baseline scale and \code{"difference"} for treatment difference.
+\item \code{group}: treatment group.
+\item \code{time}: discrete time point.
 \item \code{statistic}: type of summary statistic.
 \item \code{value}: numeric value of the estimate.
 \item \code{mcse}: Monte Carlo standard error of the estimate.
diff --git a/tests/testthat/test-brm_marginal_probabilities.R b/tests/testthat/test-brm_marginal_probabilities.R
index 80ba8584..d0f2513b 100644
--- a/tests/testthat/test-brm_marginal_probabilities.R
+++ b/tests/testthat/test-brm_marginal_probabilities.R
@@ -126,7 +126,6 @@ test_that("brm_marginal_probabilities() on change and multiple probs", {
     x$value[4L],
     mean(draws$difference[["treatment 2, visit 4"]] > -1.7)
   )
-  
   expect_equal(
     x$value[5L],
     mean(draws$difference[["treatment 2, visit 1"]] < -1.55)

From ffc5e1defeef7903ccbaaa40b09d33802f450a62 Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 13:24:11 -0400
Subject: [PATCH 15/17] fix a test

---
 R/brm_marginal_summaries.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/brm_marginal_summaries.R b/R/brm_marginal_summaries.R
index 20eb2a12..b8a87126 100644
--- a/R/brm_marginal_summaries.R
+++ b/R/brm_marginal_summaries.R
@@ -69,7 +69,7 @@ brm_marginal_summaries <- function(
   assert(level, . >= 0, . <= 1, message = "level arg must be between 0 and 1")
   table_response <- summarize_marginals(draws$response, level)
   table_change <- if_any(
-    "change" %in% names(marginals),
+    "change" %in% names(draws),
     summarize_marginals(draws$change, level),
     NULL
   )

From 69cb5d5e37916cd1858b3ca826a9026c49238130 Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 13:26:20 -0400
Subject: [PATCH 16/17] testing

---
 NAMESPACE                     | 1 -
 R/brm_marginal_summaries.R    | 2 +-
 R/brm_package.R               | 2 +-
 man/brm_marginal_summaries.Rd | 2 +-
 4 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index bd32fc12..efda81f4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,7 +13,6 @@ importFrom(brms,get_prior)
 importFrom(brms,prior)
 importFrom(brms,unstr)
 importFrom(dplyr,left_join)
-importFrom(emmeans,contrast)
 importFrom(emmeans,emmeans)
 importFrom(rlang,warn)
 importFrom(stats,as.formula)
diff --git a/R/brm_marginal_summaries.R b/R/brm_marginal_summaries.R
index b8a87126..81768ec9 100644
--- a/R/brm_marginal_summaries.R
+++ b/R/brm_marginal_summaries.R
@@ -56,7 +56,7 @@
 #'   baseline = "visit 1",
 #'   outcome = "response"
 #' )
-#' brm_marginal_summaries(draws)
+#' suppressWarnings(brm_marginal_summaries(draws))
 brm_marginal_summaries <- function(
   draws,
   level = 0.95
diff --git a/R/brm_package.R b/R/brm_package.R
index 271a7587..a13420e1 100644
--- a/R/brm_package.R
+++ b/R/brm_package.R
@@ -21,7 +21,7 @@
 #' @family help
 #' @importFrom brms brm brmsformula get_prior prior unstr
 #' @importFrom dplyr left_join
-#' @importFrom emmeans contrast emmeans
+#' @importFrom emmeans emmeans
 #' @importFrom MASS mvrnorm
 #' @importFrom rlang warn
 #' @importFrom stats as.formula model.matrix rnorm runif
diff --git a/man/brm_marginal_summaries.Rd b/man/brm_marginal_summaries.Rd
index ac70a368..7f7e0939 100644
--- a/man/brm_marginal_summaries.Rd
+++ b/man/brm_marginal_summaries.Rd
@@ -70,7 +70,7 @@ draws <- brm_marginal_draws(
   baseline = "visit 1",
   outcome = "response"
 )
-brm_marginal_summaries(draws)
+suppressWarnings(brm_marginal_summaries(draws))
 }
 \seealso{
 Other marginals: 

From 47927f5498c5f0c3ceb7141e89813a4aa332f023 Mon Sep 17 00:00:00 2001
From: wlandau-lilly <will.landau@lilly.com>
Date: Mon, 5 Jun 2023 13:36:25 -0400
Subject: [PATCH 17/17] explicit args

---
 DESCRIPTION                                      |  4 ++++
 NAMESPACE                                        | 12 ++++++++++++
 R/brm_package.R                                  |  7 ++++++-
 tests/testthat/test-brm_marginal_draws.R         | 12 ++++++++++--
 tests/testthat/test-brm_marginal_probabilities.R | 12 ++++++++++--
 tests/testthat/test-brm_marginal_summaries.R     | 12 ++++++++++--
 tests/testthat/test-brm_model.R                  |  6 +++++-
 7 files changed, 57 insertions(+), 8 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 75d75a89..b293c351 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -53,13 +53,17 @@ Depends:
   R (>= 4.0.0)
 Imports:
   brms,
+  coda,
   dplyr,
   emmeans,
   MASS,
+  posterior,
+  purrr,
   rlang,
   stats,
   tibble,
   tidyr,
+  tidyselect,
   trialr,
   utils
 Suggests:
diff --git a/NAMESPACE b/NAMESPACE
index efda81f4..6faecb9a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,15 +12,27 @@ importFrom(brms,brmsformula)
 importFrom(brms,get_prior)
 importFrom(brms,prior)
 importFrom(brms,unstr)
+importFrom(coda,as.mcmc)
 importFrom(dplyr,left_join)
 importFrom(emmeans,emmeans)
+importFrom(posterior,as_draws_df)
+importFrom(posterior,mcse_mean)
+importFrom(posterior,mcse_median)
+importFrom(posterior,mcse_quantile)
+importFrom(posterior,mcse_sd)
+importFrom(purrr,map2_df)
+importFrom(purrr,map_dbl)
+importFrom(purrr,map_df)
 importFrom(rlang,warn)
 importFrom(stats,as.formula)
+importFrom(stats,median)
 importFrom(stats,model.matrix)
 importFrom(stats,rnorm)
 importFrom(stats,runif)
+importFrom(stats,sd)
 importFrom(tibble,tibble)
 importFrom(tidyr,expand_grid)
+importFrom(tidyselect,any_of)
 importFrom(trialr,rlkjcorr)
 importFrom(utils,capture.output)
 importFrom(utils,globalVariables)
diff --git a/R/brm_package.R b/R/brm_package.R
index a13420e1..7af8a195 100644
--- a/R/brm_package.R
+++ b/R/brm_package.R
@@ -20,13 +20,18 @@
 #'     doi:10.1177/009286150804200402
 #' @family help
 #' @importFrom brms brm brmsformula get_prior prior unstr
+#' @importFrom coda as.mcmc
 #' @importFrom dplyr left_join
 #' @importFrom emmeans emmeans
 #' @importFrom MASS mvrnorm
+#' @importFrom posterior as_draws_df mcse_mean mcse_median mcse_quantile
+#'   mcse_sd
+#' @importFrom purrr map_dbl map_df map2_df
 #' @importFrom rlang warn
-#' @importFrom stats as.formula model.matrix rnorm runif
+#' @importFrom stats as.formula median model.matrix rnorm runif sd
 #' @importFrom tibble tibble
 #' @importFrom tidyr expand_grid
+#' @importFrom tidyselect any_of
 #' @importFrom trialr rlkjcorr
 #' @importFrom utils capture.output globalVariables
 NULL
diff --git a/tests/testthat/test-brm_marginal_draws.R b/tests/testthat/test-brm_marginal_draws.R
index 92927842..81441626 100644
--- a/tests/testthat/test-brm_marginal_draws.R
+++ b/tests/testthat/test-brm_marginal_draws.R
@@ -1,6 +1,10 @@
 test_that("brm_marginal_draws() on response", {
   set.seed(0L)
-  sim <- brm_simulate()
+  sim <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L
+  )
   data <- sim$data
   data$group <- paste("treatment", data$group)
   data$time <- paste("visit", data$time)
@@ -89,7 +93,11 @@ test_that("brm_marginal_draws() on response", {
 
 test_that("brm_marginal_draws() on change", {
   set.seed(0L)
-  sim <- brm_simulate()
+  sim <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L
+  )
   data <- sim$data
   data$group <- paste("treatment", data$group)
   data$time <- paste("visit", data$time)
diff --git a/tests/testthat/test-brm_marginal_probabilities.R b/tests/testthat/test-brm_marginal_probabilities.R
index d0f2513b..bdc71ac9 100644
--- a/tests/testthat/test-brm_marginal_probabilities.R
+++ b/tests/testthat/test-brm_marginal_probabilities.R
@@ -1,6 +1,10 @@
 test_that("brm_marginal_probabilities() on response", {
   set.seed(0L)
-  sim <- brm_simulate()
+  sim <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L
+  )
   data <- sim$data
   data$group <- paste("treatment", data$group)
   data$time <- paste("visit", data$time)
@@ -63,7 +67,11 @@ test_that("brm_marginal_probabilities() on response", {
 
 test_that("brm_marginal_probabilities() on change and multiple probs", {
   set.seed(0L)
-  sim <- brm_simulate()
+  sim <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L
+  )
   data <- sim$data
   data$group <- paste("treatment", data$group)
   data$time <- paste("visit", data$time)
diff --git a/tests/testthat/test-brm_marginal_summaries.R b/tests/testthat/test-brm_marginal_summaries.R
index d31f243d..7b475a93 100644
--- a/tests/testthat/test-brm_marginal_summaries.R
+++ b/tests/testthat/test-brm_marginal_summaries.R
@@ -1,6 +1,10 @@
 test_that("brm_marginal_summaries() on response", {
   set.seed(0L)
-  sim <- brm_simulate()
+  sim <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L
+  )
   data <- sim$data
   data$group <- paste("treatment", data$group)
   data$time <- paste("visit", data$time)
@@ -121,7 +125,11 @@ test_that("brm_marginal_summaries() on response", {
 
 test_that("brm_marginal_summaries() on change", {
   set.seed(0L)
-  sim <- brm_simulate()
+  sim <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L
+  )
   data <- sim$data
   data$group <- paste("treatment", data$group)
   data$time <- paste("visit", data$time)
diff --git a/tests/testthat/test-brm_model.R b/tests/testthat/test-brm_model.R
index b5e1ff2f..3ec7b84f 100644
--- a/tests/testthat/test-brm_model.R
+++ b/tests/testthat/test-brm_model.R
@@ -1,6 +1,10 @@
 test_that("brm_model() runs", {
   set.seed(0L)
-  data <- brm_simulate()$data
+  data <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L
+  )$data
   formula <- brm_formula(
     response = "response",
     group = "group",