Survival Analysis - Fitting Weibull Models for Improving Device Reliability in R

Construct Weibull model from un-censored data using fitdistrplus

To start out with, let’s take a frequentist approach and fit a 2-parameter Weibull distribution to these data. Once the parameters of the best fitting Weibull distribution of determined, they can be used to make useful inferences and predictions.

I’ll use the fitdist() function from the fitdistrplus package to identify the best fit via maximum likelihood. The parameters we care about estimating are the shape and scale.

# Fit model and extract parameters of interest
mle_wieb_nocens_fit <- fitdist(data_tbl$fatigue_duration, "weibull")
weib_shape <- mle_wieb_nocens_fit$estimate["shape"]
weib_scale <- mle_wieb_nocens_fit$estimate["scale"]

# Summarize and plot
summary(mle_wieb_nocens_fit)

## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters : 
##         estimate Std. Error
## shape   1.832201  0.1401198
## scale 100.841802  5.8068570
## Loglikelihood:  -526.2573   AIC:  1056.515   BIC:  1061.725 
## Correlation matrix:
##           shape     scale
## shape 1.0000000 0.3186071
## scale 0.3186071 1.0000000

plot(mle_wieb_nocens_fit)

The Weibull isn’t the only possible distribution we could have fit. Lognormal and gamma are both known to model time-to-failure data well. They are shown below using the denscomp() function from fitdistrplus.

# Fit gamma model, extract shape, rate
mle_gamma_nocens_fit <- fitdist(data_tbl$fatigue_duration, "gamma")
gamma_shape <- mle_gamma_nocens_fit$estimate["shape"]
gamma_rate <- mle_gamma_nocens_fit$estimate["rate"]

# Fit lognormal model, extract mean, sd
mle_lognormal_nocens_fit <- fitdist(data_tbl$fatigue_duration, "lnorm")
lnorm_meanlog <- mle_lognormal_nocens_fit$estimate["meanlog"]
lnorm_sdlog <- mle_lognormal_nocens_fit$estimate["sdlog"]

# visualize in fitdistrplus
plot.legend <- c("gamma", "lognormal", "Weibull")
denscomp(list(mle_gamma_nocens_fit, mle_lognormal_nocens_fit, mle_wieb_nocens_fit), legendtext = plot.legend)

I recreate the above in ggplot2, for fun and practice.

x_series <- seq(0, 300, by = 1)
x_tbl <- tibble(x = x_series)


wide_densities_tbl <- x_tbl %>%
  mutate(
    Weibull = dweibull(x, shape = weib_shape, scale = weib_scale),
    gamma = dgamma(x, shape = gamma_shape, rate = gamma_rate),
    lognormal = dlnorm(x, meanlog = lnorm_meanlog, sdlog = lnorm_sdlog)
  ) %>%
  gather("distribution", "value", -x)

wide_densities_tbl %>% ggplot(aes(x = x)) +
  geom_line(aes(
    x = x,
    y = value,
    color = distribution,
    linetype = distribution
  ),
  size = .8
  ) +
  geom_histogram(
    data = data_tbl,
    aes(
      x = fatigue_duration,
      y = ..density..
    ),
    binwidth = 50,
    boundary = 300,
    color = "white",
    fill = "#2c3e50",
    alpha = .6
  ) +
  labs(
    x = "Time to Fatigue Failure Event (Days)",
    y = "Density",
    title = "Histogram and theoretical densities for Fatigue Data",
    subtitle = "Parametric fits using MLE via fitdist() function"
  ) +
  theme(legend.title = element_blank())

Goodness-of-fit statistics are available and shown below for reference. If available, we would prefer to use domain knowledge and experience to identify what the true distribution is instead of these statistics which are subject to sampling variation. It is not good practice to stare at the histogram and attempt to identify the distribution of the population from which it was drawn. Nevertheless, we might look at the statistics below if we had absolutely no idea the nature of the data generating process / test.

# get goodness-of-fit
gofstat(list(mle_wieb_nocens_fit, mle_gamma_nocens_fit, mle_lognormal_nocens_fit),
  fitnames = c("weib", "gamma", "lnorm")
)

## Goodness-of-fit statistics
##                                    weib      gamma      lnorm
## Kolmogorov-Smirnov statistic 0.05066902 0.03643609 0.06948557
## Cramer-von Mises statistic   0.03408173 0.02106615 0.12076742
## Anderson-Darling statistic   0.21634543 0.17192106 0.81630880
## 
## Goodness-of-fit criteria
##                                    weib    gamma    lnorm
## Akaike's Information Criterion 1056.515 1056.408 1067.427
## Bayesian Information Criterion 1061.725 1061.618 1072.637

Using the model to infer device reliability

The model by itself isn’t what we are after. It is the vehicle from which we can infer some very important information about the reliability of the implant design. First and foremost - we would be very interested in understanding the reliability of the device at a time of interest. For instance, suppose our voice of customer research indicates that our new generation of device needs to last 10 months in vivo to be safe and competitive. Recall that each day on test represents 1 month in service. Once we fit a Weibull model to the test data for our device, we can use the reliability function to calculate the probability of survival beyond time t.³

\[\text{R} (t | \beta, \eta) = e ^ {- \bigg (\frac{t}{\eta} \bigg ) ^ {\beta}}\]

Note:

t = the time of interest (for example, 10 years)

\(\beta\) = the Weibull scale parameter

\(\eta\) = the Weibull shape parameter

This looks a little nasty but it reads something like “the probability of a device surviving beyond time t conditional on parameters \(\beta\) and \(\eta\) is [some mathy function of t, \(\beta\) and \(\eta\)]. For the model we fit above using MLE, a point estimate of the reliability at t=10 years (per the above VoC) can be calculated with a simple 1-liner:

reliability_at_10 <- exp(-(10 / weib_scale)**(weib_shape))
reliability_at_10 %>%
  scales::percent() %>%
  kable(align = "c")

In this way we infer something important about the quality of the product by fitting a model from benchtop data.

It is common to report confidence intervals about the reliability estimate but this practice suffers many limitations. The intervals change with different stopping intentions and/or additional comparisons. They also do not represent true probabilistic distributions as our intuition expects them to and cannot be propagated through complex systems or simulations. What we’d really like is the posterior distribution for each of the parameters in the Weibull model, which provides all credible pairs of \(\beta\) and \(\eta\) that are supported by the data. For that, we need Bayesian methods which happen to also be more fun.

Use grid approximation to estimate posterior

