## Coagulation Data
## March 3, 2008

###
rm(list = ls())
##
coag <- read.table("http://faculty.unlv.edu/westveld/Teaching/Sta713/Data/coag.txt", header=T)
coag[1:3,]
attach(coag)

## Make the plot on pg. 69
pdf("fig.pdf",family="Times",height=3,width=6)
par(mar=c(3,3,1,1),mgp=c(1.75,.75,0))

plot(as.integer(diet),jitter(ctime),xlab="diet",xaxt="n",
     ylab="coagulation time",pch=as.character(diet),col=as.integer(diet))  ## note: use pch=as.character(diet) for the letters

tmeans<-tapply(ctime,diet,mean)  ## apply the mean function to ctime for each diet
abline(h= tmeans, col=c(1:4))  ## add the treatment means

dev.off()

## get the ANOVA table
anova(lm(ctime~diet))

## do the randomization test
Fobs<-anova(lm(ctime~diet))$F[1]
Fsim<-NULL
for(nsim in 1:1000) {
diet.sim<-sample(diet)
Fsim<-c(Fsim, anova(lm(ctime~diet.sim))$F[1] )
}

mean(Fsim>=Fobs)

1-pf(Fobs,3,20)