Two-sample t-test and two-by-two ANOVA

Power simulation and power curves

Two-sample t-test

Application context

Listening experiment

• Task of each participant is to repeatedly adjust the frequency of a comparison tone to sound equal in pitch to a 1000-Hz standard tone
• Mean adjustment estimates the point of subjective equality \(\mu\)
• Two participants will take part in the experiment providing adjustments \(X\) and \(Y\)
• Goal is to detect a difference between their points of subjective equality \(\mu_x\) and \(\mu_y\) of 4 Hz

Model

Assumptions

• \(X_1, \ldots, X_n \sim N(\mu_x, \sigma_x^2)\) i.i.d.
• \(Y_1, \ldots, Y_m \sim N(\mu_y, \sigma_y^2)\) i.i.d.
• both samples independent
• \(\sigma_x^2 = \sigma_y^2\) but unknown

Hypothesis

• H\(_0\colon~ \mu_x - \mu_y = \delta = 0\)

Power simulation

n <- 110; m <- 110
pval <- replicate(2000, {
  x <- rnorm(n, mean = 1000 + 4, sd = 10)         # Participant 1 responses
  y <- rnorm(m, mean = 1000,     sd = 10)         # Participant 2 responses
  t.test(x, y, mu = 0, var.equal = TRUE)$p.value
})
mean(pval < 0.05)

## [1] 0.8345

Power curves

Turn into a function of n and effect size

pwrFun <- function(n = 30, d = 4, sd = 10,
                   nrep = 50) {
  n <- n; m <- n
  pval <- replicate(nrep, {
    x <- rnorm(n, mean = 1000 + d, sd = sd)
    y <- rnorm(m, mean = 1000,     sd = sd)
    t.test(x, y, mu = 0, var.equal = TRUE)$p.value
  })
  mean(pval < 0.05)
}

Set up conditions and call power function

cond <- expand.grid(d = 0:5,
                    n = c(50, 75, 100, 125))
system.time(
  cond$pwr <- mapply(pwrFun, n = cond$n, d = cond$d,
                     MoreArgs = list(nrep = 500))
)

##        User      System verstrichen 
##       1.435       0.007       1.441

## Plot results
lattice::xyplot(pwr ~ d, cond, groups = n, type = c("g", "b"),
                auto.key = list(corner = c(0, 1)))

Parameter recovery

n <- 100

out <- replicate(2000, {
  x <- rnorm(n, mean = 1000 + 4, sd = 10)
  y <- rnorm(n, mean = 1000,     sd = 10)
  t <- t.test(x, y, mu = 0, var.equal = TRUE)
  c(
    effect = -as.numeric(diff(t$estimate)),
    sd.pool = t$stderr * sqrt(2*n) / 2
  )
})

rowMeans(out)

##   effect  sd.pool 
## 3.981564 9.997941

boxplot(t(out))

Two-by-two ANOVA

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).

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)

##    (Intercept) Ahigh Bhigh Ahigh:Bhigh
## 1            1     0     0           0
## 25           1     0     1           0
## 49           1     1     0           0
## 73           1     1     1           1

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

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

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)

## [1] 0.817