This problem is simple enough that we can apply grid approximation to obtain the posterior. In this method we feed in a sequence of candidate combinations for \(\beta\) and \(\eta\) and determine which pairs were most likely to give rise to the data. The likelihood is multiplied by the prior and converted to a probability for each set of candidate \(\beta\) and \(\eta\). Flat priors are used here for simplicity - I’ll put more effort into the priors later on in this post. Since the priors are flat, the posterior estimates should agree with the maximum likelihood point estimate.

Calculate posterior via grid approximation:⁴

# function to get log-likelihood of the data for a given shape & scale pair
grid_function <- function(shape, scale) {
  dweibull(data_tbl$fatigue_duration, shape = shape, scale = scale, log = T) %>%
    sum()
}

# set up grid of possible shape, scale parameters
n <- 100
shape_grid <- seq(1, 3, length.out = n)
scale_grid <- seq(60, 130, length.out = n)
two_param_grid <- expand_grid(shape_grid, scale_grid)

# map the grid_function over all candidate parameter pairs
# multiply LL by prior and convert to probability
full_tbl <- two_param_grid %>%
  mutate(log_likelihood = map2(shape_grid, scale_grid, grid_function)) %>%
  unnest() %>%
  mutate(
    shape_prior = 1,
    scale_prior = 1
  ) %>%
  mutate(product = log_likelihood + shape_prior + scale_prior) %>%
  mutate(probability = exp(product - max(product)))

full_tbl %>%
  head(10) %>%
  kable(align = rep("c", 7))

Plot the grid approximation of the posterior. If we super-impose our point estimate from Part 1, we see the maximum likelihood estimate agrees well with the mode of the joint posterior distributions for shape and scale.

# need this data to feed to gg_label_repel to tell it where to attach label
point_tbl <- tibble(x = weib_shape, y = weib_scale)

# visualize and compare to MLE
plt_1 <- full_tbl %>%
  ggplot(aes(x = shape_grid, y = scale_grid)) +
  geom_raster(aes(fill = probability),
    interpolate = T
  ) +
  geom_point(x = weib_shape, y = weib_scale, size = 1.3) +
  geom_label_repel(
    data = point_tbl, aes(x, y),
    label = "Maximum Likelihood Estimate \n[from Part 1]",
    fill = "#8bd646ff",
    color = "black",
    segment.color = "black",
    segment.size = 1,
    #                   min.segment.length = unit(1, "lines"),
    nudge_y = -16,
    nudge_x = .5
  ) +
  scale_fill_viridis_c() +
  labs(
    title = "Posterior for Weibull Shape and Scale Parameters",
    subtitle = "Calculated with Grid Approximation",
    x = expression(eta ["shape"]),
    y = expression(beta ["scale"])
  ) +
  theme(panel.grid = element_blank())

plt_1

3d with rayshader just to flex :)

# par(mfrow = c(1, 1))
# plot_gg(plt_1, width = 5, height = 4, scale = 300, multicore = TRUE, windowsize = c(1200, 960),
#        fov = 70, zoom = 0.45, theta = 330, phi = 40)

Visualize draws from the posterior

We can sample from the grid to get the same if we weight the draws by probability.

set.seed(2020)

# take a weighted sample from the posterior
grid_approx_posterior_samples <- full_tbl %>%
  sample_n(size = 2000, replace = T, weight = probability)

# visualize
grid_approx_plot <- grid_approx_posterior_samples %>%
  ggplot(aes(x = shape_grid, y = scale_grid)) +
  geom_point(
    colour = "#e56a5dff",
    size = 2,
    alpha = 0.3
  ) +
  labs(
    x = expression(eta ["shape"]),
    y = expression(beta ["scale"])
  ) +
  geom_density_2d(color = "black", size = 1.2, alpha = .4) +
  labs(
    title = "Estimate of Joint Probabilities for Shape and Scale",
    subtitle = "Via Posterior Sampling"
  )

grid_approx_marg_plot <- ggMarginal(grid_approx_plot,
  type = "density",
  color = "white",
  alpha = 0.7,
  fill = "#e56a5dff"
)

Uncertainty in the implied time-to-failure curves

Each of the credible parameter values implies a possible Weibull distribution of time-to-failure data from which a reliability estimate can be inferred. This is a good way to visualize the uncertainty in a way that makes intuitive sense.

# sample from posterior; weighted as prob
grid_approx_posterior_samples_4plot <- full_tbl %>%
  sample_n(size = 1000, replace = T, weight = probability) %>%
  select(shape_grid, scale_grid)

# plot 1000 Weibull curves from samples parameters
weib_uncertainty_plt <- grid_approx_posterior_samples_4plot %>%
  mutate(p_y_data = map2(
    shape_grid, scale_grid,
    ~ tibble(
      x = seq(0, 200, length.out = 400),
      y = dweibull(x, .x, .y)
    )
  )) %>%
  mutate(row_id = row_number()) %>%
  unnest(p_y_data) %>%
  ggplot(aes(x = x, y = y)) +
  geom_line(aes(group = row_id), alpha = .05, color = "#440154ff") +
  labs(
    x = "Time",
    y = "Density",
    title = "Distribution of Time at Failure Event",
    subtitle = "Implied by Posterior Distribution [Sample of n=1000]"
  )

Uncertainty in the implied reliabilty of the device

Any row-wise operations performed will retain the uncertainty in the posterior distribution. This allows for a straightforward computation of the range of credible reliabilities at t=10 via the reliability function.

# use params in posterior to calculate correstponding reliabilities at t=10
grid_approx_posterior_samples_4relplot <- full_tbl %>%
  sample_n(size = 4000, replace = T, weight = probability) %>%
  select(shape_grid, scale_grid) %>%
  mutate(reliability_at_10 = exp(-(10 / scale_grid)**(shape_grid)))

# visualize reliability distribution
grid_approx_posterior_samples_4relplot %>%
  ggplot(aes(reliability_at_10)) +
  geom_histogram(aes(y = ..density..), fill = "#2c3e50", color = "white", alpha = .6) +
  labs(
    title = "Reliability Distribution",
    subtitle = "Calculated from Grid Approximation with Flat Priors"
  )

This distribution gives much richer information than the MLE point estimate of reliability. The most credible estimate of reliability is ~ 98.8%, but it could plausibly also be as low as 96%. This delta can mean the difference between a successful and a failing product and should be considered as you move through project phase gates.

Point estimate with right-censored data

