# Synthetic balanced cluster-randomized example; base R only.
# Treatment is allocated to clusters. Cluster-mean t test respects that unit.
set.seed(4004)
k <- 20; m <- 25
cluster <- rep(seq_len(k), each=m)
arm <- rep(rep(0:1, each=k/2), each=m)
u <- rnorm(k, sd=1)
y <- 0.5*arm + u[cluster] + rnorm(k*m)
d <- data.frame(cluster, arm, y)
naive <- t.test(y ~ arm, data=d, var.equal=TRUE)
means <- aggregate(y ~ cluster + arm, d, mean)
valid <- t.test(y ~ arm, data=means, var.equal=TRUE)
cat(sprintf('Participants=%d; clusters=%d; cluster size=%d\n', nrow(d),k,m))
cat(sprintf('Difference=%.3f; naive SE=%.3f; cluster SE=%.3f\n', diff(tapply(y,arm,mean)), naive$stderr, valid$stderr))
cat(sprintf('Naive CI=[%.3f, %.3f]; cluster CI=[%.3f, %.3f]\n', -naive$conf.int[2],-naive$conf.int[1],-valid$conf.int[2],-valid$conf.int[1]))
cat(sprintf('Naive p=%.4f; cluster p=%.4f; cluster df=%.0f\n',naive$p.value,valid$p.value,valid$parameter))
stopifnot(valid$stderr > naive$stderr, nrow(means)==k)

# Maintainer-only export; run from repository root. Synthetic data only.
if (Sys.getenv("CLINIC_EXPORT_DATA") == "1") {
write.csv(d,'analysis-clinic/examples/case-004-shared-data.csv',row.names=FALSE,quote=FALSE)
}