Bioequivalence and Bioavailability Forum

Hi Rocco,

❝ Where does the formula on slide 10.83 in bebac.at/lectures/Leuven2013WS2.pdf for CI for parallel design come from? I cannot seem to find reference anywhere.

Honestly, I don’t remember why I simplified the commonly used formula.
Algebra:$$s\sqrt{\tfrac{n_1+n_2}{n_1n_2}}=\sqrt{s^2(1/n_1+1/n_2)}\;\tiny{\square}$$ Comparison with the data of the example.

• Like in the presentation:
  mean.log <- function(x) mean(log(x), na.rm = TRUE)
T      <- c(100, 103, 80, 110,  78, 87, 116, 99, 122, 82, 68,  NA)
R      <- c(110, 113, 96,  90, 111, 68, 111, 93,  93, 82, 96, 137)
n1     <- sum(!is.na(T))
n2     <- sum(!is.na(R))
s1.2   <- var(log(T), na.rm = TRUE)
s2.2   <- var(log(R), na.rm = TRUE)
s0.2   <- ((n1 - 1) * s1.2 + (n2 - 1) * s2.2) / (n1 + n2 - 2)
s0     <- sqrt(s0.2)
nu.1   <- n1 + n2 - 2
nu.2   <- (s1.2 / n1 + s2.2 / n2)^2 /
          (s1.2^2 / (n1^2 * (n1 - 1)) + s2.2^2 / (n2^2 * (n2 - 1)))
t.1    <- qt(p = 1 - 0.05, df = nu.1)
t.2    <- qt(p = 1 - 0.05, df = nu.2)
PE.log <- mean.log(T) - mean.log(R)
CI.1   <- PE.log + c(-1, +1) * t.1 *
          s0 * sqrt((n1 + n2) / (n1 * n2))
CI.2   <- PE.log + c(-1, +1) * t.2 *
          sqrt(s1.2 / n1 + s2.2 / n2)
CI.t   <- 100 * exp(CI.1)
CI.w   <- 100 * exp(CI.2)
fmt    <- "%.3f %.3f %.3f %.2f    %.2f   %.2f"
cat("    method     df mean.T mean.R    PE CL.lower CL.upper",
    "\n    t-test",
    sprintf(fmt, nu.1, exp(mean.log(T)), exp(mean.log(R)),
            100 * exp(PE.log), CI.t[1], CI.t[2]),
    "\nWelch-test",
    sprintf(fmt, nu.2, exp(mean.log(T)), exp(mean.log(R)),
            100 * exp(PE.log), CI.w[1], CI.w[2]), "\n")

    method     df mean.T mean.R    PE CL.lower CL.upper
    t-test 21.000 93.554 98.551 94.93    83.28   108.20
Welch-test 20.705 93.554 98.551 94.93    83.26   108.23
  
  
• More comfortable with the t.test() function, where var.equal = TRUE gives the t-test and var.equal = FALSE the Welch-test:
  res <- data.frame(method = c("t-test", "Welch-test"), df = NA,
                             mean.T = NA, mean.R = NA, PE = NA,
                             CL.lower = NA, CL.upper = NA)
var.equal <- c(TRUE, FALSE)
for (j in 1:2) {
  x <- t.test(x = log(T), y = log(R), conf.level = 0.90,
              var.equal = var.equal[j])
  res[j, 2]   <- signif(x[[2]], 5)
  res[j, 3:4] <- signif(exp(x[[5]][1:2]), 5)
  res[j, 5]   <- round(100 * exp(diff(x[[5]][2:1])), 2)
  res[j, 6:7] <- round(100 * exp(x[[4]]), 2)
}
print(res, row.names = FALSE)

    method     df mean.T mean.R    PE CL.lower CL.upper
    t-test 21.000 93.554 98.551 94.93    83.28   108.20
Welch-test 20.705 93.554 98.551 94.93    83.26   108.23 □

The formula for Satterthwaite’s approximation of the degrees of freedom given in slide 11 contained typos (corrected in the meantime). Of course,$$\nu\approx\frac{\left(\frac{{s_{1}}^{2}}{n_1}+\frac{{s_{2}}^{2}}{n_2}\right)^2}{\frac{{s_{1}}^{4}}{n{_{1}}^{2}(n_1-1)}+\frac{{s_{2}}^{4}}{n{_{2}}^{2}(n_2-1)}}$$Satterthwaite’s approximation adjusts both for unequal variances and group sizes. The conventional t-test is fairly robust against the former but less so for the latter.
In the -function t.test() var.equal = FALSE is the default because:

• Any pre-test is bad practice (might inflate the type I error) and has low power esp. for small sample sizes.
  
• If \({s_{1}}^{2}={s_{2}}^{2}\;\wedge\;n_1=n_2\), the formula given above reduces to the simple \(\nu=n_1+n_2-2\) anyhow.
  
• In all other cases the Welch-test is conservative, which is a desirable property.