Brms survival model with interactions between a binary variable and a spline function

I would like to fit a survival model including an interaction between a binary variable 0,1 and a spline function such as Y ~ treat + treat*ns(age,3), with age as continuous variable.
After that, I would like to apply the hypothesis function to estimate credibility intervals / posterior probabilities for a sequence of patients’ ages such as seq(30,80,10).
How could I do that? Can anybody tell me the correct syntax ?
Thanks.

I think you can at least interact spline terms with other things with a nonlinear model in brms. There’s a lot more flexibility there. It might require some extra preprocessing:

library(tidyverse)
library(brms)

df = tibble(y = rnorm(50),
            x = rnorm(50),
            b = rnorm(50) > 0)

make_stancode(bf(y ~ b * spline,
                 spline ~ s(x), nl = TRUE),
              data = df)

The rest of your question I do not know the answer to unfortunately.

2 Likes

Oh I guess you’re actually looking for two splines, so maybe the model would be:

make_stancode(bf(y ~ b * spline1 + (1 - b) * spline2,
                 spline1 ~ s(x),
                 spline2 ~ s(x), nl = TRUE),
              data = df)
1 Like

And the make_stancode thing there – just ignore that. Should be brm really. make_stancode is a way to generate Stan code from the brms model so you can roughly see if what you expect to happen is happening.

1 Like

Hi, thanks for the answer. That’s for the model. But the question is then how to estimate the treatment effect and credibility intervals for each age…

For example, if I fit this model… where treat and sex are binary variables 0.1…
fit<-brm_multiple(OS | cens(1-Die) ~ 1+ treat+treat*sex, data = imp,family = brmsfamily("cox"), iter=2000)
Then I can easily estimate the effect and the intervals in men and women with…
hypothesis(fit, “treat < 0”)
hypothesis(fit, “treat +treat:sex < 0”)

The problem is when I introduce the interaction with the spline function for age as continuous variable:
fit2<-brm_multiple(OS | cens(1-Die) ~ 1+treat+treat*ns(age,3), data = imp,family = brmsfamily("cox"), iter=2000)

The question is how to use the hypothesis function here to calculate the effect in each age range, for example I would like to plot something similar to the graph I show below, which I have elaborated easily with the R rms package.
Is this possible???

Thanks…

Oof, I’m not sure honestly. Sorry for the delay.

@paul.buerkner do you know?

1 Like

The interest is to evaluate effect-modifying factors in medicine. For example in the graph above I was trying to evaluate the effect of a type of chemotherapy depending on age. But it would also be interesting to contrast effects according to year of treatment, continuous values of a biomarker, etc. This is a very important topic in medicine to individualize the therapy.
This is very simple for frequentist models but I don’t know how to do it with brms.

You can use

fi <- fitted(..., newdata = <your data for which you evaluate predictions>, 
             summary = FALSE)

and then compute the right contrasts, that is, substract the right columns (predicted response means) from each other.

mgcv has the factor-by variable smooth construct for this. The definition would be:

bf(y ~ b + s(x, by = b))

for a standard factor b. The parametric term is required because the basis created by s(x, by = b)) is centred about the model intercept, and thus we need something in the model for the group means.

An alternative formulation is to use an ordered factor for b:

bf(y ~ b + s(x) + s(x, by = b))

the parametric b effect again models the group means, the s(x) term is the spline for the reference level, and s(x, by = b) would represent l smooth differences between the relationships for the reference level and the l remaining levels of the factor b. in the example here l=1 so we again only get two smooths, but the parameterization is different to the first example. This coding is then like the standard treatment contrast coding in an ANOVA/lm().

As long as b is coded as a factor or an ordered factor, this just works in brms, and it handles these by smooths very nicely in the plotting functions etc.

2 Likes

Thanks.
I’ve tried to do this, and there are two issues.
First of all, it doesn’t work for the Cox model. It produces error (Error in bhaz_basis_matrix(data_response$y, args = x$family$baseline) : is.list(args) is not TRUE).
However, it does work with parametric models. I have therefore used a Weibull distribution, although I would like to know if Cox can be applied somehow. In the case of the Weibull AFT, the question is also how to calculate the time ratios for therapeutic effects in subjects of 50 and 70 years old. For example:

        fit_p2<-brm(SGm | cens(1-Die) ~ 1+ treat*ns(age,3), data = (aga2),family = weibull(), iter=2000)
        newdata = data.frame(expand.grid(
          treat=c(0,1), 
          Age=c(50,70)) %>%. 
            mutate(contrasts= LETTERS[1:n()])
          
          )
        fi <- as.data.frame(fitted(fit_p2, newdata=newdata, scale="linear",
                     summary = FALSE))
        names(fi) <- c("A", "B", "C", "D")
        hypothesis(fi, c("B- A < 0", "D-C < 0"))

But I don’t know if this is correct to estimate the coefficients for treat effects at the age levels.
How can I repeat that with hazard ratios in a Cox model or at least in a parametric PH model?

Unfortunately the brms experts here tend to be very busy so we can’t help you quickly here. But please consider:

This might be a bug, please try to find a minimal reproducible example and submit it to https://github.com/paul-buerkner/brms/issues

Try to find a way to simplify the question and post it as a new topic (currently the thread is a bit long and this discourages people from reading through it and answering).

Sorry I can’t help you more.

1 Like

Fitted is currently not supported for cox models but it should yield a more informative error (and it will hopefull be supported in the future). Please also keep in mind that the cox model in brms is still experimental and not officially yet supported.

3 Likes

Thanks. Can you please tell me if the above code was correct for a Weibull model?

I think it is correct yes.

1 Like

Thanks