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!
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
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!)
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.
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!
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.
Ah, this makes sense. Thank you so much for your help!
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?
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.
Thanks, this is really helpful :)