The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
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.
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
#> patient time response group biomarker1 biomarker2 biomarker3
#> <chr> <chr> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 patient_001 time_1 1.11 group_1 1.31 -0.361 1.52
#> 2 patient_001 time_2 2.15 group_1 1.31 -0.361 1.52
#> 3 patient_001 time_3 2.54 group_1 1.31 -0.361 1.52
#> 4 patient_001 time_4 -1.73 group_1 1.31 -0.361 1.52
#> 5 patient_002 time_1 1.11 group_1 0.107 -2.44 -0.139
#> 6 patient_002 time_2 2.64 group_1 0.107 -2.44 -0.139
#> 7 patient_002 time_3 1.69 group_1 0.107 -2.44 -0.139
#> 8 patient_002 time_4 0.783 group_1 0.107 -2.44 -0.139
#> 9 patient_003 time_1 0.118 group_1 1.44 -0.419 -1.54
#> 10 patient_003 time_2 2.48 group_1 1.44 -0.419 -1.54
#> # ℹ 1,190 more rows
It is good practice to convert the time variable to an ordered factor so individual discrete time points correctly represent correct chronological order.2
raw_data <- raw_data |>
brm_data_chronologize(time = "time", levels = paste0("time_", seq_len(4)))
str(raw_data$time)
#> Ord.factor w/ 4 levels "time_1"<"time_2"<..: 1 2 3 4 1 2 3 4 1 2 ...
#> - attr(*, "contrasts")= num [1:4, 1:3] 0 1 0 0 0 0 1 0 0 0 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : chr [1:4] "time_1" "time_2" "time_3" "time_4"
#> .. ..$ : chr [1:3] "2" "3" "4"
Next, 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",
group = "group",
patient = "patient",
time = "time",
covariates = c("biomarker1", "biomarker2"),
reference_group = "group_1",
reference_time = "time_1"
)
data
#> # A tibble: 1,200 × 7
#> patient time response group biomarker1 biomarker2 biomarker3
#> <chr> <ord> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 patient_001 time_1 1.11 group_1 1.31 -0.361 1.52
#> 2 patient_001 time_2 2.15 group_1 1.31 -0.361 1.52
#> 3 patient_001 time_3 2.54 group_1 1.31 -0.361 1.52
#> 4 patient_001 time_4 -1.73 group_1 1.31 -0.361 1.52
#> 5 patient_002 time_1 1.11 group_1 0.107 -2.44 -0.139
#> 6 patient_002 time_2 2.64 group_1 0.107 -2.44 -0.139
#> 7 patient_002 time_3 1.69 group_1 0.107 -2.44 -0.139
#> 8 patient_002 time_4 0.783 group_1 0.107 -2.44 -0.139
#> 9 patient_003 time_1 0.118 group_1 1.44 -0.419 -1.54
#> 10 patient_003 time_2 2.48 group_1 1.44 -0.419 -1.54
#> # ℹ 1,190 more rows
class(data)
#> [1] "brms_mmrm_data" "tbl_df" "tbl" "data.frame"
attributes <- attributes(data)
attributes$row.names <- NULL
str(attributes)
#> List of 9
#> $ names : chr [1:7] "patient" "time" "response" "group" ...
#> $ brm_outcome : 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"
#> $ class : chr [1:4] "brms_mmrm_data" "tbl_df" "tbl" "data.frame"
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 + time
For the purposes of our example, we choose a fully parameterized analysis of the raw response.
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 time2 time3 time4 biomarker1 biomarker2 groupgroup_1:time2
#> 1 1 1 0 0 0 0 1.3126508 -0.3608809 0
#> 2 1 1 0 1 0 0 1.3126508 -0.3608809 1
#> 3 1 1 0 0 1 0 1.3126508 -0.3608809 0
#> 4 1 1 0 0 0 1 1.3126508 -0.3608809 0
#> 5 1 1 0 0 0 0 0.1068624 -2.4441488 0
#> 6 1 1 0 1 0 0 0.1068624 -2.4441488 1
#> groupgroup_2:time2 groupgroup_1:time3 groupgroup_2:time3 groupgroup_1:time4 groupgroup_2:time4
#> 1 0 0 0 0 0
#> 2 0 0 0 0 0
#> 3 0 1 0 0 0
#> 4 0 0 0 1 0
#> 5 0 0 0 0 0
#> 6 0 0 0 0 0
If 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 time2 time3 time4 biomarker1 biomarker2 groupgroup_2:time2
#> 1 1 0 0 0 0 0 1.3126508 -0.3608809 0
#> 2 1 0 0 1 0 0 1.3126508 -0.3608809 0
#> 3 1 0 0 0 1 0 1.3126508 -0.3608809 0
#> 4 1 0 0 0 0 1 1.3126508 -0.3608809 0
#> 5 1 0 0 0 0 0 0.1068624 -2.4441488 0
#> 6 1 0 0 1 0 0 0.1068624 -2.4441488 0
#> groupgroup_3:time2 groupgroup_2:time3 groupgroup_3:time3 groupgroup_2:time4 groupgroup_3:time4
#> 1 0 0 0 0 0
#> 2 0 0 0 0 0
#> 3 0 0 0 0 0
#> 4 0 0 0 0 0
#> 5 0 0 0 0 0
#> 6 0 0 0 0 0
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 nlpar lb ub source
#> (flat) b default
#> (flat) b biomarker1 (vectorized)
#> (flat) b biomarker2 (vectorized)
#> (flat) b groupgroup_2 (vectorized)
#> (flat) b groupgroup_2:time2 (vectorized)
#> (flat) b groupgroup_2:time3 (vectorized)
#> (flat) b groupgroup_2:time4 (vectorized)
#> (flat) b groupgroup_3 (vectorized)
#> (flat) b groupgroup_3:time2 (vectorized)
#> (flat) b groupgroup_3:time3 (vectorized)
#> (flat) b groupgroup_3:time4 (vectorized)
#> (flat) b time2 (vectorized)
#> (flat) b time3 (vectorized)
#> (flat) b time4 (vectorized)
#> lkj(1) cortime default
#> student_t(3, 0.9, 2.5) Intercept default
#> (flat) b sigma default
#> (flat) b timetime_1 sigma (vectorized)
#> (flat) b timetime_2 sigma (vectorized)
#> (flat) b timetime_3 sigma (vectorized)
#> (flat) b timetime_4 sigma (vectorized)
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.
The result is a brms
model object with extra list
elements brms.mmrm_data
and brms.mmrm_formula
to keep track of the data and formula used to fit the model.
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: modeled_data (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 Tail_ESS
#> cortime(time_1,time_2) 0.42 0.05 0.33 0.51 1.00 3516 3008
#> cortime(time_1,time_3) -0.80 0.02 -0.84 -0.76 1.00 2361 2594
#> cortime(time_2,time_3) -0.54 0.04 -0.62 -0.46 1.00 3653 3411
#> cortime(time_1,time_4) -0.29 0.05 -0.39 -0.18 1.00 3964 3061
#> cortime(time_2,time_4) 0.03 0.06 -0.09 0.14 1.00 3444 3347
#> cortime(time_3,time_4) -0.10 0.06 -0.21 0.01 1.00 3847 3221
#>
#> Regression Coefficients:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 1.22 0.09 1.05 1.39 1.01 1290 2297
#> groupgroup_2 -1.54 0.12 -1.77 -1.30 1.01 1524 2563
#> groupgroup_3 0.18 0.13 -0.07 0.43 1.01 1248 2097
#> time2 1.31 0.10 1.12 1.51 1.00 1897 2360
#> time3 0.51 0.17 0.17 0.85 1.01 1597 2455
#> time4 -1.45 0.17 -1.78 -1.11 1.01 1393 2103
#> biomarker1 -0.01 0.01 -0.03 0.01 1.00 7162 3250
#> biomarker2 0.02 0.01 0.00 0.04 1.00 6102 3060
#> groupgroup_2:time2 0.04 0.15 -0.25 0.33 1.00 2286 2991
#> groupgroup_3:time2 0.02 0.15 -0.27 0.32 1.00 2347 2986
#> groupgroup_2:time3 -0.16 0.24 -0.63 0.31 1.01 1733 2596
#> groupgroup_3:time3 -0.19 0.25 -0.68 0.29 1.01 1541 2367
#> groupgroup_2:time4 -0.06 0.25 -0.55 0.42 1.01 1623 2689
#> groupgroup_3:time4 -0.26 0.24 -0.75 0.21 1.00 1517 2562
#> sigma_timetime_1 -0.12 0.04 -0.20 -0.04 1.00 2559 2658
#> sigma_timetime_2 0.00 0.04 -0.07 0.09 1.00 3796 3159
#> sigma_timetime_3 -0.05 0.04 -0.13 0.03 1.00 2346 2891
#> sigma_timetime_4 0.21 0.04 0.13 0.29 1.00 3982 3364
#>
#> 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).
model$brms.mmrm_data
#> # A tibble: 1,200 × 7
#> patient time response group biomarker1 biomarker2 biomarker3
#> <chr> <ord> <dbl> <chr> <dbl> <dbl> <dbl>
#> 1 patient_001 time_1 1.11 group_1 1.31 -0.361 1.52
#> 2 patient_001 time_2 2.15 group_1 1.31 -0.361 1.52
#> 3 patient_001 time_3 2.54 group_1 1.31 -0.361 1.52
#> 4 patient_001 time_4 -1.73 group_1 1.31 -0.361 1.52
#> 5 patient_002 time_1 1.11 group_1 0.107 -2.44 -0.139
#> 6 patient_002 time_2 2.64 group_1 0.107 -2.44 -0.139
#> 7 patient_002 time_3 1.69 group_1 0.107 -2.44 -0.139
#> 8 patient_002 time_4 0.783 group_1 0.107 -2.44 -0.139
#> 9 patient_003 time_1 0.118 group_1 1.44 -0.419 -1.54
#> 10 patient_003 time_2 2.48 group_1 1.44 -0.419 -1.54
#> # ℹ 1,190 more rows
Under the missing at random (MAR) assumptions, MMRMs do not require
imputation (Holzhauer and Weber (2024)).
However, if the outcomes in your data are not missing at random, or if
you are targeting an alternative estimand, then you may need to impute
missing outcomes. brms.mmrm
can leverage either of the two
alternative solutions described at https://paul-buerkner.github.io/brms/articles/brms_missings.html.
To impute missing outcomes before model fitting, first use create a
list of imputed datasets using the multiple imputation method of your
choice. The rbmi
package is uniquely suited to the multiple imputation of continuous
longitudinal clinical trial data.
variables <- rbmi::set_vars(
outcome = "response",
visit = "time",
subjid = "patient",
group = "group",
covariates = c("biomarker1", "biomarker2")
)
imputation_draws <- rbmi::draws(
data = data |>
mutate(
patient = as.factor(patient),
group = as.factor(group)
),
vars = variables,
method = rbmi::method_condmean(type = "jackknife"),
quiet = TRUE
)
imputation_run <- rbmi::impute(
draws = imputation_draws,
references = c(
group_1 = "group_1",
group_2 = "group_1",
group_3 = "group_1"
)
)
imputed_datasets <- rbmi::extract_imputed_dfs(imputation_run)
At this point, imputed_datasets
is a list of data frames
with the response variable imputed with multiple imputation. Simply
supply this list to the imputed
argument of
brm_model()
. Internally, brm_model()
calls
brms::brm_multiple(data = imputed, formula = formula)
instead of brms::brm(data = data, formula = formula)
to fit
an MMRM to each of the individual imputed datasets in the
imputed
object.
model <- brm_model(
data = data, # Yes, please supply the original non-imputed dataset too.
formula = formula,
imputed = imputed_datasets,
refresh = 0
)
Unless you set combine = FALSE
in
brm_model()
, brms
automatically combines
posterior samples across imputed datasets. This means the downstream
post-processing workflow below is exactly the same as the non-imputation
case.
Alternatively, to conduct imputation during the fitting of that
model, set model_missing_outcomes
to TRUE
in
brm_formula()
. This formula uses
response | mi()
instead of just response
on
the left-hand side to tell brms
to model each missing
outcome as a model parameter. To use this type of imputation, simply
supply the returned formula object to the formula
argument
of brm_model()
.
brm_formula(data, model_missing_outcomes = TRUE)
#> response | mi() ~ group + group:time + time + biomarker1 + biomarker2 + unstr(time = time, gr = patient)
#> sigma ~ 0 + time
Unlike imputation before model fitting, this approach requires only one fit of the model. However, that model will sample posterior draws for each missing outcome as if it were a model parameter, so the MCMC may run slower and produce a larger output object.
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:
reference_time
argument of
brm_data()
.To derive posterior draws of these marginals, use the
brm_marginal_draws()
function.
draws <- brm_marginal_draws(model = model)
draws
#> $response
#> # A draws_df: 1000 iterations, 4 chains, and 12 variables
#> group_1|time_1 group_1|time_2 group_1|time_3 group_1|time_4 group_2|time_1 group_2|time_2
#> 1 1.2 2.5 1.8 -0.25 -0.41 0.87
#> 2 1.4 2.6 1.6 -0.32 -0.32 1.11
#> 3 1.4 2.5 1.6 -0.19 -0.37 1.08
#> 4 1.4 2.6 1.6 -0.39 -0.34 0.91
#> 5 1.3 2.6 1.7 -0.36 -0.46 1.06
#> 6 1.3 2.5 1.7 -0.29 -0.52 0.89
#> 7 1.2 2.5 1.6 -0.22 -0.13 0.99
#> 8 1.2 2.6 1.8 -0.38 -0.34 0.84
#> 9 1.2 2.6 1.7 -0.36 -0.17 1.06
#> 10 1.3 2.6 1.6 -0.21 -0.16 1.00
#> group_2|time_3 group_2|time_4
#> 1 0.101 -1.8
#> 2 0.048 -1.8
#> 3 0.013 -1.7
#> 4 0.109 -1.7
#> 5 0.201 -1.8
#> 6 0.212 -1.7
#> 7 -0.196 -1.7
#> 8 0.141 -2.0
#> 9 -0.121 -1.9
#> 10 -0.086 -2.0
#> # ... 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 group_2|time_4
#> 1 1.4 0.66 -1.4 1.3 0.512 -1.4
#> 2 1.2 0.20 -1.7 1.4 0.370 -1.5
#> 3 1.1 0.21 -1.6 1.4 0.380 -1.3
#> 4 1.2 0.29 -1.7 1.2 0.451 -1.4
#> 5 1.3 0.37 -1.7 1.5 0.663 -1.3
#> 6 1.2 0.41 -1.6 1.4 0.735 -1.2
#> 7 1.3 0.42 -1.4 1.1 -0.066 -1.6
#> 8 1.4 0.64 -1.5 1.2 0.481 -1.6
#> 9 1.4 0.44 -1.6 1.2 0.049 -1.7
#> 10 1.4 0.36 -1.5 1.2 0.072 -1.8
#> group_3|time_2 group_3|time_3
#> 1 1.5 0.26
#> 2 1.4 0.28
#> 3 1.4 0.21
#> 4 1.2 0.44
#> 5 1.3 0.42
#> 6 1.3 0.33
#> 7 1.3 0.63
#> 8 1.3 0.54
#> 9 1.4 0.47
#> 10 1.4 0.31
#> # ... 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 group_3|time_4
#> 1 -0.089 -0.15 0.013 0.1075 -0.4044 -0.332
#> 2 0.224 0.17 0.212 0.1637 0.0837 -0.167
#> 3 0.311 0.17 0.272 0.2547 0.0047 0.092
#> 4 0.020 0.17 0.352 0.0071 0.1494 0.135
#> 5 0.243 0.29 0.343 0.0387 0.0480 -0.211
#> 6 0.231 0.32 0.410 0.1161 -0.0844 -0.181
#> 7 -0.180 -0.48 -0.171 0.0210 0.2130 -0.134
#> 8 -0.221 -0.16 -0.094 -0.1073 -0.0997 0.087
#> 9 -0.172 -0.39 -0.097 -0.0357 0.0301 -0.188
#> 10 -0.212 -0.29 -0.359 -0.0089 -0.0580 -0.018
#> # ... 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 group_3|time_4
#> 1 -0.086 -0.16 0.010 0.104 -0.4215 -0.252
#> 2 0.232 0.20 0.186 0.169 0.0946 -0.146
#> 3 0.316 0.19 0.221 0.259 0.0051 0.074
#> 4 0.019 0.17 0.265 0.007 0.1561 0.102
#> 5 0.248 0.32 0.263 0.040 0.0521 -0.162
#> 6 0.242 0.34 0.317 0.122 -0.0901 -0.140
#> 7 -0.164 -0.48 -0.142 0.019 0.2121 -0.111
#> 8 -0.230 -0.16 -0.078 -0.111 -0.1009 0.073
#> 9 -0.174 -0.39 -0.087 -0.036 0.0297 -0.168
#> 10 -0.213 -0.29 -0.311 -0.009 -0.0582 -0.016
#> # ... with 3990 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
#>
#> $sigma
#> # A tibble: 4,000 × 12
#> `group_1|time_1` `group_1|time_2` `group_1|time_3` `group_1|time_4` `group_2|time_1`
#> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.965 1.03 0.959 1.32 0.965
#> 2 0.810 0.967 0.885 1.14 0.810
#> 3 0.869 0.982 0.916 1.23 0.869
#> 4 0.896 1.02 0.957 1.33 0.896
#> 5 0.895 0.978 0.921 1.30 0.895
#> 6 0.861 0.955 0.937 1.29 0.861
#> 7 0.890 1.10 1.00 1.21 0.890
#> 8 0.944 0.962 0.989 1.20 0.944
#> 9 0.961 0.988 1.01 1.12 0.961
#> 10 0.932 0.996 0.996 1.15 0.932
#> # ℹ 3,990 more rows
#> # ℹ 7 more variables: `group_2|time_2` <dbl>, `group_2|time_3` <dbl>, `group_2|time_4` <dbl>,
#> # `group_3|time_1` <dbl>, `group_3|time_2` <dbl>, `group_3|time_3` <dbl>,
#> # `group_3|time_4` <dbl>
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.31 1.4
#> 2 1.3 -0.25 1.4
#> 3 1.3 -0.23 1.4
#> 4 1.3 -0.27 1.4
#> 5 1.3 -0.24 1.4
#> 6 1.3 -0.27 1.4
#> 7 1.3 -0.26 1.3
#> 8 1.3 -0.33 1.4
#> 9 1.3 -0.27 1.4
#> 10 1.3 -0.31 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.2053 0.1296 -0.0043
#> 2 -0.1020 0.1010 -0.0752
#> 3 -0.0699 0.1812 0.0470
#> 4 -0.0767 0.1023 0.0205
#> 5 -0.0034 0.2901 -0.0449
#> 6 0.0132 0.3334 -0.0367
#> 7 0.1032 -0.1752 0.1366
#> 8 0.1645 0.0078 0.1245
#> 9 0.0810 -0.1391 0.0167
#> 10 0.0895 -0.1982 0.0611
#> # ... 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.076 -0.210
#> 2 0.203 0.027
#> 3 0.251 0.117
#> 4 0.179 0.097
#> 5 0.294 -0.042
#> 6 0.320 -0.050
#> 7 -0.278 0.033
#> 8 -0.157 -0.040
#> 9 -0.220 -0.064
#> 10 -0.288 -0.028
#> # ... 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.078 -0.190
#> 2 0.204 0.039
#> 3 0.241 0.113
#> 4 0.152 0.088
#> 5 0.277 -0.023
#> 6 0.300 -0.036
#> 7 -0.263 0.040
#> 8 -0.155 -0.047
#> 9 -0.216 -0.058
#> 10 -0.272 -0.028
#> # ... with 3990 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
#>
#> $sigma
#> # A tibble: 4,000 × 3
#> `group_1|average` `group_2|average` `group_3|average`
#> <dbl> <dbl> <dbl>
#> 1 1.07 1.07 1.07
#> 2 0.951 0.951 0.951
#> 3 1.00 1.00 1.00
#> 4 1.05 1.05 1.05
#> 5 1.02 1.02 1.02
#> 6 1.01 1.01 1.01
#> 7 1.05 1.05 1.05
#> 8 1.02 1.02 1.02
#> 9 1.02 1.02 1.02
#> 10 1.02 1.02 1.02
#> # ℹ 3,990 more rows
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: 225 × 6
#> marginal statistic group time value mcse
#> <chr> <chr> <chr> <chr> <dbl> <dbl>
#> 1 difference_group lower group_2 time_2 -0.251 0.00627
#> 2 difference_group lower group_2 time_3 -0.635 0.0132
#> 3 difference_group lower group_2 time_4 -0.548 0.0113
#> 4 difference_group lower group_3 time_2 -0.268 0.0109
#> 5 difference_group lower group_3 time_3 -0.678 0.0193
#> 6 difference_group lower group_3 time_4 -0.748 0.0133
#> 7 difference_group mean group_2 time_2 0.0400 0.00310
#> 8 difference_group mean group_2 time_3 -0.160 0.00589
#> 9 difference_group mean group_2 time_4 -0.0599 0.00625
#> 10 difference_group mean group_3 time_2 0.0242 0.00307
#> # ℹ 215 more rows
The 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.83
#> 2 greater -0.1 group_2 time_3 0.405
#> 3 greater -0.1 group_2 time_4 0.572
#> 4 greater -0.1 group_3 time_2 0.809
#> 5 greater -0.1 group_3 time_3 0.348
#> 6 greater -0.1 group_3 time_4 0.249
#> 7 less 0.1 group_2 time_2 0.665
#> 8 less 0.1 group_2 time_3 0.854
#> 9 less 0.1 group_2 time_4 0.741
#> 10 less 0.1 group_3 time_2 0.709
#> 11 less 0.1 group_3 time_3 0.880
#> 12 less 0.1 group_3 time_4 0.930
Finally, 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> <ord> <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 rows
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.
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.
The axis
and facet
arguments customize the
horizontal axis and faceting variable, respectively.
Covariates can be categorical too.↩︎
Ordered factors usually have polynomial contrasts
(contr.poly()
), which makes the fixed effect
parameterization counterintuitive. However, brm_data()
manually sets the contrasts of the time variable to be treatment
contrasts (contr.treatment()
) so there are individual fixed
effects for individual discrete time points. To revert back to
polynomial contrasts, use contr.poly()
after the call to
brm_data()
.↩︎
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.