How to properly compare interacting levels

Interfaces brms
,
3

Apologies for the lack of clarity. Let me try again, this time with an example. (I am simplifying here, ignoring the possibility to have random effects, some covariates etc.)

# uniformly generated means
m = runif(27, 3, 12)

# random values, random means, same std. dev.
y = list()
for (i in 1:27) {
    y[[i]] = rnorm(30, m[i], 0.8)
}
y = unlist(y)

# factorial design
F1 = c(rep('A', 30), rep('B', 30), rep('C', 30))
F2 = c('I', 'J', 'K')
F3 = c('M', 'N', 'O')
design_matrix = expand.grid(F1=F1, F2=F2, F3=F3)

# dataset
dat = cbind(design_matrix, y)

require(brms)

summary(brm_dat_model1 <- brm(y ~
    (F1 + F2) * F3,
    data = dat,
    chains=4, iter=4000, cores=4))

This gives the summary table (it could be different in another run, given random number generators):

# Population-Level Effects: 
#           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# Intercept     4.56      0.24     4.08     5.02 1.00     3436     5065
# F1B           0.84      0.26     0.33     1.34 1.00     4201     4812
# F1C           3.81      0.26     3.32     4.33 1.00     4398     5180
# F2J           0.71      0.26     0.19     1.23 1.00     3781     5689
# F2K           1.44      0.26     0.92     1.96 1.00     4085     5495
# F3N           4.44      0.33     3.80     5.09 1.00     3615     4724
# F3O           0.89      0.34     0.23     1.55 1.00     3707     5204
# F1B:F3N       0.54      0.37    -0.17     1.26 1.00     4262     4913
# F1C:F3N      -4.21      0.36    -4.93    -3.51 1.00     4586     5286
# F1B:F3O       0.98      0.37     0.27     1.68 1.00     4321     5194
# F1C:F3O      -1.32      0.37    -2.03    -0.59 1.00     4629     5702
# F2J:F3N      -2.87      0.37    -3.60    -2.15 1.00     4271     5376
# F2K:F3N      -4.81      0.37    -5.54    -4.09 1.00     4514     5679
# F2J:F3O      -1.24      0.37    -1.97    -0.51 1.00     4478     5708
# F2K:F3O       0.72      0.37    -0.00     1.43 1.00     4540     5332
# 
# Family Specific Parameters: 
#       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# sigma     1.73      0.04     1.65     1.82 1.00     8879     5554

Now, what causes pains to grasp and work through are various comparisons when levels and combination in the table are relative to the reference: F1==A; F2==I; F3==M.

I’ve found it almost funny (and creative too) when some students approached to this by running comparable models with Intercept removed. Given the example above, something like:

brm_dat_model2 <- brm(y ~
    0 + 
    F1 : F3 + 
    F2 : F3,
    ...)

brm_dat_model3 <- brm(y ~
    0 + 
    F2 : F3 + 
    F3 : F3,
    ...)

And then they apply hypothesis(). For example:

hy = c(hy1='F1A:F3M < F1A:F3N',
    hy2='F1B:F3M == F1C:F3O')
hypothesis(brm_dat_model2, hy)

etc. etc.

This is rather clumsy and inelegant, if not altogether wrong. Yet, I see the struggle and I’m running out of ideas on how to help… In addition to this, it appears that what hypothesis() tells you “contradicts” what you observe in the plot from conditional_effects(). In fact, is there a way to get the BF directly from the conditional_effects() table?