Hello,
I’m new to Bayessian statistics and I would like to know how it is possible to test particular hypotheses in brms. The code of my model is the following:
fit <- brm(response ~ 0 + Intercept + contrast + group.IQ + (1 | subject), data = test, family = bernoulli("logit"), prior = prior, sample_prior = "yes")
Family: bernoulli
Links: mu = logit
Formula: response ~ 0 + Intercept + contrast + group.IQ + (1 | subject)
Data: test (Number of observations: 960)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~subject (Number of levels: 20)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.03 0.22 0.67 1.55 1.00 1012 1433
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.98 0.40 -1.77 -0.19 1.01 862 1376
contrasthead_heed 5.45 0.57 4.43 6.68 1.00 3055 2343
contrasthid_heed -0.49 0.30 -1.08 0.09 1.00 2541 2630
contrasthoard_hod 1.12 0.26 0.61 1.65 1.00 2345 2829
contrasthud_hard 0.72 0.26 0.21 1.23 1.00 2358 2748
contrastwhod_hood 0.00 0.27 -0.54 0.54 1.00 2406 2822
group.IQlow -0.74 0.49 -1.66 0.23 1.00 829 1192
I would like to test whether there are differences between the two levels of group.IQ, that is, group.IQlow vs groupIQhigh (for all contrasts).
Any suggestions on the code using the hypothesis function?
Would hyp1 <- c(hyp1 = " group.IQlow + Intercept = 0") be ok?
Best,
George
Assuming you have used R’s standard dummy coding for the factors, in order to test the difference between group.IQlow and group.IQhigh (the reference level) you just need to test "group.IQlow = 0" (or "group.IQlow < 0" if you had a directional hypothesis). There’s no interaction between contrast and group.IQ specified in the model, so this hypothesis test applies to all contrasts.
If you had contrast * group.IQ in the model formula, you could test the effect of group.IQ in each contrast accordingly:
c("group.IQlow = 0", # effect of group.IQ in reference contrast
"group.IQlow + contrasthead_heed:group.IQlow = 0", # effect of group.IQ in head_heed contrast
"group.IQlow + contrasthid_heed:group.IQlow = 0", # effect of group.IQ in hid_heed contrast
"group.IQlow + contrasthoard_hod:group.IQlow = 0", # effect of group.IQ in hoard_hod contrast
etc ...)
1 Like
Many thanks for your response! I have a couple of some other questions to ask (sorry for that, I’m a beginner)
- Should I use “Intercept” in hypothesis testing, e.g.,
Intercept + group.IQlow < 0?
- What would be the code for the examination of differences between contrasthead_heed and contrasthad_hard (the reference level)?
Intercept + group.IQlow < 0 would be testing whether the predicted response is less than zero for the group.IQlow group. Similarly Intercept < 0 would be testing whether the predicted response is below zero for the reference IQ group. That is why just group.IQlow tests the difference between the two groups (because Intercept + group.IQlow - Intercept = group.IQlow).
In your model, the difference between contrasthead_heed and contrasthad_hard would be tested with "contrasthead_heed > 0" (or "contrasthead_heed = 0"). This uses the same logic as for IQ: the model’s predicted response for reference IQ and reference contrast is simply Intercept; the predicted response for reference IQ and contrasthead_heed is given by Intercept + contrasthead_heed. So, subtracting one from the other, Intercept drops out. If you were comparing responses for lowIQ between the two contrasts, group.IQlow would drop out for the same reason (it is in both terms). And, of course, that would have to be the case because there is no interaction between contrast and group.IQ specified in the model.
1 Like
Dear Andy,
I’m wondering how can I find the main effect of contrast. Would the following code work?:
contrasthead_heed + contrasthid_heed + contrasthoard_hod + contrasthud_hard + contrastwhod_hood = 0
Best,
George
Hi George, there are no interactions (with “contrast”) in this model, so every population-level effect shown in the model summary is already a “main” effect. Adding the 5 terms together (as you have done) does not correspond to a main effect – if you divided by 5, this would give the average effect of all non-reference contrasts, which might be meaningful but is not a main effect of “contrast”.
Hi Andy, thanks for your prompt reply.
So, let’s assume that I include the interaction in the code:
fit <- brm(response ~ 0 + Intercept + contrast + group.IQ + contrast:group.IQ + (1 | subject), data = test, family = bernoulli("logit"), prior = prior, sample_prior = "yes")
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.90 0.40 -1.70 -0.09 1.00 1336 1674
contrasthead_heed 5.63 0.93 4.05 7.71 1.00 2465 2129
contrasthid_heed -0.27 0.36 -0.98 0.45 1.00 2379 2975
contrasthoard_hod 0.69 0.34 0.01 1.35 1.00 1979 2340
contrasthud_hard 0.49 0.34 -0.18 1.18 1.00 2239 2564
contrastwhod_hood 0.23 0.35 -0.44 0.90 1.00 2230 2644
group.IQlow -0.95 0.59 -2.12 0.18 1.00 1224 1591
contrasthead_heed:group.IQlow -0.05 1.09 -2.22 1.97 1.00 2469 2080
contrasthid_heed:group.IQlow -0.72 0.61 -1.95 0.45 1.00 2625 2979
contrasthoard_hod:group.IQlow 0.91 0.50 -0.03 1.88 1.00 2205 3067
contrasthud_hard:group.IQlow 0.50 0.50 -0.48 1.48 1.00 2309 2913
contrastwhod_hood:group.IQlow -0.58 0.56 -1.68 0.50 1.00 2506 2662
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
In this case, what code would show the main effect of contrast?
Best,
George
There are at least three ways to get a main effect when there’s an interaction. They are probably not identical but similar.
Firstly, you could just run the model without the interaction to get the main effects (i.e., your first model above) – this gives the main effect of contrast type across all participants (IQs). The model with the interaction gives you conditional effects (not main effects) – the effect of each type of contrast when groupIQ is either low or high.
Secondly, if you want just one model, you can sum code groupIQ instead of using dummy coding . Then the “raw” coefficient of each contrast type will be the main effect (midway between the effect for high groupIQ and low groupIQ).
Thirdly, you could use the marginaleffects package to calculate the average marginal effect of the contrast types.
There are complications in the latter two – interpreting the meanings of coefficients with sum-coded factors needs some care/thought – using the marginaleffects package obviously requires some additional learning. So there’s something to be said for the first approach so long as your sample of participants’ IQs corresponds reasonably well to your population of interest.
Thank you.
Let’s say I get the main effects as shown in the first model. Do these values (e.g., β, SE, CI, Rhat, etc.) in the summary(fit)really tell us something about the contribution of each variable? Do we usually use these stats when we report our results or should we use the results of hypothesis testing only?
I’m asking because I still thinking in terms of mixed-effects models where the summary(fit)is important since it gives us the p values.
Best,
George
I generally do a hypothesis test on all coefficients, or combinations of coefficients, of interest to the study. I generally prefer directional hypothesis tests, which can give the strength of evidence the effect is in the hypothesised direction.