1 The pwr_freq() Function

The pwr_freq() function takes in a set of fixed inputs pertaining to a one-sample hypothesis test and returns the corresponding statistical power of the test.

To elaborate, consider the following one-sided hypothesis test for population mean \(\theta\): \[ H_0: \theta = \theta_0 \\ H_a: \theta = \theta_1 > \theta_0, \] where \(\theta_0\) is the reference value that we will test against. Assuming the known population variance is denoted as \(\sigma^2\) and the sample mean’s distribution is approximately Gaussian, \(H_0\) is rejected if \(\bar{y} > \theta_0 + \frac{\sigma}{\sqrt{n}}Z_{1-\alpha}\), where \(\bar{y}\) is the sample mean, \(\alpha\) is the specified Type I error, \(Z_{1-\alpha}\) is the corresponding quantile of the Gaussian distribution, and \(n\) is the sample size. We can use this rejection criteria to derive the statistical power of the test, defined by \(1-\beta = P(\text{reject } H_0 | H_a \text{ is true})\). Straightforward standardization procedure leads to \[\begin{equation} \label{eq:power_func} 1-\beta = P\left(\bar{y} > \theta_0 + \frac{\sigma}{\sqrt{n}}Z_{1-\alpha} \Bigg| \theta = \theta_1 \right) = \Phi\left(\sqrt{n}\frac{\Delta}{\sigma} - Z_{1-\alpha}\right), \end{equation}\] where \(\Delta = \theta_1 - \theta_0\) denotes the critical difference and \(\Phi\) denotes the cumulative distribution function of the standard normal. Similar steps can be taken to determine the power expressions for alternative comparisons of \(H_a\). When \(H_a: \theta_1 < \theta_0\), the power is determined as \[ 1-\beta = P\left(\bar{y} < \theta_0 - \frac{\sigma}{\sqrt{n}}Z_{1-\alpha} \Bigg| \theta = \theta_1 \right) = 1 - \Phi\Big(\sqrt{n}\frac{\Delta}{\sigma} + Z_{1-\alpha}\Big). \] Finally in the two-sided case where \(H_a: \theta_1 \neq \theta_0\), the power is determined as \[\begin{align*} 1-\beta &= P\left(\bar{y} < \theta_0 - \frac{\sigma}{\sqrt{n}}Z_{1-\alpha/2} \Bigg | \theta = \theta_1 \right) + P\left(\bar{y} > \theta_0 + \frac{\sigma}{\sqrt{n}}Z_{1-\alpha/2} \Bigg | \theta = \theta_1 \right)\\ &= 1 + \Phi\left(\sqrt{n}\frac{\Delta}{\sigma} - Z_{1-\alpha/2}\right) - \Phi\left(\sqrt{n}\frac{\Delta}{\sigma} + Z_{1-\alpha/2}\right). \end{align*}\]

library(bayesassurance)

n <- seq(10, 140, 5)
pwr_vals <- bayesassurance::pwr_freq(n = n, theta_0 = 0.15, theta_1 = 0.35, sigsq = 0.3, 
                                     alt = "greater", alpha = 0.05)

head(pwr_vals$pwr_table)

pwr_vals$pwr_plot

n <- 20
pwr_freq(n = n, theta_0 = 0.15, theta_1 = 0.35, sigsq = 0.3, 
                                     alt = "greater", alpha = 0.05)
#> [1] "Power: 0.495"

2 The assurance_nd_na() Function

The assurance_nd_na() function determines the assurance of a specified outcome based on a closed-form solution that is derived under the Bayesian setting.

Suppose we seek to evaluate the tenability of \(\theta > \theta_0\) given data from a Gaussian population with mean \(\theta\) and known variance \(\sigma^2\). We assign two sets of priors for \(\theta\), one at the \(\textit{design stage}\) and the other at the \(\textit{analysis stage}\). The analysis objective specifies the condition that needs to be satisfied, which in this case, involves observing that \(P(\theta > \theta_0| \bar{y}) > 1-\alpha\). The design objective seeks a sample size that is needed to ensure that the analysis objective is met \(100\delta\%\) of the time, where \(\delta\) denotes the assurance.

Let \(\theta \sim N\big(\theta_1, \frac{\sigma^2}{n_a}\big)\) be our analysis stage prior and \(\theta \sim N\big(\theta_1, \frac{\sigma^2}{n_d}\big)\) be our design stage prior, where \(n_a\) and \(n_d\) are precision parameters that respectively quantify the degree of belief carried towards parameter \(\theta\) and the degree of belief carried towards the population from which we are drawing samples from to evaluate \(\theta\). Then, given the likelihood \(\bar{y} \sim N\big(\theta, \frac{\sigma^2}{n}\big)\), we can obtain the posterior distribution of \(\theta\) by multiplying the analysis prior and likelihood: \[\begin{equation}\label{eq: simple_posterior} N\left(\theta {\left | \theta_1, \frac{\sigma^2}{n_a} \right.}\right) \times N\left(\bar{y} {\left | \theta, \frac{\sigma^2}{n} \right.}\right) \propto N\left(\theta {\left | \frac{n_a}{n + n_a}\theta_1 + \frac{n}{n + n_a}\bar{y}, \frac{\sigma^2}{n + n_a}\right.}\right)\;. \end{equation}\]

