#read in the data
crabs=read.table("...//crab.txt",header=TRUE)
attach(crabs)
#plot all the data points
hist(count)
boxplot(count~site)

#identify observations in each group
id1=which(site==1)
id2=which(site==2)
id3=which(site==3)
id4=which(site==4)
id5=which(site==5)
id6=which(site==6)

#calculate and store sample sizes
n1=length(id1)
n2=length(id2)
n3=length(id3)
n4=length(id4)
n5=length(id5)
n6=length(id6)


#compute averages for each group
ybar1=mean(count[id1])
ybar2=mean(count[id2])
ybar3=mean(count[id3])
ybar4=mean(count[id4])
ybar5=mean(count[id5])
ybar6=mean(count[id6])

#calculate residuals
resid1=count[id1]-ybar1
resid2=count[id2]-ybar2
resid3=count[id3]-ybar3
resid4=count[id4]-ybar4
resid5=count[id5]-ybar5
resid6=count[id6]-ybar6

#concatenate all residuals
all.resid=c(resid1,resid2,resid3,resid4,resid5,resid6)

#plot boxplot of all residuals
boxplot(all.resid)
#plot boxplot of residuals by group mean
boxplot(resid1,resid2,resid3,resid4,resid5,resid6)
#residual plot (group means vs. residuals)

all.ybars=c(rep(ybar1,n1),rep(ybar2,n2),rep(ybar3,n3),rep(ybar4,n4),rep(ybar5,n5),rep(ybar6,n6))

plot(all.ybars,all.resid,pch=20)

#anova on original data
oneway.test(count~site,var.equal=TRUE)
#or, equivalently,
anova(lm(count~factor(site)))

#for multiple pairwise comparisons
TukeyHSD(aov(count~factor(site)))
#and you can plot the resulting confindence intervals
plot(TukeyHSD(aov(count~factor(site))))

#multiple comparisons using Boneferroni
pairwise.t.test(count,site,p.adjust="bonf")


#repeat this analysis for the log-transformed data
#all you need to do is to replace "count" by "log(count)"
#problem: zeroes. solution: replace count by count+1
l.count=log(count+1)
ybar1=mean(l.count[id1])
ybar2=mean(l.count[id2])
ybar3=mean(l.count[id3])
ybar4=mean(l.count[id4])
ybar5=mean(l.count[id5])
ybar6=mean(l.count[id6])

#calculate residuals
resid1=l.count[id1]-ybar1
resid2=l.count[id2]-ybar2
resid3=l.count[id3]-ybar3
resid4=l.count[id4]-ybar4
resid5=l.count[id5]-ybar5
resid6=l.count[id6]-ybar6

#concatenate all residuals
all.resid=c(resid1,resid2,resid3,resid4,resid5,resid6)

#plot boxplot of all residuals
boxplot(all.resid)
#plot boxplot of residuals by group mean
boxplot(resid1,resid2,resid3,resid4,resid5,resid6)
#residual plot (group means vs. residuals)

all.ybars=c(rep(ybar1,n1),rep(ybar2,n2),rep(ybar3,n3),rep(ybar4,n4),rep(ybar5,n5),rep(ybar6,n6))

plot(all.ybars,all.resid,pch=20)

#F test on log-transformed data
oneway.test(l.count~site,var.equal=TRUE)
#or, equivalently,
anova(lm(l.count~factor(site)))

#multiple comparisons plot for the log data
plot(TukeyHSD(aov(l.count~factor(site))))
#multiple comparisons using Boneferroni
pairwise.t.test(l.count,site,p.adjust="bonf")