Hello everyone,
I am using brms in R to check the effect of treatment and a group on some result (y). In my data (dat), there are 3 categorical variables and a continuous response variable: 1 id 1 treatment variable with 3 categories: ctrl, treat1, treat2 1 group variable with categories a, b, c And a continuous variable y.
I fit a mixed model:
library(brms)
Mod=brm(y ~ Treatment*group +(1 | id), data=dat).
At first, i want to compare the change from control for the treat1, and treat2. (treat1 vs ctrl, treat2 vs ctrl).
At second, for each treatment (treat1 and treat2) i want to compare the change from control between the groups a and b, b and c, a and c.
I want to set these contrasts in R with the hypothesis function but don’t know how to do it.
Thank you in advance for your help
1 Like
Hi, sorry your question fell through a bit.
The main thing to understand is that brms uses Dummy coding to encode factors and their interactions (you can see how the resulting model matrix - except for varying intercepts - looks like by calling model.matrix(y ~ Treatment*group, data = dat). This however puts is in a problematic situation when we model Treatment * group but want to get some sort of “overall treatment effect”: there is no such value in the model, because the model assumes a (almost) completely separate treatment effect for each group.
IMHO the only way to capture some notion of “overall treatment effect” is to assume relative frequencies of the individual groups and then compute a weighted mean of the treatment effects for individual groups.
The treatment effects for individual groups are easier (those are directly represented in the model). So to get the contrast between say treatment 1 vs. control in group b we would first look how the means for both are computed:
E(control in group b) = Intercept + b_groupb
E(treatment 1 in group b) = Intercept + b_groupb + b_Treatmenttreat1 + b_Treatmenttreat1:groupb
E(treatment 1 in group b) - E(control in group b) = b_Treatmenttreat1 + b_Treatmenttreat1:groupb
So we can run hypothesis(Mod, "Treatmenttreat1 + Treatmenttreat1:groupb > 0", class = "b") to compare those.
Does that make sense? Can you extend this to the other comparisons you care about? Feel free to ask - this stuff can be daunting. Or you can post your code for doing the other comparisons for us to review, so that you can get feedback whether you got to the bottom of it :-).
If you are unsure about the names of the corresponding parameters you can always use parnames(Mod) to see all parameters the model has.
Best of luck with your modelling work!
2 Likes