# Integrating emmeans with brms

**URL:** https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792
**Category:** brms
**Tags:** emmeans
**Created:** [September 2, 2020, 3:49am UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792 "2020-09-02T03:49:31Z")
**Posts on this page:** 9
**Page:** 1

<div class="post-metadata">

### Author: ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)
#### Post date: [September 2, 2020, 3:49am UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/1 "2020-09-02T03:49:32Z")

</div>

A recent post sparked discussion on integrating `emmeans` with `brms`. Moving it to a new topic so that it can be found later. Thx @rvlenth for bringing this up!

---

<div class="post-metadata">

### Author: ![rvlenth](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rvlenth/32/11734_2.png) [@rvlenth](https://discourse.mc-stan.org/u/rvlenth)
#### Post date: [August 30, 2020, 9:14pm UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/2 "2020-08-30T21:14:20Z")

</div>

FYI, what `emmeans` does is creates a grid of factor (and predictor) combinations, and for each node on the grid, the posterior MCMC sample of predictions at that node. Currently, the user can specify which parameter (in this case mu, sigma, or ndt. `emmeans()` may then further average selected ones of these together to obtain marginal means, back-transform them, etc.

What would be useful for a model like this is provisions to ensure that the grid involves all the predictors in all of the parts of the model; then for each node, the posterior sample of values of `ndt + exp(mu.hat + 0.5*sigma.hat^2)`. I believe it wouldn’t be too hard to do this, and I could help with the needed additions to `emm_basis.brms()`. The code would be similar to existing code for prediction. It’d help to have some understanding more generally across models handled by **brms** and how manageable is that number. I suppose we have various zero-inflated, hurdle, and truncated cases.

---

<div class="post-metadata">

### Author: ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [August 31, 2020, 8:44am UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/3 "2020-08-31T08:44:25Z")

</div>

In essence, this is what the `posterior_epred` method of brms does as it computes the mean predictions for each family ([https://github.com/paul-buerkner/brms/blob/master/R/posterior\_epred.R](https://github.com/paul-buerkner/brms/blob/master/R/posterior_epred.R)).

---

<div class="post-metadata">

### Author: ![rvlenth](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/rvlenth/32/11734_2.png) [@rvlenth](https://discourse.mc-stan.org/u/rvlenth)
#### Post date: [August 31, 2020, 8:07pm UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/4 "2020-08-31T20:07:19Z")

</div>

So… it seems like we could do some minor things to make it all work:

- We need to add an argument or enhance an existing argument to both emmeans-support methods to provide for this. Would it make sense to trap the case `dpar = "mean"`? However it is specified, I’ll refer to this as “mean” mode, subsequently.
- In `recover_data.brmsfit()`, for “mean” mode, we need to collect the terms for all fixed-effect components of the model so that we will obtain a grid of reference levels for all predictors involved in those terms. This would be similar to `emmeans:::recover_data.zeroinfl(..., mode = "response")`
- In `emm_basis.brmsfit()` in “mean” mode, it seems that all we will need to do is set the `linfct` slot to the identity matrix and `post.beta` slot to the result of `posterior_epred(object, newdata = grid, re_formula = NA)`. (And then `bhat` and `V` to its column means and covariance matrix).

This would enhance what’s possible with **brms** in the sense that **emmeans** has all kinds of built-in machinery for computing marginal averages thereof, contrasts, `by` variables, etc. In the example for `posterior_epred()`, for example, **emmeans** ’s code picks up on the fact that `carry` is nested in `period`, so in obtaining marginal averages over `carry`, it would know to average separately for each `period`.

---

<div class="post-metadata">

### Author: ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [September 1, 2020, 9:05am UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/5 "2020-09-01T09:05:41Z")

</div>

Your suggestions sound good. With regard to using `posterior_epred`, we have to be careful to disable everything except fixed effects though, since brms supports all sorts of additional terms, such as GPs and splines, which are probably another topic to discuss. Right now, it is not easily possible to (de)activate these terms but I am planning on adding such functionality in the future.

Right now, I don’t yet feel equipped to make the changes necessary to brms to be honest but would be happy to discuss this in more detail or review PRs in that direction of course.

---

<div class="post-metadata">

### Author: ![samantha.zambo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/samantha.zambo/32/8264_2.png) [@samantha.zambo](https://discourse.mc-stan.org/u/samantha.zambo)
#### Post date: [October 26, 2020, 10:25pm UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/6 "2020-10-26T22:25:48Z")

</div>

is this what the new dpar = “mean” does?

Could you clearly summarize why one would choose to use the posterior predictive distribution (dpar = “mean”) and when to use the default posterior in estimating the marginal means (emmeans)?

---

<div class="post-metadata">

### Author: ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [October 27, 2020, 7:53am UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/7 "2020-10-27T07:53:32Z")

</div>

I don’t know to be honest. There is no general answer. It depends on what you are interested it, that is, if you want to know the emmeans of a specific distributional parameters or of the implied PPD mean.

---

<div class="post-metadata">

### Author: ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [October 27, 2020, 8:46am UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/8 "2020-10-27T08:46:53Z")

</div>

Perhaps the following ideas could help. If you are interested in a specific process of your data generating distribution that is expressed by a specific distrihutional parameter(s) computer the emmeans for those. If you want an overall summary of the mean predictions, dpar = “mean” could be a good option. For most models the main parameter “mu” (default) is already the mean of the ppd so it makes literally no difference in these cases.

---

<div class="post-metadata">

### Author: ![samantha.zambo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/samantha.zambo/32/8264_2.png) [@samantha.zambo](https://discourse.mc-stan.org/u/samantha.zambo)
#### Post date: [October 27, 2020, 2:13pm UTC](https://discourse.mc-stan.org/t/integrating-emmeans-with-brms/17792/9 "2020-10-27T14:13:44Z")

</div>

So what I’m hearing is that the default of dpar = ‘mu’ is fine for general emmeans used to compare effects to see if they are different.

What about point.est = ‘mean’ vs ‘median’? I see the default in emmeans is ‘median’ while ‘mean’ is used in the brms documentation.
