pdf("C:\\Rmisc\\graphs\\kuehl_04_01.pdf")

site <- rep(1:6,each=25)
count <- c(0,0,22,3,
17,0,0,7,11,11,73,33,0,65,13,44,20,27,48,104,233,81,22,9,2,
415,466,6,14,12,0,3,1,16,55,142,10,2,145,6,4,5,124,24,204,0,0,56,0,8,
0,0,4,13,5,1,1,4,4,36,407,0,0,18,4,14,0,24,52,314,245,107,5,6,2,
0,0,0,4,2,2,5,4,2,1,0,12,1,30,0,3,28,2,21,8,82,12,10,2,0,
0,1,1,2,2,1,2,29,2,2,0,13,0,19,1,3,26,30,5,4,94,1,9,3,0,
0,0,0,2,3,0,0,4,0,5,4,22,0,64,4,4,43,3,16,19,95,6,22,0,0)

crab <- data.frame(site,count)

site <- factor(site)

tapply(count,site,mean)
tapply(count,site,md)
tapply(count,site,min)
tapply(count,site,max)

crab.mod <- lm(count~site)
summary.lm(crab.mod)
anova(crab.mod)

ehat <- residuals(crab.mod)
ybar <- predict(crab.mod)

qqnorm(ehat)
qqline(ehat)

plot(ybar,ehat,main="Residuals vs Site Mean")

ehat1 <- ehat[site==1]
ehat2 <- ehat[site==2]
ehat3 <- ehat[site==3]
ehat4 <- ehat[site==4]
ehat5 <- ehat[site==5]
ehat6 <- ehat[site==6]

levene1 <- abs(ehat1-median(ehat1))
levene2 <- abs(ehat2-median(ehat2))
levene3 <- abs(ehat3-median(ehat3))
levene4 <- abs(ehat4-median(ehat4))
levene5 <- abs(ehat5-median(ehat5))
levene6 <- abs(ehat6-median(ehat6))

levene <- c(levene1,levene2,levene3,levene4,levene5,levene6)
levene

crab.lev <- lm(levene~site)

anova(crab.lev)

dev.off()