I am trying to understand how to estimate survival curves using brms weibull regression. Using the kidney example, I fit the following model:
fit1 <- brm(time | cens(censored) ~ age + sex + disease,
data = kidney, family = weibull, inits = “0”)
Based on my understanding the brms parameterization, I need to transform the log-link linear model into a scale parameter with the following: scale = exp( linear predictors ) / gamma( 1 + 1/shape). brms provides the shape parameter.
The survival curve is then: S(t) = exp(-scale * t^shape). However, my results look totally wonky when I try to plot this curve. What am I doing wrong?
One follow-up question. I have data that are either interval censored or right censored (the event either happened sometime during the individual’s history, which varies in length, or has not yet happened). Based on my reading of the brms documentation, I think I should code this as:
censored <- ifelse(y == 1, 2, 1)
lower_bound <- ifelse(y == 1, 0.001, time) # 0.001 indicates that the lower bound of the interval is near-instantaneous, time is a vector of durations for each observation. If data are right censored then the lower bound is elapsed time.
Yes, the way you specify censoring looks correct to me, provided that actually non of your observations are not censored.
The parameterization of the weibull distribution is documented in vignette("brms_families"). Perhaps this helps you in producing the correct survival function.
For anyone else interested, I believe the following produces the correct survival function for the Weibull distribution, as parameterized in brms. We need parameters lamba and k (the shape parameter estimated in brms). The following R code converts the mean value to lambda: