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] 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
+              )
+            )
+          )
+        })
+      }
+    }
+  }
+})