Obtain the Frequentist Point Estimates

We can get the best parameter estimates for both data sets using the survreg function from the survival package.

#fit stress-in-service data using survreg from survival package
stress_in_service_fit <- survreg(Surv(stress_in_service) ~ 1,
  data = stress_in_service_tbl,
  dist = "lognormal"
)

#extract point estimates of parameters from sis-fit
sis_point_est_tbl <- tidy(stress_in_service_fit)[1, 2] %>%
  rename(meanlog = estimate) %>%
  mutate(sdlog = stress_in_service_fit$scale) %>%
  round(2) %>%
  kable(align = rep("c", 2))

sis_point_est_tbl 

#fit failure stress data using survreg from survival package
failure_stress_fit <- survreg(Surv(failure_stress, censored_surv) ~ 1,
  data = failure_stress_tbl,
  dist = "weibull"
)

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

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

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

fs_point_est_tbl

The reliability point estimate is obtained by drawing randomly from the two fitted distributions and seeing the percentage of occasions where the stress_in_service is greater than the failure_stress:

# Monte Carlo to see how often s-i-s > fs
set.seed(10)

#random draws
sis_draws <- rlnorm(n = 100000, meanlog = .88, sdlog = .50)
fs_draws <- rweibull(n = 100000, shape = 3.57, scale = 12)

#assign 1 to cases where sis_draws >= fs_draws
point_sim <- tibble(
  sis_draws = sis_draws,
  fs_draws = fs_draws
) %>%
  mutate(freq = case_when(
    sis_draws >= fs_draws ~ 1,
    TRUE ~ 0
  ))

#take freqency of 0's
reliability_pt_est <- (1 - mean(point_sim$freq))

#show as tibble
tibble(reliability = 1 - mean(point_sim$freq)) %>%
  round(3) %>% 
  kable(align = "c")

Model the Stress-in-Service with brms

For the Bayesian approach we fit the models with brms instead of survreg. The result is a posterior of plausible values for each parameter.

Before running to model, reasonable priors were established through simulation. Code and details are included in the Appendix at the end of this post so as to not derail the flow.

#Fit model to stress-in-service data. Data is known to be of lognormal form.  

# stress_in_service_model_1 <-
#  brm(
#    data = stress_in_service_tbl, family = lognormal,
#    stress_in_service ~ 1,
#   prior = c(
#     prior(normal(.5, 1), class = Intercept),
#     prior(uniform(.01, 8), class = sigma)),
#   iter = 41000, warmup = 40000, chains = 4, cores = 4,
#   seed = 4
# )

Clean up the posterior tibble and plot.

# extract posterior draws 
post_samples_stress_in_service_model_1_tbl <-
  posterior_samples(stress_in_service_model_1) %>%
  select(-lp__) %>%
  rename("mu" = b_Intercept)

#examine as tibble
post_samples_stress_in_service_model_1_tbl %>%
  head(7) %>%
  kable(align = rep("c", 2), digits = 3)

# get visual of posterior with rough idea of chain convergence
plot(stress_in_service_model_1)

Here is the summary of the stress-in-service model:

# evaluate posterior distribution with 95% CI and rhat diagnostic
summary(stress_in_service_model_1)
##  Family: lognormal 
##   Links: mu = identity; sigma = identity 
## Formula: stress_in_service ~ 1 
##    Data: stress_in_service_tbl (Number of observations: 34) 
## Samples: 4 chains, each with iter = 41000; warmup = 40000; 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     0.88      0.09     0.71     1.06 1.00     1976     2298
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.53      0.07     0.42     0.69 1.00     2399     2016
## 
## 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).

Model the Failure Stress Data with brms

The failure stress data is fit in a similar way as before. Again, prior predictive simulations are shown in the Appendix.

# failure_stress_model_1 <- brm(failure_stress | cens(censored_brms) ~ 1,
# data = failure_stress_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
# )

The following code extracts and converts the parameters from the brms default into the shape and scale that are used in the rweibull() function before displaying the summaries.

# extract posterior draws and examine as tibble
failure_stress_model_1_tbl <-
  posterior_samples(failure_stress_model_1) %>%
  select(-lp__) %>%
  rename("mu" = b_Intercept)

#compute shape and scale
post_samples_failure_stress_model_1_tbl <- posterior_samples(failure_stress_model_1) %>%
  mutate(scale = exp(b_Intercept) / (gamma(1 + 1 / shape))) %>%
  select(shape, scale)

#display as tibble
post_samples_failure_stress_model_1_tbl %>%
  head(7) %>%
  kable(align = rep("c", 2), digits = 3)

