Adding a soft sum-to-zero constraint to to a brms model

Interfaces brms
, ,
1

I would like to add a soft sum-to-zero constraint to a brms model with group-specific effects. I think I should be able to do this through stanvars. The model formula for the model I would like to fit looks something like this:

Y ~ treat + (1+treat|diag_group) + (1|cell_id)

and I would like to put sum-to-zero constraints on the diag_group random effects. I have 3 questions about this process:

Question 1:
In looking at the Stan code generated by make_stancode(), the group effects get parameters z_1 and z_2. It appears that the ordering is alphabetical, regardless of the order in which the terms appear in the formula. Can you confirm that this is the case?

Question 2:
My other two questions I think are more Stan syntax than anything. Suppose I dropped the random slope on treat, so the model formula has (1|diag_group). In this case, the group-specific effects corresponding to diag_group are defined as:

  vector[N_2] z_2[M_2];  // standardized group-level effects

where M_2=1 (random intercept only). I tried defining stanvars as follows:

stanvar(scode = "sum(z_2) ~ normal(0, N_2*.001);",
          block = "model", position = "end")

The resulting Stan code looks like it’s doing what I want, but I get an error. I think this is because z_2 is a vector of (1-element) vectors. Am I right that this is the problem, and if so, is there an easy (1-line) way to express the sum that I want?

Question 3:
When I include the random slope, z_2 is defined as follows:

  matrix[M_2, N_2] z_2;           // standardized group-level effects

where in this case M_2=2. I would like to add two sum-to-zero constraints, one for the intercepts and one for the slopes. How can I express this?

Thanks,
Jonathan

  • Operating System: Windows
  • brms Version: 2.16.3
3 Likes
2

Update: I have been successful in adding the sum-to-zero constraint on the intercepts with the following code:

stanvar(scode = "sum(z_2[1]) ~ normal(0, N_2*.001);",
          block = "model", position = "end")

For groups that have a random slope in addition to the random intercept, the following works for the slopes:

stanvar(scode = "sum(z_2[2]) ~ normal(0, N_2*.001);",
          block = "model", position = "end")

When there is more than one group of random effects, you still need to know which index corresponds to which variable (i.e., z_1 vs. z_2 vs. z_3) - it appears that they are ordered alphabetically by the group name, but just to be sure I wrote a little function to compare the Stan data (generated from make_standata) to the original data frame to identify each group.

2 Likes