First, I’ll set up a function to generate simulated data from a Weibull distribution and censor any observations greater than 100. I set the function up in anticipation of using the survreg() function from the survival package in R. The syntax is a little funky so some additional detail is provided below.

# function to generate random Weibull data and censor data > 100
rweibull_cens_mod_fcn <- function(n, shape, scale) {
  raw_times <- rweibull(n, shape = shape, scale = scale)
  tibble(failure_time_raw = raw_times) %>%
    mutate(time = case_when(
      failure_time_raw < 100 ~ failure_time_raw,
      TRUE ~ 100
    )) %>%
    mutate(censor = case_when(
      time == 100 ~ 0,
      TRUE ~ 1
    )) %>%
    mutate(censor = censor == 1) %>%
    mutate(time = time %>% round(digits = 2)) %>%
    select(-failure_time_raw)
}

# set seed for repeatability
set.seed(54)
first_test_fit_tbl <- rweibull_cens_mod_fcn(30, 3, 100)
first_test_fit_tbl %>%
  head(10) %>%
  kable(align = rep("c", 2))

I admit this looks a little strange because the data that were just described as censored (duration greater than 100) show as “FALSE” in the censored column. This is due to the default syntax of the survreg() function in the survival package that we intend to fit the model with:⁵

survival package defaults for censoring:

• 0 or FALSE for censoring, 1 or TRUE for observed event

To further throw us off the trail, the survreg() function returns “scale”" and “intercept”" that must be converted to recover the shape and scale parameters that align with the rweibull() function used to create the data. Don’t fall for these tricks - just extract the desired information as follows:

survival package defaults for parameterizing the Weibull distribution:

• survreg’s scale parameter = 1/(rweibull shape parameter)
• survreg’s intercept = log(rweibull scale parameter)

Ok let’s see if the model can recover the parameters when we providing survreg() the tibble with n=30 data points (some censored):

# fit model to simulated data (n=30)
test_fit <- survival::survreg(Surv(time, censor) ~ 1,
  data = first_test_fit_tbl,
  dist = "weibull"
)

# evaluate the fit
summary(test_fit)

## 
## Call:
## survival::survreg(formula = Surv(time, censor) ~ 1, data = first_test_fit_tbl, 
##     dist = "weibull")
##               Value Std. Error    z       p
## (Intercept)  4.5469     0.0896 50.7 < 2e-16
## Log(scale)  -0.9266     0.1972 -4.7 2.6e-06
## 
## Scale= 0.396 
## 
## Weibull distribution
## Loglik(model)= -106.7   Loglik(intercept only)= -106.7
## Number of Newton-Raphson Iterations: 5 
## n= 30

Extract and covert shape and scale with broom::tidy() and dplyr:

# extract scale parameter
scale <- tidy(test_fit)[1, 2] %>%
  rename(scale = estimate) %>%
  exp() %>%
  round(2)

# extract shape parameter
shape <- tidy(test_fit)[2, 2] %>%
  rename(shape = estimate) %>%
  exp() %>%
  .^-1 %>%
  round(2)

# summarize
point_estimates <- bind_cols(shape, scale)
point_estimates %>% kable(align = rep("c", 2))

What has happened here? We know the true parameters are shape = 3, scale = 100 because that’s how the data were generated. These point estimates are pretty far off. Is the survreg() fitting function broken? Are there too few data and we are just seeing sampling variation? Is it confused by the censored data? I honestly don’t know. To answer these questions, we need a new function that fits a model using survreg() for any provided sample size. The data to make the fit are generated internal to the function. The function returns a tibble with estimates of shape and scale for that particular trial:

Simulation to understand point estimate sensitivity to sample size

# seed for repeatability
set.seed(2025)

# fcn takes sample size n, simulates weibull(3, 100) data, fits model, estimates params
test_map_fcn <- function(n) {
  holder <- rweibull_cens_mod_fcn(n, 3, 100)

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

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

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

  point_estimates <- bind_cols(shape, scale)
  point_estimates
}

Now that we have a function that takes a sample size n and returns fitted shape and scale values, we want to apply the function across many values of n. Let’s look at what happens to our point estimates of shape and scale as the sample size n increases from 10 to 1000 by 1.

# create sequence of n's
n_sim_mle <- seq(10, 1000, by = 1) %>%
  tibble() %>%
  rename(n = ".")

# map the fitting function across vector of n's.
results_tbl <- n_sim_mle %>%
  mutate(results = map(n, test_map_fcn)) %>%
  unnest()

# peek at format of results
results_tbl %>%
  head(5) %>%
  kable(align = rep("c", 3))

Visualize results of the simulation:

results_tbl %>%
  gather(key = key, value = value, -n) %>%
  ggplot(aes(x = n, y = value)) +
  geom_point(color = "#5C126EFF", alpha = 0.3) +
  facet_wrap(~key, scales = "free_y") +
  geom_smooth(color = "#F17020FF", alpha = 0.9) +
  labs(
    x = "Sample Size",
    y = "",
    title = "Estimated Shape and Scale Parameters for Different Data Set Sizes",
    subtitle = "Weibull Regressions using survival package.  Data generated from Weibull(3, 100)"
  )

This simulation is illuminating. It’s apparent that there is sampling variability effecting the estimates. On average, the true parameters of shape = 3 and scale = 100 are correctly estimated. But on any given experimental run, the estimate might be off by quite a bit. The precision increases with sample size as expected but the variation is still relevant even at large n.

Based on this simulation we can conclude that our initial point estimate of 2.5, 94.3 fit from n=30 is within the range of what is to be expected and not a software bug or coding error.

Estimates for product reliability at 15, 30, 45, and 60 months are shown below.

results_tbl %>%
  mutate(reliability_at_15 = exp(-(15 / scale)**(shape))) %>%
  mutate(reliability_at_30 = exp(-(30 / scale)**(shape))) %>%
  mutate(reliability_at_45 = exp(-(45 / scale)**(shape))) %>%
  mutate(reliability_at_60 = exp(-(60 / scale)**(shape))) %>%
  select(-c(shape, scale)) %>%
  gather(key = "key", value = "value", -n) %>%
  ggplot(aes(x = n, y = value)) +
  geom_point(aes(color = key), alpha = .4) +
  #  geom_smooth(aes(group = key), color = "black", size = .4, alpha = .5) +
  ylim(c(.73, 1)) +
  scale_color_viridis_d(option = "C", begin = 0, end = .8, direction = -1) +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  labs(
    title = "Reliability Point Estimates",
    subtitle = "Effect of Sample Size and Sampling Variability",
    x = "Sample Size",
    y = "Reliability Estimate"
  ) +
  theme(legend.title = element_blank())

The above gives a nice sense of the uncertainty in the reliability estimate as sample size increases, but you can’t actually simulate a confidence interval from those data because there aren’t enough data points at any one sample size. To do that, we need many runs at the same sample size.

Simulation of 95% confidence intervals on reliability

In the code below, I generate n=1000 simulations of n=30 samples drawn from a Weibull distribution with shape = 3 and scale = 100. For each set of 30 I fit a model and record the MLE for the parameters.

# each simulation will be like a benchtop test of n=30 parts
ci_reps_tbl <- tibble(sample_size = rep(30, 1000))

set.seed(98)

# loop: draw 30 from rweibull(30, 100); fit model, estimate parameters
ci_results_tbl <- ci_reps_tbl %>%
  mutate(results = map(sample_size, test_map_fcn)) %>%
  unnest()

# peek at results
ci_results_tbl %>%
  head(10) %>%
  kable(align = rep("c", 3))

95% of the reliability estimates lik above the .05 quantile. This means the .05 quantile is the analogous boundary for a simulated 95% confidence interval. In the code below, the .05 quantile of reliability is estimated for each time requirement of interest where we have 1000 simulation at each.

# function will take a time of interest and calculate the implied reliability
# only works with the params from the ci_results_tbl above
reliability95_fcn <- function(t_requirement) {
  ci_results2tbl <- ci_results_tbl %>%
    mutate(reliability_at_t = exp(-(t_requirement / scale)**(shape))) %>%
    mutate(q_05 = quantile(reliability_at_t, .05))
  ci_results2tbl
}

# map the function across a sequence of candidate time requirements
reliability_seq_tbl <- tibble(time_requirement = seq(5, 100, by = 5)) %>% mutate(rel_tbl = map(time_requirement, reliability95_fcn))

# calculate the .05 quantile for reliability in each set of 1000 (each candidate time requirement)
reliability_summary_tbl <- reliability_seq_tbl %>%
  unnest(rel_tbl) %>%
  group_by(time_requirement) %>%
  summarize(rel_q05 = mean(q_05)) %>%
  ungroup()

# Peek
reliability_summary_tbl %>% kable(align = rep("c", 2))

# reference tibbles for labeling w/ gg_label_repel
label_ci_tbl <- tibble(x = 45, y = 0.8429503)
label_95rel_tbl <- tibble(x = 10, y = 0.95)

# visualize
reliability_seq_tbl %>%
  unnest(rel_tbl) %>%
  ggplot(aes(x = time_requirement, y = reliability_at_t)) +
  geom_jitter(alpha = 0.03) +
  #  geom_point(data = reliability_summary_tbl, aes(x = time_requirement, y = rel_q05), color = "#c44c74ff", size = 1) +
  geom_smooth(data = reliability_summary_tbl, aes(x = time_requirement, y = rel_q05), color = "#c44c74ff", size = 1, alpha = .2) +
  geom_hline(aes(yintercept = .95), color = "#240691ff", size = 1) +
  geom_label_repel(
    data = label_ci_tbl, aes(x, y),
    label = "1-sided 95% Confidence\n bound on Reliability",
    fill = "#c44c74ff",
    color = "#f0f921ff",
    segment.color = "#c44c74ff",
    segment.size = 1,
    #                   min.segment.length = unit(1, "lines"),
    nudge_y = -.5,
    nudge_x = -5
  ) +
  geom_label_repel(
    data = label_95rel_tbl, aes(x, y),
    label = "Product Reliability\n Requirement: 95%",
    fill = "#240691ff",
    color = "#f0f921ff",
    segment.color = "#240691ff",
    segment.size = 1,
    #                   min.segment.length = unit(1, "lines"),
    nudge_y = -.25,
    nudge_x = 5
  ) +
  labs(
    title = "Approximate 95% Confidence Interval for Various Durability Requirements",
    subtitle = "Each Sim: 1000 models fit from 1000 draws of n=30 from Weibull(3, 100)",
    x = "Device Requirement (Service Life without Failure)",
    y = "Reliability"
  )

There’s a lot going on here so it’s worth it to pause for a minute. If I was to try to communicate this in words, I would say:

• Assume we have designed a medical device that fails according to a Weibull distribution with shape = 3 and scale = 100.
• Assume the service life requirement for the device is known and specified within the product’s requirements
• Assume we can only test n=30 units in 1 test run and that testing is expensive and resource intensive
• The n=30 failure/censor times will be subject to sampling variability and the model fit from the data will likely not be Weibull(3, 100)
• The variability in the parameter estimates is propagated to the reliability estimates - a distribution of reliability is generated for each potential service life requirement (in practice we would only have 1 requirement)
• The .05 quantile of the reliability distribution at each requirement approximates the 1-sided lower bound of the 95% confidence interval. This threshold changes for each candidate service life requirement.

Why does any of this even matter? Here’s the TLDR of this whole section:

Suppose the service life requirement for our device is 24 months (2 years). Our boss asks us to set up an experiment to verify with 95% confidence that 95% of our product will meet the 24 month service requirement without failing. The industry standard way to do this is to test n=59 parts for 24 days (each day on test representing 1 month in service). If all n=59 pass then we can claim 95% reliability with 95% confidence. However, if we are willing to test a bit longer then the above figure indicates we can run the test to failure with only n=30 parts instead of n=59. If it cost a lot to obtain and prep test articles (which it often does), then we just saved a ton of money and test resources by treating the data as variable instead of attribute. ⁶ We also get information about the failure mode for free.

Use brm() to generate a posterior distribution for shape and scale

The formula for asking brms to fit a model looks relatively the same as with survival. To start, we fit a simple model with default priors.

# fit model with brm() [use default priors]
# mtf_weib_fit <- brm(time | cens(censored) ~ 1,
# data = months_to_failure_tbl_4.3.3, family = weibull())

Just like with the survival package, the default parameterization in brms can easily trip you up. We are fitting an intercept-only model meaning there are no predictor variables. 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\). There is no doubt that this is a rambling post - even so, it is not within scope to try to explain link functions and GLM’s (I’m not expert enough to do it anyways, refer to Statistical Rethinking by McElreath). In short, to convert to scale we need to both undo the link function by taking the exponent and then refer to the brms documentation to understand how the mean \(\mu\) relates to the scale \(\beta\). The operation looks like this:⁷

#scale = exp(Intercept)/(gamma(1 + 1/shape))

Examine the results:

# saveRDS(mtf_weib_fit, file = "censored_data_mtf_weib_fit.rds")
mtf_weib_fit <- readRDS(file = "censored_data_mtf_weib_fit.rds")

summary(mtf_weib_fit)

##  Family: weibull 
##   Links: mu = log; shape = identity 
## Formula: time | cens(censored) ~ 1 
##    Data: months_to_failure_tbl_4.3.3 (Number of observations: 30) 
## 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 Rhat Bulk_ESS Tail_ESS
## Intercept     4.41      0.08     4.25     4.58 1.00     2349     1782
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     2.74      0.52     1.82     3.84 1.00     2530     2281
## 
## 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).

Gut-check on convergence of chains. Things look good visually and Rhat = 1 (also good).

plot(mtf_weib_fit)

Within the tibble of posterior draws we convert the intercept to scale using the formula previously stated.

# extract scale and shape
post_mtf_weib_samples <- posterior_samples(mtf_weib_fit) %>%
  mutate(scale = exp(b_Intercept) / (gamma(1 + 1 / shape))) %>%
  select(shape, scale)

# peek
post_mtf_weib_samples %>%
  head(10) %>%
  kable(align = rep("c", 2), digits = 2)

Here is our first look at the posterior drawn from a model fit with censored data. We know the data were simulated by drawing randomly from a Weibull(3, 100) so the true data generating process is marked with lines.

mtf_plot <- post_mtf_weib_samples %>%
  ggplot(aes(x = shape, y = scale)) +
  geom_point(
    colour = "#453781FF",
    size = 2,
    alpha = 0.1
  ) +
  geom_hline(aes(yintercept = 100), size = .5, alpha = .3) +
  geom_vline(aes(xintercept = 3), size = .5, alpha = .3) +
  labs(
    title = "Credible Parameters for Shape and Scale",
    subtitle = "Run-Out Data Treated as Right-Censored",
    x = expression(eta ["shape"]),
    y = expression(beta ["scale"])
  )


mtf_marg_plot <- ggMarginal(mtf_plot,
  type = "density",
  color = "white",
  alpha = 0.7,
  fill = "#453781FF"
)

It looks like we did catch the true parameters of the data generating process within the credible range of our posterior. However, it is certainly not centered. Once again we should question: is the software working properly? Is the sample size a problem? Are the priors appropriate? Was the censoring specified and treated appropriately?

Investigation of how to treat censored data points

Let’s start with the question about the censoring. One question that I’d like to know is: What would happen if we omitted the censored data completely or treated it like the device failed at the last observed time point? This hypothetical should be straightforward to simulate. Let’s fit a model to the same data set, but we’ll just treat the last time point as if the device failed there (i.e. we’ll have lots of failures at t=100).

# This model treats the data as though the devices failed at the last observed timepoint
# mtf_weib_nocens_fit <- brm(time ~ 1,
# data = months_to_failure_tbl_4.3.3, family = weibull())

# saveRDS(mtf_weib_nocens_fit, file = "mtf_weib_nocens_fit.rds")
mtf_weib_nocens_fit <- readRDS(file = "mtf_weib_nocens_fit.rds")

Now another model where we just omit the censored data completely (i.e. remove any units that don’t fail from the data set completely and fit a model to the rest).

# This model just omits the censored data completely
# mtf_weib_omit_fit <- brm(time ~ 1,
# data = months_to_failure_tbl_4.3.3 %>% filter(censored == 0), family = weibull())

 plot(mtf_weib_nocens_fit)

 plot(mtf_weib_omit_fit)

Create tibble of posterior draws from partially censored, un-censored, and censor-omitted models with identifier column.

cens_model_post_draws <- posterior_samples(mtf_weib_fit) %>%
  mutate(scale = exp(b_Intercept) / (gamma(1 + 1 / shape))) %>%
  select(shape, scale) %>%
  mutate(model = "Weibull with Censoring")

uncens_model_post_draws <- posterior_samples(mtf_weib_nocens_fit) %>%
  mutate(scale = exp(b_Intercept) / (gamma(1 + 1 / shape))) %>%
  select(shape, scale) %>%
  mutate(model = "Weibull (Treat Last Timepoint as Failure)")

omit_model_post_draws <- posterior_samples(mtf_weib_omit_fit) %>%
  mutate(scale = exp(b_Intercept) / (gamma(1 + 1 / shape))) %>%
  select(shape, scale) %>%
  mutate(model = "Weibull Omitting Censored Data")

combined_weib_tbl <- cens_model_post_draws %>%
  bind_rows(uncens_model_post_draws) %>%
  bind_rows(omit_model_post_draws) %>%
  mutate(model = as_factor(model))

combined_weib_tbl %>%
  head(5) %>%
  kable(align = rep("c", 3), digits = 2)

Here we compare the effect of the different treatments of censored data on the parameter estimates. Intervals are 95% HDI.

d_shp_1 <- combined_weib_tbl %>% ggplot(aes(x = shape)) +
  geom_density(aes(fill = model), size = 0, alpha = 0.4) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(x = expression(eta["shape"])) +
  stat_pointintervalh(aes(y = 0),
    point_interval = mode_hdi, .width = .95, show.legend = FALSE
  ) +
  facet_wrap(~model) +
  theme_classic() +
  theme(
    legend.position = "none",
    axis.line.y = element_blank(),
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    legend.title = element_blank()
  )


d_shp_2 <- combined_weib_tbl %>% ggplot(aes(x = shape)) +
  geom_density(aes(fill = model), size = 0, alpha = 0.4) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(
    x = expression(eta["shape"]),
    title = "Parameter Estimates",
    subtitle = "Effect of Including Censored Data in Model"
  ) +
  theme_classic() +
  theme(
    axis.line.y = element_blank(),
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    legend.title = element_blank()
  )

d_scale_1 <- combined_weib_tbl %>% ggplot(aes(x = scale, fill = model)) +
  geom_density(size = 0, alpha = 0.4) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(x = expression(beta["scale"])) +
  stat_pointintervalh(aes(y = 0),
    point_interval = mode_hdi, .width = .95, show.legend = FALSE
  ) +
  facet_wrap(~model) +
  theme_classic() +
  theme(
    legend.position = "none",
    axis.line.y = element_blank(),
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    legend.title = element_blank()
  )

