Is there a standard way to calculate effect size from a brms model?

Interfaces brms
,
1

Hi folks, is there a reasonably straightforward way to get posterior distributions over effect size (Cohen’s d, or similar) from a fitted model in brms? I’ve seen a few posts for how to do it using generated quantities in raw stancode, but felt like there must be some more user-friendly method that I’ve just missed. Thanks very much!

2

Perhaps something like:

library(brms)
mtcars$am <- factor(mtcars$am)

m <- brm(hp ~ am, mtcars)
draws <- as.data.frame(m)
draws$cohens_d <- draws$b_am1 / draws$sigma
3

Thanks for this - however this gives me negative estimates (β = -0.485 [-1.19, 0.225]). Aren’t effect sizes supposed to be positive? (or have I deeply misunderstood something - if so, my bad!)

4

That depends on what effect size measure you’re interested in. For a difference in means, it just depends on which group you subtract from which group. https://en.wikipedia.org/wiki/Effect_size#Difference_family:_Effect_sizes_based_on_differences_between_means

If you would have multiple studies comparing treatment vs control, and perform a meta-analysis, you’d really like to know the signs since they tell you whether the mean of the treatment was smaller or larger than the control.

1 Like
5

Ah, makes sense. Another silly question, then: is the sigma you used above able to stand in for the pooled SD in the formulation of Cohen’s d? Or is what you suggested a different measure of effect size that isn’t Cohen’s d? Thanks!

6

You need a measure of the within group variance, which is what the residual standard deviation supplies in this case. See e.g. https://stats.stackexchange.com/questions/171941/cohens-d-from-regression-coefficient

That said, there are lots of ways to express effect sizes, I was just using this an example since you didn’t give any details of what you wanted to do.

1 Like
7

Ah, this makes sense. Thank you so much for your help!

8

Hi again @Ax3man - sorry to bother you, but I’m trying to use the method above on a large hierarchical model, and I’m not exactly sure what parameter to extract. For the sake of example, if I have a model like

DV ~ Param1*Param2 + (1+Param1*Param2|Subject)

with a binary DV with about 30 subjects with 30 datapoints each.

If I call the code you mentioned, I don’t see a “sigma”, I see a bunch of sd_ parameters. Which one is the pooled residual standard error, to divide the coefficient by?

9

For a binary outcome you will need to use a different approach. As Wouter mentioned above, his approach was for estimating the standardised difference in means between two groups (i.e., Cohen’s d). When you have a binary outcome, the mean and variance aren’t as meaningful.

If you’re calling brms with family = bernoulli(link = "logit") then you can exponentiate the coefficients to give Odds Ratios, which can be used as measures of effect size for binary outcomes.

2 Likes
10

Thanks, this is really helpful :)

1 Like
11

There is an effect size δ_t (delta total) taking into account all possible variances. This is especially useful in hierarchical models. See section A.5 here

12

Cohen’s d is an attempt at a unitless index that is in units of subject to subject variation. So I think the denominator needs to be the standard deviation of Y, if the residuals have a symmetric distribution. I don’t recommend Cohen’s d in general, preferring to stick to the more interpretable original data units.

1 Like
13

I’m probably missing something, but I divided by the residual standard deviation, which I think is equal to the within-group standard deviation of Y, and analogous to Cohen’s within-group pooled SD.

1 Like