A simulation evaluating the confidence interval coverage for a difference in means, allowing for unequal variances, based on Welch's approximation for the degrees of freedom. The code demonstrates the basic structure of a simulation study in R.

Simulate Welch t-test

#----------------------------------------------
# data-generating model 
#----------------------------------------------

two_group_data <- function(iterations, n, p, var_ratio, delta) {
  Group <- c(rep("C", n * p), rep("T", n * (1 - p)))
  Y_C <- matrix(rnorm(iterations * n * p, mean = 0, sd = 1), n * p, iterations)
  Y_T <- matrix(rnorm(iterations * n * (1 - p), mean = delta, sd = sqrt(var_ratio)), n * (1 - p), iterations)
  dat <- data.frame(Group, rbind(Y_C, Y_T))
  return(dat)
}


#----------------------------------------------
# estimation procedures
#----------------------------------------------

CI_Welch <- function(dat, alpha = 0.05) {
  CI <- apply(dat[,-1], 2, function(X) t.test(X ~ dat$Group, var.equal=FALSE, conf.level = 1 - alpha)$conf.int)
  return(t(CI))
}


#----------------------------------------------
# performance statistics
#----------------------------------------------

coverage <- function(CI, delta) {
  covered <- (CI[,1] < delta) & (delta < CI[,2])
  return(mean(covered))
}


#----------------------------------------------
# Simulation driver
#----------------------------------------------

run_sim <- function(iterations, n, p, var_ratio, delta) {
  dat <- two_group_data(iterations, n, p, var_ratio, delta)
  Welch <- coverage(CI_Welch(dat), delta)
  return(c(Welch = Welch))
}


#------------------------------------
# Design parameters
#------------------------------------

iterations <- 1000
n <- 12 * (1:5)
p <- 1 / c(2,3,4,6)
R <- 2^(-2:2)
delta <- 0

parms <- expand.grid(iterations=iterations, n=n, p=p, var_ratio=R, delta=delta)


#--------------------------------------
# run simulations in serial
#--------------------------------------

library(plyr)
set.seed(20110325)
system.time(results <- maply(parms, .fun = run_sim, .drop=FALSE))
save(results, file="BF results.Rdata")


#--------------------------------------
# plot results
#--------------------------------------

library(reshape)
library(ggplot2)

load("BF results.Rdata")
dimnames(results)$p <- c("1/6","1/4","1/3","1/2")
names(dimnames(results))[4] <- "R"
results_long <- melt(results)

qplot(data = results_long,
      x = n, y = value,
      geom = "line") +
  facet_grid(R ~ p, labeller = label_both) + 
  labs(y = "Coverage") + 
  theme_bw()