94 lines
2.6 KiB
R
94 lines
2.6 KiB
R
#' ---
|
|
#' title: "Two-by-two ANOVA: Power simulation of the interaction test"
|
|
#' author: ""
|
|
#' date: "Last modified: 2026-01-09"
|
|
#' ---
|
|
|
|
#' ## Application context
|
|
|
|
#'
|
|
#' Effect of fertilizers
|
|
#'
|
|
#' In an experiment, two fertilizers (A and B, each either low or high dose)
|
|
#' will be combined and the yield of peas (Y) in kg be observed. Goal is to
|
|
#' detect an increase of the Fertilizer-A effect by an additional 12 kg when
|
|
#' combined with a high dose of Fertilizer B (interaction effect).
|
|
#'
|
|
|
|
#+ echo = FALSE
|
|
dat <- data.frame(
|
|
A = rep(1:2, each = 2),
|
|
B = rep(1:2, times = 2),
|
|
y = c(30, 30 + 5, 30 + 30, 30 + 30 + 5 + 12)
|
|
)
|
|
par(mai = c(.6, .6, .1, .1), mgp = c(2, .7, 0))
|
|
plot(y ~ A, dat, type = "n", xlim = c(0.8, 2.2), ylim = c(20, 80),
|
|
xlab = "Fertilizer A", ylab = "Yield (kg)", xaxt = "n")
|
|
lines(y ~ A, dat[dat$B == 1, ], col = "darkblue")
|
|
lines(y ~ A, dat[dat$B == 2, ], col = "darkblue")
|
|
lines(1:2, c(30 + 5, 30 + 30 + 5), lty = 2, col = "darkblue")
|
|
axis(1, 1:2, c("low", "high"))
|
|
arrows(2, 30 + 30 + 6, 2, 30 + 30 + 5 + 11, code = 3, length = 0.1,
|
|
col = "darkgray")
|
|
text(c(1, 2, 1, 2, 2.07), c(27, 55, 40, 80, 65 + 6),
|
|
c(expression(mu), expression(mu + alpha[2]),
|
|
expression(mu + beta[2]),
|
|
expression(mu + alpha[2] + beta[2] + (alpha * beta)[22]),
|
|
"12 kg")
|
|
)
|
|
text(1.5, 27 + 30/2, "Fertilizer B: low", srt = 21, col = "darkgray")
|
|
text(1.5, 38 + (30 + 12)/2, "Fertilizer B: high", srt = 32, col = "darkgray")
|
|
|
|
#' ## Model
|
|
|
|
#'
|
|
#' Assumptions
|
|
#'
|
|
#' - $Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} +
|
|
#' \varepsilon_{ijk}$
|
|
#' - $\varepsilon_{ijk} \sim N(0, \sigma^2) \text{ i.i.d.}$
|
|
#' - $i = 1, \dots, I$; $j = 1, \dots, J$; $k = 1, \dots, K$
|
|
#' - $\alpha_1 = \beta_1 := 0$
|
|
#'
|
|
#' Hypothesis
|
|
#'
|
|
#' - $H_0^{AB}\colon~ (\alpha\beta)_{ij} = 0 \text{ for all } i,j$
|
|
#'
|
|
|
|
#' ## Setup
|
|
|
|
set.seed(1704)
|
|
n <- 96
|
|
dat <- data.frame(
|
|
A = factor(rep(1:2, each = n/2), labels = c("low", "high")),
|
|
B = factor(rep(rep(1:2, each = n/4), 2), labels = c("low", "high"))
|
|
)
|
|
X <- model.matrix(~ A * B, dat)
|
|
unique(X)
|
|
beta <- c(mu = 30, a2 = 30, b2 = 5, ab22 = 12)
|
|
means <- X %*% beta
|
|
|
|
lattice::xyplot(I(means + rnorm(n, sd = 10)) ~ A, dat, groups = B,
|
|
type = c("g", "p", "a"), auto.key = TRUE, ylab = "Yield (kg)")
|
|
|
|
#' ## Parameter recovery
|
|
#+ cache = TRUE
|
|
|
|
out <- replicate(2000, {
|
|
y <- means + rnorm(n, sd = 10) # y = mu + a + b + ab + e
|
|
m <- aov(y ~ A * B, dat)
|
|
c(coef(m), sigma = sigma(m))
|
|
})
|
|
boxplot(t(out))
|
|
|
|
#' ## Power simulation
|
|
#+ cache = TRUE
|
|
|
|
pval <- replicate(2000, {
|
|
y <- means + rnorm(n, sd = 10)
|
|
m <- aov(y ~ A * B, dat)
|
|
summary(m)[[1]]$"Pr(>F)"[3] # test of interaction
|
|
})
|
|
mean(pval < 0.05)
|
|
|