## @knitr lecture9prep library(mvtnorm) library(ggplot2) ## @knitr ztest pz <- function(sample, mu = 0, tails=2) { #calculate z z <- (mean(sample) - mu)/(sd(sample)/sqrt(length(sample))) #return the two-tailed answer return(tails * pnorm(abs(z), lower.tail = F)) } ## @knitr ztest2 set.seed(697) pz(rnorm(5000)) pz(rnorm(5, mean=1)) ## @knitr powZOneN pzPower <- function(n, effect=0, n.sims=500, alpha=0.05){ #a vector of p-values p <- numeric(n.sims) #just your run of the mill power by #simulation for one choice of n and effect for(i in 1:n.sims){ samp <- rnorm(n, effect) #Hey - it's our pz function! p[i] <- pz(samp) } #calculate power and return it 1-sum(p > alpha)/n.sims } ## @knitr powZManyN pzPowerSample <- function(nVec, effect=0, n.sims=500, alpha=0.05, ...){ #Create a vector for the output we're going to get powVec <- numeric(length(nVec)) #loop over each sample size for(i in 1:length(nVec)){ #What's this? Another new function! powVec[i] <- pzPower(nVec[i], effect, n.sims, alpha) } #make the plot & return output plot(nVec, powVec, ...) powVec } ## @knitr powZdemo pzPowerSample(3:10, effect=1, type="b")