Using regression random effects as outcomes or predictors in other models in BRMS

Interfaces brms
, ,
1

Sorry if this has been posted before, been googling to try and find an answer.

Let’s say I have a model with a random effect (in this case, a random slope for a signal detection theory model). In this case, each subject (indexed by subno) has their own slope for cond.

Now let’s say i want to expand this model so that I can test how other predictors (e.g. age) predict the random participant slope for cond.

Or, perhaps, i would like to use the random slope for cond as a predictor for some other variable y.

 # devtools::install_github("cran/sdtalt") (not on CRAN)
# Using example data set
library(sdtalt)
library(tidyverse)
library(brms)
library(tidybayes)
data(confcontr)

dat_long = confcontr %>%
           mutate(cond = isold - 1/2)

m1 = brm(
  sayold  ~ 1 + cond + (0 + cond | subno),
  data   = dat_long,
  family = bernoulli(link = "probit"),
  warmup = 1000, 
  iter   = 3000, 
  chains = 4, 
  init  = "0", 
  cores  = 4,
  seed   = 11
)

What is commonly done, is to extract the samples from this model and calculate each participant’s mean/median cond score (pps_cond). Then, use these estimates in a second, seperate model (e.g., y ~ pps_cond) or (pps_cond ~ age).

However, this approach wouldn’t account for uncertainty in the pps_cond scores. Alternatively, I was thinking i could randomly draw pps_cond scores for each participant from the posterior, and use brms multiple imputation function to run lots of models across different random of pps_cond scores, but I’m sure there’s a more elegant solution.

I’ve been looking into brms latent variable syntax, and I’m a little confused how it could apply it in the above situation.

2

Did you ever find a satisfying answer to this? I’d like to do something similar.

3

I too am in a similar situation, wanting to combine an expectation from one (logistic) model with an expectation from a second model that is conditional on the first. Any success for you?

4

I’d say the problem is rather different if you want to implement predictors of the group-level parameters, or use the group-level parameters as predictors of something else. In the case where you want to implement predictors of the group-level parameters; you could view this as a nonlinear modelling problem where the intercept and slope are nonlinear parameters, and they are described by linear submodels that implement the between-subject variation and the predictors. This is a typical form of the problem as applied in pharmacometrics for example. It corresponds to the group-level values including both systematic and random variation.
I’m not sure if there is a way to specify this in the typical linear brms syntax, but you could do it with the nonlinear syntax.
Here’s a toy example (excuse the potentially poor priors, they would need to be set with care in the real case):

example_nlin_data <- data.frame(ID = rep(seq(1,20), each = 10), x1 = runif(200,-5,5), x2 = rep(runif(10,-5,5), each = 20))
example_nlin_data$ID_intercept <- rep(rnorm(20,0,1), each = 10)
example_nlin_data$ID_slope <- 1+(example_nlin_data$x2*2)
example_nlin_data$y1 <- example_nlin_data$ID_intercept+(example_nlin_data$x1*example_nlin_data$ID_slope) + rnorm(200,0,5)

ggplot(data = example_nlin_data, aes(x = x1, y = y1)) + geom_point(alpha = 0.6) + facet_wrap(~ID, ncol = 4) + theme_bw()

example_nlin_mod1 <- brm(bf(y1 ~ nlinintercept + (x1*nlinslope), nlinintercept ~ 1 + (1|ID), nlinslope ~ 1 + x2 + (1|ID), nl = TRUE), data = example_nlin_data, control = list(adapt_delta = 0.95), prior = set_prior('normal(0,5)', class = 'b', nlpar = 'nlinintercept') + set_prior('normal(0,5)', class = 'b', nlpar = 'nlinslope'))
summary(example_nlin_mod1)

So in this example the subject-level slope in the linear relationship between y1 and x1 is itself dependent on a linear function of the predictor x2 (which has one value per subject), with intercept=1 and slope=2.
Those are reasonably well-recovered:

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: y1 ~ nlinintercept + (x1 * nlinslope) 
         nlinintercept ~ 1 + (1 | ID)
         nlinslope ~ 1 + x2 + (1 | ID)
   Data: example_nlin_data (Number of observations: 200) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~ID (Number of levels: 20) 
                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(nlinintercept_Intercept)     0.78      0.50     0.04     1.88 1.00     1604     1989
sd(nlinslope_Intercept)         0.19      0.14     0.01     0.52 1.00     1527     2055

Population-Level Effects: 
                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nlinintercept_Intercept     0.15      0.40    -0.62     0.93 1.00     5021     2964
nlinslope_Intercept         1.05      0.13     0.79     1.32 1.00     5581     2738
nlinslope_x2                2.00      0.05     1.91     2.09 1.00     5112     2912

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     4.96      0.27     4.48     5.53 1.00     6341     2784

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

For the case where you want to use the subject-level parameters as predictors for something else, I think that could be viewed as a joint modeling problem (might not be very workable in brms, but readily implemented in Stan).

1 Like