R code for data simulations with both large and small effect sizes to see the variation in p-values and Bayes factors. Inspiration comes from work done by Geoff Cumming and Daniel Lakens.

############################################

########Simulation of data with large effect size#######

############################################

##set up

library(BayesFactor)

library(MBESS)

library(pwr)

library(ggplot2)

options(scipen=20)

N.Sim = 1000 #number of simulated experiments

n<-30 #sample size in each group

sd = 12.5 #average cohen’s d of .8

#Cohen’s D independent samples

d.indt = function (m1 = m1, m2 = m2, sd1 = sd1, sd2 = sd2, n1 = n, n2 = n) {

spooled = sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n1+n2 – 2))

d = (m1 – m2) / spooled

d

}

#z test

ztest = function(m1 = m1, m2 = m2, sd1 = sd, sd2 = sd, n1 = n, n2 = n){

denom = sqrt(((sd^2)/n1)+((sd^2)/n2))

z = (m1-m2)/denom

p = 2*pnorm(-abs(z))

list(Z = z, p = p)

}

mydata = data.frame(sim = 1:100, P.Value = 1:100, PosthocPower = 1:100,

CohensD = 1:100,

CILo = 1:100, CIhi = 1:100, BayesFactor = 1:100)

round = 0

#Simulation Loop

for(i in 1:N.Sim){

#data

x<-rnorm(n = n, mean = 20, sd = sd) #Cohen’s d of .8

y<-rnorm(n = n, mean = 10, sd = sd) #Cohen’s d of .8

#descriptives

m1 = mean(x)

m2 = mean(y)

sd1 = sd(x)

sd2 = sd(y)

n1 = n

n2 = n

#z-test

z=ztest(m1 = m1, m2 = m2, sd1 = sd, sd2 = sd, n1 = n, n2 = n)

round = round+1

mydata$P.Value[round] = z$p

CI = ci.smd(ncp = z$Z, n.1=n, n.2=n, conf.level=0.95) #CI of Cohen’s d .8

mydata$CILo[round] = CI$Lower.Conf.Limit.smd

mydata$CIhi[round] = CI$Upper.Conf.Limit.smd

D = d.indt(m1 = m1, m2 = m2, sd1 = sd1, sd2 = sd2, n1 = n, n2 = n)

mydata$CohensD[round] = D

power = pwr.t.test(n=n, sig.level = .05, d=D,

type = “two.sample”,

alternative = “two.sided”)

mydata$PosthocPower[round] = power$power

#BayesFactor

BF10.current<-exp(ttest.tstat(z$Z,n,n,rscale=sqrt(2)/2)$bf)

mydata$BayesFactor[round] = BF10.current #add value to the vector

}

summary(mydata)

hist(exp(mydata$P.Value), breaks = 100)

hist(mydata$BayesFactor, breaks = 100)

hist(mydata$PosthocPower, breaks = 100)

############################################

########Simulation of data with small effect size#######

############################################

#set up

N.Sim = 100 #number of simulated experiments

n<-30 #sample size in each group

sd = 25 #cohen’s d of .2

mydata2 = data.frame(sim = 1:100, P.Value = 1:100, PosthocPower = 1:100,

CohensD = 1:100,

CILo = 1:100, CIhi = 1:100, BayesFactor = 1:100)

round = 0

#Simulation Loop

for(i in 1:N.Sim){

#data

x<-rnorm(n = n, mean = 15, sd = sd)

y<-rnorm(n = n, mean = 10, sd = sd)

#descriptives

m1 = mean(x)

m2 = mean(y)

sd1 = sd(x)

sd2 = sd(y)

n1 = n

n2 = n

#z-test

z=ztest(m1 = m1, m2 = m2, sd1 = sd, sd2 = sd, n1 = n, n2 = n)

round = round+1

mydata2$P.Value[round] = z$p

CI = ci.smd(smd = z$Z, n.1=n, n.2=n, conf.level=0.95)

mydata2$CILo[round] = CI$Lower.Conf.Limit.smd

mydata2$CIhi[round] = CI$Upper.Conf.Limit.smd

D = d.indt(m1 = m1, m2 = m2, sd1 = sd1, sd2 = sd2, n1 = n, n2 = n)

mydata2$CohensD[round] = D

power = pwr.t.test(n=n, sig.level = .05, d=D,

type = “two.sample”,

alternative = “two.sided”)

mydata2$PosthocPower[round] = power$power

#BayesFactor

BF10.current<-exp(ttest.tstat(z$Z,n,n,rscale=sqrt(2)/2)$bf)

mydata2$BayesFactor[round] = BF10.current #add value to the vector

}

summary(mydata2)

hist(exp(mydata2$P.Value), breaks = 100)

hist(mydata2$BayesFactor, breaks = 100)

hist(mydata2$PosthocPower, breaks = 100)

################################################################

## Recent Comments