power.anova.testLet’s now revisit the binding fraction example for sample size calculation. Suppose that we want to test equal mean binding fractions among antibiotics against the alternative \[ H_1: \mu_p = \mu + 3, \mu_t = \mu+3, \mu_s = \mu - 6, \mu_E = \mu, \mu_C = \mu \] so that \[ \mu_ 1 = \mu_2 = 3, \mu_3 = -6, \mu_4 = \mu_5 = 0. \]
Assume \(\sigma = 3\) and we need to use \(\alpha = 0.05\). The noncentrality parameter is given by \[ \gamma = n[(\frac{3}{3})^2 +(\frac{3}{3})^2 + (\frac{-6}{3})^2] \] and the \(\alpha = 0.05\) critical value for \(H_0\) is given by \[ F^* = F(5 - 1, 5n - 5, 0.05) \] We need the area to the right of \(F^*\) for the noncentral \(F\) distribution with degrees of freedom \(4\) and \(5(n - 1)\) and noncentrality parameter \(\gamma = 6n\) to be greater or equal to the desired power level of \(1 - \beta = 0.8\).
Now we will do the “real” calculations to finish this example.
We proceed by filling out the columns of the following table step by step.
| \(n\) | ndf | ddf | NCP | QF | PF | POWER |
|---|---|---|---|---|---|---|
| 2 | ||||||
| 3 | ||||||
| 4 | ||||||
| 5 | ||||||
| 6 | ||||||
| 7 | ||||||
| 8 | ||||||
| 9 | ||||||
| 10 |
rm(list=ls()) # clean up workspace
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.8
## ✓ tidyr 1.2.0 ✓ stringr 1.4.0
## ✓ readr 2.1.2 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()binding.fraction.data <- matrix(
c(29.6, 24.3, 28.5, 32.0,
27.3, 32.6, 30.8, 34.8,
5.8, 6.2, 11.0, 8.3,
21.6, 17.4, 18.3, 19,
29.2, 32.8, 25.0, 24.2), nrow = 5, ncol = 4, byrow = TRUE
)
rownames(binding.fraction.data) <- c("Penicillin G", "Tetracyclin", "Streptomycin", "Erythromycin", "Chloramphenicol")
binding.df <- as_tibble(binding.fraction.data, rownames = NA) %>%
rownames_to_column(var = "antibiotic") %>%
pivot_longer(cols = 2:5, values_to = "bf", names_to = NULL)## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.Let’s consider a sequence of \(n\) values from \(2\) to \(10\).
n <- 2:10\[ \Pr(MS[Trt]/MS[E] < F^*; H_1 \mbox{ is true}) \]
\[ \Pr(MS[Trt]/MS[E] > F^*; H_1 \mbox{ is true}) \]
power.anova.testRead the documentation of power.anova.test function in R. You may use this webpage for reference webpage link.
Repeat the power calculations above using power.anova.test function. Do you obtain the same values?