Post-hoc test brms

Interfaces brms
1

To begin, I just mention that I am totally new to bayesian statistic (and I don’t have a very solid statistical background).

I ran the following model

prior1 <- c(
  prior(normal(0, 50), class = b),
  prior(exponential(0.1), class = sd),
  prior(exponential(0.1), class = sigma))

BMvpa <- brm (
  RT ~ 1 + GroupC*PrimeC*CongC*EmoC*SexC
  + (1 + CongC *PrimeC*EmoC || ID) 
  + (1 + GroupC*PrimeC*CongC*EmoC*SexC || Target),
  data = df1,
  family = exgaussian(),
  prior = prior1,
  warmup = 2000, iter = 5000,
  chains = 3,
  cores= parallel::detectCores(),
  sample_prior = TRUE
)

Her is part of my output:

                            Predictors Estimate     SE   Lower   Upper  Rhat   BF01
1                                   Intercept  700.175 13.470 674.330 726.617 1.000   <NA>
2                                       Group   49.027 33.542 -17.261 115.122 1.000   0.51
3                                       Prime  -12.197  2.816 -17.799  -6.655 1.000  0.017
4                                  Congruency  -15.879  2.798 -21.507 -10.435 1.001 <0.001
5                                     Emotion   17.092  6.860   3.740  30.373 1.000   0.32
6                                         Sex    5.362 24.381 -42.347  52.871 1.001  2.031
...
21                    Group x Prime x Emotion   22.339  8.509   5.464  39.136 1.001   0.24
...

Each variable (Group, Sex, Prime, Congruency , Emotion) is dichotomous and coded -0.5 (TD, M, LSF, ICG, Joy), +0.5 (ASD, F, HSF, CG, Anger).

I would like to go more in details in the interaction Group x Prime x Emotion (represented on the figure that i posted here) and would like to know the posterior distribution regarding the effect of Prime for each group and each emotion.

I thought about 2 strategies.

1/ First using emmeans :

BMGPE_emm <- emmeans(BMvpa, ~ GroupC:PrimeC:EmoC)
BMGPE_fac <- update(BMGPE_emm, levels =list(GroupC= c("TD","ASD"), PrimeC= c("LSF","HSF"), EmoC=c("joy","anger")))
contRT1 <- as.data.frame(contrast(BMGPE_fac, method = "pairwise", by = c("GroupC","EmoC")))

output:

1 LSF - HSF     TD   joy  8.706978 -0.3446188  18.33661
2 LSF - HSF    ASD   joy 20.487280 10.2622944  30.46115
3 LSF - HSF     TD anger 15.029082  6.2702713  24.67623
4 LSF - HSF    ASD anger  4.412052 -5.6749680  14.60393

I am not sure about this because I would have expected only negative estimates (HSF primes reducing response time). Additionally, is there a possibility to compute a Bayes factor here?

2/ Second, using the hypothesis function (I read the blogpost of Matti Vuorre). I think it would be the best but the result I got are even more strange and I think I probably made a mistake (I was expected only negative estimates) :

> hypothesis(BMvpa,c(qAJ = "PrimeC + 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
  Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1        qAJ     3.77     17.35   -30.51    37.79       3.56      0.78     
---

> hypothesis(BMvpa,c(qAA = "PrimeC + 0.5*GroupC + 0.5*EmoC = 0"))
Hypothesis Tests for class b:
  Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1        qAA    20.86     17.46   -13.54    55.14       1.63      0.62     
---

> hypothesis(BMvpa,c(qTJ = "PrimeC - 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
  Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1        qTJ   -45.26     17.34   -78.93   -10.78       0.15      0.13    *
---

> hypothesis(BMvpa,c(qTA = "PrimeC - 0.5*GroupC - 0.5*EmoC = 0"))
Hypothesis Tests for class b:
  Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1        qTA   -45.26     17.34   -78.93   -10.78       0.15      0.13    *
---

I am sure I am wrong on what I did but I don’t manage to find what is wrong (I check the coding of my variables but I think they are right).

Could somebody help me?

(I add a small part of my data)

df1.csv (598.1 KB)

2

Hi @AdeLac74 and welcome to the Stan forum.

You can get the simple effects of Prime for each combination of Group and Emotion with both emmeans and hypothesis(). You just need to change your emmeans call a little bit, and make sure the hypothesis strings are coded correctly. The example below assumes that you are using the default effect coding of factor variables Group, Emotion and Prime. If you use numeric contrast coding, emmeans will probably behave differently, and you will need to adjust the hypothesis strings to include the appropriate multipliers (-.5 etc.). I recommend including categorical variables as R factors (either with contrast coding (see ?contrasts; e.g. contrasts(data$variable) <- matrix(c(-.5, .5))) or default effect coding so that the functions (brm(), emmeans(), hypothesis()) can deal with them appropriately.

library(knitr)
library(brms)
library(emmeans)
library(tidyverse)

dat <- read_csv("df1.csv")

dat <- select(dat, RT, Group, Prime, Emotion)
# Group only has 1 level
set.seed(1)
dat$Group <- sample(c("TD", "ASD"), size = nrow(dat), replace = TRUE)

fit <- brm(
  RT ~ Group*Prime*Emotion, 
  data = dat, 
  file = "GroupPrimeEmotion"
  )

summary(fit)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: RT ~ Group * Prime * Emotion 
##    Data: dat (Number of observations: 3000) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                     774.93     15.65   743.76   805.03 1.00     1676
## GroupTD                         9.43     22.29   -33.69    52.89 1.00     1414
## PrimeLSF                       -3.96     22.20   -46.00    39.39 1.00     1296
## Emotionjoy                    -24.22     22.18   -67.89    19.56 1.00     1569
## GroupTD:PrimeLSF                7.58     31.11   -53.97    68.98 1.00     1037
## GroupTD:Emotionjoy            -17.91     31.53   -80.79    43.09 1.00     1378
## PrimeLSF:Emotionjoy            31.85     31.29   -30.73    91.44 1.00     1409
## GroupTD:PrimeLSF:Emotionjoy   -12.80     44.73   -98.84    75.24 1.00     1162
##                             Tail_ESS
## Intercept                       2774
## GroupTD                         2126
## PrimeLSF                        2394
## Emotionjoy                      2447
## GroupTD:PrimeLSF                2244
## GroupTD:Emotionjoy              2118
## PrimeLSF:Emotionjoy             2503
## GroupTD:PrimeLSF:Emotionjoy     2325
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma   302.62      3.89   295.03   310.06 1.00     4095     2186
## 
## 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).

For some reason, emmeans gives the opposite signs (note the direction of the subtractions).

ems <- emmeans(fit, ~Prime, by = c("Group", "Emotion"))
summary(pairs(ems), point.est = mean)

## Group = ASD, Emotion = anger:
##  contrast  estimate lower.HPD upper.HPD
##  HSF - LSF     3.96     -39.5      45.7
## 
## Group = TD, Emotion = anger:
##  contrast  estimate lower.HPD upper.HPD
##  HSF - LSF    -3.62     -48.4      38.2
## 
## Group = ASD, Emotion = joy:
##  contrast  estimate lower.HPD upper.HPD
##  HSF - LSF   -27.89     -71.9      12.5
## 
## Group = TD, Emotion = joy:
##  contrast  estimate lower.HPD upper.HPD
##  HSF - LSF   -22.67     -66.8      18.4
## 
## Point estimate displayed: mean 
## HPD interval probability: 0.95

With hypothesis:

h <- c(
  "Prime; Group=ASD; Emotion=anger" = "PrimeLSF = 0",
  "Prime; Group=ASD; Emotion=joy" = "PrimeLSF + PrimeLSF:Emotionjoy = 0",
  "Prime; Group=TD; Emotion=anger" = "PrimeLSF + GroupTD:PrimeLSF = 0",
  "Prime; Group=TD; Emotion=joy" = "PrimeLSF + GroupTD:PrimeLSF + PrimeLSF:Emotionjoy + GroupTD:PrimeLSF:Emotionjoy = 0"
)
hypothesis(fit, h)

## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 Prime; Group=ASD;...nger    -3.96     22.20   -46.00    39.39         NA
## 2 Prime; Group=ASD;...=joy    27.89     21.59   -14.69    69.98         NA
## 3 Prime; Group=TD; ...nger     3.62     22.12   -39.25    47.70         NA
## 4 Prime; Group=TD; ...=joy    22.67     21.96   -20.22    65.06         NA
##   Post.Prob Star
## 1        NA     
## 2        NA     
## 3        NA     
## 4        NA     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
3

Hi Matti,

First, I thank you very much for your quick and detailed answer which helps me to understand better, but I am still struggling, probably because of the numeric contrast coding I used… (I recoded my variables manually, which might be inappropriate?)
Each variable (Group, Sex, Prime, Congruency , Emotion) is dichotomous and coded -0.5 (TD, LSF, Joy), +0.5 (ASD,HSF, Anger) and I really struggle to adjust the hypothesis string to find similar results than those with emmeans…
For instance, what would be the correct writting for this hypothesis:
“Prime; Group=ASD; Emotion=joy” = “PrimeLSF + PrimeLSF:Emotionjoy = 0” given my contrast code?

(I wonder if emmeans gives the opposite signs because it takes the alphabetical order in your example? )

4

If you use numerically coded factors, useful functions such as conditional_effects() won’t do their job properly. However there is nothing wrong in manually coding factors and it can sometimes be very useful. You just have to know exactly how to multiply parameters to get the effects you need.

To write that hypothesis with your contrast codes, you’d do something like

"Prime + Prime:Emotion*0.5 = 0"

To get the effect of prime when emotion is 0.5 (or whatever “joy” was).

You can also reverse the emmeans comparisons: https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html. It looks like emmeans does a lot of useful stuff for figuring out complex interactions so probably you can just go with that, instead of hypothesis(). It also helps to keep a plot open showing the conditional effects of your model (e.g. conditional_effects()) while you try to figure out the interactions, so you can visually confirm that you understand the resulting coefficients.

5

Thank you again Matti for all your very usefull advices! :-)
I will look at this today!
Cheers

6

Hi, matti, your answer was very useful, but I still have a question that, if Prime was encoded to (-.5, .5) and emotion also encoded to (-.5, .5). And you enter hypothesis (“Prime + Prime:Emotion*0.5 = 0”), how does this function know 0.5 is refer to, because I think that 0.5 can refer to prime like primeLSF and 0.5 can also can refer to emotion like joy. It is very thankful if you could help me figure it out ~