Two-sample normal sample size

Introduction

Limited support is provided for 2-sample design with a normally distributed random variable as the outcome. Users are encouraged to look at guidance such as in Jennison and Turnbull (2000). We provide a tool where for a large sample case where a reasonable estimate of standard deviation is available, a reasonable sample size can be computed based straightforward distribution theory outlined below.

The problem considered

The overall sample size notation used for gsDesign is to consider a standardized effect size parameter which is referred to as θ in Jennison and Turnbull (2000). We begin with the 2-sample normal problem where we assume a possibly different standard deviation in each treatment group. For j = 1, 2, we let X_j, i, i = 1, 2, …n_j represent independent and identically distributed observations following a normal distribution with mean μ_j and standard deviation σ_j. The natural parameter for comparing the two distributions is

δ = μ₂ − μ₁

and we wish to test if δ > 0 in a one-sided testing scenario to test for superiority of treatment 2 over treatment 1. We could also consider testing, say, δ > δ₀ for a non-inferiority scenario with δ₀ < 0 or super superiority if δ₀ > 0. While normally a t-test would be used for this, for large sample sizes this would be nearly equivalent to a Z-test defined by:

$$Z=\frac{\bar X_2 - \bar X_1-\delta_0}{\sqrt{\sigma^2_2/n_2 + \sigma_1^2/n_1}}\approx \frac{\bar X_2 - \bar X_1}{\sqrt{s^2_2/n_2 + s_1^2/n_1}}=t$$ where X̄_j is the sample mean and s_j² is the sample variance for group j = 1, 2. The far right hand side of this is Welch’s t-test. For our examples we use this t-test and show that the sample size computation based on the Z-test above works well for the chosen problems.

Sample size

Thus, n₂ = rn/(1 + r), n₁ = n/(1 + r) and when r = 1 we have n₁ = n₂ = n/2. Now that we have completed needed notation, those not interested in the theory behind the sample size and power calculation used may skip the rest of this section.

We let σ² = (1 + r)(σ₁² + σ₂²/r) and define θ = (δ − δ₀)/σ. Under the given assumptions, $$Z \sim \text{Normal}\left(\sqrt n\theta,1\right).$$ Under the null hypothesis that δ = δ₀, we have Z ∼ Normal(0, 1). Thus, regardless of n we have P₀[Z ≥ Φ⁻¹(1 − α)] = α. Under the alternate hypothesis that δ = δ₁ and we denote a corresponding θ₁. We define the type II error β and power 1 − β by

$$ \begin{align} 1-\beta =& P_1[Z\ge \Phi^{-1}(1-\alpha)]\\ =& P[Z-\sqrt n\theta_1\ge \Phi^{-1}(1-\alpha)-\sqrt n\theta_1]\\ =&\Phi(\Phi^{-1}(1-\alpha)-\sqrt n\theta_1)). \end{align}$$

If the power 1 − β is fixed, we can invert this formula to compute sample size with:

$$n= \left(\frac{\Phi^{-1}(1-\beta)+\Phi^{-1}(1-\alpha)}{\theta_1}\right)^2.$$

For 2-sided testing, we simply substitute α/2 for α in the above two formulas.

Examples

We consider two examples to check the above formulas vs. nNormal(). We then confirm that the approximation is working well by simulating and confirming that the power and Type I error approximations are useful. Finally, we provide a simple group sequential design example.

r <- 2
sigma <- sqrt((1 + r) * (1.6^2 + 1.25^2 / r))
theta <- 0.8 / sigma
((qnorm(.9) + qnorm(.975)) / theta)^2
#> [1] 164.5684

nNormal(delta1 = 0.8, sd = 1.6, sd2 = 1.25, alpha = 0.025, n = 200, ratio = 2)
#> [1] 0.9466825

pwr <- pnorm(qnorm(.975) - sqrt(200) * theta, lower.tail = FALSE)
pwr
#> [1] 0.9466825

n <- 100:200
pwrn <- nNormal(delta1 = 0.8, sd = 1.6, sd2 = 1.25, alpha = 0.025, n = n, ratio = 2)
plot(n, pwrn, type = "l")

delta1 <- seq(.5, 1, .025)
pwrdelta1 <- nNormal(delta1 = delta1, sd = 1.6, sd2 = 1.25, alpha = 0.025, n = 200, ratio = 2)
plot(delta1, pwrdelta1, type = "l")

nsim <- 1000000
delta <- 0.8
sd1 <- 1.6
sd2 <- 1.25
n1 <- 67
n2 <- 133
deltahat <- rnorm(n = nsim, mean = delta, sd = sd1 / sqrt(n1)) -
  rnorm(n = nsim, mean = 0, sd = sd2 / sqrt(n2))
s <- sqrt(
  sd1^2 * rchisq(n = nsim, df = n1 - 1) / (n1 - 1) / n1 +
    sd2^2 * rchisq(n = nsim, df = n2 - 1) / (n2 - 1) / n2
)
z <- deltahat / s
mean(z >= qnorm(.975))
#> [1] 0.946545

sqrt(pwr * (1 - pwr) / nsim)
#> [1] 0.0002246659

d <- gsDesign(
  k = 2,
  n.fix = nNormal(delta1 = 0.8, sd = 1.6, sd2 = 1.25, alpha = 0.025, beta = .1, ratio = 2),
  delta1 = 0.8
)
d %>%
  gsBoundSummary(deltaname = "Mean difference") %>%
  kable(row.names = FALSE)

cat(summary(d))

gsDesign(
  k = 2,
  delta = theta,
  delta1 = 0.8
) %>%
  gsBoundSummary(deltaname = "Mean difference") %>%
  kable(row.names = FALSE)