How to draw posterior samples with a different covariate value other than the intercept?

Interfaces brms
1

Use the dataset kidney that comes with brms as an example. Suppose I have the following model:

library(brms)
fit <- brm(time | cens(censored) ~ age + sex + (1+age||patient), data = kidney, family = "exponential")

When I draw posterior samples with

pe <- fixef(fit, summary = FALSE) # population-level samples
ge <- ranef(fit, summary = FALSE) # group-level samples

I assume that the resulting samples for sex (males and females) are associated with the variable age being at the intercept (e.g., 0 or mean). If I want to obtain posterior samples at a different value of age (e.g., 37.5), is there a way to achieve that?

2

Hi, check out the fitted() method. Try the newdata argument with the covariate value that you’d like to use. Something like

fitted(fit, newdata = data.frame(age = 37.5))
3

I do realize that both fitted and predict are indeed two inferential tools available in brms. However, these two functions are meant to make inferences or predictions for the response variable, not for an explanatory variable (e.g., sex in the example I gave above). And that is why I would like to draw posterior samples using fixef and ranef with a particular value of a covariate (e.g., age in this case).

4

Do you want to predict sex as a function of age? Or time as a function of sex with age = 37.5?

In the former case you might want to revise your model. In the latter case, try:

fitted(fit, newdata = data.frame(age = 37.5, sex = 0:1))
5

Thanks for the suggestion! I guess you meant

fitted(fit, newdata = data.frame(age = 37.5, sex = c('female', 'male')))

which gives

**Error in get(g, data) : object 'patient' not found**

Again, my interest is not about predictions at the individual level. For my particular dataset, I would like to draw posterior samples that are best captured by the following two lines except that I want to control age at a value other than the intercept:

pe <- fixef(fit, summary = FALSE) # population-level samples
ge <- ranef(fit, summary = FALSE) # group-level samples
6

“Fitted” values (linear predictions) for a new patient, using just the fixed effects.

fitted(
  fit, 
  newdata = data.frame(age = 37.5, sex = c('female', 'male'), patient = "new"), 
  re_formula = NA,
  allow_new_levels = TRUE
)

“Fitted” values (linear predictions) for a new patient, using the random effects too.

fitted(
  fit, 
  newdata = data.frame(age = 37.5, sex = c('female', 'male'), patient = "new"), 
  re_formula = NULL,
  allow_new_levels = TRUE
)
2 Likes
7

Yes, I agree if I want to make inference about a random patient. However, my goal is different: I want to draw posterior samples that are better captured by the following two lines except that I want to control age at a value other than the intercept:

pe <- fixef(fit, summary = FALSE) # population-level samples
ge <- ranef(fit, summary = FALSE) # group-level samples
8 
pe <- fixef(fit, summary = FALSE) # population-level samples

@tjmahr’s first answer is correct re that one.

For the group level samples, you want to get patient-level predictions for people with the counterfactual that they are 37.5 years old? You again want to use the fitted() method, setting all predictors to their desired values, borrowing from @tjmahr:

library(tidyverse)  # I don't remember where distinct() is from
my_predictor_values <- distinct(kidney, patient, sex)
my_predictor$age <- 37.5
fitted(
 fit, 
 newdata = my_predictor_values, 
 re_formula = NULL
)

I am not 100% sure on what you want, but I am near 100% sure that the answer is fitted().

3 Likes
9

It does seem the suggestions by you and @tjmahr are very close to what I’m looking for! By default predict() offers summarized results, which is not what I want. However, now I just realized that, if I set summary = FALSE in predict(), I may be able to achieve my goal.

Thanks a lot, @matti and @tjmahr!

10

Still struggling with this. First, we have the following model

library(brms)
fit <- brm(time | cens(censored) ~ age + sex + (1+age||patient), data = kidney, chains=4, iter = 2000, family = "exponential")

Let’s focus on the population effect of sex for the moment:

fixef(fit, summary = TRUE)

which gives

            Estimate  Est.Error       Q2.5      Q97.5
Intercept  3.827167116 0.63288949  2.6049586 5.09145851
age       -0.003162866 0.01220938 -0.0268914 0.02037655
sexfemale  1.453764900 0.42597116  0.6152186 2.30179523

If I interpret the above results correctly, the following two rows in the result above

           Estimate  Est.Error       Q2.5      Q97.5
Intercept  3.827167116 0.63288949  2.6049586 5.09145851
sexfemale  1.453764900 0.42597116  0.6152186 2.30179523

correspond to the male effect and the difference between female and male, respectively, while age is controlled at 0 (intercept). However, I couldn’t reconcile the above result with what predict renders:

fitted(fit, newdata = data.frame(age = 0, sex = c('female', 'male'), patient = "new"), re_formula = NA, allow_new_levels = TRUE)
      Estimate Est.Error     Q2.5    Q97.5
[1,] 228.28978 140.74064 69.64470 577.8507
[2,]  56.33436  40.93747 13.53066 162.6269

What am I missing?

Also, if I want to predict the population effect while sex is controlled at the average between the two levels, how would I set it up in predict. The following does not seem to work:

fitted(fit, newdata = data.frame(age = 0, sex = 'xx', patient = "new"), re_formula = NA, allow_new_levels = TRUE)

11

The model uses an exponential family. fitted() returns values on the response scale by default. You can use the response argument to control that.

fitted(
  fit, 
  newdata = data.frame(age = 0, sex = c('female', 'male'), patient = "new"), 
  re_formula = NA, 
  allow_new_levels = TRUE, 
  scale = "response"
)
#>       Estimate Est.Error     Q2.5    Q97.5
#> [1,] 226.51666 135.03088 70.22007 596.8289
#> [2,]  55.57147  39.66607 13.16078 161.0541

fitted(
  fit, 
  newdata = data.frame(age = 0, sex = c('female', 'male'), patient = "new"), 
  re_formula = NA, 
  allow_new_levels = TRUE, 
  scale = "linear"
)
#>      Estimate Est.Error     Q2.5    Q97.5
#> [1,] 5.276526 0.5359080 4.251634 6.391630
#> [2,] 3.814932 0.6328713 2.577241 5.081741
2 Likes
12

Perfect! You solved my first puzzle.

OK, I think I’ve found a solution for my second question: If I want to predict the population effect while sex is controlled at the average between the two levels, how would I set it up in predict?

First, draw posterior samples for the two sexes:

samples <- fitted(fit, newdata = data.frame(age = 0, sex = c('female', 'male'), patient = "new"), re_formula = NA, allow_new_levels = TRUE, summary=FALSE, scale = "linear")

Then make inferences about sex effects based on the above samples.