Brms: How to keep it consistent for mu and shape in posterior_linpred (Weibull time-to-event)?

Interfaces brms
1

For models with a second parameter (e.g. shape for Weibull time-to-event), for which we do distributional regression (and where we have random effects on both the main mean model and the distributional model), how do we keep things consistent when predicting with posterior_linpred? What I’m interested in doing is looking at the probability at time T=t for time to event predicted for a new random effects level. For that, I want the linear predictor for both the mean and the shape.

What I worry about is that one can trade-off mean vs. shape to some degree to fit the data, so unless I get the linear predictors for the same sampled random effects on these, they would not “fit together”. I assume that posterior_linpred preserves the order of MCMC samples, so I should be fine on the fixed effects, but I’m not sure about the random sampling of the random effects. I.e. how can I ensure the same (correlated) random effects are draw when I separately request them for “mu” and “sigma” (posterior_linpred does not seem to allow me to ask for them at the same time, so I have to call the function twice) and then paste the samples together (as shown below)?

PS (updates): The case below is substantially simplified, but this really makes a difference in practice, because mixing mean and shape parameters that “don’t belong together” (as far as the random effects “are concerned”), leads to a much wider distribution than is truly the case and could via the non-linear link function even lead to wrong estimates.

library(tidyverse)
library(brms)

set.seed(1234)

# Simulate example data
# There's 10 observations per subject (subject ID = subject)
# event = 1 (event observed, status="none"), 0 (censored, status="right")
# time = event time or censoring time

mydata <- expand_grid(subject=1L:10L, repeatid=1L:10L) %>%
  mutate(subject_effect = rnorm(n=10)[subject],
         x = rexp(n=n(), rate=exp(1+subject_effect)),
         c = rexp(n=n(), rate=1),
        event = 1L*(x<=c),
        status=ifelse(event==1L, "none", "right"),
        time = pmin(x,c))

# Fit brms model with intercept and random subject effect
brmfit1 <- brm( formula = bf( time | cens(status) ~ 1 + (1|b|subject), 
                             shape ~ 1 + (1|b|subject) ), 
    family = weibull(),
    data = mydata,
    prior = prior(class=Intercept, normal(1, 2.5)) +
        prior(class=sd, normal(0,2.5)) +
        prior(class=Intercept, normal(0, 1), dpar=shape) +
        prior(class=sd, normal(0,0.5), dpar=shape)) 

# This is the tricky bit: I predict for a hypothetical new subject (subject=11)
# When the random effects for that new subject are drawn, how do I correctly preserve their correlaiton
tibble(
    # First predict mean parameter
    logmu = posterior_linpred(object=brmfit1,
                              newdata=tibble(subject=11L),
                              re_formula=NULL,
                              dpar="mu",
                              allow_new_levels=T, 
                              sample_new_levels="gaussian")[,1],
    # Second predict shape parameter
    logshape = posterior_linpred(object=brmfit1,
                                 newdata=tibble(subject=11L),
                                 re_formula=NULL,
                                 dpar="shape",
                                 allow_new_levels=T, 
                                 sample_new_levels="gaussian")[,1]) %>%
    mutate(prob_event_before_t_1 = pweibull(q=1,
                                            shape=exp(logshape), 
                                            scale=exp(logmu - lgamma(1+1/exp(logshape))))) %>%
ggplot(aes(x=prob_event_before_t_1)) +
geom_density()
  • Operating System: Windows 10
  • brms Version: 2.18.5
2

One option could be to call prepare_predictions directly, which will then internally create the samples for the new random effects levels. Then, the resulting object can be passed to posterior_epred(., dpar = "<dpar>", scale = "linear").