Issue with passing variable to brms formula

Interfaces brms
1

I am fitting a large number of models with smooth terms. For each model, I set a different value of k and pass it as a variable to the s() term in the brms formula. I noticed that if the variable changes in the global environment, it will cause undesired behavior when predicting from the model. See errors below. Is there a better recommended way to pass variables to a brms formula to avoid this problem?

my_k <- 5

fit <- brm(
  mpg ~ s(hp, k = my_k),
  data = mtcars
)

predict(fit)

my_k <- 1

predict(fit)

my_k <- 1000

predict(fit)

The first call to predict(fit) works as intended. The second two calls to predict(fit), after the value of my_k has changed in the global environment, do not. The first one gives incorrect predictions with the warning:

Warning message:
In smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) :
  basis dimension, k, increased to minimum possible

The second gives this error:

Error in smooth.construct.tp.smooth.spec(object, dk$data, dk$knots) : 
  A term has fewer unique covariate combinations than specified maximum degrees of freedom

This is clearly because the “wrong” k value is being assumed when predicting from the fit. How can I pass a variable to the k argument of s() without it later depending on values in the environment?

Please also provide the following information in addition to your question:

  • Operating System: Windows 11
  • brms Version: 2.21.0
2

I think this is the solution:

myform <- paste('mpg ~ s(hp, k =', my_k, ')')
fit <- brm(bf(myform), data = mtcars)

If it is done this way, later calls to predict(fit) are unaffected by whatever value my_k currently has in the global environment.

1 Like
3

I think what was going wrong in your initial code is that even though the actual model fit uses the current value of my_k (5 in this case) to fit the model, one or more elements of the fit object returned by brm contain an unevaluated version of my_k. The predict function then uses an unevaluated version of my_k when setting up the predictions. If my_k has changed in the interim, a problem then occurs when my_k is eventually evaluated (and is no longer equal to 5) while creating the predictions.

To make this more concrete, run the following code:

my_k <- 5

fit <- brm(
  mpg ~ s(hp, k = my_k),
  data = mtcars
)

fit is a list with several elements (run names(fit) to see the names of these elements), one of which is called "formula". Let’s look what’s in the "formula" element:

fit$formula # you could also run: fit[["formula"]]

mpg ~ s(hp, k = my_k)

Note the my_k rather than 5 in the formula. The predict function (for brmsfit objects, which is what fit is) calls a few other functions, including brms:::prepare_predictions.brmsfit to get the terms that go into the model formula. This function uses fit$formula to get the terms in the model formula, and I think this is where things eventually go wrong if my_k has changed between when the model is fit and when the predictions are extracted.

Your solution works because it forces evaluation of my_k before the model is fit, so the value 5 appears in fit$formula, rather than the unevaluated my_k. For example:

my_k = 5
myform <- paste('mpg ~ s(hp, k =', my_k, ')')

myform
[1] "mpg ~ s(hp, k = 5 )"

fit2 <- brm(bf(myform), data = mtcars)
fit2$formula

mpg ~ s(hp, k = 5)

You could also create a helper function to reduce the amount of coding needed if you want to fit several different models.

fit_fnc = function(kval) {
  brm(paste('mpg ~ s(hp, k =', kval, ')'), data = mtcars)
}

The code below will return a list of three model fits, one for each of three k values (the set_names() is just so that each list element will be named for the value of k that was used for that model fit):

library(tidyverse) # For the map function

fits = c(3,5,7) %<% set_names() %>% map(fit_fnc)

map(fits, ~ .x$formula)

$`3`
mpg ~ s(hp, k = 3) 

$`5`
mpg ~ s(hp, k = 5) 

$`7`
mpg ~ s(hp, k = 7)
1 Like
4

Since k affects the parameters block, you might want to try the 'update" function from brms to update the model with alternate values of k

1 Like