n1<-35                      # Sample size for x
n2<-25                      # Sample size for y
sigma1<-50                  # Population sd for x
sigma2<-50                  # Population sd for y
mu1<-40                     # Population mean for x
mu2<-0                      # Population mean for y
nsim<-10000                 # Number of trials
pval1<-numeric(nsim)        # Storage for p-value of pooled-t
pval2<-numeric(nsim)        # Storage for p-value of Welch-t
for(i in 1:nsim){
xsim<-rnorm(n1,mu1,sigma1)  # Generate x-data
ysim<-rnorm(n2,mu2,sigma2)  # Generate y-data
pval1[i]<-t.test(xsim,ysim,alternative="two.sided",var.equal=TRUE)$p.value
pval2[i]<-t.test(xsim,ysim,alternative="two.sided",var.equal=FALSE)$p.value
}
cbind(mean(pval1<.05),mean(pval2<.05))

# Write a function 
welch_sim<-function(nsim=10000, n1=30, n2=30, mu1=0, mu2=0, sigma1=1,sigma2=1){
pval1<-numeric(nsim)        # Storage for p-value of pooled-t
pval2<-numeric(nsim)        # Storage for p-value of Welch-t
for(i in 1:nsim){
xsim<-rnorm(n1,mu1,sigma1)  # Generate x-data
ysim<-rnorm(n2,mu2,sigma2)  # Generate y-data
pval1[i]<-t.test(xsim,ysim,alternative="two.sided",var.equal=TRUE)$p.value
pval2[i]<-t.test(xsim,ysim,alternative="two.sided",var.equal=FALSE)$p.value
}                           # End of for() loop
return(c(mean(pval1<.05),mean(pval2<.05)))
}                           # End of function()

# Call the function
welch_sim(10000,35,25,0,0,50,50)
welch_sim(10000,35,25,10,0,50,50)
welch_sim(10000,35,25,20,0,50,50)
welch_sim(10000,35,25,30,0,50,50)
welch_sim(10000,35,25,40,0,50,50)
# Use unequal variance
welch_sim(10000,35,25,mu1=0,0,20,50)
welch_sim(10000,35,25,mu1=0,0,50,20)

# Use for() loop to automate
muvec<-c(-40,-35,-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,40)
outmat<-matrix(rep(0,2*length(muvec)), ncol=2)
for(j in 1:length(muvec)){
  outmat[j,]<-welch_sim(10000,35,25,muvec[j],0,20,50)
}
outmat
colnames(outmat)<-c("pooled","welch")
cbind(muvec,outmat)
power1<-outmat[,1]
power2<-outmat[,2]
plot(power1~muvec,ylim=c(0,1),type="l",xlab="mu1-mu2",ylab="Power",
     main="Power curves")
lines(power2~muvec,col="red",type="l",lty=2)
legend("topleft",legend=c("pooled","welch"),lty=c(1,2))
abline(h=.05,col="blue",lty=6)