require(msm)
require(viopoints)

groups=c("F1.acetophenone","F1.propanol")
treatments=c("acetophenone","propanol")
days=c(1,2)
n=12
startle.baseline=10
startle.sd=1.5
day.effect=-5
scent.effect=2


out=NULL
for(d in days){
for(g in groups){

if(d==1 & g==groups[1]){treatment=treatments[2]}
if(d==1 & g==groups[2]){treatment=treatments[1]}
if(d==2 & g==groups[1]){treatment=treatments[1]}
if(d==2 & g==groups[2]){treatment=treatments[2]}

if(d==1){adj=0}
if(d==2){adj=day.effect}

leader<-rtnorm(n,startle.baseline+adj,startle.sd, lower=0)
test<-rtnorm(n,startle.baseline+adj+scent.effect,startle.sd, lower=0)

OPS.scores<-100*(test-leader)/leader

out<-rbind(out,cbind(d,treatment,g,OPS.scores))
}
}

out=as.data.frame(out, stringsAsFactors = F)
out[,4]<-as.numeric(out[,4])


par(mfrow=c(1,2))
for(j in 1:length(treatments)){
sub<-out[which(out[,2]==treatments[j]),]
viopoints(sub[,4]~sub[,3], pch=c(15,16), col="Black",
ylim=c(min(out[,4]),max(out[,4])),
ylab=paste("Percentage OPS to", treatments[j]))
boxplot(sub[,4]~sub[,3], add=T, outline=F)
abline(h=0, lty=3)
}