d_scale_2 <- combined_weib_tbl %>% ggplot(aes(x = scale)) +
  geom_density(aes(fill = model), size = 0, alpha = 0.4) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(
    x = expression(beta["scale"]),
    title = "Parameter Estimates",
    subtitle = "Effect of Including Censored Data in Model"
  ) +
  theme_classic() +
  theme(
    axis.line.y = element_blank(),
    strip.background = element_blank(),
    strip.text.x = element_blank(),
    legend.title = element_blank()
  )

d_shp_2 + d_shp_1 + plot_layout(
  ncol = 1,
  guides = "collect"
)

d_scale_2 + d_scale_1 + plot_layout(
  ncol = 1,
  guides = "collect"
)

When we omit the censored data or treat it as a failure, the shape parameter shifts up and the scale parameter shifts down. In both cases, it moves farther away from true. This should give is confidence that we are treating the censored points appropriately and have specified them correctly in the brm() syntax.

Plotting the joint distributions for the three groups:

combined_weib_plt <- combined_weib_tbl %>%
  ggplot(aes(x = shape, y = scale)) +
  geom_point(aes(color = model, fill = model), size = 2, alpha = 0.08) +
  geom_hline(aes(yintercept = 100), size = .5, alpha = .3) +
  geom_vline(aes(xintercept = 3), size = .5, alpha = .3) +
  scale_color_manual(values = c("#453781FF", "#EA4F88", "#FDE725FF")) +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) + # force legend icons to be alpha = 1 instead of .08
  labs(
    title = "Credible Parameters for Shape and Scale",
    subtitle = "Effect of Treating Censored Data Points in Different Ways",
    x = expression(eta ["shape"]),
    y = expression(beta ["scale"])
  ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  )

cens_marg_plot <- ggMarginal(combined_weib_plt,
  groupColour = TRUE,
  groupFill = TRUE,
  type = "density",
  alpha = 0.7
)

Our censored data set (purple) is closest to true. But we still don’t know why the highest density region of our posterior isn’t centered on the true value.

Evaluation of priors

We haven’t looked closely at our priors yet (shame on me) so let’s do that now. The default priors are viewed with prior_summary().

# get default priors
prior_summary(mtf_weib_fit) %>%
  select(prior, class) %>%
  kable(align = rep("c", 2))

I was taught to visualize what the model thinks before seeing the data via prior predictive simulation. I made a good-faith effort to do that, but the results are funky for brms default priors. I an not an expert here, but I believe this is because very vague default Gamma priors aren’t good for prior predictive simulations but quickly adapt to the first few data points they see.⁸. The prior must be placed on the intercept when must be then propagated to the scale which further muddies things. All in all there isn’t much to see. A lot of the weight is at zero but there are long tails for the defaults. I have all the code for this simulation for the defaults in the Appendix.

If you take this at face value, the model thinks the reliability is always zero before seeing the model. Again, I think this is a special case for vague gamma priors but it doesn’t give us much confidence that we are setting things up correctly.

After viewing the default predictions, I did my best to iterate on the priors to generate something more realistic. My process was manual and my general plan was to force some crdibility over higher values of shape using a uniform distribution. The intercept values were more diffcult to iterate. Here are the revised priors I tried:

tibble(prior = c("student_t(3, 5, 5)", "uniform(0, 10)"),
       class = c("Intercept", "shape")) %>% kable(align = rep("c", 2))

As can be seen, the revised priors were able to spread some credibility up across the middle reliability values but ended up a lot of mass on either end, which wasn’t to goal.

Regardless, I refit the model with the (potentially) improved more realistic (but still not great) priors and found minimal difference in the model fit as shown below.

Note: all models throughout the remainder of this post use the “better” priors (even though there is minimal difference in the model fits relative to brms default).

The above analysis, while not comprehensive, was enough to convince me that the default brms priors are not the problem with initial model fit (recall above where the mode of the posterior was not centered at the true data generating process and we wondered why). I do need to get better at doing these prior predictive simulations but it’s a deep, dark rabbit hole to go down on an already long post. Given the low model sensitivity across the range of priors I tried, I’m comfortable moving on to investigate sample size.

The original model was fit from n=30. We need a simulation that lets us adjust n.

Evaluate sensitivity of posterior to sample size

Here we write a function to generate censored data of different shape, scale, and sample size. The syntax of the censoring column is brms (1 = censored).

rweibull_cens_gen_fcn <- function(n, shape, scale) {
  raw_times <- rweibull(n, shape = shape, scale = scale)
  tibble(failure_time_raw = raw_times) %>%
    mutate(time = case_when(
      failure_time_raw < 100 ~ failure_time_raw,
      TRUE ~ 100
    )) %>%
    mutate(censor = case_when(
      time == 100 ~ 1,
      TRUE ~ 0
    )) %>%
    select(-failure_time_raw)
}

Now the function above is used to create simulated data sets for different sample sizes (all have shape 3, scale = 100)

set.seed(3980)

fit_weib_45_tbl <- rweibull_cens_gen_fcn(45, 3, 100)
fit_weib_60_tbl <- rweibull_cens_gen_fcn(60, 3, 100)
fit_weib_100_tbl <- rweibull_cens_gen_fcn(100, 3, 100)
fit_weib_200_tbl <- rweibull_cens_gen_fcn(200, 3, 100)
fit_weib_300_tbl <- rweibull_cens_gen_fcn(300, 3, 100)
fit_weib_400_tbl <- rweibull_cens_gen_fcn(400, 3, 100)
fit_weib_500_tbl <- rweibull_cens_gen_fcn(500, 3, 100)
fit_weib_600_tbl <- rweibull_cens_gen_fcn(600, 3, 100)
fit_weib_700_tbl <- rweibull_cens_gen_fcn(700, 3, 100)
fit_weib_800_tbl <- rweibull_cens_gen_fcn(800, 3, 100)

Fit and save a model to each of the above data sets.

# mtf_700_weib_fit <- brm(time | cens(censor) ~ 1,
# data = fit_weib_700_tbl, family = weibull(),
#    prior = c(
#     prior(student_t(3, 5, 5), class = Intercept),
#     prior(uniform(0, 10), class = shape)
#    ),
#    iter = 41000, warmup = 40000, chains = 4, cores = 4,
#    seed = 4
#  )

Draw from the posterior of each model and combine into one tibble along with the original fit from n=30.

