diff --git a/DESCRIPTION b/DESCRIPTION
index 1101fbab..b84e1397 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -53,6 +53,8 @@ Depends:
   R (>= 4.0.0)
 Imports:
   brms,
+  dplyr,
+  emmeans,
   MASS,
   rlang,
   stats,
@@ -65,6 +67,13 @@ Suggests:
   markdown (>= 1.1),
   rmarkdown (>= 2.4),
   testthat (>= 3.0.0)
+LinkingTo: 
+  BH,
+  Rcpp,
+  RcppEigen,
+  RcppParallel,
+  rstan,
+  StanHeaders
 Encoding: UTF-8
 Language: en-US
 VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
index 22b1991e..c166e393 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,12 +3,16 @@
 export(brm_formula)
 export(brm_model)
 export(brm_simulate)
+export(brm_summary)
 importFrom(MASS,mvrnorm)
 importFrom(brms,brm)
 importFrom(brms,brmsformula)
 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)
 importFrom(stats,model.matrix)
diff --git a/R/brm_formula.R b/R/brm_formula.R
index c1948ec7..232eff6c 100644
--- a/R/brm_formula.R
+++ b/R/brm_formula.R
@@ -1,6 +1,6 @@
 #' @title Model formula
 #' @export
-#' @family model
+#' @family models
 #' @description Build a model formula for an MMRM.
 #' @return An object of class `"brmsformula"` returned from
 #'   `brms::brmsformula()`. It contains the fixed effect parameterization,
@@ -58,7 +58,7 @@ brm_formula <- function(
   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_vec(covariates, "covariates arg must be a nonempty chr vector")
+  assert_chr_vec(covariates, "covariates arg must be a character vector")
   assert_lgl(intercept)
   assert_lgl(effect_group)
   assert_lgl(effect_time)
diff --git a/R/brm_model.R b/R/brm_model.R
index e5489faf..1dfae85f 100644
--- a/R/brm_model.R
+++ b/R/brm_model.R
@@ -1,6 +1,6 @@
 #' @title Basic MMRM
 #' @export
-#' @family model
+#' @family models
 #' @description Fit a basic MMRM model using `brms`.
 #' @return A fitted model object from `brms`.
 #' @param data A tidy data frame with one row per patient per discrete
diff --git a/R/brm_package.R b/R/brm_package.R
index f23e5e92..271a7587 100644
--- a/R/brm_package.R
+++ b/R/brm_package.R
@@ -20,6 +20,8 @@
 #'     doi:10.1177/009286150804200402
 #' @family help
 #' @importFrom brms brm brmsformula get_prior prior unstr
+#' @importFrom dplyr left_join
+#' @importFrom emmeans contrast emmeans
 #' @importFrom MASS mvrnorm
 #' @importFrom rlang warn
 #' @importFrom stats as.formula model.matrix rnorm runif
diff --git a/R/brm_summary.R b/R/brm_summary.R
new file mode 100644
index 00000000..2d1bb4d5
--- /dev/null
+++ b/R/brm_summary.R
@@ -0,0 +1,166 @@
+#' @title Summarize an MMRM.
+#' @export
+#' @family results
+#' @description Summarize a basic MMRM model fit.
+#' @details Currently assumes the response variable is `CHG`
+#'   (change from baseline) and not `AVAL` (raw response).
+#' @return A `tibble` with summary statistics of the marginal posterior.
+#' @inheritParams brm_formula
+#' @param model Fitted `brms` model object from [brm_model()].
+#' @param control Character of length 1, name of the control arm
+#'   in the `group` column in the data.
+#' @examples
+#' set.seed(0L)
+#' sim <- brm_simulate()
+#' data <- sim$data
+#' formula <- brm_formula(
+#'   response = "response",
+#'   group = "group",
+#'   time = "time",
+#'   patient = "patient",
+#'   effect_base = FALSE,
+#'   interaction_base = FALSE
+#' )
+#' tmp <- utils::capture.output(
+#'   suppressMessages(
+#'     suppressWarnings(
+#'       model <- brm_model(
+#'         data = data,
+#'         formula = formula,
+#'         chains = 1,
+#'         iter = 100,
+#'         refresh = 0
+#'       )
+#'     )
+#'   )
+#' )
+#' brm_summary(
+#'   model = model,
+#'   group = "group",
+#'   time = "time",
+#'   patient = "patient",
+#'   control = 1
+#' )
+brm_summary <- function(
+  model,
+  base = "BASE",
+  group = "TRT01P",
+  time = "AVISIT",
+  patient = "USUBJID",
+  covariates = character(0),
+  control = "Placebo"
+) {
+  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_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(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_response <- emmeans::emmeans(
+    object = model,
+    specs = as.formula(sprintf("~%s:%s", time, group)),
+    weights = "proportional",
+    nuisance = nuisance
+  )
+  table_response <- brm_summary_response(
+    data = data,
+    emmeans_response = emmeans_response
+  )
+  table_diff <- brm_summary_diff(
+    data = data,
+    emmeans_response = emmeans_response,
+    group = group,
+    time = time,
+    nuisance = nuisance,
+    control = control
+  )
+  dplyr::left_join(
+    x = table_response,
+    y = table_diff,
+    by = c(group, time)
+  )
+}
+
+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
+  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
+    }
+  }
+  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
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index e09bbfba..98d2c2af 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -17,5 +17,5 @@ reference:
 - title: Results
   contents:
   - '`brm_convergence`'
