treatment<-c(rnormp2(larger.bimodal.proportion*100000,0+effect,treatment.sd1,treatment.shape1),
rnormp2(smaller.bimodal.proportion*100000,0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
if(batchmode!=T){
if(using.RStudio!=T){
dev.new(width=600, height=500)
}
}
if(using.RStudio==T){
live.plot.in.own.window=F
live.plot=F
}
layout(matrix(c(1,2,3,4,5,6), 3, 2, byrow = TRUE))
hist(treatment, col=rgb(0,0,1,.5), freq=F,
breaks=seq(min(c(treatment,control))-1,max(c(treatment,control))+1,by=.1),
xlab="Value", ylim=c(0,1.2),
main=c("Population Distributions",
paste("Difference in Means=", round(effect,2)),
paste("Bimodal Proportions=", round(larger.bimodal.proportion,2), ",", round(smaller.bimodal.proportion,2),
" Bimodal Offset=", round(bimodal.offset,2))
),
sub=paste("Bimodal Effect on Treatment Only=",bimodal.treatment.only)
)
hist(control, col=rgb(1,0,0,.5), freq=F,
breaks=seq(min(c(treatment,control))-1,max(c(treatment,control))+1,by=.1),
add=T)
legend(x="topleft", legend=c(paste("Treatment:",
"mu=", 0+effect,
", sd1=", round(treatment.sd1,2),
", p1=", round(treatment.shape1,2),
", sd2=", round(treatment.sd2,2),
", p2=", round(treatment.shape2,2)
),
paste("Control:",
"mu=",0,
", sd1=", round(control.sd1,2),
", p1=", round(control.shape1,2),
", sd2=", round(control.sd2,2),
", p2=", round(control.shape2,2))
),
col=c(rgb(0,0,1,.5),rgb(1,0,0,.5)), lwd=6, bty="n")
plot(initial.samp.size,1, type="l", lwd=2, col="White",
xlim=c(initial.samp.size,(max.samps+samp.interval)), xlab="Sample Size (per group)",
ylim=c(0,1), ylab="P-value",
main=c("Representative P-value Traces (Significant: P < 0.05)",
paste("Outlier Cutoff(",outlier.method,") =",round(outlier.cutoff,2),
", Outlier Bias (High,Low) =", "(", round(outlier.bias.high,2), ",", round(outlier.bias.low,2), ")",
sep=""),
paste("# of Outcomes=", number.of.outcomes)
),
sub=paste("Sequential Sampling=", sequential, ", U-test=", also.perform.Utest)
)
abline(h=.05, lwd=2)
###Start Simulation###
all.results=NULL
pb <- txtProgressBar(min = 0, max = runs*max.samps, style = 3)
for(r in 1:runs){
###Start Run###
run.result=NULL
experiment.results=vector("list", number.of.outcomes)
for(t in 1:number.of.outcomes){
if(bimodal.treatment.only==T){
random.sample<-rep(1,initial.samp.size)
}else{
random.sample<-bimodal.pop[runif(initial.samp.size,1,1000)]
}
control<-c(rnormp2(length(which(random.sample==1)),0,control.sd1,control.shape1),
rnormp2(length(which(random.sample==0)),0+bimodal.offset,control.sd2,control.shape2))
random.sample<-bimodal.pop[runif(initial.samp.size,1,1000)]
treatment<-c(rnormp2(length(which(random.sample==1)),0+effect,treatment.sd1,treatment.shape1),
rnormp2(length(which(random.sample==0)),0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
experiment.result<-rbind(control,treatment)
experiment.results[[t]]<-experiment.result
}
x=initial.samp.size
while(x <(max.samps+samp.interval)){
sample.result=NULL
n.cutoff=min(n.cutoff,x)
test.results=NULL
for(t in 1:number.of.outcomes){
control<-experiment.results[[t]][1,]
treatment<-experiment.results[[t]][2,]
if(also.perform.Utest==T){
pval1<-t.test(treatment,control)$p.value
pval2<-wilcox.test(treatment,control)$p.value
pval<-min(pval1,pval2)
#Drop Outliers
if(drop.outliers==T){
if(pval>.05){
control.temp<-control
treatment.temp<-treatment
if(outlier.method=="SD"){
high.cutoff<-mean(control.temp, na.rm=T)+outlier.cutoff*sd(control.temp, na.rm=T)
low.cutoff<-mean(control.temp, na.rm=T)-outlier.cutoff*sd(control.temp, na.rm=T)
high.outliers<-which(control.temp>high.cutoff)
low.outliers<-which(control.temp
control.temp[high.outliers]<-NA
control.temp[low.outliers]<-NA
high.cutoff<-mean(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*sd(treatment.temp, na.rm=T)
low.cutoff<-mean(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*sd(treatment.temp, na.rm=T)
high.outliers<-which(treatment.temp>high.cutoff)
low.outliers<-which(treatment.temp
treatment.temp[high.outliers]<-NA
treatment.temp[low.outliers]<-NA
}
if(outlier.method=="MAD"){
high.cutoff<-median(control.temp, na.rm=T)+outlier.cutoff*mad(control.temp, na.rm=T)
low.cutoff<-median(control.temp, na.rm=T)-outlier.cutoff*mad(control.temp, na.rm=T)
high.outliers<-which(control.temp>high.cutoff)
low.outliers<-which(control.temp
control.temp[high.outliers]<-NA
control.temp[low.outliers]<-NA
high.cutoff<-median(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*mad(treatment.temp, na.rm=T)
low.cutoff<-median(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*mad(treatment.temp, na.rm=T)
high.outliers<-which(treatment.temp>high.cutoff)
low.outliers<-which(treatment.temp
treatment.temp[high.outliers]<-NA
treatment.temp[low.outliers]<-NA
}
if(
length(treatment.temp[which(!is.na(treatment.temp))]) > (n.cutoff - 1) &&
length(control.temp[which(!is.na(control.temp))]) > (n.cutoff - 1)
){
pval1<-t.test(treatment.temp,control.temp)$p.value
pval2<-wilcox.test(treatment.temp,control.temp)$p.value
pval<-min(pval1,pval2)
control<-control.temp
treatment<-treatment.temp
#Replace Outliers for Small N
}else{
while(
length(treatment.temp[which(!is.na(treatment.temp))]) < n.cutoff &&
length(control.temp[which(!is.na(control.temp))]) < n.cutoff
){
if(bimodal.treatment.only==T){
random.sample<-rep(1,length(which(is.na(control.temp))))
}else{
random.sample<-bimodal.pop[runif(length(which(is.na(control.temp))),1,1000)]
}
control.temp[which(is.na(control.temp))]<-c(rnormp2(length(which(random.sample==1)),
0,control.sd1,control.shape1),
rnormp2(length(which(random.sample==0)),
0+bimodal.offset,control.sd2,control.shape2))
random.sample<-bimodal.pop[runif(length(which(is.na(treatment.temp))),1,1000)]
treatment.temp[which(is.na(treatment.temp))]<-c(rnormp2(length(which(random.sample==1)),
0+effect,treatment.sd1,treatment.shape1),
rnormp2(length(which(random.sample==0)),
0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
if(outlier.method=="SD"){
high.cutoff<-mean(control.temp, na.rm=T)+outlier.cutoff*sd(control.temp, na.rm=T)
low.cutoff<-mean(control.temp, na.rm=T)-outlier.cutoff*sd(control.temp, na.rm=T)
high.outliers<-which(control.temp>high.cutoff)
low.outliers<-which(control.temp
control.temp[high.outliers]<-NA
control.temp[low.outliers]<-NA
high.cutoff<-mean(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*sd(treatment.temp, na.rm=T)
low.cutoff<-mean(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*sd(treatment.temp, na.rm=T)
high.outliers<-which(treatment.temp>high.cutoff)
low.outliers<-which(treatment.temp
treatment.temp[high.outliers]<-NA
treatment.temp[low.outliers]<-NA
}
if(outlier.method=="MAD"){
high.cutoff<-median(control.temp, na.rm=T)+outlier.cutoff*mad(control.temp, na.rm=T)
low.cutoff<-median(control.temp, na.rm=T)-outlier.cutoff*mad(control.temp, na.rm=T)
high.outliers<-which(control.temp>high.cutoff)
low.outliers<-which(control.temp
control.temp[high.outliers]<-NA
control.temp[low.outliers]<-NA
high.cutoff<-median(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*mad(treatment.temp, na.rm=T)
low.cutoff<-median(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*mad(treatment.temp, na.rm=T)
high.outliers<-which(treatment.temp>high.cutoff)
low.outliers<-which(treatment.temp
treatment.temp[high.outliers]<-NA
treatment.temp[low.outliers]<-NA
}
} #end while sample size too small loop
pval1<-t.test(treatment.temp,control.temp)$p.value
pval2<-wilcox.test(treatment.temp,control.temp)$p.value
pval<-min(pval1,pval2)
control<-control.temp
treatment<-treatment.temp
} # end replace pvals
} # end if pval...
} #end drop outliers
# Do if U-test==F
}else{
pval<-t.test(treatment,control)$p.value
#Drop Outliers
if(drop.outliers==T){
if(pval>.05){
control.temp<-control
treatment.temp<-treatment
if(outlier.method=="SD"){
high.cutoff<-mean(control.temp, na.rm=T)+outlier.cutoff*sd(control.temp, na.rm=T)
low.cutoff<-mean(control.temp, na.rm=T)-outlier.cutoff*sd(control.temp, na.rm=T)
high.outliers<-which(control.temp>high.cutoff)
low.outliers<-which(control.temp
control.temp[high.outliers]<-NA
control.temp[low.outliers]<-NA
high.cutoff<-mean(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*sd(treatment.temp, na.rm=T)
low.cutoff<-mean(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*sd(treatment.temp, na.rm=T)
high.outliers<-which(treatment.temp>high.cutoff)
low.outliers<-which(treatment.temp
treatment.temp[high.outliers]<-NA
treatment.temp[low.outliers]<-NA
}
if(outlier.method=="MAD"){
high.cutoff<-median(control.temp, na.rm=T)+outlier.cutoff*mad(control.temp, na.rm=T)
low.cutoff<-median(control.temp, na.rm=T)-outlier.cutoff*mad(control.temp, na.rm=T)
high.outliers<-which(control.temp>high.cutoff)
low.outliers<-which(control.temp
control.temp[high.outliers]<-NA
control.temp[low.outliers]<-NA
high.cutoff<-median(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*mad(treatment.temp, na.rm=T)
low.cutoff<-median(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*mad(treatment.temp, na.rm=T)
high.outliers<-which(treatment.temp>high.cutoff)
low.outliers<-which(treatment.temp
treatment.temp[high.outliers]<-NA
treatment.temp[low.outliers]<-NA
}
if(
length(treatment.temp[which(!is.na(treatment.temp))]) > (n.cutoff - 1) &&
length(control.temp[which(!is.na(control.temp))]) > (n.cutoff - 1)
){
pval<-t.test(control.temp,treatment.temp)$p.value
control<-control.temp
treatment<-treatment.temp
#Replace Outliers for Small N
}else{
while(
length(treatment.temp[which(!is.na(treatment.temp))]) < n.cutoff &&
length(control.temp[which(!is.na(control.temp))]) < n.cutoff
){
if(bimodal.treatment.only==T){
random.sample<-rep(1,length(which(is.na(control.temp))))
}else{
random.sample<-bimodal.pop[runif(length(which(is.na(control.temp))),1,1000)]
}
control.temp[which(is.na(control.temp))]<-c(rnormp2(length(which(random.sample==1)),
0,control.sd1,control.shape1),
rnormp2(length(which(random.sample==0)),
0+bimodal.offset,control.sd2,control.shape2))
random.sample<-bimodal.pop[runif(length(which(is.na(treatment.temp))),1,1000)]
treatment.temp[which(is.na(treatment.temp))]<-c(rnormp2(length(which(random.sample==1)),
0+effect,treatment.sd1,treatment.shape1),
rnormp2(length(which(random.sample==0)),
0+effect+bimodal.offset,treatment.sd2,treatment.shape2))
if(outlier.method=="SD"){
high.cutoff<-mean(control.temp, na.rm=T)+outlier.cutoff*sd(control.temp, na.rm=T)
low.cutoff<-mean(control.temp, na.rm=T)-outlier.cutoff*sd(control.temp, na.rm=T)
high.outliers<-which(control.temp>high.cutoff)
low.outliers<-which(control.temp
control.temp[high.outliers]<-NA
control.temp[low.outliers]<-NA
high.cutoff<-mean(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*sd(treatment.temp, na.rm=T)
low.cutoff<-mean(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*sd(treatment.temp, na.rm=T)
high.outliers<-which(treatment.temp>high.cutoff)
low.outliers<-which(treatment.temp
treatment.temp[high.outliers]<-NA
treatment.temp[low.outliers]<-NA
}
if(outlier.method=="MAD"){
high.cutoff<-median(control.temp, na.rm=T)+outlier.cutoff*mad(control.temp, na.rm=T)
low.cutoff<-median(control.temp, na.rm=T)-outlier.cutoff*mad(control.temp, na.rm=T)
high.outliers<-which(control.temp>high.cutoff)
low.outliers<-which(control.temp
control.temp[high.outliers]<-NA
control.temp[low.outliers]<-NA
high.cutoff<-median(treatment.temp, na.rm=T)+outlier.bias.high*outlier.cutoff*mad(treatment.temp, na.rm=T)
low.cutoff<-median(treatment.temp, na.rm=T)-outlier.bias.low*outlier.cutoff*mad(treatment.temp, na.rm=T)
high.outliers<-which(treatment.temp>high.cutoff)
low.outliers<-which(treatment.temp
treatment.temp[high.outliers]<-NA
treatment.temp[low.outliers]<-NA
}
} #end while sample size too small loop
pval<-t.test(treatment.temp,control.temp)$p.value
control<-control.temp
treatment<-treatment.temp
} # end replace pvals
}# end if pval>.05
} #end drop outliers
}#end If U-test==F
measured.effect=mean(treatment, na.rm=T)-mean(control, na.rm=T)
test.result<-cbind(pval,measured.effect)
test.results<-rbind(test.results,test.result)
}
pval<-min(test.results[,1])
measured.effect<-test.results[which(test.results[,1]==min(test.results[,1]))[1],2]
sample.result<-cbind(sample.result,r, x,pval,measured.effect)
run.result<-rbind(run.result,sample.result)
##Update Live Plot
if(live.plot==T){
lines(run.result[,2],run.result[,3], type="l", lwd=3, col=rainbow(runs)[r],
xlim=c(initial.samp.size,max.samps+initial.samp.size),ylim=c(0,1)
)
if(live.plot.in.own.window==T && x==initial.samp.size && r==1){
dev.new()
plot(initial.samp.size,1, type="l", lwd=2, col="White",
xlim=c(initial.samp.size,(max.samps+samp.interval)), xlab="Sample Size (per group)",
ylim=c(0,1), ylab="P-value",
main=c("Representative P-value Traces (Significant: P < 0.05)",
paste("Outlier Cutoff(SDs) =",round(outlier.cutoff,2),
", Outlier Bias (High,Low) =", "(", round(outlier.bias.high,2), ",", round(outlier.bias.low,2), ")"),
paste("# of Outcomes=", number.of.outcomes)
),
sub=paste("Sequential Sampling=", sequential, ", U-test=", also.perform.Utest)
)
abline(h=.05, lwd=2)
lines(run.result[,2],run.result[,3], type="l", lwd=3, col=rainbow(runs)[r],
xlim=c(initial.samp.size,max.samps+initial.samp.size),ylim=c(0,1)
)
dev.set(which=dev.prev())
}
if(live.plot.in.own.window==T && x*r>initial.samp.size){
dev.set(which=dev.next())
lines(run.result[,2],run.result[,3], type="l", lwd=3, col=rainbow(runs)[r],
xlim=c(initial.samp.size,max.samps+initial.samp.size),ylim=c(0,1)
)
dev.set(which=dev.prev())
}
# End Update Live Plot
Sys.sleep(sleep.time)
}
#Generate New Samples
if(sequential==T){
##Take Sequential Samples##
for(t in 1:number.of.outcomes){
#Control
if(bimodal.treatment.only==T){
random.sample<-rep(1,samp.interval)
}else{
random.sample<-bimodal.pop[runif(samp.interval,1,1000)]
}
control<-append(experiment.results[[t]][1,],c(rnormp2(length(which(random.sample==1)),
0,control.sd1,control.shape1),
rnormp2(length(which(random.sample==0)),
0+bimodal.offset,control.sd2,control.shape2)
))
#Treatment
random.sample<-bimodal.pop[runif(samp.interval,1,1000)]
treatment<-append(experiment.results[[t]][2,],c(rnormp2(length(which(random.sample==1)),
0+effect,treatment.sd1,treatment.shape1),
rnormp2(length(which(random.sample==0)),
0+effect+bimodal.offset,treatment.sd2,treatment.shape2)
))
experiment.result<-rbind(control,treatment)
experiment.results[[t]]<-experiment.result
}#End Sequential
}else{
##Take Independant Samples (non sequential)##
experiment.results=vector("list", number.of.outcomes)
for(t in 1:number.of.outcomes){
#Control
if(bimodal.treatment.only==T){
random.sample<-rep(1,x)
}else{
random.sample<-bimodal.pop[runif(x,1,1000)]
}
control<-c(rnormp2(length(which(random.sample==1)),
0,control.sd1,control.shape1),
rnormp2(length(which(random.sample==0)),
0+bimodal.offset,control.sd2,control.shape2)
)
#Treatment
random.sample<-bimodal.pop[runif(x,1,1000)]
treatment<-c(rnormp2(length(which(random.sample==1)),
0+effect,treatment.sd1,treatment.shape1),
rnormp2(length(which(random.sample==0)),
0+effect+bimodal.offset,treatment.sd2,treatment.shape2)
)
experiment.result<-rbind(control,treatment)
experiment.results[[t]]<-experiment.result
} #End independant
} #End sampling
setTxtProgressBar(pb, x + max.samps*(r-1))
x<-x+samp.interval
} # End Run
all.results<-rbind(all.results,run.result)
if(live.plot==F && r==runs){
colnum<-1
for(k in (runs-9):runs){
xvar<-all.results[which(all.results[,1]==k),2]
yvar<-all.results[which(all.results[,1]==k),3]
lines(xvar,yvar, type="l", lwd=3, col=rainbow(10)[colnum],
xlim=c(initial.samp.size,max.samps+initial.samp.size),ylim=c(0,1)
)
colnum<-colnum+1
}
}
} # End Simulation
close(pb)
###Extract Significant Results###
rownames(all.results)<-seq(1,nrow(all.results))
sigs<-which(all.results[,3]<.05)
percent.sigs<-100*hist(all.results[sigs,2], plot=F,
breaks=seq(initial.samp.size-1, max.samps+2*samp.interval,by=(samp.interval))
)$counts/runs
###Extract First Significant Result for each Simulation Run###
first.sigs=NULL
first.sigs.filtered=NULL
for(i in 1:runs){
temp<-all.results[which(all.results[,1]==i),]
if(length(temp[which(temp[,3]<.05),2])==0){
first.sig<-Inf
}else{
first.sig<-min(temp[which(temp[,3]<.05),2])
}
first.sigs<-rbind(first.sigs,first.sig)
}
first.sigs.filtered<-first.sigs[which(first.sigs!=Inf),]
percent.atLeastOne.sig<-100*cumsum(
hist(first.sigs.filtered,
breaks=seq(initial.samp.size-1, max.samps+2*samp.interval,by=(samp.interval)),
plot=F
)$counts
)/runs
), and I might have made copy paste mistakes. 
you don't need to be a scientist to question everything, how do you think Socrates got killed?