mtf_30_weib_fit_post_draws <- post_mwnp_samples # this is original data fit to milely informed priors, see Appendix
mtf_60_weib_fit_post_draws <- posterior_samples(mtf_60_weib_fit)
mtf_100_weib_fit_post_draws <- posterior_samples(mtf_100_weib_fit)
mtf_200_weib_fit_post_draws <- posterior_samples(mtf_200_weib_fit)
mtf_300_weib_fit_post_draws <- posterior_samples(mtf_300_weib_fit)
mtf_400_weib_fit_post_draws <- posterior_samples(mtf_400_weib_fit)
mtf_500_weib_fit_post_draws <- posterior_samples(mtf_500_weib_fit)
mtf_600_weib_fit_post_draws <- posterior_samples(mtf_600_weib_fit)
mtf_700_weib_fit_post_draws <- posterior_samples(mtf_700_weib_fit)
mtf_800_weib_fit_post_draws <- posterior_samples(mtf_800_weib_fit)


combined_n_tbl <- bind_rows(
  "60"  = mtf_60_weib_fit_post_draws,
  "100" = mtf_100_weib_fit_post_draws,
  "200" = mtf_200_weib_fit_post_draws,
  "300" = mtf_300_weib_fit_post_draws,
  "400" = mtf_400_weib_fit_post_draws,
  "500" = mtf_500_weib_fit_post_draws,
  "600" = mtf_600_weib_fit_post_draws,
  "700" = mtf_700_weib_fit_post_draws,
  "800" = mtf_800_weib_fit_post_draws,
  .id = "model"
) %>%
  mutate(model = as_factor(model)) %>%
  mutate(scale = exp(b_Intercept) / (gamma(1 + 1 / shape))) %>%
  select(shape, scale, model)

# have to add in the n=30 set seperately since it is already converted from Intercept to shape
combined_n_tbl <- mtf_30_weib_fit_post_draws %>%
  mutate(model = "30") %>%
  bind_rows(combined_n_tbl) %>%
  mutate(model = as_factor(model))

Finally we can visualize the effect of sample size on precision of posterior estimates.

combined_n_plt <- combined_n_tbl %>%
  ggplot(aes(x = shape, y = scale)) +
  geom_point(aes(color = model, fill = model), size = 2, alpha = 0.01) +
  scale_color_viridis_d() +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) + # force legend icons to be alpha = 1 instead of .05
  geom_hline(yintercept = 100, colour = "black", alpha = 0.3, size = .5) +
  geom_vline(xintercept = 3, colour = "black", alpha = 0.3, size = .5) +
  labs(
    title = "Effect of Sample Size on Parameter Estimation",
    subtitle = "Model: Weibull with Censoring; True ~ Weibull(shape = 3, scale = 100)",
    x = expression(eta ["shape"]),
    y = expression(beta ["scale"])
  ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  )

cens_marg_2_plot <- ggMarginal(combined_n_plt,
  groupColour = TRUE,
  groupFill = TRUE,
  type = "density",
  alpha = 0.7
)

This figure tells a lot. We simply needed more data points to zero in on the true data generating process. At n=30, there’s just a lot of uncertainty due to the randomness of sampling.

This plot looks really cool, but the marginal distributions are bit cluttered. This is a perfect use case for ggridges which will let us see the same type of figure but without overlap.

cn1 <- combined_n_tbl %>%
  gather(key = "key", value = "value", -model) %>%
  mutate(true_p = case_when(
    key == "scale" ~ 100,
    TRUE ~ 3
  )) %>%
  ggplot(aes(x = value, y = model, group = model, fill = model)) +
  geom_density_ridges(size = 1 / 3, rel_min_height = .005, alpha = .65) +
  geom_vline(aes(xintercept = true_p), color = "#711a6eff") +
  stat_pointintervalh(
    point_interval = mode_hdi,
    .width = .95,
    show.legend = FALSE
  ) +
  coord_flip() +
  scale_fill_viridis_d() +
  facet_wrap(~key, scales = "free_y") +
  labs(
    x = "Parameter Value",
    y = "Number of Data Points Used to Fit Model",
    title = "Model Fitting Variability due to Random Sampling [New Draw Each Fit]",
    subtitle = "True Distribution ~ Weibull(3, 100)"
  ) +
  theme(legend.position = "")

cn1

Set of 800 to demonstrate Bayesian updating.

set.seed(125)
fit_weib_800_bu_tbl <- rweibull_cens_gen_fcn(800, 3, 100)

fit_weib_30_bu_tbl <- fit_weib_800_bu_tbl[1:30, ]
fit_weib_60_bu_tbl <- fit_weib_800_bu_tbl[1:60, ]
fit_weib_100_bu_tbl <- fit_weib_800_bu_tbl[1:100, ]
fit_weib_200_bu_tbl <- fit_weib_800_bu_tbl[1:200, ]
fit_weib_300_bu_tbl <- fit_weib_800_bu_tbl[1:300, ]
fit_weib_400_bu_tbl <- fit_weib_800_bu_tbl[1:400, ]
fit_weib_500_bu_tbl <- fit_weib_800_bu_tbl[1:500, ]
fit_weib_600_bu_tbl <- fit_weib_800_bu_tbl[1:600, ]
fit_weib_700_bu_tbl <- fit_weib_800_bu_tbl[1:700, ]
fit_weib_800_bu_tbl <- fit_weib_800_bu_tbl[1:800, ]

Fit a model the first set of n=30

# mtf_30_bu_weib_fit <- brm(time | cens(censor) ~ 1,
# data = fit_weib_30_bu_tbl, family = weibull(),
#    prior = c(
#     prior(student_t(3, 5, 5), class = Intercept),
#     prior(uniform(0, 10), class = shape)
#    ),
#    iter = 41000, warmup = 40000, chains = 4, cores = 4,
#    seed = 4
#  )

We use the update() function in brms to update and save each model with additional data.

# mtf_60_bu_weib_fit <- update(mtf_30_bu_weib_fit, newdata = fit_weib_60_bu_tbl)
# mtf_100_bu_weib_fit <- update(mtf_60_bu_weib_fit, newdata = fit_weib_100_bu_tbl)
# mtf_200_bu_weib_fit <- update(mtf_100_bu_weib_fit, newdata = fit_weib_200_bu_tbl)
# mtf_300_bu_weib_fit <- update(mtf_200_bu_weib_fit, newdata = fit_weib_300_bu_tbl)
# mtf_400_bu_weib_fit <- update(mtf_300_bu_weib_fit, newdata = fit_weib_400_bu_tbl)
# mtf_500_bu_weib_fit <- update(mtf_400_bu_weib_fit, newdata = fit_weib_500_bu_tbl)
# mtf_600_bu_weib_fit <- update(mtf_500_bu_weib_fit, newdata = fit_weib_600_bu_tbl)
# mtf_700_bu_weib_fit <- update(mtf_600_bu_weib_fit, newdata = fit_weib_700_bu_tbl)
# mtf_800_bu_weib_fit <- update(mtf_700_bu_weib_fit, newdata = fit_weib_800_bu_tbl)

