I was wondering if it is possible to calculate the effect size in terms of Cohen’s d from the output of a mixed-effect model fitted in brms. I have seen different posts about this but the answers were not very consistent for me.
Suppose I had the following output:
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: popular ~ 1 + sex + extrav + (1 | class)
## Data: popular2data (Number of observations: 2000)
## Samples: 2 chains, each with iter = 3000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~class (Number of levels: 100)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 0.80 0.06 0.69 0.93 389 1.01
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 2.14 0.12 1.90 2.37 585 1.00
## sex 1.25 0.04 1.18 1.33 4125 1.00
## extrav 0.44 0.02 0.41 0.47 3382 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 0.77 0.01 0.75 0.79 4470 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Would it be possible to calculate the effect size for the predictor sex?
The background is that I am writing a thesis and my supervisor is not very familiar with Bayesian methods, so he would like to have Cohen’s d as he understands this better.
If Cohen’s d really doesn’t make any sense in this context, I would be happy if someone has recommendations on which measure of effect size one could introduce to more frequentist-minded people. I am very new to this topic myself but find it incredibly exciting.
Cohen’s D is very possible in Bayesian frameworks. In fact, in Kruschke’s influential 2013 paper on Bayesian t-tests (Bayesian estimation supersedes the t test - PubMed), one of the advantages of Bayesian computations of effect sizes is that we obtain the posterior distribution of the effect size, rather than a point estimate. While standard confidence intervals can be computed for the typical point-estimated Cohen’s D, access to the full posterior allows much more interesting inference (e.g., probability that the effect size is within certain effect ranges like “small” or “large”). There is a separate package in R for fitting the Bayesian t-test (https://cran.r-project.org/web/packages/BEST/vignettes/BEST.pdf), and there is a helpful blog post from Matti Vuorre on how to do this in brms (Sometimes I R: How to Compare Two Groups with Robust Bayesian Estimation in R).
Ultimately, you are not doing a typical t-test as you have other variables of interest in the model. This resource describes computing Cohen’s D in brms from a model kind of similar to what you have: A An introduction to Bayesian multilevel models using brms | Understanding rumination as a form of inner speech. There is also a convenience function in the report package from the easystats group that converts rstanarm models to effect sizes: Report the effect size(s) of a model or a test — report_effectsize • report. You might check out the source code (report/report_effectsize.R at master · easystats/report · GitHub) to see how they get the estimate, or you may even just consider estimating the model in rstanarm (should get very similar results, and might even be able to pass the same inits from the brms fit to the rstanarm fit) just so that you can get the estimates for your advisor (because at the end of the day, clean code is great but a defended thesis is best). Your model seems pretty simple, so I’d suspect that running it again in rstanarm is a trivial investment in time if you just want to get an estimate and move on with your life.
Edit to include a last minute thought:
You can also just consider standardizing your model and getting the standardized regression coefficients. Those should be reasonably understandable effect sizes in their own right and don’t require any fooling around with the posteriors. You can either standardize your variables prior to fitting the model (recoding Introduction to Mediation, Moderation, and Conditional Process Analysis) or you can utilize the convenience function standardize_parameters() from the effectsize package (also from the easystats group). You can read more about the latter option here: Parameters standardization — standardize_parameters • effectsize
In hierarchical models the variance is divided into within- and between-level variance. There is an effect size δ_t (delta total) taking into account all possible variances. See section A.5 here