Sensitivity Analysis

Sensitivity analysis involves fixing the sample size, desired power, and alpha level, and determining the effect size that can be detected with the desired power. We demonstrate this functionality using a t-test.

The fixed parameters are the sample size N, the alpha level, and the goal power. We want to determine the corresponding effect size:

N <- 100
alpha <- 0.01
goal_power <- 0.95

We define our simulation function:

simfun_sensitivity <- function(esize) {

    # Generate a data set
    dat <- rnorm(n = N, mean = esize)
    # Test the hypothesis
    res <- t.test(dat)
    res$p.value < alpha
}

We can find the corresponding effect size using:

# Effect Size Finding res1 <- find.design(simfun
# = simfun_sensitivity, boundaries =
# c(0,1),integer=FALSE, power =
# goal_power,surrogate = 'gpr',evaluations=8000)
res1 <- extensions_results[[1]]
summary(res1)
#> 
#> Call:
#> find.design(simfun = simfun_sensitivity, boundaries = c(0, 1), 
#>     power = goal_power, evaluations = 8000, surrogate = "gpr", 
#>     integer = FALSE)
#> 
#> Design: esize = 0.431758284218821
#> 
#> Power: 0.95127,  SE: 0.0028
#> Evaluations: 8000,  Time: 26.02,  Updates: 16
#> Surrogate: Gaussian process regression

This leads us to an effect size of about 0.43, as measured by Cohen’s d. To confirm the correctness of this result, we check it analytically using the pwr package.

library(pwr)
esize_correct <- pwr.t.test(n = N, sig.level = alpha,
    power = goal_power, type = "one.sample")$d
esize_correct
#> [1] 0.4293178

The analytical result is very close to the result obtained with mlpwr. With a small simulation study, we further confirm that the analytical result is correct.

a <- replicate(10000, simfun_sensitivity(esize_correct))
mean(a)
#> [1] 0.9521

Compromise Analysis

In compromise analysis, the sample size and effect size are fixed in the DGF, and the goal is to find a decision criterion that optimizes the desired ratio of alpha and beta errors. Again, we demonstrate this functionality using a t-test scenario.

The fixed parameters are the sample size N, the effect size (given by Cohen’s d) and the desired ratio of the alpha and beta errors:

N <- 100
esize <- 0.3
desired_ratio <- 1

We now define our simulation function. Its argument is a given threshold for the test statistic. This function then generates two data sets - one under the null, one under the alternative hypothesis - and evaluates whether the use of this threshold leads to a Type I error (for the data set generated under the null hypothesis) or a Type II error (for the data set generated under the alternative hypothesis).

simfun_compromise <- function(crit) {

    # Generate a data set
    dat <- rnorm(n = N, mean = 0)

    # Test the hypothesis
    res <- t.test(dat)
    a <- res$statistic > crit

    # Generate a data set
    dat <- rnorm(n = N, mean = esize)

    # Test the hypothesis
    res <- t.test(dat)
    b <- res$statistic < crit

    return(c(a, b))
}

We run this function one to see how the output looks - TRUE indicates that a Type I (1st entry) or Type II error (2nd entry) was observed:

simfun_compromise(0.1)
#>     t     t 
#> FALSE FALSE

However, we are not interested in the Type I and Type II errors themselves, but in their ratio. We therefore need to define a function that organizes how the individual simulations are aggregated:

aggregate_fun <- function(x) {
    y <- rowMeans(matrix(x, nrow = 2))
    y[2]/y[1]
}

We can find the optimal decision criterion for reaching the desired ratio of alpha and beta errors using:

# res2 <- find.design(simfun = simfun_compromise,
# boundaries = c(1,2), power =
# desired_ratio,integer = FALSE, aggregate_fun =
# aggregate_fun,surrogate='svr',use_noise=FALSE,evaluations=8000)
res2 <- extensions_results[[2]]
summary(res2)
#> 
#> Call:
#> find.design(simfun = simfun_compromise, boundaries = c(1, 2), 
#>     power = desired_ratio, evaluations = 8000, surrogate = "svr", 
#>     aggregate_fun = aggregate_fun, integer = FALSE, use_noise = FALSE)
#> 
#> Design: crit = 1.48455055391577
#> 
#> Power: 1.00086,  SE: NA
#> Evaluations: 8000,  Time: 6.89,  Updates: 6
#> Surrogate: Support vector regression

The analysis with mlpwr leads to a desired threshold of about 1.48. The ratio of Type I and Type II errors is given as “power” in this output and close to the desired value of 1.

Note that we used some additional arguments in this case. The ‘integer’ argument tells the algorithm whether to look for integer design parameters - this is set to ‘FALSE’ for our case of effect sizes. The ‘use_noise=FALSE’ indicates that we do not account for the Monte Carlo error at each design parameter. This functionality can be added in the future. Note also that we use ‘power’ as an argument to the find.design function and in the output of summary(res2) to indicate the desired ratio of alpha and beta error. This is a makeshift solution and we consider rewording this argument to ‘goal’ in the future.

To check the correctness of our result, we compare it with an analytical solution obtained from the t-distribution and an numerical optimizer:

fun <- \(alpha) {
    beta <- alpha * desired_ratio
    abs(qt(1 - alpha, ncp = 0, df = N - 1)  # crit for alpha
 -
        qt(beta, ncp = esize * sqrt(N), df = N - 1)  # crit for beta
)
}
alpha <- optim(0.1, fun, lower = 1e-09, upper = 1 -
    1e-09, method = "L-BFGS-B")$par
crit <- qt(1 - alpha, ncp = 0, df = N - 1)
crit
#> [1] 1.503731

This leads to a threshold of about 1.50, which is very close to the solution of mlpwr. Using simulations, we can show that the analytical solution is indeed correct:

a <- replicate(10000, simfun_compromise(crit))
aggregate_fun(a)
#> [1] 0.992515