From a2515e7acadff0df10f5b63a2475dde9c8e78ccf Mon Sep 17 00:00:00 2001 From: wlandau Date: Tue, 25 Apr 2023 16:03:09 -0400 Subject: [PATCH 1/6] Add brm_summary() --- NAMESPACE | 1 + R/brm_formula.R | 4 +- R/brm_model.R | 2 +- R/brm_summary.R | 160 +++++++++++++++++++++++++++++++++++++++++++++ _pkgdown.yml | 2 +- man/brm_formula.Rd | 4 +- man/brm_model.Rd | 4 +- man/brm_summary.Rd | 74 +++++++++++++++++++++ 8 files changed, 243 insertions(+), 8 deletions(-) create mode 100644 R/brm_summary.R create mode 100644 man/brm_summary.Rd diff --git a/NAMESPACE b/NAMESPACE index 22b1991e..b69d2afd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ export(brm_formula) export(brm_model) export(brm_simulate) +export(brm_summary) importFrom(MASS,mvrnorm) importFrom(brms,brm) importFrom(brms,brmsformula) 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_summary.R b/R/brm_summary.R new file mode 100644 index 00000000..29ce12cd --- /dev/null +++ b/R/brm_summary.R @@ -0,0 +1,160 @@ +#' @title Summarize an MMRM. +#' @export +#' @family results +#' @description Summarize a basic MMRM model fit. +#' @return A `tibble` with summary statistics of the marginal posterior. +#' @inheritParams brm_formula +#' @param model Fitted `brms` model object from [brm_model()] +#' @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_change( + 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_change <- function( + data, + emmeans_response, + group, + time, + nuisance, + control +) { + contrasts_diff <- list() + reference <- tibble::as_tibble(emmeans_response) + for (level_group in setdiff(sort(unique(reference[[group]])), control)) { + for (level_time in sort(unique(reference[[time]]))) { + contrast_treatment <- as.integer( + reference[[group]] == group & reference[[time]] == level_time + ) + 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..3f1d30f9 --- /dev/null +++ b/man/brm_summary.Rd @@ -0,0 +1,74 @@ +% 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.} +} +\value{ +A \code{tibble} with summary statistics of the marginal posterior. +} +\description{ +Summarize a basic MMRM model fit. +} +\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} From 3e9054f0a40fa6ad4c7a8ea2685d746e2dcd958f Mon Sep 17 00:00:00 2001 From: wlandau Date: Mon, 15 May 2023 13:36:30 -0400 Subject: [PATCH 2/6] Fix and test brm_summary() --- R/brm_summary.R | 14 ++++--- man/brm_summary.Rd | 4 ++ tests/testthat/test-brm_summary.R | 62 +++++++++++++++++++++++++++++++ 3 files changed, 75 insertions(+), 5 deletions(-) create mode 100644 tests/testthat/test-brm_summary.R diff --git a/R/brm_summary.R b/R/brm_summary.R index 29ce12cd..c97ba2d7 100644 --- a/R/brm_summary.R +++ b/R/brm_summary.R @@ -2,6 +2,8 @@ #' @export #' @family results #' @description Summarize a basic MMRM model fit. +#' @details Currently assumes the response variable is `CHG` +#' (change from baseline) and not `AVAL` (raw response). #' @return A `tibble` with summary statistics of the marginal posterior. #' @inheritParams brm_formula #' @param model Fitted `brms` model object from [brm_model()] @@ -91,7 +93,7 @@ brm_summary <- function( data = data, emmeans_response = emmeans_response ) - table_diff <- brm_summary_diff_change( + table_diff <- brm_summary_diff( data = data, emmeans_response = emmeans_response, group = group, @@ -117,7 +119,7 @@ brm_summary_response <- function(data, emmeans_response) { out } -brm_summary_diff_change <- function( +brm_summary_diff <- function( data, emmeans_response, group, @@ -125,15 +127,17 @@ brm_summary_diff_change <- function( nuisance, control ) { - contrasts_diff <- list() reference <- tibble::as_tibble(emmeans_response) + contrasts_diff <- list() for (level_group in setdiff(sort(unique(reference[[group]])), control)) { for (level_time in sort(unique(reference[[time]]))) { contrast_treatment <- as.integer( - reference[[group]] == group & reference[[time]] == level_time + (reference[[group]] == level_group) & + (reference[[time]] == level_time) ) contrast_control <- as.integer( - reference[[group]] == control & reference[[time]] == level_time + (reference[[group]] == control) & + (reference[[time]] == level_time) ) contrast <- contrast_treatment - contrast_control contrasts_diff[[length(contrasts_diff) + 1L]] <- contrast diff --git a/man/brm_summary.Rd b/man/brm_summary.Rd index 3f1d30f9..c895a1b3 100644 --- a/man/brm_summary.Rd +++ b/man/brm_summary.Rd @@ -38,6 +38,10 @@ 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() diff --git a/tests/testthat/test-brm_summary.R b/tests/testthat/test-brm_summary.R new file mode 100644 index 00000000..c412918a --- /dev/null +++ b/tests/testthat/test-brm_summary.R @@ -0,0 +1,62 @@ +test_that("brm_summary()", { + set.seed(0L) + sim <- brm_simulate( + n_group = 2L, + n_patient = 100L, + n_time = 4L, + hyper_beta = 1, + hyper_sigma = 1, + hyper_correlation = 1 + ) + data <- sim$data + formula <- brm_formula( + response = "response", + group = "group", + time = "time", + patient = "patient", + effect_base = FALSE, + interaction_base = FALSE + ) + tmp <- utils::capture.output( + suppressMessages( + suppressWarnings( + model <- brm_model( + data = data, + formula = formula, + chains = 1, + iter = 100, + refresh = 0 + ) + ) + ) + ) + out <- brm_summary( + model = model, + group = "group", + time = "time", + patient = "patient", + control = 1 + ) + cols <- c( + "time", + "group", + "response_mean", + "response_lower", + "response_upper", + "diff_mean", + "diff_lower", + "diff_upper" + ) + expect_equal(sort(colnames(out)), sort(cols)) + expect_equal(nrow(out), 8L) + expect_equal(as.integer(out$time), rep(seq_len(4), each = 2)) + expect_equal(as.integer(out$group), rep(seq_len(2), times = 4)) + expect_true(all(out$response_mean > out$response_lower)) + expect_true(all(out$response_mean > out$response_upper)) + expect_true(all(out$diff_mean > out$diff_lower)) + expect_true(all(out$diff_mean > out$diff_upper)) + expect_equal( + out$diff_mean[seq_len(4)], + out$response_mean[out$group == 2] - out$response_mean[out$group == 1] + ) +}) From 9de49ab26a1cbf31207c4dfbcc06490dde5453db Mon Sep 17 00:00:00 2001 From: wlandau Date: Mon, 15 May 2023 13:39:51 -0400 Subject: [PATCH 3/6] checks --- DESCRIPTION | 2 ++ NAMESPACE | 3 +++ R/brm_package.R | 2 ++ R/brm_summary.R | 2 +- 4 files changed, 8 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1101fbab..7693d923 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -53,6 +53,8 @@ Depends: R (>= 4.0.0) Imports: brms, + dplyr, + emmeans, MASS, rlang, stats, diff --git a/NAMESPACE b/NAMESPACE index b69d2afd..c166e393 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,6 +10,9 @@ 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_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 index c97ba2d7..c68483cf 100644 --- a/R/brm_summary.R +++ b/R/brm_summary.R @@ -136,7 +136,7 @@ brm_summary_diff <- function( (reference[[time]] == level_time) ) contrast_control <- as.integer( - (reference[[group]] == control) & + (reference[[group]] == control) & (reference[[time]] == level_time) ) contrast <- contrast_treatment - contrast_control From 3c299a70b9867380f875933c200b4232284ac695 Mon Sep 17 00:00:00 2001 From: wlandau Date: Mon, 15 May 2023 13:55:37 -0400 Subject: [PATCH 4/6] Fix tests --- R/brm_summary.R | 4 +++- man/brm_summary.Rd | 5 ++++- tests/testthat/test-brm_summary.R | 15 ++++++++------- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/R/brm_summary.R b/R/brm_summary.R index c68483cf..2d1bb4d5 100644 --- a/R/brm_summary.R +++ b/R/brm_summary.R @@ -6,7 +6,9 @@ #' (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 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() diff --git a/man/brm_summary.Rd b/man/brm_summary.Rd index c895a1b3..fa464928 100644 --- a/man/brm_summary.Rd +++ b/man/brm_summary.Rd @@ -15,7 +15,7 @@ brm_summary( ) } \arguments{ -\item{model}{Fitted \code{brms} model object from \code{\link[=brm_model]{brm_model()}}} +\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.} @@ -31,6 +31,9 @@ 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. diff --git a/tests/testthat/test-brm_summary.R b/tests/testthat/test-brm_summary.R index c412918a..ce7aa04c 100644 --- a/tests/testthat/test-brm_summary.R +++ b/tests/testthat/test-brm_summary.R @@ -49,14 +49,15 @@ test_that("brm_summary()", { ) expect_equal(sort(colnames(out)), sort(cols)) expect_equal(nrow(out), 8L) - expect_equal(as.integer(out$time), rep(seq_len(4), each = 2)) - expect_equal(as.integer(out$group), rep(seq_len(2), times = 4)) + expect_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)) - expect_true(all(out$diff_mean > out$diff_upper)) + 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_len(4)], - out$response_mean[out$group == 2] - out$response_mean[out$group == 1] + out$diff_mean[seq(5, 8)], + out$response_mean[out$group == 2] - out$response_mean[out$group == 1], + tolerance = 0.01 ) }) From badf4a556f54e565bd33fb41ddac87775bc89754 Mon Sep 17 00:00:00 2001 From: wlandau Date: Mon, 15 May 2023 14:03:38 -0400 Subject: [PATCH 5/6] BH --- DESCRIPTION | 1 + 1 file changed, 1 insertion(+) diff --git a/DESCRIPTION b/DESCRIPTION index 7693d923..7123cddf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -63,6 +63,7 @@ Imports: trialr, utils Suggests: + BH, knitr (>= 1.30), markdown (>= 1.1), rmarkdown (>= 2.4), From c34a65a26c0e96d203dff0b14e5cb395b25265b0 Mon Sep 17 00:00:00 2001 From: wlandau Date: Mon, 15 May 2023 14:09:38 -0400 Subject: [PATCH 6/6] pkgs --- DESCRIPTION | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 7123cddf..b84e1397 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -63,11 +63,17 @@ Imports: trialr, utils Suggests: - BH, knitr (>= 1.30), 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