# Getting brms to not use a reference category (one coefficient per category level)

**URL:** <https://discourse.mc-stan.org/t/getting-brms-to-not-use-a-reference-category-one-coefficient-per-category-level/25109>\
**Category:** brms\
**Tags:** specification, meta-analysis\
**Created:** [November 3, 2021, 5:07pm UTC](https://discourse.mc-stan.org/t/getting-brms-to-not-use-a-reference-category-one-coefficient-per-category-level/25109 "2021-11-03T17:07:53Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Bjoern](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bjoern/32/5429_2.png) [@Bjoern](https://discourse.mc-stan.org/u/Bjoern)\
**Post date:** [November 3, 2021, 5:07pm UTC](https://discourse.mc-stan.org/t/getting-brms-to-not-use-a-reference-category-one-coefficient-per-category-level/25109/1 "2021-11-03T17:07:53Z")

</div>

I’m trying out a simple network meta-analysis with `brms` following the approach in [Piepho et al. (2012; “The Use of Two-Way Linear Mixed Models in Multitreatment Meta-Analysis”)](https://doi.org/10.1111/j.1541-0420.2012.01786.x). Ideally, I’d want to specify a model like `est | se(stderr) ~ 0 + treatment + study` and want to have a coefficient for each unique value of treatment and study. I realize that that’s in a sense over-parameterized and that in a frequentist approach I would only have number of categories - 1 coefficients. However, by having priors, I presumably can presumably deal with the over-parametrization (and yes, I realize sampling will be less efficient). Why do I want to? I feel it’s easier for me to think about what priors are reasonable. Otherwise, if I set independent priors, the a-priori uncertainty is different for the reference category than the other categories etc., which seems a bit unnatural.

I tried to declare my factors differently by first setting some `options(contrasts=...)` options etc., but I did not figure out how to tell `brms` what I want here.

```
library(tidyverse)
library(brms)

set.seed(123)
example = tibble(
  study = rep(1:20, each=2), 
  treatment = c( rep(c(0L,1L), 4), rep(c(0L,2L), 4), rep(c(1L,2L), 4), rep(c(2L,3L), 4), rep(c(1L,3L), 4))) %>%
  mutate(est=rnorm(n=n(), treatment, 1) + rep(rnorm(n=20, 0, 1), each=2), # outcome per trial arm based on treatment and random study effect
         stderr=1) %>% 
  mutate(treatment=factor(treatment), 
         study=factor(study))

brmfit1 = brm(data=example, 
              est | se(stderr) ~ 0 + treatment + study, 
              family = gaussian(link="identity"),
              prior=prior(student_t(4, 0, 2.5), class = "b"))
brmfit1

```

As you can see, study 1 is used as a reference level for studies:

```
Term Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
treatment0 -0.50 0.55 -1.55 0.59 1.00 1296 2196
treatment1 0.61 0.50 -0.38 1.60 1.00 872 1773
treatment2 1.50 0.53 0.47 2.56 1.00 857 1734
treatment3 2.48 0.61 1.25 3.69 1.00 944 1800
study2 0.97 0.80 -0.56 2.56 1.00 1957 2710
study3 0.09 0.78 -1.47 1.65 1.00 2063 2600
study4 2.05 0.81 0.41 3.66 1.00 1906 2510
study5 1.03 0.81 -0.54 2.66 1.00 1889 2567
study6 0.15 0.79 -1.39 1.72 1.00 1763 2607
study7 0.31 0.81 -1.27 1.89 1.00 1761 2548
study8 0.60 0.79 -0.95 2.20 1.00 1936 2984
study9 0.46 0.80 -1.10 2.05 1.00 1676 2329
study10 0.43 0.80 -1.12 1.97 1.00 1525 2170
study11 0.06 0.79 -1.49 1.57 1.00 1691 2817
study12 -0.40 0.82 -1.98 1.21 1.00 1840 2308
study13 -0.64 0.82 -2.23 0.95 1.00 1764 2539
study14 2.20 0.84 0.52 3.83 1.00 1668 2502
study15 0.31 0.83 -1.27 1.98 1.00 1684 2122
study16 1.93 0.84 0.26 3.55 1.00 1617 2772
study17 -0.20 0.80 -1.71 1.34 1.00 1697 2408
study18 1.67 0.84 0.11 3.30 1.00 1623 2542
study19 0.75 0.83 -0.87 2.36 1.00 1805 2746
study20 0.31 0.81 -1.32 1.91 1.00 1621 2540

```

This was under Windows 10 using R version 4.1.1 and `brms` 2.16.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:** [November 10, 2021, 1:18pm UTC](https://discourse.mc-stan.org/t/getting-brms-to-not-use-a-reference-category-one-coefficient-per-category-level/25109/2 "2021-11-10T13:18:18Z")

</div>

I think you are missing some important background on meta-analyses. Yoy may want to check out [Bayesian meta-analysis in brms | A. Solomon Kurz](https://solomonkurz.netlify.app/post/bayesian-meta-analysis/) or other general resource.

Including `study` as a “fixed” effect means there will be no sharing of information across studies. If you treat study as “varying” effect (i.e. `est | se(stderr) ~ 0 + treatment + (1 | study)`, you’ll get partial pooling of the results and also a separate coefficient for each study. Your prior will then be not on the study coefficients, but on the standard deviation of effects in individual studies from the overall mean (which is IMHO a reasonably meaningful quantity).

Note however that in the model as you wrote it (and also in the “varying effect” version) you assume the for all pairs of studies `treatment_study_a - control_study_a = treatment_study_b - control_study_b`, which appears dubious and you may in fact want to have a `(1 + treatment | study)` term.

Best of luck with your metaanalysis.