Implementing R2D2 prior in a mutlivariate betabinomial regression

Modeling
1

I have been having trouble getting a betabinomial mutlvariate regression with splines to converge using the default priors in brms. It was recommended to me to use an R2D2 prior with main=TRUE. However I am having trouble getting the prior to work.

Briefly I am conducting a longitudinal multivariate regression in brms of three count outcomes - days of heroin, cannabis, and alcohol use in the previous 28-day measurement period - based on days of amphetamine use at start of treatment. This is the dataset

atsUseInOTP_alloutcomes_noMissing_di.RData (8.4 KB)

And this is the model I am trying to run

bform_spline2_mvr <- bf(mvbind(cann28, alc28, heroin28) | trials(set) ~ ats_baseline + s(yearsFromStart) + s(yearsFromStart, ats_baseline) + (1|p|encID))


# 
fit_mvr_spline2 <- brm(formula = bform_spline2_mvr,
                       family = beta_binomial(),
                       data = d,
                       seed = 1234,
                       refresh = 0,
                       chains = 4,
                       cores = 4,
                       warmup = 1e3,
                       iter = 3e3)

Now this works but leads to huge numbers of divergent transitions and ESS of 10-20. I suggest not running it as it takes a long time.

Based on this case study I tried a version of the model with an R2D2 prior, as follows

fit_mvr_spline <- brm(formula = bform_spline2_mvr,
                      family = beta_binomial(),
                      data = d,
                       prior = prior(R2D2(mean_R2 = 1/3, 
                                          prec_R2 = 3, 
                                          cons_D2 = 1/3,
                                          main = TRUE),
                                     class = b),
                       save_pars = save_pars(all=TRUE),
                       seed = 1234,
                       refresh = 0,
                       chains = 4,
                       cores = 4,
                       warmup = 1e3,
                       iter = 3e3)

But I get the following error message

Error: The following priors do not correspond to any model parameter: 
b ~ R2D2(mean_R2 = 1/3, prec_R2 = 3, cons_D2 = 1/3, main = TRUE)
Function 'default_prior' might be helpful to you.

I’ve treid removing arguments and different values, and different classes but to no avail. Can anyone tell me what I’m doing wrong?

2

Here too, resp = "..." is needed (not just class). ;)

3

Thank you so much @amynang. So would you have multiple priors within the same brms model, one for each response variable?

I have three outcomes - cann28, alc28 and heroin28- so would it be something like below?

fit_mvr_spline2 <- brm(formula = bform_spline2_mvr,
                       family = beta_binomial(),
                       data = d,
                       prior = prior(R2D2(mean_R2 = 1/3, 
                                          prec_R2 = 3, 
                                          cons_D2 = 1/3,
                                          main = TRUE),
                                     class = b,
                                     resp = "cann28"),
                               prior(R2D2(mean_R2 = 1/3, 
                                          prec_R2 = 3, 
                                          cons_D2 = 1/3,
                                          main = TRUE),
                                     class = b,
                                     resp = "heroin28"),
                               prior(R2D2(mean_R2 = 1/3, 
                                          prec_R2 = 3, 
                                          cons_D2 = 1/3,
                                          main = TRUE),
                                     class = b,
                                     resp = "alc28"),
                       save_pars = save_pars(all=TRUE),
                       seed = 1234,
                       refresh = 0,
                       chains = 4,
                       cores = 4,
                       warmup = 1e3,
                       iter = 3e3)
4

Check the R2D2 prior documentation at R2D2 Priors in brms — R2D2 • brms. Specifically, check what is said about option main and what is said in section Details. The documentation states it, but could be more clear that only on R2D2 prior should have main=TRUE, and only that prior defines the R2D2 parameters. It seems that if you have only one R2D2 prior, you should use main=FALSE as main=TRUE implies there would be other R2D2 priors. Also, note that you should put R2D2 prior also for sds (SDs of smoothing splines) as in the case study you linked to. I don’t think you need to define resp for class = b, or at least I didn’t do that.

5

Hi @llewmills. Saw your post and have a few thoughts.

I note you said

So my first thought is that you are not going to simply fix this with R2D2 priors. It sounds like there is some pathology in the model and you need to explore this. Did you make some pairs plots? That would be a good starting point. Don’t get me wrong I think the R2D2 priors are great, but they won’t automatically fix a hidden pathology in the model.

