Independent t-test: Power simulation and power curves
================
Last modified: 2026-01-09

## 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

``` r
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.831

## Power curves

Turn into a function of n and effect size

``` r
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

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

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

![](powersim-ttest_files/figure-gfm/unnamed-chunk-3-1.png)<!-- -->

## Parameter recovery

``` r
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.978400 9.993142

``` r
boxplot(t(out))
```

![](powersim-ttest_files/figure-gfm/unnamed-chunk-4-1.png)<!-- -->