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…
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 Issues · paul-buerkner/brms · GitHub
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
A question about this.
Can you update the status of this?
At present the implementation of the Cox model is mature enough to be used in papers?
We already have the complete analysis done but I don’t know if I should repeat it whole with another model.
1 Like
I would hope it is a reasonable implementation of the cox model but brms of course comes without any warranty so you always use it at your.own risk.
1 Like
Thanks for the comment.
Our article is already written.
As the author of the package, I simply wanted to ask for your advice. Do you advise us to repeat everything with an AFT or should we leave it as it is and submit the paper?
I used brms because I don’t know any other software that has Bayesian Cox model, with multiple imputation, and non-linear effects (natural splines) of the variables.
The alternative is to repeat the whole analysis with an AFT model.
What would you do here?
I would leave it as is and submit the paper. :-)
1 Like
The tedious but safe way would be to run SBC on the model. But that might not be the best investment of your time - depends on how much you are willing to “pay” to further reduce the particular risk of the model being problematic.
Best of luck with your paper!
1 Like