Simulate the Strength of Each Part

set.seed(805)
stent_strength_tbl <- tibble(part_strength_mpa = rnorm(n = 25, mean = 600, sd = 13))

stent_strength_tbl %>%
  gt_preview()

See Test Levels and Input Test Results

The gonogo() function will initiate the test plan for sequential testing that prescribes test levels based on the outcome of each run. This operation happens in the console as shown in the screenshot below.

We’ll use the Neyer option by setting the “test” argument equal to 2. One key feature of sequential testing is that high and low initial guesses for the mean must be provided at the onset of testing along with an initial guess for the standard deviation. These guesses help set the scale. We’ll provide vague guesses that would be reasonable given our domain knowledge of cobalt-based alloys and their fatigue properties.

If we were really running the test on physical parts we’d want to set up a test at the prescribed stimulus level, then report the result to the console by inputing 2 numbers: the level you actually ran the test at (should be the same as what was recommended) and the outcome (1 for failure or 0 for survival) separated by a single space. After the results from the first test are entered, the recommended level of the next test is presented. When all runs are complete, the results are stored in the assigned object.

# z <- gonogo(mlo = 500, mhi = 700, sg = 10, test = 2, reso = 1) #input results in console

The resulting object z is a list that has a lot of information. The following code will pull the run history and clean it up nicely with tools and formatting from gt package.

z$d0 %>%
  gt() %>%
  cols_align("center") %>%
  opt_row_striping(row_striping = TRUE) %>%
  tab_header(
    title = "Run History of Sensitivity Test of n=25 Stents",
    subtitle = "X and EX are Stimulus Levels, Y = Outcome (1 = part failed, 0 = part survived)"
  ) %>%
  cols_width(everything() ~ px(100)) %>%
  tab_style(
    style = list(
      cell_text(weight = "bold")
    ),
    locations = cells_column_labels(columns = everything())
  )

Gonogo has some really convenient built-in visualization tools. The run history can be plotted using gonogo’s ptest function with argument plt=1.

ptest(z, plt = 1)

## NULL

The values we really care about are the estimates of the mean and sd of the latent failure strength distribution. Numerical MLE’s for these parameters can be quickly extracted as follows:

z$musig
## [1] 600.383131   7.894673

We see that these are quite close to the true parameters of mean = 600 MPa, sd = 13 MPa! That really impressive for only 25 parts and vague initial guesses.

I like to occasionally look under the hood and see how the outputs are generated. I don’t always understand it but it helps me learn. In this case, the calculations to extract an estimate are pretty simple. gonogo is using a glm with probit link function. The results are replicated here with manual calculations.

xglm <- glm(Y ~ X, family = binomial(link = probit), data = z$d0)

ab <- as.vector(xglm$coef)
muhat <- -ab[1] / ab[2]
sighat <- 1 / ab[2]

tibble(
  true_parameter = c(600, 13),
  gonogo_estimate = c(muhat, sighat)
) %>%
  mutate(gonogo_estimate = gonogo_estimate %>% round(1)) %>%
  gt()

Extract a tibble of 2-sided, pointwise confidence limits about all quantiles

This first function takes the output from the gonogo function and calculates 2-sided .90 confidence bounds on the various quantiles. The lims() function is from gonogo.

# helper function

conf_limits_fcn <- function(tbl, conf = .9) {
  conf_limits_tbl_full <- lims(ctyp = 2, dat = tbl, conf = conf, P = seq(from = .01, to = .99, by = .01)) %>%
    as_tibble() %>%
    rename(
      ciq_lower = V1,
      estimate_q = V2,
      ciq_upper = V3,
      ip_lower = V4,
      estimate_p = V5,
      cip_upper = V6
    )
}

rileysim_conf_tbl <- conf_limits_fcn(z$d0)

rileysim_conf_tbl %>%
  gt_preview()

Extract a point estimate for the lower 2-sided confidence limit for the .10 quantile.

We present results to FDA using the terminology “confidence and reliability”. By extracting the lower end of the 2-sided 90% confidence band on the .10 quantile, we obtain an estimate of the 95/90 tolerance bound i.e. a bound above which we would expect 90% or more of the population to lie, with 95% confidence. This would be a conservative estimate of the endurance limit for this stent design.

lower_conf_fcn <- function(tbl, percentile) {
  lower_scalar <- tbl %>%
    mutate(qprobs = qnorm(estimate_p)) %>%
    arrange(qprobs) %>%
    filter(estimate_p == percentile) %>%
    mutate(ciq_lower = ciq_lower %>% round(1)) %>%
    pluck(1)

  lower_qnorm <- tbl %>%
    mutate(qprobs = qnorm(estimate_p)) %>%
    arrange(qprobs) %>%
    filter(estimate_p == percentile) %>%
    pluck(7)

  lower_tbl <- tibble(
    estimate_q = lower_scalar,
    estimate_p = percentile,
    qprobs = lower_qnorm
  )

  return(lower_tbl)
}

lower_conf_fcn(rileysim_conf_tbl, percentile = .1) %>% pluck(1)
## [1] 580.2

