require(msm) #used for truncated normal distribution "rtnorm()"
require(viopoints) #used for viopoints plot"


# Group Names
groups=c("F1.acetophenone","F1.propanol")

# Treatment Names
treatments=c("acetophenone","propanol")

# Days experiments were run
days=c(1,2)

# Sample size (both groups)
n=12

# Mean of baseline "startle" amount (Day 1, no scent, "leader trial")
startle.baseline=10

# SD of baseline "startle" amount
startle.sd=1.5

# Habituation effect (Add to startle.baseline if day= 2)
day.effect=-5

# Scent effect (Add to startle.baseline for trials including scent)
scent.effect=2


out=NULL
# Cycle through Days
for(d in days){
  # Cycle through Groups
  for(g in groups){

    # "Counterbalanced design"
    #On day 1: Group 1 <-treatment 2;  Group 2 <-treatment 1
    #On day 2: Group 1 <-treatment 1;  Group 2 <-treatment 2
    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 day 1 add zero to baseline mean, if day 2 then add day effect
    if(d==1){adj=0}
    if(d==2){adj=day.effect}

    # Simulate "leader" results (no scent)
    leader<-rtnorm(n,startle.baseline+adj,startle.sd, lower=0)

    # Simulate "test" results (with scent)
    test<-rtnorm(n,startle.baseline+adj+scent.effect,startle.sd, lower=0)

    # Calculate "OPS.scores" using normalization strategy described in paper
    OPS.scores<-100*(test-leader)/leader

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

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


# Plot by treatment and Group
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)
}