# save each one
# saveRDS(mtf_800_bu_weib_fit, file = "mtf_800_bu_weib_fit.rds")

mtf_800_bu_weib_fit <- readRDS(file = "mtf_800_bu_weib_fit.rds")
mtf_700_bu_weib_fit <- readRDS(file = "mtf_700_bu_weib_fit.rds")
mtf_600_bu_weib_fit <- readRDS(file = "mtf_600_bu_weib_fit.rds")
mtf_500_bu_weib_fit <- readRDS(file = "mtf_500_bu_weib_fit.rds")
mtf_400_bu_weib_fit <- readRDS(file = "mtf_400_bu_weib_fit.rds")
mtf_300_bu_weib_fit <- readRDS(file = "mtf_300_bu_weib_fit.rds")
mtf_200_bu_weib_fit <- readRDS(file = "mtf_200_bu_weib_fit.rds")
mtf_100_bu_weib_fit <- readRDS(file = "mtf_100_bu_weib_fit.rds")
mtf_60_bu_weib_fit <- readRDS(file = "mtf_60_bu_weib_fit.rds")
mtf_30_bu_weib_fit <- readRDS(file = "mtf_30_bu_weib_fit.rds")

Extract posterior draws for each one.

mtf_30_bu_weib_fit_post_draws <- posterior_samples(mtf_30_bu_weib_fit)
mtf_60_bu_weib_fit_post_draws <- posterior_samples(mtf_60_bu_weib_fit)
mtf_100_bu_weib_fit_post_draws <- posterior_samples(mtf_100_bu_weib_fit)
mtf_200_bu_weib_fit_post_draws <- posterior_samples(mtf_200_bu_weib_fit)
mtf_300_bu_weib_fit_post_draws <- posterior_samples(mtf_300_bu_weib_fit)
mtf_400_bu_weib_fit_post_draws <- posterior_samples(mtf_400_bu_weib_fit)
mtf_500_bu_weib_fit_post_draws <- posterior_samples(mtf_500_bu_weib_fit)
mtf_600_bu_weib_fit_post_draws <- posterior_samples(mtf_600_bu_weib_fit)
mtf_700_bu_weib_fit_post_draws <- posterior_samples(mtf_700_bu_weib_fit)
mtf_800_bu_weib_fit_post_draws <- posterior_samples(mtf_800_bu_weib_fit)

Combine into single tibble and convert intercept to scale. Some data wrangling is in anticipation for ggplot().

combined_nbu_tbl <- bind_rows(
  "30"  = mtf_30_bu_weib_fit_post_draws, # note: this one is the original dataset
  "60"  = mtf_60_bu_weib_fit_post_draws,
  "100" = mtf_100_bu_weib_fit_post_draws,
  "200" = mtf_200_bu_weib_fit_post_draws,
  "300" = mtf_300_bu_weib_fit_post_draws,
  "400" = mtf_400_bu_weib_fit_post_draws,
  "500" = mtf_500_bu_weib_fit_post_draws,
  "600" = mtf_600_bu_weib_fit_post_draws,
  "700" = mtf_700_bu_weib_fit_post_draws,
  "800" = mtf_800_bu_weib_fit_post_draws,
  .id = "model"
) %>%
  mutate(model = as_factor(model)) %>%
  mutate(scale = exp(b_Intercept) / (gamma(1 + 1 / shape))) %>%
  select(shape, scale, model)

cn2 <- combined_nbu_tbl %>%
  gather(key = "key", value = "value", -model) %>%
  mutate(true_p = case_when(
    key == "scale" ~ 100,
    TRUE ~ 3
  )) %>%
  ggplot(aes(x = value, y = model, group = model, fill = model)) +
  geom_density_ridges(size = 1 / 3, rel_min_height = .005, alpha = .65) +
  geom_vline(aes(xintercept = true_p), color = "#711a6eff") +
  stat_pointintervalh(
    point_interval = mode_hdi,
    .width = .95,
    show.legend = FALSE
  ) +
  coord_flip() +
  scale_fill_viridis_d() +
  facet_wrap(~key, scales = "free_y") +
  labs(
    x = "Parameter Value",
    y = "Number of Data Points Used to Fit Model",
    title = "Model Fitting Variability due to Random Sampling [Updated Model Each Fit]",
    subtitle = "True Distribution ~ Weibull(3, 100)"
  ) +
  theme(legend.position = "")

cn1

#cn2

The precision increase here is more smooth since supplemental data is added to the original set instead of just drawing completely randomly for each sample size. This is Bayesian updating.

Evaluate Sensitivity of Reliability Estimate to Sample Size.

To wrap things up, we should should translate the above figures into a reliability metric because that is the prediction we care about at the end of the day. I chose an arbitrary time point of t=40 to evaluate the reliability.

combined_nbu_rr_tbl <- combined_nbu_tbl %>% mutate(reliability_at_t40 = exp(-(40 / scale)**(shape)))

rq05_tbl <- combined_nbu_rr_tbl %>% group_by(model) %>%
  summarize(rel_q05 = quantile(reliability_at_t40, .05)) %>%
  mutate(model = as.numeric(model)) %>%
  ungroup()

reliability_ridge_plt <- combined_nbu_rr_tbl %>%
  gather(key = "key", value = "value", -c(model, reliability_at_t40)) %>%
  ggplot(aes(x = reliability_at_t40, y = model, group = model, fill = model)) +
  geom_density_ridges(size = 1 / 3, rel_min_height = .005, alpha = .8) +
  stat_pointintervalh(
    point_interval = mode_hdi,
    .width = .95,
    show.legend = FALSE
  ) +
  coord_flip() +
  scale_fill_viridis_d() +
  labs(
    x = "Reliability at t=40",
    y = "Number of Data Points Used to Fit Model",
    title = "Uncertainty in Reliability Posterior",
    subtitle = "Reliability Assessed at t = 40"
  ) +
  theme(legend.position = "") +
  xlim(c(.85, 1))



reliability_ridge_plt