Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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.start+adj,startle.sd, lower=0)
- test<-rtnorm(n,startle.start+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)
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement