This R Markdown document contains the code for performing simulation studies related to the present method for setting acceptance criteria.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.1.2     ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(furrr)
## Loading required package: future
## 
## Attaching package: 'future'
## 
## The following object is masked from 'package:rmarkdown':
## 
##     run
library(cmstatr)
library(tikzDevice)

plan(multisession, workers = 12)
factors <- read_csv(
  "factors.csv",
  col_types = cols(
    alpha = col_double(),
    n = col_double(),
    m = col_double(),
    k1 = col_double(),
    k2 = col_double()
  ))
mu_pop <- 0
sd_pop <- 1

n_qual_1 <- 5000
n_equiv_1 <- 5000

n_qual_2 <- 2500
n_equiv_2 <- 2500

sim_alpha <- 0.05

delta_mean <- seq(0, 2, length.out = 9)
delta_sd <- seq(1, 5, length.out = 9)
delta_min_indiv <- seq(0, 4, length.out = 9)
delta_mixture <- seq(0, 4, length.out = 9)
pi.mixture <- 0.4
sim_n_values <- unique(factors$n[factors$n <= 100])
sim_m_values <- unique(factors$m[factors$m %% 2 == 0])
k_factors <- unique(factors$m) %>% 
  map_dfr(function(m) {
    k <- k_equiv(sim_alpha, m)
    tibble(
      m = m,
      k1 = k[1],
      k2 = k[2]
    )
  })
sim_equiv <- map_dfr(
  sim_m_values,
  function(cur_m) {
    k1 <- k_factors$k1[k_factors$m == cur_m]
    k2 <- k_factors$k2[k_factors$m == cur_m]
    
    future_map_dfr(
      sim_n_values,
      .options = furrr_options(seed = 1234),
      function(cur_n) {
        r1 <- factors$k1[factors$n == cur_n & factors$m == cur_m
                         & factors$alpha == sim_alpha]
        r2 <- factors$k2[factors$n == cur_n & factors$m == cur_m
                         & factors$alpha == sim_alpha]
        
        equiv_samples <- future_map_dfr(
          1:n_equiv_1,
          .options = furrr_options(seed = 2345),
          function(i_equiv) {
            x <- rnorm(cur_m, mu_pop, sd_pop)
            data.frame(
              min = min(x),
              avg = mean(x)
            )
          }
        )
        
        future_map_dfr(
          1:n_qual_1,
          .options = furrr_options(seed = 3456),
          function(i_qual) {
            x_qual <- rnorm(cur_n, mu_pop, sd_pop)
            avg_qual <- mean(x_qual)
            sd_qual <- sd(x_qual)
            
            rbind(
              data.frame(
                method = factor("Vangel", c("Two-Sample", "Vangel")),
                accept_min = equiv_samples$min > avg_qual - k1 * sd_qual,
                accept_avg = equiv_samples$avg > avg_qual - k2 * sd_qual
              ),
              data.frame(
                method = factor("Two-Sample", c("Two-Sample", "Vangel")),
                accept_min = equiv_samples$min > avg_qual - r1 * sd_qual,
                accept_avg = equiv_samples$avg > avg_qual - r2 * sd_qual
              )
            )
          }
        ) %>% 
          mutate(accept = accept_min & accept_avg) %>% 
          group_by(method) %>%
          summarise(`Rejection Rate` = sum(!accept) / n(),
                    .groups = "drop") %>%
          mutate(n = cur_n, m = cur_m)
      }
    )
  }
)

Power Simulation – Generic