This posterior distribution gives us \(P(\theta > \theta_0 | \bar{y})\) and the assurance is then defined as \[\begin{equation}\label{eq:assurance} \delta = P_{\bar{y}}\left\{\bar{y}: P(\theta > \theta_0 | \bar{y}) > 1 - \alpha\right\}. \end{equation}\] The assurance expression can be expanded out further by using the marginal distribution of \(\bar{y}\), which is obtained by \[ \int{N\left(\theta {\left|\theta_1, \frac{\sigma^2}{n_d}\right.}\right) \times N\left(\bar{y} {\left | \theta, \frac{\sigma^2}{n} \right.}\right) d\theta} = N\left(\bar{y} {\left |\theta_1, \Big(\frac{1}{n} + \frac{1}{n_d}\Big) \sigma^2 \right.}\right) . \] Since the assurance definition is conditioned on \(\bar{y}\), we use this to standardize the assurance expression to obtain the following closed-form solution: \[ \delta(\Delta, n) = \Phi\Bigg(\sqrt{\frac{nn_d}{n+n_d}} \Bigg[\frac{n+n_a}{n}\frac{\Delta}{\sigma} + Z_{\alpha}\frac{\sqrt{n+n_a}}{n}\Bigg]\Bigg), \] where \(\Delta = \theta_1 - \theta_0\). Similar steps can be taken to derive the closed-form expressions of assurance for \(\theta < \theta_0\) and \(\theta \neq \theta_0\). If we wish to seek the tenability of \(\theta < \theta_0\), the assurance is given by \[ \delta(\Delta, n) = 1-\Phi\Bigg(\sqrt{\frac{nn_d}{n+n_d}} \Bigg[\frac{n+n_a}{n}\frac{\Delta}{\sigma} + Z_{1-\alpha}\frac{\sqrt{n+n_a}}{n}\Bigg]\Bigg). \] Finally, for evaluating the tenability of \(\theta \neq \theta_0\), the assurance is given by \[ \delta(\Delta, n) = 1-\Phi\Bigg(\sqrt{\frac{nn_d}{n+n_d}} \Bigg[\frac{n+n_a}{n}\frac{\Delta}{\sigma} + Z_{1-\alpha/2}\frac{\sqrt{n+n_a}}{n}\Bigg]\Bigg) + \Bigg[\frac{n+n_a}{n}\frac{\Delta}{\sigma} + Z_{\alpha/2}\frac{\sqrt{n+n_a}}{n}\Bigg]\Bigg). \]

n <- seq(10, 500, 5)
n_a <- 1e-8
n_d <- 1e+8
theta_0 <- 0.15
theta_1 <- 0.25
sigsq <- 0.104
assur_vals <- assurance_nd_na(n = n, n_a = n_a, n_d = n_d, theta_0 = theta_0, 
                              theta_1 = theta_1, sigsq = sigsq, alt = "greater", 
                              alpha = 0.05)

head(assur_vals$assurance_table)

assur_vals$assur_plot

3 The pwr_curves() Function

The pwr_curves() function constructs a single plot based on the combined set of inputs used to compute the frequentist power and Bayesian assurance. The plot includes three curves that include the power curve, the assurance curve obtained under a closed-form solution, and optionally, the assurance curve obtained under a simulation setting. The optional simulated assurance points are obtained using bayes_sim, which is discussed in a later vignette.

Load in the bayesassurance package and specify the following parameters:

1. n: sample size (either scalar or vector)
2. n_a: Precision parameter that quantifies the degree of belief carried towards parameter \(\theta\)
3. n_d: Precision parameter that quantifies the degree of belief carried towards the population from which we are drawing samples from to evaluate \(\theta\).
4. theta_0: parameter value that is known a priori (typically provided by the client)
5. theta_1: alternative parameter value that will be tested in comparison to #’ theta_0. See alt for specification options.
6. sigsq: known variance \(\sigma^2\)
7. alt: specifies alternative test case, where alt = "greater" tests if \(\theta > \theta_0\), alt = "less" tests if \(\theta < \theta_0\), and alt = "two.sided" performs a two-sided test. By default, alt = "greater".
8. alpha: significance level
9. bayes_sim: when set to TRUE, this indicates that the user wants to include simulated assurance results in the outputted plot. Default setting is FALSE.
10. mc_iter: If bayes_sim = TRUE, specifies the number of MC samples to evaluate. Default set to 5000 when not specified.

n <- seq(100, 300, 10)
n_a <- 1e-8
n_d <- 1e+8
theta_0 <- 0.15
theta_1 <- 0.25
sigsq <- 0.104
alpha <- 0.05

out1 <- bayesassurance::pwr_curve(n, n_a, n_d, theta_0, theta_1, sigsq,  
                                  alt = "greater", alpha, bayes_sim = FALSE)

head(out1$power_table)

head(out1$assurance_table)

n <- seq(100, 300, 10)
n_a <- 1e-8
n_d <- 1e+8
theta_0 <- 0.15
theta_1 <- 0.25
sigsq <- 0.104
alpha <- 0.05

set.seed(10)
out2 <- bayesassurance::pwr_curve(n, n_a, n_d, theta_0, theta_1, sigsq, alt = "greater", 
                                  alpha, bayes_sim = TRUE)

head(out2$bayes_sim_table)

out2$plot