A mixed model of repeated measures (MMRM) analyzes longitudinal clinical trial data. In a longitudinal dataset, there are multiple patients, and each patient has multiple observations at a common set of discrete points in time.
Data
To use the brms.mmrm package, begin with a longitudinal
dataset with one row per patient observation and columns for the
response variable, treatment group indicator, discrete time point
indicator, patient ID variable, and optional baseline covariates such as
age and region. If you do not have a real dataset of your own, you can
simulate one from the package. The following dataset has the raw
response variable, the essential factor variables, and continuous
baseline covariates.1 In general, the outcome variable can either
be the raw response or change from baseline.
library(brms.mmrm)
library(dplyr)
library(magrittr)
set.seed(0L)
raw_data <- brm_simulate_simple(
n_group = 3,
n_patient = 100,
n_time = 4
) %>%
extract2("data") %>%
brm_simulate_continuous(c("biomarker1", "biomarker2", "biomarker3"))
raw_data
#> # A tibble: 1,200 × 7
#> response group time patient biomarker1 biomarker2 biomarker3
#> <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 1.11 group_1 time_1 patient_001 1.31 -0.361 1.52
#> 2 2.15 group_1 time_2 patient_001 1.31 -0.361 1.52
#> 3 2.54 group_1 time_3 patient_001 1.31 -0.361 1.52
#> 4 -1.73 group_1 time_4 patient_001 1.31 -0.361 1.52
#> 5 1.11 group_1 time_1 patient_002 0.107 -2.44 -0.139
#> 6 2.64 group_1 time_2 patient_002 0.107 -2.44 -0.139
#> 7 1.69 group_1 time_3 patient_002 0.107 -2.44 -0.139
#> 8 0.783 group_1 time_4 patient_002 0.107 -2.44 -0.139
#> 9 0.118 group_1 time_1 patient_003 1.44 -0.419 -1.54
#> 10 2.48 group_1 time_2 patient_003 1.44 -0.419 -1.54
#> # ℹ 1,190 more rowsNext, create a special classed dataset that the package will recognize. The classed data object contains a pre-processed version of the data, along with attributes to declare the outcome variable, whether the outcome is response or change from baseline, the treatment group variable, the discrete time point variable, control group, baseline time point, and the covariates selected for analysis.
data <- brm_data(
data = raw_data,
outcome = "response",
role = "response",
group = "group",
patient = "patient",
time = "time",
covariates = c("biomarker1", "biomarker2"),
reference_group = "group_1",
reference_time = "time_1"
)
data
#> # A tibble: 1,200 × 6
#> response group time patient biomarker1 biomarker2
#> <dbl> <chr> <chr> <chr> <dbl> <dbl>
#> 1 1.11 group_1 time_1 patient_001 1.31 -0.361
#> 2 2.15 group_1 time_2 patient_001 1.31 -0.361
#> 3 2.54 group_1 time_3 patient_001 1.31 -0.361
#> 4 -1.73 group_1 time_4 patient_001 1.31 -0.361
#> 5 1.11 group_1 time_1 patient_002 0.107 -2.44
#> 6 2.64 group_1 time_2 patient_002 0.107 -2.44
#> 7 1.69 group_1 time_3 patient_002 0.107 -2.44
#> 8 0.783 group_1 time_4 patient_002 0.107 -2.44
#> 9 0.118 group_1 time_1 patient_003 1.44 -0.419
#> 10 2.48 group_1 time_2 patient_003 1.44 -0.419
#> # ℹ 1,190 more rows
class(data)
#> [1] "brms_mmrm_data" "tbl_df" "tbl" "data.frame"
roles <- attributes(data)
roles$row.names <- NULL
str(roles)
#> List of 14
#> $ names : chr [1:6] "response" "group" "time" "patient" ...
#> $ class : chr [1:4] "brms_mmrm_data" "tbl_df" "tbl" "data.frame"
#> $ brm_outcome : chr "response"
#> $ brm_role : chr "response"
#> $ brm_group : chr "group"
#> $ brm_time : chr "time"
#> $ brm_patient : chr "patient"
#> $ brm_covariates : chr [1:2] "biomarker1" "biomarker2"
#> $ brm_reference_group: chr "group_1"
#> $ brm_reference_time : chr "time_1"
#> $ brm_levels_group : chr [1:3] "group_1" "group_2" "group_3"
#> $ brm_levels_time : chr [1:4] "time_1" "time_2" "time_3" "time_4"
#> $ brm_labels_group : chr [1:3] "group_1" "group_2" "group_3"
#> $ brm_labels_time : chr [1:4] "time_1" "time_2" "time_3" "time_4"Formula
Next, choose a brms model formula for the fixed effect
and variance parameters. The brm_formula() function from
brms.mmrm makes this process easier. A cell means
parameterization for this particular model can be expressed as follows.
It specifies one fixed effect parameter for each combination of
treatment group and time point, and it makes the specification of
informative priors straightforward through the prior
argument of brm_model().
brm_formula(
data = data,
intercept = FALSE,
baseline = FALSE,
group = FALSE,
time = FALSE,
baseline_time = FALSE,
group_time = TRUE
)
#> response ~ 0 + group:time + biomarker1 + biomarker2 + unstr(time = time, gr = patient)
#> sigma ~ 0 + timeFor the purposes of our example, we choose a fully parameterized analysis of the raw response.
formula <- brm_formula(
data = data,
intercept = TRUE,
baseline = FALSE,
group = TRUE,
time = TRUE,
baseline_time = FALSE,
group_time = TRUE
)
formula
#> response ~ group + group:time + time + biomarker1 + biomarker2 + unstr(time = time, gr = patient)
#> sigma ~ 0 + timeParameterization
The formula is not the only factor that ultimately determines the
fixed effect parameterization. The ordering of the categorical variables
in the data, as well as the contrast option in R, affect
the construction of the model matrix. To see the model matrix that will
ultimately be used in brm_model(), run
brms::make_standata() and examine the X
element of the returned list.
The contrast option accepts a named vector of two
character vectors which govern model.matrix() contrasts for
unordered and ordered variables, respectively.
The make_standata() function lets you see the data that
brms will generate for Stan. This includes the fixed
effects model matrix X. Note the differences in the
groupgroup_* additive terms between the matrix below and
the one above.
head(brms::make_standata(formula = formula, data = data)$X)
#> Intercept groupgroup_1 groupgroup_2 timetime_1 timetime_2 timetime_3
#> 1 1 1 0 1 0 0
#> 2 1 1 0 0 1 0
#> 3 1 1 0 0 0 1
#> 4 1 1 0 0 0 0
#> 5 1 1 0 1 0 0
#> 6 1 1 0 0 1 0
#> biomarker1 biomarker2 groupgroup_1:timetime_1 groupgroup_2:timetime_1
#> 1 1.3126508 -0.3608809 1 0
#> 2 1.3126508 -0.3608809 0 0
#> 3 1.3126508 -0.3608809 0 0
#> 4 1.3126508 -0.3608809 0 0
#> 5 0.1068624 -2.4441488 1 0
#> 6 0.1068624 -2.4441488 0 0
#> groupgroup_1:timetime_2 groupgroup_2:timetime_2 groupgroup_1:timetime_3
#> 1 0 0 0
#> 2 1 0 0
#> 3 0 0 1
#> 4 0 0 0
#> 5 0 0 0
#> 6 1 0 0
#> groupgroup_2:timetime_3
#> 1 0
#> 2 0
#> 3 0
#> 4 0
#> 5 0
#> 6 0If you choose a different contrast method, a different model matrix may result.
options(
contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")
)
# different model matrix than before:
head(brms::make_standata(formula = formula, data = data)$X)
#> Intercept groupgroup_2 groupgroup_3 timetime_2 timetime_3 timetime_4
#> 1 1 0 0 0 0 0
#> 2 1 0 0 1 0 0
#> 3 1 0 0 0 1 0
#> 4 1 0 0 0 0 1
#> 5 1 0 0 0 0 0
#> 6 1 0 0 1 0 0
#> biomarker1 biomarker2 groupgroup_2:timetime_2 groupgroup_3:timetime_2
#> 1 1.3126508 -0.3608809 0 0
#> 2 1.3126508 -0.3608809 0 0
#> 3 1.3126508 -0.3608809 0 0
#> 4 1.3126508 -0.3608809 0 0
#> 5 0.1068624 -2.4441488 0 0
#> 6 0.1068624 -2.4441488 0 0
#> groupgroup_2:timetime_3 groupgroup_3:timetime_3 groupgroup_2:timetime_4
#> 1 0 0 0
#> 2 0 0 0
#> 3 0 0 0
#> 4 0 0 0
#> 5 0 0 0
#> 6 0 0 0
#> groupgroup_3:timetime_4
#> 1 0
#> 2 0
#> 3 0
#> 4 0
#> 5 0
#> 6 0Priors
Some analyses require informative priors, others require
non-informative ones. Please use brms to
construct a prior suitable for your analysis. The brms
package has documentation on how its default priors are constructed and
how to set your own priors. Once you have an R object that represents
the joint prior distribution of your model, you can pass it to the
brm_model() function described below. The
get_prior() function shows the default priors for a given
dataset and model formula.
brms::get_prior(data = data, formula = formula)
#> prior class coef group resp dpar
#> (flat) b
#> (flat) b biomarker1
#> (flat) b biomarker2
#> (flat) b groupgroup_2
#> (flat) b groupgroup_2:timetime_2
#> (flat) b groupgroup_2:timetime_3
#> (flat) b groupgroup_2:timetime_4
#> (flat) b groupgroup_3
#> (flat) b groupgroup_3:timetime_2
#> (flat) b groupgroup_3:timetime_3
#> (flat) b groupgroup_3:timetime_4
#> (flat) b timetime_2
#> (flat) b timetime_3
#> (flat) b timetime_4
#> lkj(1) cortime
#> student_t(3, 0.9, 2.5) Intercept
#> (flat) b sigma
#> (flat) b timetime_1 sigma
#> (flat) b timetime_2 sigma
#> (flat) b timetime_3 sigma
#> (flat) b timetime_4 sigma
#> nlpar lb ub source
#> default
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> default
#> default
#> default
#> (vectorized)
#> (vectorized)
#> (vectorized)
#> (vectorized)Model
To run an MMRM, use the brm_model() function. This
function calls brms::brm() behind the scenes, using the
formula and prior you set in the formula and
prior arguments.
model <- brm_model(data = data, formula = formula, refresh = 0)The result is a brms model object.
model
#> Family: gaussian
#> Links: mu = identity; sigma = log
#> Formula: response ~ group + group:time + time + biomarker1 + biomarker2 + unstr(time = time, gr = patient)
#> sigma ~ 0 + time
#> Data: data[!is.na(data[[attr(data, "brm_outcome")]]), ] (Number of observations: 1200)
#> Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup draws = 4000
#>
#> Correlation Structures:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> cortime(time_1,time_2) 0.42 0.05 0.33 0.51 1.00 3325
#> cortime(time_1,time_3) -0.80 0.02 -0.84 -0.76 1.00 2161
#> cortime(time_2,time_3) -0.54 0.04 -0.62 -0.46 1.00 3218
#> cortime(time_1,time_4) -0.29 0.05 -0.39 -0.18 1.00 4082
#> cortime(time_2,time_4) 0.03 0.06 -0.09 0.14 1.00 3227
#> cortime(time_3,time_4) -0.10 0.06 -0.21 0.01 1.00 3697
#> Tail_ESS
#> cortime(time_1,time_2) 2911
#> cortime(time_1,time_3) 2816
#> cortime(time_2,time_3) 3092
#> cortime(time_1,time_4) 2939
#> cortime(time_2,time_4) 2861
#> cortime(time_3,time_4) 2853
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 1.21 0.09 1.04 1.39 1.00 1319
#> groupgroup_2 -1.53 0.12 -1.77 -1.29 1.00 1394
#> groupgroup_3 0.19 0.13 -0.06 0.44 1.00 1386
#> timetime_2 1.32 0.10 1.11 1.51 1.00 1994
#> timetime_3 0.52 0.17 0.17 0.87 1.00 1479
#> timetime_4 -1.44 0.17 -1.79 -1.11 1.00 1484
#> biomarker1 -0.01 0.01 -0.03 0.01 1.00 6061
#> biomarker2 0.02 0.01 0.00 0.04 1.00 4965
#> groupgroup_2:timetime_2 0.04 0.14 -0.24 0.32 1.00 2310
#> groupgroup_3:timetime_2 0.02 0.14 -0.26 0.31 1.00 2186
#> groupgroup_2:timetime_3 -0.17 0.25 -0.65 0.31 1.00 1589
#> groupgroup_3:timetime_3 -0.21 0.25 -0.69 0.28 1.00 1596
#> groupgroup_2:timetime_4 -0.08 0.25 -0.55 0.42 1.00 1544
#> groupgroup_3:timetime_4 -0.27 0.25 -0.75 0.21 1.00 1609
#> sigma_timetime_1 -0.12 0.04 -0.20 -0.04 1.00 2469
#> sigma_timetime_2 0.01 0.04 -0.08 0.09 1.00 3347
#> sigma_timetime_3 -0.05 0.04 -0.13 0.03 1.00 2168
#> sigma_timetime_4 0.21 0.04 0.13 0.29 1.00 3483
#> Tail_ESS
#> Intercept 2071
#> groupgroup_2 1918
#> groupgroup_3 2158
#> timetime_2 2641
#> timetime_3 2113
#> timetime_4 1840
#> biomarker1 2949
#> biomarker2 2885
#> groupgroup_2:timetime_2 2727
#> groupgroup_3:timetime_2 2725
#> groupgroup_2:timetime_3 2149
#> groupgroup_3:timetime_3 2363
#> groupgroup_2:timetime_4 2053
#> groupgroup_3:timetime_4 2129
#> sigma_timetime_1 2791
#> sigma_timetime_2 3135
#> sigma_timetime_3 2622
#> sigma_timetime_4 2850
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).Marginals
Regardless of the choice of fixed effects formula,
brms.mmrm performs inference on the marginal distributions
at each treatment group and time point of the mean of the following
quantities:
- Response.
- Change from baseline, if you set
roleto"change"inbrm_data(). - Treatment difference, in terms of change from baseline.
- Effect size: treatment difference divided by the residual standard deviation.
To derive posterior draws of these marginals, use the
brm_marginal_draws() function.
draws <- brm_marginal_draws(data = data, formula = formula, model = model)
draws
#> $response
#> # A draws_df: 1000 iterations, 4 chains, and 12 variables
#> group_1|time_1 group_2|time_1 group_3|time_1 group_1|time_2 group_2|time_2
#> 1 1.2 -0.25 1.5 2.5 1.05
#> 2 1.2 -0.23 1.4 2.6 1.16
#> 3 1.2 -0.31 1.4 2.6 0.99
#> 4 1.2 -0.32 1.4 2.6 1.12
#> 5 1.2 -0.26 1.3 2.5 1.02
#> 6 1.1 -0.19 1.4 2.5 1.00
#> 7 1.1 -0.30 1.4 2.4 1.08
#> 8 1.2 -0.33 1.5 2.6 1.01
#> 9 1.1 -0.35 1.3 2.5 1.16
#> 10 1.2 -0.44 1.4 2.6 0.91
#> group_3|time_2 group_1|time_3 group_2|time_3
#> 1 2.7 1.8 -0.0050
#> 2 2.8 1.7 0.0115
#> 3 2.7 1.8 0.0939
#> 4 2.8 1.7 -0.0090
#> 5 2.6 1.8 -0.0371
#> 6 2.7 1.8 -0.0146
#> 7 2.8 1.9 -0.1218
#> 8 2.8 1.8 0.0043
#> 9 2.8 1.7 -0.0098
#> 10 2.8 1.6 0.1413
#> # ... with 3990 more draws, and 4 more variables
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
#>
#> $difference_time
#> # A draws_df: 1000 iterations, 4 chains, and 9 variables
#> group_1|time_2 group_1|time_3 group_1|time_4 group_2|time_2 group_2|time_3
#> 1 1.3 0.55 -1.6 1.3 0.25
#> 2 1.3 0.47 -1.3 1.4 0.24
#> 3 1.4 0.59 -1.5 1.3 0.41
#> 4 1.4 0.52 -1.3 1.4 0.31
#> 5 1.4 0.60 -1.3 1.3 0.23
#> 6 1.4 0.68 -1.4 1.2 0.17
#> 7 1.3 0.77 -1.3 1.4 0.18
#> 8 1.4 0.57 -1.4 1.3 0.33
#> 9 1.3 0.61 -1.2 1.5 0.34
#> 10 1.4 0.44 -1.2 1.4 0.58
#> group_2|time_4 group_3|time_2 group_3|time_3
#> 1 -1.7 1.2 0.107
#> 2 -1.8 1.3 0.227
#> 3 -1.7 1.3 0.244
#> 4 -1.5 1.4 0.410
#> 5 -1.5 1.3 0.574
#> 6 -1.9 1.3 0.350
#> 7 -1.4 1.4 0.391
#> 8 -1.6 1.2 0.043
#> 9 -1.4 1.5 0.377
#> 10 -1.4 1.5 0.297
#> # ... with 3990 more draws, and 1 more variables
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
#>
#> $difference_group
#> # A draws_df: 1000 iterations, 4 chains, and 6 variables
#> group_2|time_2 group_2|time_3 group_2|time_4 group_3|time_2 group_3|time_3
#> 1 0.0069 -0.30 -0.028 -0.089 -0.447
#> 2 0.0631 -0.23 -0.498 0.016 -0.244
#> 3 -0.1212 -0.19 -0.173 -0.123 -0.347
#> 4 0.0026 -0.21 -0.211 -0.081 -0.115
#> 5 -0.1005 -0.37 -0.212 -0.037 -0.027
#> 6 -0.2115 -0.51 -0.460 -0.044 -0.329
#> 7 0.0383 -0.59 -0.153 0.064 -0.377
#> 8 -0.1074 -0.24 -0.215 -0.229 -0.529
#> 9 0.1657 -0.27 -0.214 0.148 -0.238
#> 10 -0.0431 0.14 -0.233 0.060 -0.147
#> group_3|time_4
#> 1 -0.19
#> 2 -0.40
#> 3 -0.27
#> 4 -0.53
#> 5 -0.37
#> 6 -0.19
#> 7 -0.42
#> 8 -0.47
#> 9 -0.38
#> 10 -0.40
#> # ... with 3990 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
#>
#> $effect
#> # A draws_df: 1000 iterations, 4 chains, and 6 variables
#> group_2|time_2 group_2|time_3 group_2|time_4 group_3|time_2 group_3|time_3
#> 1 0.0070 -0.32 -0.022 -0.090 -0.465
#> 2 0.0624 -0.25 -0.426 0.016 -0.258
#> 3 -0.1264 -0.20 -0.138 -0.128 -0.373
#> 4 0.0026 -0.22 -0.178 -0.082 -0.118
#> 5 -0.0994 -0.39 -0.183 -0.036 -0.028
#> 6 -0.2190 -0.58 -0.370 -0.046 -0.375
#> 7 0.0403 -0.59 -0.128 0.067 -0.379
#> 8 -0.1078 -0.25 -0.164 -0.230 -0.544
#> 9 0.1702 -0.26 -0.192 0.152 -0.226
#> 10 -0.0462 0.15 -0.195 0.064 -0.158
#> group_3|time_4
#> 1 -0.15
#> 2 -0.34
#> 3 -0.22
#> 4 -0.45
#> 5 -0.32
#> 6 -0.15
#> 7 -0.35
#> 8 -0.36
#> 9 -0.34
#> 10 -0.33
#> # ... with 3990 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}If you need samples from these marginals averaged across time points,
e.g. an “overall effect size”, brm_marginal_draws_average()
can average the draws above across discrete time points (either all or a
user-defined subset).
brm_marginal_draws_average(draws = draws, data = data)
#> $response
#> # A draws_df: 1000 iterations, 4 chains, and 3 variables
#> group_1|average group_2|average group_3|average
#> 1 1.3 -0.28 1.4
#> 2 1.4 -0.27 1.4
#> 3 1.3 -0.30 1.4
#> 4 1.3 -0.27 1.4
#> 5 1.3 -0.27 1.3
#> 6 1.3 -0.31 1.4
#> 7 1.3 -0.27 1.4
#> 8 1.3 -0.31 1.4
#> 9 1.3 -0.24 1.4
#> 10 1.3 -0.32 1.4
#> # ... with 3990 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
#>
#> $difference_time
#> # A draws_df: 1000 iterations, 4 chains, and 3 variables
#> group_1|average group_2|average group_3|average
#> 1 0.076 -0.033 -0.165
#> 2 0.158 -0.065 -0.052
#> 3 0.175 0.015 -0.072
#> 4 0.213 0.073 -0.030
#> 5 0.216 -0.013 0.073
#> 6 0.225 -0.167 0.039
#> 7 0.279 0.044 0.033
#> 8 0.211 0.024 -0.199
#> 9 0.258 0.152 0.103
#> 10 0.210 0.163 0.048
#> # ... with 3990 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
#>
#> $difference_group
#> # A draws_df: 1000 iterations, 4 chains, and 2 variables
#> group_2|average group_3|average
#> 1 -0.109 -0.24
#> 2 -0.223 -0.21
#> 3 -0.160 -0.25
#> 4 -0.139 -0.24
#> 5 -0.229 -0.14
#> 6 -0.393 -0.19
#> 7 -0.235 -0.25
#> 8 -0.187 -0.41
#> 9 -0.106 -0.16
#> 10 -0.047 -0.16
#> # ... with 3990 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
#>
#> $effect
#> # A draws_df: 1000 iterations, 4 chains, and 2 variables
#> group_2|average group_3|average
#> 1 -0.110 -0.23
#> 2 -0.203 -0.20
#> 3 -0.154 -0.24
#> 4 -0.130 -0.22
#> 5 -0.223 -0.13
#> 6 -0.389 -0.19
#> 7 -0.227 -0.22
#> 8 -0.173 -0.38
#> 9 -0.093 -0.14
#> 10 -0.032 -0.14
#> # ... with 3990 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}The brm_marginal_summaries() function produces posterior
summaries of these marginals, and it includes the Monte Carlo standard
error (MCSE) of each estimate.
summaries <- brm_marginal_summaries(draws, level = 0.95)
summaries
#> # A tibble: 165 × 6
#> marginal statistic group time value mcse
#> <chr> <chr> <chr> <chr> <dbl> <dbl>
#> 1 difference_group lower group_2 time_2 -0.240 0.00686
#> 2 difference_group lower group_2 time_3 -0.651 0.0115
#> 3 difference_group lower group_2 time_4 -0.550 0.00994
#> 4 difference_group lower group_3 time_2 -0.261 0.00608
#> 5 difference_group lower group_3 time_3 -0.691 0.0166
#> 6 difference_group lower group_3 time_4 -0.751 0.0114
#> 7 difference_group mean group_2 time_2 0.0372 0.00303
#> 8 difference_group mean group_2 time_3 -0.174 0.00625
#> 9 difference_group mean group_2 time_4 -0.0783 0.00634
#> 10 difference_group mean group_3 time_2 0.0213 0.00307
#> # ℹ 155 more rowsThe brm_marginal_probabilities() function shows
posterior probabilities of the form,
\[ \begin{aligned} \text{Prob}(\text{treatment effect} > \text{threshold}) \end{aligned} \]
or
\[ \begin{aligned} \text{Prob}(\text{treatment effect} < \text{threshold}) \end{aligned} \]
brm_marginal_probabilities(
draws = draws,
threshold = c(-0.1, 0.1),
direction = c("greater", "less")
)
#> # A tibble: 12 × 5
#> direction threshold group time value
#> <chr> <dbl> <chr> <chr> <dbl>
#> 1 greater -0.1 group_2 time_2 0.828
#> 2 greater -0.1 group_2 time_3 0.375
#> 3 greater -0.1 group_2 time_4 0.532
#> 4 greater -0.1 group_3 time_2 0.814
#> 5 greater -0.1 group_3 time_3 0.324
#> 6 greater -0.1 group_3 time_4 0.242
#> 7 less 0.1 group_2 time_2 0.670
#> 8 less 0.1 group_2 time_3 0.866
#> 9 less 0.1 group_2 time_4 0.761
#> 10 less 0.1 group_3 time_2 0.712
#> 11 less 0.1 group_3 time_3 0.890
#> 12 less 0.1 group_3 time_4 0.936Finally, brm_marignal_data() computes marginal means and
confidence intervals on the response variable in the data, along with
other summary statistics.
summaries_data <- brm_marginal_data(data = data, level = 0.95)
summaries_data
#> # A tibble: 84 × 4
#> statistic group time value
#> <chr> <chr> <chr> <dbl>
#> 1 lower group_1 time_1 1.40
#> 2 lower group_1 time_2 2.71
#> 3 lower group_1 time_3 1.91
#> 4 lower group_1 time_4 0.0151
#> 5 lower group_2 time_1 -0.150
#> 6 lower group_2 time_2 1.23
#> 7 lower group_2 time_3 0.218
#> 8 lower group_2 time_4 -1.59
#> 9 lower group_3 time_1 1.57
#> 10 lower group_3 time_2 2.95
#> # ℹ 74 more rowsVisualization
The brm_plot_compare() function compares means and
intervals from many different models and data sources in the same plot.
First, we need the marginals of the data.
brm_plot_compare(
data = summaries_data,
model1 = summaries,
model2 = summaries
)
If you omit the marginals of the data, you can show inference on change from baseline or the treatment effect.
brm_plot_compare(
model1 = summaries,
model2 = summaries,
marginal = "difference_group" # treatment effect
)
Additional arguments let you control the primary comparison of interest (the color aesthetic), the horizontal axis, and the faceting variable.
brm_plot_compare(
model1 = summaries,
model2 = summaries,
marginal = "difference_group",
compare = "group",
axis = "time",
facet = "source" # model1 vs model2
)
Finally, brm_plot_draws() can plot the posterior draws
of the response, change from baseline, or treatment difference.
brm_plot_draws(draws = draws$difference_group)
The axis and facet arguments customize the
horizontal axis and faceting variable, respectively.
brm_plot_draws(
draws = draws$difference_group,
axis = "group",
facet = "time"
)