plot(failure_stress_model_1)

summary(failure_stress_model_1)
##  Family: weibull 
##   Links: mu = log; shape = identity 
## Formula: failure_stress | cens(censored_brms) ~ 1 
##    Data: failure_stress_tbl (Number of observations: 34) 
## Samples: 4 chains, each with iter = 41000; warmup = 40000; 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     2.38      0.06     2.28     2.49 1.00     2257     2279
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     3.56      0.55     2.56     4.67 1.00     1963     2169
## 
## 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).

Visualization of Uncertainty - Credible Curves for Stress-in-Service and Failure Stress

I haven’t ever used the gganimate package and this seems like a nice opportunity. The code below draws a small handful of parameters from the posterior and plots them to visualize the uncertainty in both distributions.

#take 25 sets of parameters, convert to lnorm curves
lnorm_stress_curve_tbl <- post_samples_stress_in_service_model_1_tbl[1:25, ] %>%
  mutate(plotted_y_data = map2(
    mu, sigma,
    ~ tibble(
      x = seq(0, 20, length.out = 100),
      y = dlnorm(x, .x, .y)
    )
  )) %>%
  unnest(plotted_y_data) %>%
  mutate(model = "Stress in Service [lnorm]") %>%
  rename(
    param_1 = mu,
    param_2 = sigma
  )

#take 25 sets of parameters, convert to Weib curves
weib_stress_curve_tbl <- post_samples_failure_stress_model_1_tbl[1:25, ] %>%
  mutate(plotted_y_data = map2(
    shape, scale,
    ~ tibble(
      x = seq(0, 20, length.out = 100),
      y = dweibull(x, .x, .y)
    )
  )) %>%
  unnest(plotted_y_data) %>%
  mutate(model = "Failure Stress [weib]") %>%
  rename(
    param_1 = shape,
    param_2 = scale
  )

#combine
a <- bind_rows(lnorm_stress_curve_tbl, weib_stress_curve_tbl) %>% mutate(param_1_fct = as_factor(param_1))

#visualize
p <- a %>%
  ggplot(aes(x, y)) +
  geom_line(aes(x, y, group = param_1_fct, color = model), alpha = 1, size = 1) +
  labs(
    x = "Stress",
    y = "Density",
    title = "Credible Failure Stress and Service Stress Distributions",
    subtitle = "n=25 curves sampled from the posterior"
  ) +
  scale_color_viridis_d(end = .8) +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  theme(legend.title = element_blank()) +
  theme(legend.position = "bottom") +
  transition_states(param_1_fct, 0, 1) +
  shadow_mark(past = TRUE, future = TRUE, alpha = .3, color = "gray50", size = .4)

#gganimate 
animate(p, nframes = 50, fps = 2.5, width = 900, height = 600, res = 120, dev = "png", type = "cairo")

Building the Credible Reliability Distribution

Having now obtained the posterior distributions for both stress-in-service and failure stress, we can select random sets of parameters and compare a random (but credible) pair of distributions. By simulating from each random pair of distributions and calculating a single value of reliability as before, we can build out a credible reliability distribution. The blue vertical line indicates the frequentist point estimate we obtained at the beginning of the analysis.

#set number of simulations
n_sims <- 10000
set.seed(1001)

#stress-in-service (lognormal) simulations
labeled_post_ln_tbl <- post_samples_stress_in_service_model_1_tbl %>%
  mutate(
    model = "lognormal"
  ) %>%
  rename(
    param1 = mu,
    param2 = sigma
  ) %>%
  mutate(nested_data_ln = map2(param1, param2, ~ rlnorm(n_sims, .x, .y)))

#failure stress (Weibull) simulations
labeled_post_wb_tbl <- post_samples_failure_stress_model_1_tbl %>%
  mutate(
    model = "weibull"
  ) %>%
  rename(
    param1 = shape,
    param2 = scale
  ) %>%
  mutate(nested_data_wb = map2(param1, param2, ~ rweibull(n_sims, .x, .y)))

#combine and calculate reliability for each pair to build reliability distribution
all_post_samples_tbl <- bind_cols(labeled_post_ln_tbl, labeled_post_wb_tbl) %>%
  select(nested_data_ln, nested_data_wb) %>%
  mutate(reliability = map2_dbl(nested_data_ln, nested_data_wb, ~ mean(.x < .y)))

Visualize the results with some help from tidybayes::geom_halfeyeh()

