An adaptation of the Welch t-test simulation, for running in parallel on the Stampede server of the Texas Advanced Computing Cluster.

Simulate Welch t-test in parallel

#----------------------------------------------
# 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, seed=NULL) {
  if (!is.null(seed)) set.seed(seed)
  dat <- two_group_data(iterations, n, p, var_ratio, delta)
  Welch <- coverage(CI_welch(dat), delta)
  return(c(Welch = Welch))
}

source_func <- ls()

#------------------------------------
# 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)
parms$seed <- round(runif(nrow(parms)) * 2^30)

#--------------------------------------
# run simulations in parallel
#--------------------------------------

library(Rmpi)
library(snow)
library(foreach)
library(iterators)
library(doSNOW)
library(plyr)

# set up parallel processing
cluster <- getMPIcluster()
registerDoSNOW(cluster)

# export source functions
clusterExport(cluster, source_func)

# execute simulations
BFtime <- system.time(BFresults <- mdply(parms, .fun = run_sim, .parallel=TRUE))
stopCluster(cluster)

print(BFtime)

save(BFresults, BFtime, file="BFresults.Rdata")