Power Analysis for a DV Test – Frequentist and Bayesian Estimation in R

Simulate the Data

reps <- 1000
set.seed(127)
n <- 9


single_sim_tbl <- 
  tibble(
  n = n,
  mean = sum_data$mean, #mean from predicate data
  sd = sum_data$sd) %>% #ds from predicate data
  slice(rep(1:n(), each = reps)) %>%
  mutate(replicate = seq(from = 1, to = reps)) %>%
  mutate(sim_data = pmap(.l = list(n = n, mean = mean, sd = sd), .f = rnorm)) %>% #create the random data points, n=9 per sim
  unnest() %>%
  group_by(mean, sd, replicate) %>%
  summarize(ltl = normtol.int(sim_data, P = .9, alpha = .1)) %>% #calculate 90/90 tolerance interval
  unnest() %>%
  ungroup() %>%
  rename(ltl = "1-sided.lower") %>%
  select(1:7) %>%
  mutate(fail = case_when( 
    ltl <= 300 ~ 1,
    TRUE ~ 0)) %>% #flag failures
  mutate(across(c(x.bar, ltl), round, digits = 1))


single_sim_tbl %>%
  gt_preview()

Calculate Frequency of Success

Having generated and stored the above simulation, it is straightforward to calculate the predicted success rate. We see that the long run frequency of success is 84%, conditional on all of our strict assumptions about population parameters and fixed sample size:

single_sim_tbl %>%
  summarize(total_simulated_experiments = n(),
            sum_pass = n() - sum(fail)) %>%
  mutate(prop_pass = sum_pass/total_simulated_experiments) %>%
  gt_preview()

Visualize

Here we look at the means and lower tolerance limits (LTLs) from the simulation to help us gut check the conclusions and further understand what is going on.

single_sim_tbl %>% 
  pivot_longer(cols = c(x.bar, ltl), names_to = "param", values_to = "value") %>%
  mutate(param = case_when(param == "x.bar" ~ "mean",
                           TRUE ~ "lower_tolerance_limit")) %>%
  ggplot(aes(x = param, y = value)) +
  geom_line(aes(group = replicate), size = .1, alpha = .1) +
  geom_jitter(aes(color = param), size = .4, alpha = .4, width = .01) +
  geom_hline(yintercept = 300, color = "firebrick", linetype = 2) +
  theme_classic() +
  theme(legend.title=element_blank()) +
  labs(title = "Means and Lower Tolerance Limits for 1000 Simulated Experiments",
       subtitle = "Each sim: n=9 data points draws from a normal distribution",
       x = "",
       y = "Eb - Er [mV SCE]",
       caption = "red dotted line indicates lower spec limit")

And here we look by simulation run number and add the marginal distributions for mean and sd.

int_plot <- single_sim_tbl %>%
  ggplot(aes(group = replicate)) +
  geom_segment(aes(x = replicate, y = x.bar, xend = replicate, yend = ltl), size = .2, alpha = .2) +
  geom_point(data = single_sim_tbl %>% 
               pivot_longer(cols = c(x.bar, ltl), names_to = "param", values_to = "value") %>%
               mutate(param = case_when(param == "x.bar" ~ "Sample Mean", TRUE ~ "1-Sided Lower Tolerance Limit")), 
               aes(x = replicate, y = value, color = param), size = .6) +
  geom_hline(yintercept = 300, linetype = 2, color = "firebrick") +
  theme_classic() +
  theme(legend.position = "bottom") +
  labs(y = "Eb - Er (mV SCE)",
       x = "Simulated Experiment Replicate Number",
       title = "Means and Lower Tolerance Limits for 1000 Simulated Experiments",
       subtitle = "Each sim: n=9 data points draws from a normal distribution") +
  theme(legend.title=element_blank())

ggMarginal(int_plot,
  groupColour = TRUE,
  groupFill = TRUE,
  type = "density",
  alpha = 0.7,
  margins = "y"
)

Simulated the Data

