Fitting time-to-event data with Weibull Hazard using brm function

Good afternoon, colleagues!

I have recently been the user of R package ‘brms’.
I’m trying to use it to fit time-to-event data about the survival of patients.
I planned to use a Hazard function of this kind: the Weibull Hazard multiplied by the factor responsible for the effect of covariates.
My question is this: did I understand correctly that in this package this fitting in brm will be recorded as follows?

fit2 ← brm(time (time) | cens (CENSOR) ~ WHOSTATN + BASESLD + NLine,
data = df1, family = weibull(), inits = “0”,
control = list(adapt_delta = 0.95, max_treedepth = 12))

Is it enough, that the type of hazard is “contained” in the likelihood function?


  • Operating System: Windows 7
  • brms Version: 2.3.1

I don’t understand the time (time) specification. I would rather write:

fit2 <- brm(time | cens(CENSOR) ~ WHOSTATN + BASESLD + NLine,
            data = df1, family = weibull(), inits = “0”,
            control = list(adapt_delta = 0.95, max_treedepth = 12))

The rest seems sensible to me.

Thank you, Paul!
With such time specification it works well and reasonably. Compilation and calculations — all is ok!
But without this specification I obtained such error message:

Error: Family ‘weibull’ requires response greater than 0.

Maybe, such specification is the bad way, but without this trick compilation can not start.

Depending on the input, the time() transformation in R does something I am not sure you want to do. Check the output of time(df1$time).

Rather, you should take the error message of brms seriously. Perhaps there is something in the time variable which is incorrect or you have 0 in it, which cannot be handled by such models. Alternatively, you may need to transform time to a numeric variable first via as.numeric(time).

1 Like

Paul, you are absolutely right!
I didn’t think, that time = 0 can not be used.
Much grateful!

But I still have a question: why brm doesn’t give the results of estimating the two parameters in the Weibull risk function?
Since there is not only the shape parameter k, but also the “lambda” is a scale parameter:
(k/lambda)*(x/lambda)^(k-1).

Please see https://cran.r-project.org/web/packages/brms/vignettes/brms_families.html for the parameterization of the weibull family (and other families) in brms.

1 Like

Paul, this is precisely the parametrization of the Weibul distribution, to which I am accustomed. And in the summary we can see:

Family: weibull
Links: mu = log; shape = identity
Formula: time | cens(CENSOR) ~ AGE
Data: df1 (Number of observations: 250)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 14.99 3.37 8.92 22.22 2021 1.00
AGE -0.04 0.05 -0.14 0.06 2143 1.00

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
shape 0.28 0.03 0.22 0.35 1813 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

We can see value of shape parameter.
But, if I’m interested in looking at the resulting value of the scale parameter, how can it be derived?

The parameterization is in terms of mean and shape not in terms of scale and shape.

You can extract the mean fitted values via fitted(fit2). If you know the transformation formula to compute the scale based on mean and shape, you can use

ps_mean <- fitted(fit2, summary = FALSE)
ps_shape <- as.vector(as.matrix(fit2, pars = "^shape$"))
ps_scale <- trans(ps_mean, ps_shape)
posterior_summary(ps_scale)

where trans() is a vectorized function computing the scale based on mean and shape.

1 Like