I also notice your spline structure is the following:

I am struggling to understand what you are doing here - why include the spline on yearsFromStart and then again on yearsFromStart by ats_baseline? I’m wondering if this could give you an identifiability problem. Again the pairs plot might be informative here.

6

Thank you @jroon I think you are right that the prior alone will not fix this model. I will examine the pairs plots. But in the meantime, do you think the spline on the s(yearsFromStart) term unnecessary, given the inclusion of the spline in the interaction term s(yearsFromStart, ats_baseline) ? Perhaps should I try leaving the main effect term yearsFromStart without the spline, just as a main effect?

I should point out that a version of the model bf(mvbind(cann28, alc28, heroin28) | trials(set) ~ s(yearsFromStart, by = ats_baseline) + (1|p|encID)) while not without some problems (some observations with bad pareto-k, semi-resolved after re-looing) at least converges with healthy ESS and with good diagnostics (traceplots, histograms of posterior predictive replicates, and pit-ecdf calibration plots). So this model is just a further exploration.

7

Ok so a few things to unpack there.

First s(yearsFromStart, by = ats_baseline) and s(yearsFromStart, ats_baseline) are not the same thing. If you leave out by = you are smoothing across both variables at once (as best I can tell). If you use by it does different things depending on varaible type. Best to look at the helpfile for ?mgcv::s() for the details (this is what brms is using under the hood to set up the splines):

Screenshot 2025-10-15 at 18.03.59
Screenshot 2025-10-15 at 18.03.591682×700 161 KB

The bits in red might partly explain the issues you are having by the way and why I recommend to look at the pairs files for relevant variables.

On another front - judging by the variable name yearsFromStart you have used the spline to transform the time variable. I actually wouldnt’ do this. Mathematically, it’s fine of course…. but it’s hard to interpret and hard to explain the results to non-statisticians - how do you explain spline transformed non-linear time to drug dependency clinicians for example? However, its typically easier to explain non-linear effects of some baseline characteristic or drug treatment over time to the subject matter experts. It also maybe easier to reason about when making modelling choices too. For example, does it make sense clinically forats_baseline (whatever that is) to have a varying effect over time with respect to your outcome? I would think about such things before comparing the model metrics. Just my 2 cents, I hope it’s helpful!

8

Thank you again for the feedback @jroon. This is a case of, for the sake of brevity, you only seeing the 10% of the iceberg above the surface (in this case that 10% is the post you replied to). Professor Vehtari has been giving me advice on the modelling process for almost a year and has explained to me the implications of the different interaction syntaxes for s(). What started as a gaussian model became binomial, then betabinomial, then went from three simple regressions to a multivariate model. I’m not actually convinced multivariate is conceptually a great idea for multiple drug classes. It’s quite rare for people in opioid treatment (our sample) to be using all three of these drugs at the same levels at the same time: people tend to use one or two at a time and not the others. The observations with the bad pareto-k levels tend to be these cases, where they might have, for example, used heroin and cannabis on 28 out of the previous 28 days but not used alcohol at all.

As to your point about clinical knowledge: I am not really a statistician, in the way Professor Vehtari is (few are). I am an addiction researcher who does a lot of stats so I have a lot of experience in clinical addiction research and work closely with clinicians. One of the reasons I am conducting this analysis is that there is not a lot of evidence from naturalistic data like this, about long-term trajectories of primary and secondary substance use among people being treated for opioid dependence. So while clinicians might have hunches about what happens there is little empirical data. But all are in agreement that, like all human behaviour, these relationships are never linear. That is why I am persisting with this very challenging project. What gets published will form justification for requesting a much larger data extract from our Ministry of Health, an extract large enough to make the inferences that can be drawn from the data much more certain. I agree with you that the coefficients from spline models are very hard to interpret, but I have scripts that make this possible, for example estimating the difference at various timepoints between people who used 0 days amphetamine at start of treatment and people who used 14 days (see here)

9

Hi @llewmills. Oh sorry I did not realise you had already gotten such advice. My comments about not wanting to explain non-linear time to clinicians are informed by the fact that I used to be one before going into research/stats and I am therefore very aware of the typical level of statistics they hold → I would not want to explain non-linear time variables to my former colleagues 😅, but I know they are used to the idea of drugs having non-linear effects over time.