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