# equiv_deviation has the arguments:
# i_equiv, delta, cur_m
power_simulation <- function(delta_list, qual_generation, equiv_generation) {
  map_dfr(
    sim_m_values,
    function(cur_m) {
      k1 <- k_factors$k1[k_factors$m == cur_m]
      k2 <- k_factors$k2[k_factors$m == cur_m]
      
      future_map_dfr(
        sim_n_values,
        .options = furrr_options(seed = 12345),
        function(cur_n) {
          r1 <- factors$k1[factors$n == cur_n & factors$m == cur_m
                           & factors$alpha == sim_alpha]
          r2 <- factors$k2[factors$n == cur_n & factors$m == cur_m
                           & factors$alpha == sim_alpha]
          
          t_crit <- qt(1 - sim_alpha / 2, cur_n - 1)
          f_crit <- qf(1 - sim_alpha / 2, cur_m - 1, cur_n - 1)
          
          future_map_dfr(
            delta_list,
            .options = furrr_options(seed = 23456),
            function(delta) {
              equiv_samples <- future_map_dfr(
                1:n_equiv_2,
                .options = furrr_options(seed = 34567),
                function(i_equiv) {
                  x <- equiv_generation(i_equiv, delta, cur_m)
                  data.frame(
                    min = min(x),
                    avg = mean(x),
                    sd_eq = sd(x)
                  )
                }
              )
              
              future_map_dfr(
                1:n_qual_2,
                .options = furrr_options(seed = 45678),
                function(i_qual) {
                  x_qual <- qual_generation(cur_n)
                  avg_qual <- mean(x_qual)
                  sd_qual <- sd(x_qual)
                  
                  # Ref. Escobar, pp. 55
                  t_min <- avg_qual - t_crit * sqrt(1 / cur_n + 1 / cur_m) * sd_qual
                  
                  # Ref. Escobar, pp. 58
                  f_max <- sd_qual * sqrt(f_crit)
                  
                  accept_min <- equiv_samples$min > avg_qual - r1 * sd_qual
                  accept_avg <- equiv_samples$avg > avg_qual - r2 * sd_qual
                  accept_t <- equiv_samples$avg >= t_min
                  accept_f <- equiv_samples$sd_eq <= f_max
                  accept_vangel_min <- equiv_samples$min > avg_qual - k1 * sd_qual
                  accept_vangel_avg <- equiv_samples$avg > avg_qual - k2 * sd_qual
                  
                  data.frame(
                    accept = c(accept_min & accept_avg,
                               accept_t & accept_f,
                               accept_vangel_min & accept_vangel_avg),
                    method = as.factor(c(rep_len("Two-Sample", n_equiv_2),
                                         rep_len("MSD", n_equiv_2),
                                         rep_len("Vangel", n_equiv_2)))
                  )
                }
              ) %>% 
                group_by(method) %>% 
                summarise(`Rejection Rate` = sum(!accept) / n(),
                          .groups = "drop") %>%
                mutate(n = cur_n, m = cur_m, delta = delta)
            }
          )
        }
      )
    }
  )
}

Power Simulation (Reduction in Mean)

sim_power_mean <- power_simulation(
  delta_mean,
  function(cur_n) {
    rnorm(cur_n, mu_pop, sd_pop)
  },
  function(i_equiv, delta, cur_m) {
    rnorm(cur_m, mu_pop - delta * sd_pop, sd_pop)
  }
)

Power Simulation (Increase in SD)

sim_power_sd <- power_simulation(
  delta_sd,
  function(cur_n) {
    rnorm(cur_n, mu_pop, sd_pop)
  },
  function(i_equiv, delta, cur_m) {
    rnorm(cur_m, mu_pop, sd_pop * delta)
  }
)

Power Simulation (Reduction in Minimum Individual)

sim_power_min_indiv <- power_simulation(
  delta_min_indiv,
  function(cur_n) {
    rnorm(cur_n, mu_pop, sd_pop)
  },
  function(i_equiv, delta, cur_m) {
    c(
      rnorm(cur_m - 1, mu_pop, sd_pop),
      rnorm(1, mu_pop - delta * sd_pop, sd_pop)
    )
  }
)

Power Simulation (Mixture Distribution)

sim_power_mixture <- power_simulation(
  delta_mixture,
  function(cur_n) {
    rnorm(cur_n, mu_pop, sd_pop)
  },
  function(i_equiv, delta, cur_m) {
    m.1 <- rbinom(1, cur_m, pi.mixture)
    c(
      rnorm(cur_m - m.1, mu_pop, sd_pop),
      rnorm(m.1, mu_pop - delta * sd_pop, sd_pop)
    )
  }
)

Power Simulation (Lognormal Data)

For a lognormal distribution, the coefficient of variation (CV) is given by:

\[ CV = \sqrt{\exp\left(\sigma^2\right) - 1} \]

This can be rearranged to solve for \(\sigma\) as follows. We’ll be setting a CV that is typical of composite materials:

\[ \sigma = \sqrt{\log\left(CV^2 + 1\right)} \]

The mean for a lognormal distribution is given by:

\[ mean = m \exp\left(\frac{1}{2}\sigma^2\right) \]

We’ll set the mean equal to 100.

The standard deviation is given by:

\[ SD = m \exp\left(\frac{1}{2}\sigma^2\right)\sqrt{\exp\left(\sigma^2\right)-1} \]

