Three-factor chi-squared test of independence

3d_chisq.R

# Demonstrate that this is chi-squared distributed
# with xyz - (x-1) - (y-1) - (z-1) -1 d.f.
stat <- numeric(1000)
for (i in seq_along(stat)) {
    N <- 1000
    x <- rmultinom(1, size=N, prob=rep(1/27, 27))
    dim(x) <- c(3,3,3)

    p1 <- vapply(1:3, function(i) sum(x[i,,]), 1) / N
    p2 <- vapply(1:3, function(i) sum(x[,i,]), 1) / N
    p3 <- vapply(1:3, function(i) sum(x[,,i]), 1) / N

    e <- outer(outer(p1, p2), p3) * N
    stat[i] <- sum((x - e)^2/e)
}