Difference in fitting to right censored data between MLE and Bayesian method

Question

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

Answer 1

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$.