We’ll set the acceptance mean so that it’s \(\delta\) times the SD lower than the mean used for the qualification data.

\[ mean_{accept} = mean_{pop} - \delta SD \\ m_{accept} \exp\left(\frac{1}{2}\sigma^2\right) = m_{pop} \exp\left(\frac{1}{2}\sigma^2\right) - \delta m_{pop} \exp\left(\frac{1}{2}\sigma^2\right)\sqrt{\exp\left(\sigma^2\right)-1} \\ m_{accept} = m_{pop} - \delta m_{pop} \sqrt{\exp\left(\sigma^2\right)-1} \]

sigma_pop <- sqrt(log(0.06 ^ 2 + 1))
m_pop <- 100 / exp(1/2 * sigma_pop ^ 2)
delta_lognorm <- seq(0, 2, length.out = 9)

The PDFs of the acceptance samples in this simulation for \(\delta=0\), \(\delta=1\) and \(\delta=2\)

ggplot(data.frame(x = c(80 , 120)), aes(x = x)) +
  stat_function(fun = dlnorm, args = list(
    meanlog = log(m_pop),
    sdlog = sigma_pop)) +
  stat_function(fun = dlnorm, args = list(
    meanlog = log(m_pop - 1.0 * m_pop * sqrt(exp(sigma_pop^2) - 1)),
    sdlog = sigma_pop)) +
  stat_function(fun = dlnorm, args = list(
    meanlog = log(m_pop - 2.0 * m_pop * sqrt(exp(sigma_pop^2) - 1)),
    sdlog = sigma_pop))

sim_power_lognormal <- power_simulation(
  delta_lognorm,
  function(cur_n) {
    rlnorm(cur_n, log(m_pop), sigma_pop)
  },
  function(i_equiv, delta, cur_m) {
    rlnorm(cur_n, log(m_pop - delta * m_pop * sqrt(exp(sigma_pop^2) - 1)), sigma_pop)
  }
)

Power Simulation (Weibull Data)

For a Weibull distribution with scale parameter \(\eta\) and shape parameter \(\beta\), the CV is given by:

\[ CV = \left\{ \frac{ \Gamma\left[\left(\beta + 2\right) / \beta\right] }{ \left\{\Gamma\left[\left(\beta + 1\right) / \beta\right]\right\}^2 } -1 \right\}^\frac{1}{2} \]

shape_pop <- uniroot(
  function(b) {
    cv <- 0.06
    sqrt(
      gamma((b + 2) / b) /
        (gamma((b + 1) / b)) ^ 2
      - 1
    ) - cv
  },
  c(0.1, 1000)
)$root

The mean is given by:

\[ mean = \eta \Gamma\left[\left(\beta + 1\right) / \beta\right] \]

scale_pop <- 100 / gamma((shape_pop + 1) / shape_pop)

The standard deviation for a Weibull distribution is given by:

\[ SD = \eta \sqrt{ \Gamma\left[\left(\beta + 2\right) / \beta\right] - \left\{\Gamma\left[\left(\beta + 1\right) / \beta\right]\right\}^2 } \]

The mean of the acceptance sample is set as follows:

\[ mean_{accept} = mean_{pop} - \delta SD \\ \eta_{accept} \Gamma\left[\left(\beta + 1\right) / \beta\right] = \eta_{pop} \Gamma\left[\left(\beta + 1\right) / \beta\right] - \delta \eta_{pop} \sqrt{ \Gamma\left[\left(\beta + 2\right) / \beta\right] - \left\{\Gamma\left[\left(\beta + 1\right) / \beta\right]\right\}^2 } \\ \eta_{accept} = \eta_{pop} \left(1 - \delta \frac{ \sqrt{ \Gamma\left[\left(\beta + 2\right) / \beta\right] - \left\{\Gamma\left[\left(\beta + 1\right) / \beta\right]\right\}^2 } } { \Gamma\left[\left(\beta + 1\right) / \beta\right] } \right) \]

delta_weibull <- seq(0, 2, length.out = 9)

The PDFs of the acceptance samples in this simulation for \(\delta=0\), \(\delta=1\) and \(\delta=2\)

ggplot(data.frame(x = c(80 , 120)), aes(x = x)) +
  stat_function(fun = dweibull, args = list(
    shape = shape_pop,
    scale = scale_pop)) +
  stat_function(fun = dweibull, args = list(
    shape = shape_pop,
    scale = scale_pop * (1 - 1.0 *(
      sqrt(
        gamma((shape_pop + 2) / shape_pop) -
          (gamma((shape_pop + 1) / shape_pop)) ^ 2
      ) / gamma((shape_pop + 1) / shape_pop)
    )))) +
  stat_function(fun = dweibull, args = list(
    shape = shape_pop,
    scale = scale_pop * (1 - 2.0 *(
      sqrt(
        gamma((shape_pop + 2) / shape_pop) -
          (gamma((shape_pop + 1) / shape_pop)) ^ 2
      ) / gamma((shape_pop + 1) / shape_pop)
    ))))

sim_power_weibull <- power_simulation(
  delta_weibull,
  function(cur_n) {
    rweibull(cur_n, shape = shape_pop, scale = scale_pop)
  },
  function(i_equiv, delta, cur_m) {
    rweibull(
      cur_n,
      shape = shape_pop,
      scale = scale_pop * (1 - delta *(
        sqrt(
          gamma((shape_pop + 2) / shape_pop) -
            (gamma((shape_pop + 1) / shape_pop)) ^ 2
        ) / gamma((shape_pop + 1) / shape_pop)
      )))
  }
)
save(
  mu_pop,
  sd_pop,
  n_qual_1,
  n_equiv_1,
  n_qual_2,
  n_equiv_2,
  sim_alpha,
  delta_mean,
  delta_sd,
  delta_min_indiv,
  sim_equiv,
  sim_power_mean,
  sim_power_sd,
  sim_power_min_indiv,
  sim_power_mixture,
  sim_power_lognormal,
  sim_power_weibull,
  file = "sim_results.RData"
)
sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## Random number generation:
##  RNG:     L'Ecuyer-CMRG 
##  Normal:  Inversion 
##  Sample:  Rejection 
##  
## locale:
##  [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
##  [5] LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
##  [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tikzDevice_0.12.3.1 cmstatr_0.9.1       furrr_0.3.0        
##  [4] future_1.27.0       forcats_0.5.1       stringr_1.4.0      
##  [7] dplyr_1.0.9         purrr_0.3.4         readr_2.1.2        
## [10] tidyr_1.2.0         tibble_3.1.8        ggplot2_3.3.6      
## [13] tidyverse_1.3.2     rmarkdown_2.14     
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.3          sass_0.4.2          bit64_4.0.5        
##  [4] vroom_1.5.7         jsonlite_1.8.0      modelr_0.1.8       
##  [7] bslib_0.4.0         assertthat_0.2.1    highr_0.9          
## [10] googlesheets4_1.0.0 cellranger_1.1.0    yaml_2.3.5         
## [13] globals_0.15.1      pillar_1.8.0        backports_1.4.1    
## [16] glue_1.6.2          digest_0.6.29       rvest_1.0.2        
## [19] colorspace_2.0-3    htmltools_0.5.3     pkgconfig_2.0.3    
## [22] broom_1.0.0         listenv_0.8.0       haven_2.5.0        
## [25] scales_1.2.0        tzdb_0.3.0          googledrive_2.0.0  
## [28] generics_0.1.3      farver_2.1.1        ellipsis_0.3.2     
## [31] cachem_1.0.6        withr_2.5.0         cli_3.3.0          
## [34] magrittr_2.0.3      crayon_1.5.1        readxl_1.4.0       
## [37] evaluate_0.15       fs_1.5.2            fansi_1.0.3        
## [40] parallelly_1.32.1   MASS_7.3-58         xml2_1.3.3         
## [43] SuppDists_1.1-9.7   tools_4.2.1         hms_1.1.1          
## [46] gargle_1.2.0        lifecycle_1.0.1     munsell_0.5.0      
## [49] reprex_2.0.1        kSamples_1.2-9      compiler_4.2.1     
## [52] jquerylib_0.1.4     rlang_1.0.4         grid_4.2.1         
## [55] filehash_2.4-3      labeling_0.4.2      gtable_0.3.0       
## [58] codetools_0.2-18    DBI_1.1.3           R6_2.5.1           
## [61] lubridate_1.8.0     knitr_1.39          fastmap_1.1.0      
## [64] bit_4.0.4           utf8_1.2.2          stringi_1.7.8      
## [67] parallel_4.2.1      vctrs_0.4.1         dbplyr_2.2.1       
## [70] tidyselect_1.1.2    xfun_0.31