Hypothesis test brms

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

1. Should I use “Intercept” in hypothesis testing, e.g., `Intercept + group.IQlow < 0`?
2. 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

Thank you so much!

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

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.