# Multivariate analysis

**URL:** <https://discourse.mc-stan.org/t/multivariate-analysis/34556>\
**Category:** brms\
**Tags:** brms\
**Created:** [March 21, 2024, 2:19pm UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556 "2024-03-21T14:19:39Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![adriaancampo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/adriaancampo/32/17926_2.png) [@adriaancampo](https://discourse.mc-stan.org/u/adriaancampo)\
**Post date:** [March 21, 2024, 2:19pm UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/1 "2024-03-21T14:19:39Z")

</div>

Hi everybody,

I am new to brms, and bayesian statistics. I would like to dig a bit deeper though, and I tried applying it to my own data.

I am running a multivariate model. I have 3 metrics, which I follow over Time (factor), and which I want to compare between 2 Groups (factor).  
The metrics are heavily correlated, and the distributions are definitely not normal (rather skewed). It is not an option to apply a transformation on the data.

I found in the online material how to run a multivariate model in brms (here: [Estimating Multivariate Models with brms](https://cran.r-project.org/web/packages/brms/vignettes/brms_multivariate.html)). However, this approach doesn’t support non normal underlying distributions.

Would it be an option to construct the formula like this?

Form1 ← bf(Response1 ~ 1 + Group \* Time + (1|Participant)) + skew\_normal()  
Form2 ← bf(Response2 ~ 1 + Group \* Time + (1|Participant)) + skew\_normal()  
Form3 ← bf(Response3 ~ 1 + Group \* Time + (1|Participant)) + skew\_normal()

ModelFormula ← Form1 + Form2 + Form3

But how do I account for the fact that the data are correlated? And how can I correct for multiple comparisons?

Thanks for your reply!  
A

---

<div class="post-metadata">

**Author:** ![wpetry](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wpetry/32/4852_2.png) [@wpetry](https://discourse.mc-stan.org/u/wpetry)\
**Post date:** [March 21, 2024, 4:54pm UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/2 "2024-03-21T16:54:29Z")

</div>

> [@adriaancampo](#):
>
> Would it be an option to construct the formula like this?
> 
> Form1 ← bf(Response1 ~ 1 + Group \* Time + (1|Participant)) + skew\_normal()  
> Form2 ← bf(Response2 ~ 1 + Group \* Time + (1|Participant)) + skew\_normal()  
> Form3 ← bf(Response3 ~ 1 + Group \* Time + (1|Participant)) + skew\_normal()
> 
> ModelFormula ← Form1 + Form2 + Form3

This will fit a multivariate model, but it won’t include residual correlations between the responses (`rescor = FALSE`) because `brms` doesn’t have a multivariate skew-normal.

> [@adriaancampo](#):
>
> But how do I account for the fact that the data are correlated?

Depending on your data and goals, you might be able to exploit correlations in the effects of `Participant` by specifying `(1 | p | Participant)`. The letter `p` doesn’t mean anything special, it’s just used as a grouping marker that let’s `brms` know which group-level effects should be correlated. I think the argument here might be that participants have consistent, correlated responses to 1-3 that mean that they differ in the set of intercepts their specific time series will have.

> [@adriaancampo](#):
>
> And how can I correct for multiple comparisons?

I’m not sure where the multiple comparisons issue would arise with your problem. Andrew Gelman has a [blog post](https://statmodeling.stat.columbia.edu/2016/08/22/bayesian-inference-completely-solves-the-multiple-comparisons-problem/) that summarizes thinking on multiple comparisons in Bayesian analyses.

---

<div class="post-metadata">

**Author:** ![Solomon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/solomon/32/2495_2.png) [@Solomon](https://discourse.mc-stan.org/u/Solomon)\
**Post date:** [March 22, 2024, 2:27pm UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/3 "2024-03-22T14:27:12Z")

</div>

As a half-measure, you can allow correlations in your upper-level deviation terms with syntax like this:

```no-highlight
Form1 <- bf(Response1 ~ 1 + Group * Time + (1 |p| Participant)) + skew_normal()
Form2 <- bf(Response2 ~ 1 + Group * Time + (1 |p| Participant)) + skew_normal()
Form3 <- bf(Response3 ~ 1 + Group * Time + (1 |p| Participant)) + skew_normal()

```

Notice the double `||` with some key in between. `p` does not need to be in your data, and you can use other characters like `X` or whatever. The trick is that by placing that `p` in all 3 of your formulas, you will allow the random intercepts for `Response1`, `Response2`, and `Response3` to correlate. You can learn more about this in the `brmsformula` section of the brms reference manual ([link](https://cran.r-project.org/web/packages/brms/brms.pdf)).

---

<div class="post-metadata">

**Author:** ![adriaancampo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/adriaancampo/32/17926_2.png) [@adriaancampo](https://discourse.mc-stan.org/u/adriaancampo)\
**Post date:** [April 30, 2024, 10:32am UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/4 "2024-04-30T10:32:02Z")

</div>

Hi, @Solomon and @wpetry, thanks for your helpful advice!

---

<div class="post-metadata">

**Author:** ![Johannes\_Schwenke](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/johannes_schwenke/32/20381_2.png) [@Johannes\_Schwenke](https://discourse.mc-stan.org/u/Johannes_Schwenke)\
**Post date:** [April 30, 2025, 3:41pm UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/5 "2025-04-30T15:41:47Z")

</div>

I have a question in a similar vein, so replying here instead of opening a new thread.  
I’d like to jointly estimate the survival and the rate of recurrent events (hospitalizations) for participants (needed for a while-alive estimand). They are likely negative correlated.

```no-highlight
formula_hosp_negbinomial <- bf(
    n_re | rate(survival_time) ~ trt + ecog_fstcnt,
    family = negbinomial()
    )

    formula_death_weibull <- bf(
    survival_time | cens(censoring) ~ trt + ecog_fstcnt,
    family = weibull()
    )

```

Could I use brms to somehow account for this correlation beyond what is captured by the covariates? I don’t think a random intercept at the patient level can work, given that I only have 1 observation per patient.

---

<div class="post-metadata">

**Author:** ![Solomon](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/solomon/32/2495_2.png) [@Solomon](https://discourse.mc-stan.org/u/Solomon)\
**Post date:** [April 30, 2025, 4:13pm UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/6 "2025-04-30T16:13:35Z")

</div>

I’m not aware that brms can combine different likelihoods like this.

---

<div class="post-metadata">

**Author:** ![Johannes\_Schwenke](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/johannes_schwenke/32/20381_2.png) [@Johannes\_Schwenke](https://discourse.mc-stan.org/u/Johannes_Schwenke)\
**Post date:** [May 1, 2025, 5:03am UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/7 "2025-05-01T05:03:51Z")

</div>

@Solomon thanks for the swift reply!

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [May 1, 2025, 6:33am UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/8 "2025-05-01T06:33:40Z")

</div>

This sounds like the kind of problem a multi-state survival model would be suited to. Alas while I know they exist, I’ve no experience with them myself, and I don’t know if they are implemented in Stan. There’s an explanation here: [flexsurv: A Platform for Parametric Survival Modeling in R - PMC](https://pmc.ncbi.nlm.nih.gov/articles/PMC5868723/#S27)

---

<div class="post-metadata">

**Author:** ![Johannes\_Schwenke](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/johannes_schwenke/32/20381_2.png) [@Johannes\_Schwenke](https://discourse.mc-stan.org/u/Johannes_Schwenke)\
**Post date:** [May 1, 2025, 8:05am UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/9 "2025-05-01T08:05:14Z")

</div>

Interesting. I know that [Frank Harrell](https://onlinelibrary.wiley.com/doi/full/10.1002/sim.10133) has developed packages for multi-state models which IIRC are based on stan, so I will look into that.

I also found [this publication](https://academic.oup.com/biometrics/article/79/4/3792/7585592) in which the authors address the correlation between survival and the recurrent event rate with a shared frailty term. The authors also provide their rstan code.

---

<div class="post-metadata">

**Author:** ![lldb14](https://avatars.discourse-cdn.com/v4/letter/l/f19dbf/32.png) [@lldb14](https://discourse.mc-stan.org/u/lldb14)\
**Post date:** [May 2, 2025, 8:08pm UTC](https://discourse.mc-stan.org/t/multivariate-analysis/34556/10 "2025-05-02T20:08:08Z")

</div>

Who is the likelihood of brms?  
What is the marginal correlation between responses for brms?  
Who is the posterior of brms?  
How does brms write when it has a random effect a posteriori?  
How does brms induce correlation between responses?