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

Interfaces brms
,
1

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”). 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.

2

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 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.