Skip to content

Instantly share code, notes, and snippets.

@jepusto
Last active December 29, 2015 14:49
Show Gist options
  • Save jepusto/7686463 to your computer and use it in GitHub Desktop.
Save jepusto/7686463 to your computer and use it in GitHub Desktop.
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.
#----------------------------------------------
# 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment