2
$\begingroup$

I am fitting a Weibull curve to right censored data. I am doing it by general MLE method using Survival::survreg() as well as Bayesian method using brms::brm.

I am pretty sure that I am getting the model right.

In the results, I am getting similar shape parameter for the Weibull curve for both methods, but the Scale parameter is always (in different datasets) smaller in the Bayesian method.

Take this code for instance, that creates a distribution based on a shape and scale and then fits to them with both methods:

rweibull_cens <- function(n, shape, scale) {
  a_random_death_time <- rweibull(n, shape = shape, scale = scale) 
  a_random_censor_time <- rweibull(n, shape = shape, scale = scale)
  observed_time <- pmin(a_random_censor_time, a_random_death_time)
  censor <- observed_time == a_random_death_time
  tibble(time = observed_time, censor = censor)
}

n = 1e3

rweibull_cens(n, shape = 1.5, scale = 200) -> d


df_fit <- survival::survreg(Surv(time, censor) ~ 1,
  data = d,
  dist = "weibull"
)

scale <- tidy(df_fit )[1, 2] %>%
  rename(scale = estimate) %>%
  exp() %>%
  round(3)

shape <- tidy(df_fit )[2, 2] %>%
  rename(shape = estimate) %>%
  exp() %>%
  .^-1 %>%
  round(3)

d %>%
  mutate(censor = if_else(censor == 0, 1, 0)) %>%
  brm(time | cens(censor) ~ 1, data = ., family = "weibull",  cores = 4) -> bfit

print(bfit, digits = 3)
report = summary(bfit)
bayes_shape = round(report$spec_pars$Estimate,3)
bayes_intercept = round(report$fixed$Estimate,3)
bayes_scale = round(exp(bayes_intercept),3)

print(paste("n =", n))
print(paste("Bayes Shape = ",bayes_shape))
print(paste("Shape = ",shape))
print(paste("Bayes Scale = ",bayes_scale))
print(paste("Scale = ",scale))

I was wondering if anyone can help me why this behavior happens. Thanks

$\endgroup$
5
  • 1
    $\begingroup$ What type of priors are you using and how large is your data set? Both are important when comparing MLE and Bayesian estimation. $\endgroup$ Sep 1, 2021 at 0:11
  • $\begingroup$ The database is large, about 1000 observations. The prior for the intercept is studnetT, and I did play with it but it converges to the smaller value than of MLE either way. My only guess is that there might be inherent difference in the way these methods are accounting for censoring. $\endgroup$ Sep 1, 2021 at 2:35
  • 1
    $\begingroup$ @PrestonBotter I did as well provide a sample code that replicates this issue. $\endgroup$ Sep 1, 2021 at 3:11
  • $\begingroup$ @PrestonBotter were you able to run the code snipet? $\endgroup$ Sep 5, 2021 at 22:52
  • 1
    $\begingroup$ I'm puzzled by this, too. You don't need censored data to reproduce this. It doesn't seem to be sample-size related. Even putting a strong prior onto the known scale value leads to an underestimate with brms. If you have a chance, see what happens with uncensored data and alternate Bayesian Weibull-fitting software. I'm not familiar enough with this type of modeling to know if this behavior is inevitable or represents a problem specific to brms. I can't come back to this for a week or so. If the issue isn't resolved by then, I'll put a bounty on this question. $\endgroup$
    – EdM
    Sep 7, 2021 at 13:50

1 Answer 1

3
$\begingroup$

This is yet another case where you must pay very close attention to the parameterization of a Weibull model.* This doesn't depend on having censored data, as the same issue arises with uncensored data.

The brms package uses a shape, scale parameterization as in the standard R rweibull() function. In the brms vignette on parameterization, the authors use $\alpha$ to represent the shape and $s$ to represent the scale of a Weibull model.

The "Intercept" in a brm() Weibull model, however, is not $\log(s)$, as assumed by the OP in the statement bayes_scale = round(exp(bayes_intercept),3). As the vignette explains, the package authors

use $\mu$ to refer to the main model parameter, which is usually the mean of the response distribution or some closely related quantity. In a regression framework, $\mu$ is not estimated directly but computed as $\mu=g(\eta)$, where $\eta$ is a predictor term ... and $g$ is the response function (i.e., inverse of the link function).

The "Intercept" term returned by brm() is thus $\log (\mu)$, the log of the mean of the distribution. For a given shape $\alpha$ the Weibull scale parameter $s$ is then:

$$ s = \mu/\Gamma(1+1/\alpha),$$

as described in the vignette.

When I followed the code in this question (using set.seed(101) before generating the random Weibull sample), brm() gave an Intercept of 5.238432 and a shape of 1.461953. With the above formula, that gives a scale $s$ of 207.99.

With the same random Weibull sample, survreg() gave an Intercept (here, the log of the scale in the rweibull() parameterization) of 5.335683, for a scale of 207.61. The apparent discrepancy between the maximum-likelihood and Bayesian estimates of scale thus disappears once you recognize that the Intercept in the latter case, from brm(), is $\log(\mu)$, not $\log(s)$.

As I searched the web for similar cases, I found that this definition of the Intercept for a brm() Weibull model isn't always appreciated, as for example on this page. If you generate samples based on a Weibull shape value near 1, as on that page, you wouldn't notice the discrepancy caught by the OP here, as $\Gamma(2)=1$.


*I got the crucial hint for this solution from this page by Riley King, that illustrates frequentist and Bayesian Weibull modeling of both uncensored and censored data. The author used a grid approximation directly on shape and scale for a Bayesian posterior estimate on uncensored data, then moved to brm() for censored data. He read the manual closely enough to recognize that:

The parameters that get estimated by brm() are the Intercept and shape. We can use the shape estimate as-is, but it’s a bit tricky to recover the scale. The key is that brm() uses a log-link function on the mean $\mu$.

$\endgroup$
5
  • $\begingroup$ Wonderful work @EdM. I think that the decision of the brms developer to use such a non-standard parameterization needs to be rethought. $\endgroup$ Sep 8, 2021 at 13:07
  • $\begingroup$ @FrankHarrell thanks. This might also become an issue with the rstanarm package, whose development version is introducing survival analysis functions. (I haven't succeeded in installing the development version.) I don't know if this choice of "Intercept" definition is inherent in stan or was a separate choice for the brms front end. $\endgroup$
    – EdM
    Sep 8, 2021 at 15:36
  • $\begingroup$ My quick read of the rstanarm survival functions is that they use standard notation for the Weibull distribution. Someone please correct me if otherwise. I wish this were in CRAN. rstanarm has not been updated on CRAN in a long time. $\endgroup$ Sep 8, 2021 at 16:52
  • $\begingroup$ Thanks a lot @EdM! $\endgroup$ Sep 13, 2021 at 16:45
  • 1
    $\begingroup$ @FrankHarrell, here is one of the authors of the stan_surv paper: it is indeed a pity that it hasn’t been added to CRAN, yet. This is not due to instability or whatsoever of the survival functionality, but rather due to infrastructure problems (early 2020 it was due to compile time limits on windows machines). I just contacted the main authors of the rstanarm package again to flag this. $\endgroup$
    – ermeel
    Oct 14, 2022 at 3:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.