usage.Rmd

---
title: "Usage"
output:
  rmarkdown::html_vignette:
    toc: true
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Usage}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



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.^[Covariates can be categorical too.] In general, the outcome variable can either be the raw response or change from baseline.


```r
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 rows
```

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.


```r
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_labels_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_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()`.


```r
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.


```r
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 + time
```

# Parameterization

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.


```r
options(contrasts = c(unordered = "contr.SAS", ordered = "contr.poly"))
```

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.


```r
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                       0
```

If you choose a different contrast method, a different model matrix may result.


```r
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                       0
```

# Priors

Some analyses require informative priors, others require non-informative ones. Please use [`brms`](https://paul-buerkner.github.io/brms/) to construct a prior suitable for your analysis. The [`brms`](https://paul-buerkner.github.io/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.


```r
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.


```r
model <- brm_model(data = data, formula = formula, refresh = 0)
```



The result is a `brms` model object.


```r
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
#> 
#> Regression Coefficients:
#>                         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:

1. Response.
2. Change from baseline, if you set `role` to `"change"` in `brm_data()`.
3. Treatment difference, in terms of change from baseline.
4. Effect size: treatment difference divided by the residual standard deviation.

To derive posterior draws of these marginals, use the `brm_marginal_draws()` function.


```r
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_1|time_2 group_1|time_3 group_1|time_4 group_2|time_1
#> 1             1.2            2.5            1.8         -0.386          -0.25
#> 2             1.2            2.6            1.7         -0.083          -0.23
#> 3             1.2            2.6            1.8         -0.278          -0.31
#> 4             1.2            2.6            1.7         -0.150          -0.32
#> 5             1.2            2.5            1.8         -0.160          -0.26
#> 6             1.1            2.5            1.8         -0.288          -0.19
#> 7             1.1            2.4            1.9         -0.181          -0.30
#> 8             1.2            2.6            1.8         -0.207          -0.33
#> 9             1.1            2.5            1.7         -0.070          -0.35
#> 10            1.2            2.6            1.6         -0.022          -0.44
#>    group_2|time_2 group_2|time_3 group_2|time_4
#> 1            1.05        -0.0050           -1.9
#> 2            1.16         0.0115           -2.0
#> 3            0.99         0.0939           -2.0
#> 4            1.12        -0.0090           -1.9
#> 5            1.02        -0.0371           -1.8
#> 6            1.00        -0.0146           -2.0
#> 7            1.08        -0.1218           -1.7
#> 8            1.01         0.0043           -1.9
#> 9            1.16        -0.0098           -1.8
#> 10           0.91         0.1413           -1.9
#> # ... 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).


```r
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.


```r
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 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}
$$


```r
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.936
```

Finally, `brm_marignal_data()` computes marginal means and confidence intervals on the response variable in the data, along with other summary statistics.


```r
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 rows
```

# Visualization

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.


```r
brm_plot_compare(
  data = summaries_data,
  model1 = summaries,
  model2 = summaries
)
```

![plot of chunk response](usage/response-1.png)

If you omit the marginals of the data, you can show inference on change from baseline or the treatment effect.


```r
brm_plot_compare(
  model1 = summaries,
  model2 = summaries,
  marginal = "difference_group" # treatment effect
)
```

![plot of chunk difference](usage/difference-1.png)

Additional arguments let you control the primary comparison of interest (the color aesthetic), the horizontal axis, and the faceting variable.


```r
brm_plot_compare(
  model1 = summaries,
  model2 = summaries,
  marginal = "difference_group",
  compare = "group",
  axis = "time",
  facet = "source" # model1 vs model2
)
```

![plot of chunk differencecustom](usage/differencecustom-1.png)

Finally, `brm_plot_draws()` can plot the posterior draws of the response, change from baseline, or treatment difference.


```r
brm_plot_draws(draws = draws$difference_group)
```

![plot of chunk draws](usage/draws-1.png)

The `axis` and `facet` arguments customize the horizontal axis and faceting variable, respectively.


```r
brm_plot_draws(
  draws = draws$difference_group,
  axis = "group",
  facet = "time"
)
```

![plot of chunk drawscustom](usage/drawscustom-1.png)