-  - '`brm_summary`'
   - '`brm_plot`'
+  - '`brm_summary`'
diff --git a/man/brm_formula.Rd b/man/brm_formula.Rd
index 646459f2..8917419e 100644
--- a/man/brm_formula.Rd
+++ b/man/brm_formula.Rd
@@ -77,7 +77,7 @@ brm_formula(
 )
 }
 \seealso{
-Other model: 
+Other models: 
 \code{\link{brm_model}()}
 }
-\concept{model}
+\concept{models}
diff --git a/man/brm_model.Rd b/man/brm_model.Rd
index 6e4e44d7..86b15267 100644
--- a/man/brm_model.Rd
+++ b/man/brm_model.Rd
@@ -75,7 +75,7 @@ model
 brms::prior_summary(model)
 }
 \seealso{
-Other model: 
+Other models: 
 \code{\link{brm_formula}()}
 }
-\concept{model}
+\concept{models}
diff --git a/man/brm_summary.Rd b/man/brm_summary.Rd
new file mode 100644
index 00000000..fa464928
--- /dev/null
+++ b/man/brm_summary.Rd
@@ -0,0 +1,81 @@
+% 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(
+  model,
+  base = "BASE",
+  group = "TRT01P",
+  time = "AVISIT",
+  patient = "USUBJID",
+  covariates = character(0),
+  control = "Placebo"
+)
+}
+\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{control}{Character of length 1, name of the control arm
+in the \code{group} column in the data.}
+}
+\value{
+A \code{tibble} with summary statistics of the marginal posterior.
+}
+\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()
+data <- sim$data
+formula <- brm_formula(
+  response = "response",
+  group = "group",
+  time = "time",
+  patient = "patient",
+  effect_base = FALSE,
+  interaction_base = FALSE
+)
+tmp <- utils::capture.output(
+  suppressMessages(
+    suppressWarnings(
+      model <- brm_model(
+        data = data,
+        formula = formula,
+        chains = 1,
+        iter = 100,
+        refresh = 0
+      )
+    )
+  )
+)
+brm_summary(
+  model = model,
+  group = "group",
+  time = "time",
+  patient = "patient",
+  control = 1
+)
+}
+\concept{results}
diff --git a/tests/testthat/test-brm_summary.R b/tests/testthat/test-brm_summary.R
new file mode 100644
index 00000000..ce7aa04c
--- /dev/null
+++ b/tests/testthat/test-brm_summary.R
@@ -0,0 +1,63 @@
+test_that("brm_summary()", {
+  set.seed(0L)
+  sim <- brm_simulate(
+    n_group = 2L,
+    n_patient = 100L,
+    n_time = 4L,
+    hyper_beta = 1,
+    hyper_sigma = 1,
+    hyper_correlation = 1
+  )
+  data <- sim$data
+  formula <- brm_formula(
+    response = "response",
+    group = "group",
+    time = "time",
+    patient = "patient",
+    effect_base = FALSE,
+    interaction_base = FALSE
+  )
+  tmp <- utils::capture.output(
+    suppressMessages(
+      suppressWarnings(
+        model <- brm_model(
+          data = data,
+          formula = formula,
+          chains = 1,
+          iter = 100,
+          refresh = 0
+        )
+      )
+    )
+  )
+  out <- brm_summary(
+    model = model,
+    group = "group",
+    time = "time",
+    patient = "patient",
+    control = 1
+  )
+  cols <- c(
+    "time",
+    "group",
+    "response_mean",
+    "response_lower",
+    "response_upper",
+    "diff_mean",
+    "diff_lower",
+    "diff_upper"
+  )
+  expect_equal(sort(colnames(out)), sort(cols))
+  expect_equal(nrow(out), 8L)
+  expect_equal(as.integer(out$time), rep(seq_len(4), 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))
+  expect_equal(
+    out$diff_mean[seq(5, 8)],
+    out$response_mean[out$group == 2] - out$response_mean[out$group == 1],
+    tolerance = 0.01
+  )
+})