Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Simulate Power Distributions
library(MBESS)
library(pwr)
nSims <- 100000 #number of simulated experiments
p <-numeric(nSims) #set up empty container for all simulated p-values
obs_pwr <-numeric(nSims) #set up empty container
t <-numeric(nSims) #set up empty container
d_all<-numeric(nSims)
N<-33 #number of participants
for(i in 1:nSims){ #for each simulated experiment
x<-rnorm(n = N, mean = 100, sd = 20) #produce 100 simulated participants
y<-rnorm(n = N, mean = 110, sd = 20) #produce 100 simulated participants
z<-t.test(x,y) #perform the t-test
d<-smd(Mean.1= mean(x), Mean.2=mean(y), s.1=sd(x), s.2=sd(y), n.1=N, n.2=N, Unbiased=TRUE)
d_all[i]<-d
obs_pwr[i]<-pwr.t.test(n = N, d = d, sig.level = 0.05, type = "two.sample", alternative = "two.sided")$power
p[i]<-z$p.value #get the p-value and store it
t[i]<-t.test(x, y, alternative = "two.sided", paired = FALSE, var.equal = TRUE, conf.level = 0.95)$statistic
}
#Calculate power in simulation
cat ("The power (in %) is")
sum(p < 0.05)/nSims*100
#now plot histograms of p-values, t-values, d, observed power
hist(p, main="Histogram of p-values", xlab=("Observed p-value"))
hist(t, main="Histogram of t-values", xlab=("Observed t-value"))
hist(d_all, main="Histogram of d", xlab=("Observed d"))
hist(obs_pwr, main="Histogram of observed power", xlab=("Observed power (%)"))
plot(p,obs_pwr)
abline(v = 0.05, lwd = 2, col = "red", lty = 2)
abline(h = 0.5, lwd = 2, col = "red", lty = 2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment