Assessing influence of population-level predictors on one specific trial

Interfaces brms
,
1

I have an experimental design with one between-subjects factor (gender, 2 levels) and one within-subjects factor (cond, 2 levels). Each subject encounters multiple trials. The data is ordinal and I am using brms to construct a linear mixed model with a probit cumulative function.

In addition to assessing the influence of the group and condition overall, I would like to assess their influence on the behaviour on one specific trial. Instead of constructing a second model using the dataset only including this trial, I wanted to use the full model and then use posterior_linpred() to assess the posterior distributions in a similar manner as hypothesis(). I am computing the posterior draws for the group-level effect trial only, including the slopes and intercept, but dropping the group-level effect subject. Does this approach make sense or am I missing something?

Here is a minimal example using data from the psych package.

# loading libraries
library(brms)
library(tidyverse)

# using parallel processing
options(mc.cores = parallel::detectCores())

# get some data
df = psych::bfi %>% 
  drop_na() %>% 
  slice_sample(n = 400) %>% 
  mutate(subID = 1:n()) %>% 
  select(subID, gender, N1:N5, A1:A5) %>% 
  pivot_longer(c(N1:N5, A1:A5)) %>%
  mutate(
    gen   = if_else(gender == 2, "female", "male"),
    cond  = substr(name, 1, 1),
    trial = as.numeric(substr(name, 2, 2))
  ) %>% mutate_if(is.character, as.factor)

# set the contrasts 
contrasts(df$gen)  = contr.sum(2)
contrasts(df$cond) = contr.sum(2)

# run the model
m = brm(value ~ gen * cond + ( cond | subID ) + ( gen * cond | trial ), 
        data = df, seed = 4444, 
        prior = c(prior(normal(-0.97, 1), class = Intercept, coef = 1),
                  prior(normal(-0.43, 1), class = Intercept, coef = 2),
                  prior(normal( 0.00, 1), class = Intercept, coef = 3),
                  prior(normal( 0.43, 1), class = Intercept, coef = 4),
                  prior(normal( 0.97, 1), class = Intercept, coef = 5),
                  prior(normal( 0.00, 1), class = b),
                  prior(exponential(1),      class = sd),
                  prior(lkj(2),              class = cor)),
        family = cumulative("probit"))

# compute posterior draws of linear predictor
mt.post = posterior_linpred(m, newdata = df %>% select(gen, cond, trial) %>% distinct(), 
                             re_formula = ~ ( gen * cond | trial ) )
colnames(mt.post) = df %>% select(gen, cond, trial) %>% distinct() %>% 
  transmute(comb = paste(gen, cond, trial, sep = '_')) %>% pull()

# summarise posterior distributions
as.data.frame(mt.post) %>% 
  mutate(
    h.gender      = (female_N_5 + female_A_5)/2 - (male_N_5 + male_A_5)/2, 
    h.condition   = (female_A_5 + male_A_5)/2 - (female_N_5 + male_N_5)/2, 
    h.interaction = (female_N_5 - female_A_5) - (male_N_5 - male_A_5) 
  ) %>% select(starts_with("h")) %>%
  pivot_longer(cols = everything(), names_to = "Effect") %>%
  group_by(Effect) %>%
  summarise(
    Estimate = mean(value), 
    CI.Lower = quantile(value, 0.025, names = F), 
    CI.Upper = quantile(value, 0.975, names = F),
    Post.Prob = mean(value > 0)*100
  )
  • Operating System: macOS Tahoe 26.5.1
  • brms Version: ‘2.23.0’
2

Hi,
this looks roughly correct to me. The only potential “problem” is that you don’t just assess the effect of one specific trial, but you allow the group and condition to have different effect in each individual trial + you share a bit of information across all trials due to the random effect structure. This is desirable if you believe trials do generally matter and each trial is unique, you just for some reasone are interested in the effects of one trial in particular.

If you instead believe there is just one trial that is special and all of the other trials are identical, you may be better served by explicitly coding something like:

mutate(
   ....,
   is_trial_5 = as.numeric(trial == 5) 
)
...

m = brm(value ~ gen * cond * is_trial_5 + ( cond | subID )

Your code could also likely be a little bit more elegant using tidybayes::add_linpred_draws (but the code looks good as is)

Hope that makes sense and best of luck modelling.