Scenario

Suppose we have binomial data in a balanced one-way layout with \(a\) samples (e.g., treatment groups) and \(n\) trials per sample so that the total sample size equals \(N = an\). We are to determine the power to reject the global null hypothesis that the success probabilities \(\pi_i\) of all treatments \(i = 1, \dots, a\) are equal to some common success probability \(\pi_0\).

Numerical example

Assume \(a=6\) treatments are tested with \(n=500\) subjects and the outcome is binary. The success probability for treatments 1 through 5 is \(\pi = 0.5\) whereas it is \(\pi + b = 0.5 + 0.05 = 0.55\) for treatment \(6\), and we want to assess \(H_0: \pi_1 = \dots = \pi_a = \pi_0\).

a <- 6
n <- 500
pi <- 0.5 
b <- 0.05

The alternative we are about to investigate is:

(p <- c(rep(pi, a - 1), pi + b))
## [1] 0.50 0.50 0.50 0.50 0.50 0.55

The hypothesized common success probability \(\pi_0\) is:

(p0 <- rep(mean(c(rep(pi, a - 1), pi + b)), a))
## [1] 0.5083 0.5083 0.5083 0.5083 0.5083 0.5083

Product binomial sampling

Under product binomial sampling the noncentrality parameter for \(\chi^2\) is computed as \(\lambda = \delta^T diag(\pi_0)^{-1} \delta\) where \(\delta = \sqrt{n}(\pi - \pi_0)\):

delta <- sqrt(n) * (p - p0)
(t(matrix(delta)) %*% solve(diag(p0)) %*% matrix(delta))
##       [,1]
## [1,] 2.049

We may formulate \(\lambda\) also similar to equation (6.11) of Agresti (2013):

(lambdaChi2 <- n * sum((p - p0)^2 / p0))
## [1] 2.049

We approximate the power by \(P[\chi^2_1(\lambda_{\chi^2}) > c]\) with \(c\) the 95% quantile of \(\chi^2_1\):

crit <- qchisq(0.95, 1)
1 - pchisq(crit, 1, ncp=lambdaChi2)
## [1] 0.2989

The noncentrality parameter for \(G^2\) resembles equation (6.12) of Agresti (2013):

(lambdaG2 <- 2 * n * sum(p * log(p / p0)))
## [1] 2.006

We approximate the power by \(P[\chi^2_1(\lambda_{G^2}) > c]\) with \(c\) as above:

1 - pchisq(crit, 1, ncp=lambdaG2)
## [1] 0.2937

Check via simulation

We simulate the powers of \(\chi^2\) and \(G^2\) for the above setting under product binomial sampling with 10,000 replications:

g.test <- function(x) {
    expected <- rowSums(x) %o% colSums(x)/sum(x)
    g <- 2 * sum(x * log(x/expected))
    df <- (nrow(x) - 1) * (ncol(x) - 1)
    p <- 1 - pchisq(g, df)
    return(p)
}
runs <- 10000
Chi2 <- G2 <- 0
for(i in 1:runs){
  MyData <- data.frame(Drug = LETTERS[1:a],
                       Success = c(rbinom(1, n, pi), rbinom(1, n, pi), rbinom(1, n, pi),
                                   rbinom(1, n, pi), rbinom(1, n, pi), rbinom(1, n, pi + b)),
                       Total = n)
  succ <- xtabs(Success ~ Drug, MyData)
  fail <- xtabs(Total ~ Drug, MyData) - succ
  MyTable <- rbind(succ, fail)
  if(chisq.test(MyTable)$p.value < 0.05){
    Chi2 <- Chi2 + 1
  }
  if(g.test(MyTable) < 0.05){
    G2 <- G2 + 1
  }
}
Chi2 / runs
## [1] 0.3013
G2 / runs
## [1] 0.303

This is very close to the power calculated unter product binomial sampling.

Multinomial sampling

Under multinomial sampling (which we do not have here) the noncentrality parameter \(\lambda_M\) for \(\chi^2\) is computed with \(\delta_M = \sqrt{N}(\pi - \pi_0)\):

deltaM <- sqrt(n * a) * (p - p0)
t(matrix(deltaM)) %*% solve(diag(p0)) %*% matrix(deltaM)
##      [,1]
## [1,] 12.3

Or we can write it like in equation (6.11) of Agresti (2013):

(lambdaM <- n * a * sum((p - p0)^2 / p0))
## [1] 12.3

Now we approximate the power by \(P[\chi^2_{a - 1}(\lambda_M) > c_M]\) with \(c_M\) the 95% quantile of \(\chi^2_{a - 1}\):

critM <- qchisq(0.95, a - 1)
1 - pchisq(critM, a - 1, ncp=lambdaM)
## [1] 0.7802

This is radically different from the power calculated under product binomial sampling!

Similarly, the noncentrality parameter for \(G^2\) under multinomial sampling is that of equation (6.12) of Agresti (2013):

(lambdaM <- 2 * n * a * sum(p * log(p / p0)))
## [1] 12.04

And the power is approximated as:

1 - pchisq(critM, a - 1, ncp=lambdaM)
## [1] 0.7701