How to code post-hoc comparisons with interactions when using brms and the "hypothesis" function

Interfaces brms
1

I’m confused about how to correctly code comparisons between specific factor levels when interactions are involved.

There is a good description of how to run comparisons with the hypothesis function here without interations, where you simply add the levels of interest to the intercept.

In this conversation, @paul.buerkner describes how to do this when an interaction is involved, by adding both the interactions and main effects to the intercept. But when I try and apply this to my example I get strange results. I must have interpreted Paul’s example incorrectly, and am not doing this correctly. Here is my example:

My model is a multilevel proportional binomial model of seedling survival, with a set of predictors that includes a couple of interactions. This is the model output:

 Family: binomial 
  Links: mu = logit 
Formula: seedling.count | trials(seedling.count.old) ~ Fertiliser + Shelter + Grass2 + mo(pasture.cover) + light.st + (1 | Block) + Shelter:mo(pasture.cover) + Shelter:Grass2 
   Data: d2 (Number of observations: 281) 
Samples: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
         total post-warmup samples = 12000

Group-Level Effects: 
~Block (Number of levels: 20) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.29      0.53     1.48     3.52 1.00     3033     5253

Population-Level Effects: 
                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                     -4.14      0.97    -6.11    -2.25 1.00     3111     5605
FertiliserF_2                 -0.25      0.15    -0.55     0.04 1.00    17872     9622
ShelterS_2                     1.76      0.72     0.38     3.18 1.00     4992     7110
Grass2G_1b                    -0.15      0.38    -0.90     0.58 1.00    10311     9896
Grass2G_2a                    -0.08      0.59    -1.27     1.06 1.00     9308     8573
Grass2G_2b                    -0.37      0.63    -1.64     0.83 1.00     8163     8270
Grass2G_3a                     0.26      0.61    -0.97     1.45 1.00     9432     8791
Grass2G_3b                     0.26      0.64    -1.02     1.49 1.00    10170     8649
light.st                       0.45      0.10     0.26     0.64 1.00    16633     9492
ShelterS_2:Grass2G_1b          0.77      0.55    -0.32     1.86 1.00     8568     8489
ShelterS_2:Grass2G_2a          1.24      0.64    -0.00     2.52 1.00     8634     8774
ShelterS_2:Grass2G_2b          0.67      0.65    -0.60     1.96 1.00     7991     8489
ShelterS_2:Grass2G_3a          1.46      0.65     0.20     2.73 1.00     8973     9042
ShelterS_2:Grass2G_3b          0.65      0.68    -0.67     1.99 1.00     9366     9585
mopasture.cover                0.11      0.14    -0.17     0.38 1.00     4939     7238
mopasture.cover:ShelterS_2    -0.44      0.13    -0.70    -0.18 1.00     4893     6971

Simplex Parameters: 
                               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
mopasture.cover1[1]                0.22      0.17     0.01     0.60 1.00     9145     6511
mopasture.cover1[2]                0.16      0.14     0.00     0.51 1.00    13373     6151
mopasture.cover1[3]                0.15      0.13     0.00     0.47 1.00    13097     6842
mopasture.cover1[4]                0.17      0.14     0.01     0.51 1.00    13337     6790
mopasture.cover1[5]                0.16      0.14     0.00     0.51 1.00    13833     7836
mopasture.cover1[6]                0.15      0.13     0.00     0.48 1.00    12070     8217
mopasture.cover:ShelterS_21[1]     0.06      0.06     0.00     0.21 1.00    12720     6682
mopasture.cover:ShelterS_21[2]     0.09      0.07     0.00     0.27 1.00    11174     5866
mopasture.cover:ShelterS_21[3]     0.17      0.10     0.01     0.39 1.00     8816     5009
mopasture.cover:ShelterS_21[4]     0.17      0.11     0.01     0.43 1.00    10516     5837
mopasture.cover:ShelterS_21[5]     0.29      0.16     0.02     0.61 1.00     8361     6314
mopasture.cover:ShelterS_21[6]     0.23      0.11     0.03     0.47 1.00    12202     6691

Samples were drawn 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 interested in examining various contrasts of the Shelter:Grass2 interaction.

  • Shelter has two levels: S_1 (reference), S_2
  • Grass2 has 6 levels: G_1a (reference), G_1b, G_2a, G_2b, G_3a, G_3b

Here is the conditional_effects plot for the Shelter:Grass2 interaction:

image
image736×454 5.01 KB

For the purposes of this example, these are the comparisons I want to make:

    1. Compare G_3a and G_3b (pink and dark blue) within Shelter treatment S_2
    1. Compare S_1 and S_2 within Grass2 treatment G_3b (pink)

