t1err-sim.R

### Utility functions

# Do a single t test simulation
# report the p value
ttest.sim = function(n, func, true.mean = 0, alpha = 0.05){
  x = func(n) - true.mean
  t.test(x)$p.value
}

# Do a sequence of M t tests, report significance 
# as a proportion of M
Msim = function(n, M, func, true.mean, alpha = 0.05)
  mean(replicate(M, ttest.sim(n, func, true.mean, alpha))<alpha)


### Our populations that violate normality

# beta with values stuck near 0 and 1
beta.func1 = function(n){
  rbeta(n, .1, .1)
}

# t distribution with 3 degrees of freedom
# (volatile sd)
t.func = function(n){
  rt(n,3)
}


### simulation setup

# sequence of sample sizes at which to simulate
Ns = seq(5,50,5)
# number of simulations per sample size
M = 100000
# nominal alpha
alpha = 0.05

### Perform simulations
beta.sim1 = sapply(Ns, Msim, M = M, func = beta.func1, true.mean = .5, alpha = alpha)
t.sim = sapply(Ns, Msim, M = M, func = t.func, true.mean = 0, alpha = alpha)
norm.sim = sapply(Ns, Msim, M = M, func = rnorm, true.mean = 0, alpha = alpha)


### plot the results
par(las=1, lwd=2)
plot(Ns,norm.sim, pch=19, ty='l', ylim = c(0,.1), xlab="Sample size", ylab="Type I error rate", log = "x")
abline(h = alpha, col="red", lty=2)
lines(Ns, beta.sim1, col="blue")
lines(Ns, t.sim, col="green")
legend("topright", legend = c("normal","beta", "t","Nominal alpha"), lty=c(1,1,1,2), col = c("black","blue","green","red"))