# Goal: Simulation to study size and power in a simple problem. set.seed(101) # The data generating process: a simple uniform distribution with stated mean dgp <- function(N,mu) {runif(N)-0.5+mu} # Simulate one FIXED hypothesis test for H0:mu=0, given a true mu for a sample size N one.test <- function(N, truemu) { x <- dgp(N,truemu) muhat <- mean(x) s <- sd(x)/sqrt(N) # Under the null, the distribution of the mean has standard error s threshold <- 1.96*s (muhat < -threshold) || (muhat > threshold) } # Return of TRUE means reject the null # Do one experiment, where the fixed H0:mu=0 is run Nexperiments times with a sample size N. # We return only one number: the fraction of the time that H0 is rejected. experiment <- function(Nexperiments, N, truemu) { sum(replicate(Nexperiments, one.test(N, truemu)))/Nexperiments } # Measure the size of a test, i.e. rejections when H0 is true experiment(10000, 50, 0) # Measurement with sample size of 50, and true mu of 0. # Power study: I.e. Pr(rejection) when H0 is false # (one special case in here is when the H0 is actually true) muvalues <- seq(-.15,.15,.01) # When true mu < -0.15 and when true mu > 0.15, # the Pr(rejection) veers to 1 (full power) and it's not interesting. # First do this with sample size of 50 results <- NULL for (truth in muvalues) { results <- c(results, experiment(10000, 50, truth)) } par(mai=c(.8,.8,.2,.2)) plot(muvalues, results, type="l", lwd=2, ylim=c(0,1), xlab="True mu", ylab="Pr(Rejection of H0:mu=0)") abline(h=0.05, lty=2) # Now repeat this with sample size of 100 (should yield a higher power) results <- NULL for (truth in muvalues) { results <- c(results, experiment(10000, 100, truth)) } lines(muvalues, results, lwd=2, col="blue") legend(x=-0.15, y=.2, lwd=c(2,1,2), lty=c(1,2,1), cex=.8, col=c("black","black","blue"), bty="n", legend=c("N=50", "Size, 0.05", "N=100"))