The workflow below is very similar to the single simulation section with a few exceptions. The expand_grid() function is what generated all combinations of the specified sequences for mean, sd, and n. slice() is used to copy and append rows, and a “scenario” column is added to track each permutation of the 3 variables.

set.seed(2015)
reps <-  1000

mean = c(seq(from = 300, to = 500, by = 50), sum_data$mean)
sd = c(seq(from = 75, to = 175, by = 25), sum_data$sd)
n = seq(from = 5, to = 17, by = 2)

full_even_sim_tbl <-  expand_grid(mean, sd, n) %>%
  mutate(mean = mean %>% round(0),
         sd = sd %>% round(0)) %>%
  mutate(scenario = seq(from = 1, to = nrow(.))) %>%
  slice(rep(1:n(), each = reps)) %>%
  mutate(replicate_number = rep(seq(from = 1, to = reps), nrow(.)/reps)) %>%
  mutate(sim_data = pmap(.l = list(n = n, mean = mean, sd=sd), .f = rnorm)) %>%
  unnest(cols = everything()) %>%
  group_by(mean, sd, scenario, replicate_number, n) %>%
  summarize(ltl = normtol.int(sim_data, P = .9, alpha = .1)) %>%
  unnest(cols = everything()) %>%
  ungroup() %>%
  rename(ltl = '1-sided.lower',
         utl = '1-sided.upper') %>%
  mutate(fail = case_when(ltl <= 300 ~ 1,
         TRUE ~ 0))

full_even_sim_tbl %>% 
  mutate(across(c(x.bar, ltl, utl), round, digits = 1)) %>%
  gt_preview()

Calculate Frequency of Success

Now, each scenario (combination of assumed mean, sd, and sample size) has its own probability (long run frequency) of passing:

summary_tbl <- full_even_sim_tbl %>%
  group_by(mean, sd, n, max(replicate_number)) %>%
  summarize(number_fail = sum(fail)) %>%
  rename(reps = 'max(replicate_number)') %>%
  mutate(avg = stringr::str_glue("mean Eb-Er = {as_factor(mean)} mV [SCE]")) %>%
  mutate(prop_pass = 1-(number_fail / reps)) %>%
  ungroup()

summary_tbl %>%
  gt_preview()

Visualize – Power Curves

Plotting these data produce so-called power curves from which we can easily see trends for probability of passing based on the parameter assumptions and our choice of sample size. Constructing these curves provide a principled way to select a reasonable sample size.

summary_tbl %>%
  ggplot(aes(x = sd, y = prop_pass)) +
  geom_point(aes(color = as_factor(n))) +
  geom_line(aes(color = as_factor(n))) +
  geom_vline(xintercept = sum_data$sd, color = "firebrick", linetype = 2) +
  facet_wrap(~ avg) +
  lims(x = c(75, 175)) +
  scale_color_discrete(name = "Sample Size") +
#scale_color_viridis_d(option = "C", end = .8) +
  labs(x = "Standard Deviation (mV [SCE])",
       y = str_glue("Passing Proportion\nBased on n={reps} simulations each")) +
  theme_bw() +
  labs(title = "Predicted Power Curves for Simulated Corrosion Experiments",
       subtitle = "as a function of mean, standard deviation, and sample size",
       caption = "513 mV facet represents predicate mean, red dotted line represents predicate sd")

We note that at the lower means and higher sds we would have very low chance of passing.

Specify Priors

We put these into this table form because we’ll want to plot these curves to confirm visually that our assumptions make sense.

prior_pred_tbl <- tibble(
  mu = rstudent_t(n = 300, df = 3, mu = 500, sigma = 100),
  sig = runif(300, min = 20, max = 200) %>% abs()
) %>%
  mutate(row_id = row_number()) %>%
  select(row_id, everything()) %>%
  mutate(plotted_y_data = map2(mu, sig, ~ tibble(
    x = seq(-200, 1000, length.out = 10000),
    y = dnorm(x, .x, .y)
  ))) %>%
#  slice(1) %>%
  unnest() %>%
  mutate(model = "mu ~ student(3, 500, 100), sigma ~ unif(20,200)")

prior_pred_tbl %>%
  gt_preview()

Visualize Prior Pred

This vis shows a variety of possible outcomes for Eb-Er before seeing the data.

prior_pred_tbl %>%
  ggplot(aes(x = x, y = y, group = row_id)) +
  geom_line(aes(x, y), size = .5, alpha = .2, color = "#2c3e50") +
  labs(
    x = "mV [SCE]",
    y = "",
    title = "Prior Predictive Simulations"
  ) +
  theme_minimal() + theme(
  axis.text.y = element_blank(),
  axis.ticks.y = element_blank())

Specify Data

We re-enter our legacy data.

p_tbl <- tibble(p = c(424, 543, 571))

Fit Model with brms

Having specified our data and priors, we can assume a likelihood and brms can do the rest. Here we assume Gaussian (normal) for the likelihood of the data and copy our priors from the previous section. Once fit, we’ll have access to a posterior in table form comprised of credible values for the mean and sd for many possible normal distributions. This is similar to the grid of values we build manually in the frequentist section to generate power curve, but now now they are weighted by plausibility.

pa_mod_1 <-
 brm(
   data = p_tbl, family = gaussian,
   p ~ 1,
  prior = c(
    prior(student_t(3, 500, 100), class = Intercept),
    prior(uniform(1, 200), class = sigma)),
  iter = 41000, warmup = 40000, chains = 4, cores = 4,
  seed = 4
)

summary(pa_mod_1)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: p ~ 1 
##    Data: p_tbl (Number of observations: 3) 
##   Draws: 4 chains, each with iter = 41000; warmup = 40000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   508.58     50.22   402.45   608.44 1.00     1314      988
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma   105.34     40.47    44.68   190.13 1.00     1198     1360
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Visualize Chains and Posterior

This section provided diagnostic plots for the Markov Chain Monte Carlo. Chains look well dispersed and relatively similar (good) with a “fuzzy caterpillar” look that we want.

plot(pa_mod_1)

Examine Posterior Draws

Let’s take a peek at the table of samples from the posterior. They look reasonable.

pa_mod_1_post_tbl <- posterior_samples(pa_mod_1) %>%
  select(-lp__) %>%
  rename("mu" = b_Intercept)

pa_mod_1_post_tbl %>%
  gt_preview()

Simulate the Data

Knowing that the posterior draws represent credible values population mean and sd, our goal is to produce a new simulation that conducts and a simulated experiment over every credible combination. This is very similar to our frequentist version above, but instead of just copying arbitrary combinations of parameters over a grid, here the posterior draws have already selected credible values for the parameters. This is sort of like a “weighted” version of the frequentist simulation that is compatible with our pre-existing beliefs as established by the priors and then updated though the predicate data.

set.seed(2015)
reps <-  1000

n = seq(from = 5, to = 17, by = 2)

b_sim_tbl <- pa_mod_1_post_tbl %>%
  rename(mean = mu,
         sd = sigma) %>%
  mutate(scenario = seq(from = 1, to = nrow(.))) %>%
  slice(rep(1:n(), each = length(n))) %>%
  mutate(num_draws = rep(x = n, times = nrow(.)/length(n))) %>%
  mutate(dat = pmap(.l = list(n = num_draws, mean = mean, sd=sd), .f = rnorm)) %>%
  unnest(cols = everything()) %>%
  group_by(mean, sd, scenario, num_draws) %>%
  summarize(ltl = normtol.int(dat, P = .9, alpha = .1)) %>%
  unnest(cols = everything()) %>%
  ungroup() %>%
  arrange(scenario) %>%
  rename(ltl = '1-sided.lower',
         utl = '1-sided.upper') %>%
  mutate(fail = case_when(ltl <= 300 ~ 1,
         TRUE ~ 0)) %>%
  select(-utl) %>%
  mutate(across(is.numeric, round, digits = 1))
  
b_sim_tbl %>%
  gt_preview()