head(Cone11.12.1)
names(Cone11.12.1)<-c("Site", "Habitat", "Species")
boxplot(Species~Site, data=Cone11.12.1, main="Species Richness by Site", ylab="Number of Species")
barplot(Species~Site, data=Cone11.12.1)

library(nlme)
model2=lme(Species~Site, random=~1|Habitat, data=Cone11.12.1, method="REML")
anova.lme(model2, type="sequential", adjustSigma = FALSE)

library(multcomp)
posthoc =glht(model2, linfct=mcp(Site="Tukey"))
mcs2=summary(posthoc, test=adjusted("single-step"))
mcs2
cld(mcs2, level=0.05, decreasing=TRUE)

library(multcompView)
library(lsmeans)
lsm2 = lsmeans(model2, pairwise ~ Habitat, adjust="tukey")
lsm2
cld(lsm2, alpha=0.05, Letters=letters, adjust="tukey")

model.fixed2 = gls (Species~Habitat, data=Cone11.12.1,method="REML")
anova(model2, model.fixed2)
head(Cone11.12.2)
names(Cone11.12.2)<-c("Site","Species")
cats<-barplot(Cone11.12.2$Species, names.arg=c("Site 1", "Site 2", "Site 3", "Site 4", "Site 5"), ylim = c(0, 8), ylab="Number of Species", cex.lab="2", cex.axis = "2", cex.names="2")
error.bar<- function(plot, means, upper, lower=upper, length=0.1) {
   arrows(plot,means+upper, plot, means-lower, angle=90, code=3, length=length)}
error.bar(cats, Cone11.12.2$Species, c(.25, 2.25, NA, .25, .25))