#visualize
all_post_samples_tbl %>%
  ggplot(aes(x = reliability, y = 0)) +
  geom_halfeyeh(
    fill = "firebrick",
    point_interval = median_qi, .width = .95, alpha = .9
  ) +
  geom_vline(xintercept = reliability_pt_est, color = "dodgerblue", size = 1.1, alpha = .7) +
  #  stat_dotsh(quantiles = 100, size = .5) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(
    title = "Distribution of Predicted Reliability",
    subtitle = "Marked by median and 95% probability interval. Vertical line is the point estimate"
  ) +
  theme_bw() +
  theme(panel.grid = element_blank())

Behold, a full reliability distribution supported by the data! So much better for decision making than the point estimate!

Thank you for reading.

Appendix - Prior Predictive Simulation

set.seed(45)
mu_prior <- rlnorm(100000, meanlog = .5, sdlog = 1)
mu_prior_tbl <- mu_prior %>%
  as_tibble() %>%
  filter(value > 0)

muuuuu <- mu_prior_tbl %>% ggplot(aes(x = mu_prior)) +
  geom_histogram(aes(y = ..density..), fill = "#2c3e50", color = "white", alpha = .6) +
  scale_x_continuous(trans = "log10")

mu_prior_tbl %>%
  mutate(orig = log(value)) %>%
  pull(orig) %>%
  mean()
## [1] 0.5033306
mu_prior_tbl %>%
  mutate(orig = log(value)) %>%
  pull(orig) %>%
  sd()
## [1] 1.001799
set.seed(45)
sigma_prior <- runif(100000, .01, 8)

p0_priors_tbl <- sigma_prior %>%
  as_tibble() %>%
  bind_cols(mu_prior_tbl) %>%
  rename(sigma = value, mu = value1)


sigmaaa <- p0_priors_tbl %>% ggplot(aes(x = sigma_prior)) +
  geom_histogram(aes(y = ..density..), fill = "#2c3e50", color = "white", alpha = .6)

muuuuu + sigmaaa + plot_annotation(title = "Prior Predicitve Simulations for mu and sigma")

p0 <- p0_priors_tbl[1:1000, ] %>%
  mutate(row_id = row_number()) %>%
  mutate(plotted_y_data = pmap(
    list(sigma, mu, row_id),
    ~ tibble(
      x = seq(.1, 100, length.out = 1000),
      y = dlnorm(x, .x, .y),
      z = row_id
    )
  )) %>%
  unnest(plotted_y_data) %>%
  filter(x > 1) %>%
  ggplot(aes(x, y)) +
  geom_line(aes(group = row_id), alpha = .15, color = "#2c3e50") +
  labs(
    x = "Stress-in-Service",
    y = "Density",
    title = "Implied Stress-in-Service Possibilities",
    subtitle = "Generated from Priors Only"
  ) +
  scale_x_continuous(trans = "log10") +
  ylim(c(0, 1))

p0

# seed for reproducibility
set.seed(12)

# Evaluate Mildly Informed Priors
shape_prior <- runif(100000, 0, 10)
shape_prior_tbl <- shape_prior %>% as_tibble()
shaaaape <- shape_prior_tbl %>% ggplot(aes(x = shape_prior)) +
  geom_histogram(aes(y = ..density..), binwidth = 1, boundary = 10, fill = "#2c3e50", color = "white", alpha = .6)


intercept_prior <- rstudent_t(100000, 3, 5, 5)

priors_tbl <- intercept_prior %>%
  as_tibble() %>%
  bind_cols(shape_prior_tbl) %>%
  rename(intercept = value, shape = value1) %>%
  mutate(scale_prior = exp(intercept) / (gamma(1 + 1 / shape))) %>%
  filter(scale_prior < 1000) %>%
  select(-intercept)

scaaaale <- priors_tbl %>% ggplot(aes(x = scale_prior)) +
  geom_histogram(aes(y = ..density..), binwidth = 10, boundary = 100, fill = "#2c3e50", color = "white", alpha = .6) +
  ylim(c(0, .005))

shaaaape + scaaaale + plot_annotation(title = "Prior Predicitve Simulations for Shape and Scale")

p1 <- priors_tbl[1:500, ] %>%
  mutate(plotted_y_data = map2(
    shape, scale_prior,
    ~ tibble(
      x = seq(0, 200, length.out = 400),
      y = dweibull(x, .x, .y)
    )
  )) %>%
  unnest(plotted_y_data) %>%
  ggplot(aes(x, y)) +
  geom_line(aes(group = shape), alpha = .2, color = "#2c3e50") +
  xlim(c(0, 50)) +
  ylim(c(0, .5)) +
  labs(
    x = "Failure Stress Distributions",
    y = "Density",
    title = "Implied Failure Stress Possibilities",
    subtitle = "Generated from Priors Only"
  )

p1