Plot CDF

Now some plotting functions to help visualize things. First, we use our confidence limit dataframe from above to construct a cdf of failure strength. With some clever use of geom_label_repel() and str_glue() we can embed our data right onto the plot for easy communication. The estimated endurance limit from our test sequence is 580 MPa.

plot_cdf_fcn <- function(tbl, percentile = .1, units = "MPa") {
  tbl %>%
    ggplot(aes(x = ciq_lower, y = estimate_p)) +
    geom_line(size = 1.3) +
    geom_line(aes(x = ciq_upper, y = estimate_p), size = 1.3) +
    geom_line(aes(x = estimate_q, y = estimate_p), linetype = 2) +
    geom_ribbon(aes(xmin = ciq_lower, xmax = ciq_upper), alpha = .1, fill = "purple") +
    geom_point(data = lower_conf_fcn(tbl, percentile = percentile), aes(estimate_q, estimate_p), color = "purple", size = 2) +
    geom_label_repel(
      data = lower_conf_fcn(tbl, percentile = percentile),
      aes(estimate_q, estimate_p),
      label = str_glue("1-sided, 95% Confidence Bound:\n{lower_conf_fcn(rileysim_conf_tbl, percentile = .1) %>% pluck(1)} {units}"),
      fill = "#8bd646ff",
      color = "black",
      segment.color = "black",
      segment.size = .5,
      nudge_y = .13,
      nudge_x = -25
    ) +
    geom_segment(data = tbl %>%
      filter(estimate_p == percentile), aes(x = ciq_lower, xend = ciq_upper, y = estimate_p, yend = estimate_p), linetype = 3) +
    theme_bw() +
    labs(
      x = "Stress (MPa)",
      y = "Percentile",
      caption = "Dotted line is estimate of failure strength distribution\nConfidence bands are 2-sided 90% (i.e. 1-sided 95 for lower)",
      title = "Sensitivity Testing of Stent Durability",
      subtitle = str_glue("placeholder")
    )
}

a <- plot_cdf_fcn(tbl = rileysim_conf_tbl, percentile = .1)

a +
  labs(subtitle = str_glue("Estimates: Mean = {muhat %>% round(1)} MPa, Sigma = {sighat %>% round(1)} MPa"))

Normal Probability Plot

Engineers in my industry have an unhealthy obsession to normal probability plots. Personally I don’t like them, but Minitab presents them by default and engineers in my industry will look for data in this format. As part of this project I taught myself how to convert distributional data from a CFD into a normal probability plot. I show the helper function here and the result.

norm_probability_plt_fcn <- function(tbl, percentile, units = "MPa") {
  conf_tbl <- tbl %>%
    mutate(qprobs = qnorm(estimate_p)) %>%
    arrange(qprobs)

  probs_of_interest <- c(0.01, 0.05, seq(0.1, 0.9, by = 0.1), 0.95, 0.99)

  probs_of_interest_tbl <- conf_tbl %>%
    filter(estimate_p %in% probs_of_interest)

  conf_tbl %>%
    ggplot() +
    geom_line(aes(x = ciq_lower, y = qprobs), size = 1.3) +
    geom_line(aes(x = ciq_upper, y = qprobs), size = 1.3) +
    geom_line(aes(x = estimate_q, y = qprobs), linetype = 2) +
    geom_ribbon(aes(xmin = ciq_lower, xmax = ciq_upper, y = qprobs), alpha = .1, fill = "purple") +
    geom_point(data = lower_conf_fcn(tbl, percentile = percentile), aes(estimate_q, qprobs), color = "purple", size = 2) +
    geom_label_repel(
      data = lower_conf_fcn(tbl, percentile = percentile),
      aes(estimate_q, qprobs),
      label = str_glue("1-sided, 95% Confidence Bound:\n{lower_conf_fcn(rileysim_conf_tbl, percentile = .1) %>% pluck(1)} {units}"),
      fill = "#8bd646ff",
      color = "black",
      segment.color = "black",
      segment.size = .5,
      nudge_y = 1,
      nudge_x = -20
    ) +
    geom_segment(data = conf_tbl %>%
      filter(estimate_p == percentile), aes(x = ciq_lower, xend = ciq_upper, y = qprobs, yend = qprobs), linetype = 3) +
    scale_y_continuous(limits = range(conf_tbl$qprobs), breaks = probs_of_interest_tbl$qprobs, labels = 100 * probs_of_interest_tbl$estimate_p) +
    theme_bw() +
    labs(
      x = "Stress (MPa)",
      y = "Percentile",
      caption = "Dotted line is estimate of failure strength distribution\nConfidence bands are 2-sided 90%",
      title = "Sensitivity Testing of Stent Durability",
      subtitle = str_glue("TBD_2")
    )
}

b <- norm_probability_plt_fcn(rileysim_conf_tbl, percentile = .1)

b +
  labs(subtitle = str_glue("Estimates: Mean = {muhat %>% round(1)} MPa, Sigma = {sighat %>% round(1)} MPa"))