Here is the code I’ve used for Comparison 1 and the output. The estimate is negative, which indicates that the G_3b treatment has lower seedling survival than the G_3a treatment (and this matches what is seen in the conditional effects plot) :

hypothesis(seedling2_R5_a, "Intercept + ShelterS_2:Grass2G_3b + ShelterS_2 + Grass2G_3b = Intercept + ShelterS_2:Grass2G_3a + ShelterS_2 + Grass2G_3a")

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (Intercept+Shelte... = 0    -0.81      0.61    -2.04     0.32         NA        NA

Comparison 2, however, gives a positive estimate, indicating that survival is higher for S_2 than S_1. But the conditional effects plot shows the opposite.

hypothesis(seedling2_R5_a, "Intercept + ShelterS_2:Grass2G_3b + ShelterS_2 + Grass2G_3b = Intercept + Grass2G_3b")

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (Intercept+Shelte... = 0     2.41      0.78     0.94     4.02         NA        NA    *

Please could I have some help understanding what has gone wrong with my comparisons, and how I can code them correctly.

Thanks very much

  • Operating System: Windows
  • brms Version: 2.15.0
2

I don’t fully understand what the impact of monotonic predictors are in this scenario. But is the problem here caused because – in the second contrast – you are comparing ShelterS_1 with ShelterS_2 but not including the mopasture.cover:ShelterS_2 interaction?

3

Aha, thanks @andymilne, that must be it. I’ve re-run the model without the monotonic interaction, and now all the comparisons make sense.

Although I still can’t figure out how to adjust the code to include the monotonic effect, because no matter what format I try an error comes up saying that the parameter cannot be found. I must be coding it wrong but I feel like I’ve tried every possible way that mopasture.cover is coded in the output.

There is one other thing I don’t understand. Why is it necessary to include the main effects along with the interaction term? When I try to re-create the numbers for the summary table, I only need to include the interaction term. For example, in the summary table of my original question, the comparison between ShelterS_2:Grass2G_3a and the Intercept (which is ShelterS_1:Grass2G_1a) has an estimate of 1.46 with CIs 0.2 and 2.73. This code gives me the correct numbers: hypothesis(seedling2_R5_a, "Intercept + ShelterS_2:Grass2G_3b = Intercept"). When I add the main effects as well, like this, the numbers are wrong: hypothesis(seedling2_R5_a, "Intercept + ShelterS_2:Grass2G_3b + ShelterS_2 + Grass2G_3a = Intercept"). But only including the interaction doesn’t work when I’m trying Comparison 1 (from the original question), I definitely need to include the main effects to get numbers that look right.

It would be excellent if someone could explain the theory behind what the hypothesis function is doing and why both the main effects and interactions need to be included (and why that doesn’t work when doing a comparison from the summary output).

4

It seems this is a known bug in brms: Hypothesis with monotonic effects · Issue #778 · paul-buerkner/brms · GitHub There is a workaround suggested there but maybe easier, for now, just to treat pasture.cover as a normal (unordered) factor?

Those seem to be quite different comparisons and neither lines up with your original comparisons of interest. If there is no monotonic effect, your previous hypothesis tests correctly lined up with your two comparisons.

Note that, by cancelling out terms,
"Intercept + ShelterS_2 + Grass2G_3b + ShelterS_2:Grass2G_3b = Intercept + ShelterS_2 + Grass2G_3a + ShelterS_2:Grass2G_3a"
is the same as
"Grass2G_3b + ShelterS_2:Grass2G_3b = Grass2G_3a + ShelterS_2:Grass2G_3a",
and
"Intercept + ShelterS_2 + Grass2G_3b + ShelterS_2:Grass2G_3b = Intercept + Grass2G_3b"
is the same as
ShelterS_2 + ShelterS_2:Grass2G_3b = 0"

5

Ok thanks, good to know.

I don’t think I explained my question very well. I was trying to do a new third comparison, one where the ‘answer’ is already provided by the summary, to test whether I could get to the answer and therefore whether my coding method was correct. Originally I ran hypothesis(seedling2_R5_a, "Intercept + ShelterS_2:Grass2G_3b + ShelterS_2 + Grass2G_3a = Intercept") , which I thought was they same method I’d used to run my two comparisons from the original question (adding both the interaction and two main effects, like Paul says in the linked discussion). But this gave the wrong numbers, the correct numbers were achieved by hypothesis(seedling2_R5_a, "Intercept + ShelterS_2:Grass2G_3b = Intercept"), which just has the interaction added and not the two main effects. I’m just not fully understanding exactly what each parameter is doing. But your examples of cancelling out terms makes it slightly clearer to see what’s going on.