welch_student.R

welch_student.R

require(car) #Car package required for Levene's test

n1<-38 #size condition x
n2<-22 #size condition y
sd1<-1.11 #sd condition x
sd2<-1.84 #sd condition y 
m1<-0
m2<-0
trueD<-(m2-m1)/(sqrt((((n1 - 1)*((sd1^2))) + (n2 - 1)*((sd2^2)))/((n1+n2)-2)))
trueD

nSims <- 5000 #number of simulated experiments (large numbers might take a while)
p1 <-numeric(nSims) #set up empty container for all simulated Student's t-test p-values 
p2 <-numeric(nSims) #set up empty container for all simulated Welch's t-test p-values
pvalueLevene<-numeric(nSims) #set up empty container for all Levene's t-test p-values
#create variables for dataframe
catx<-rep("x",n1)
caty<-rep("y",n2)
condition<- c(catx,caty)

#run simulations
for(i in 1:nSims){ #for each simulated experiment
  sim_x<-rnorm(n = n1, mean = m1, sd = sd1) #simulate participants condition x
  sim_y<-rnorm(n = n2, mean = m2, sd = sd2) #simulate participants condition y
  #perform Student and Welch t-test
  p1[i]<-t.test(sim_x,sim_y, alternative = "two.sided", var.equal = TRUE)$p.value #perform the t-test and store p-value
  p2[i]<-t.test(sim_x,sim_y, alternative = "two.sided", var.equal = FALSE)$p.value #perform the t-test and store p-value
  #create dataframe for levene's test
  xy<- c(sim_x,sim_y)
  alldata<-data.frame(xy,condition)
  #perform Levene's test
  pvalueLevene[i]<-leveneTest(alldata$xy ~ alldata$condition, data = alldata)$"Pr(>F)"[1:1]
}

#empty dataframe
alldata$xy<-NULL
alldata$condition<-NULL

#now plot the histogram for p-value distributions
hist(p1, main="Histogram of p-values from Student's t-test under the null ", xlab=("Observed p-value"))
hist(p2, main="Histogram of p-values from Welch's t-test under the null", xlab=("Observed p-value"))

#create plot of p-values (just section p < .015 plotted)
plot(p1,p2,ylim=c(0,0.15),xlim=c(0,0.15),xlab="p-values from Student's t-test", ylab="p-values from Welch's t-test")
abline(h=0.05) #0.05 line for Welch's test
abline(v=0.05) #0.05 line for Student's t 
abline(0,1, col="red", lwd=4) #draw diagonal red line

#Calculate Type 1 error in simulation for Student t-test
var_ratio<-sd1^2/sd2^2
errorrate<-sum(p1 < 0.05)/nSims*100 
cat ("The Type 1 error rate (if the true effect is 0) or the power using Student's t-test is ",errorrate," with a ratio between the variances of ",var_ratio,". Note the nominal Type 1 error rate should be 5%, and equal variances imply a variance ratio of 1")

#Calculate Type 1 error in simulation for Welch's t-test
errorrate2<-sum(p2 < 0.05)/nSims*100 
cat ("The Type 1 error rate (if the true effect is 0) or the power using Welch's t-test is ",errorrate2," with a ratio between the variances of ",var_ratio,". Note the nominal Type 1 error rate should be 5%, and equal variances imply a variance ratio of 1")

observedpowerLevene<-sum(pvalueLevene < 0.05)/nSims*100 
cat ("The observed power for the Levene test is ",observedpowerLevene)