bray_curtis.R

# This function will be used later in the script
# It generates a plot for 2dimensional MDS data that is colored by treatment
# It superimposes on that plot an ellipse which captures a proprotion of the variance of the MDS points
ggplot.NMDS<-function(XX,ZZ,COLORS){
	library(ggplot2)
	MDS1<-data.frame(scores(XX))$NMDS1
	MDS2<-data.frame(scores(XX))$NMDS2
	Aggregate.Fraction<-ZZ
	
	NMDS<-data.frame(MDS1,MDS2,Aggregate.Fraction)
	
	NMDS.mean=aggregate(NMDS[,1:2],list(group=Aggregate.Fraction),mean)
	
	veganCovEllipse<-function (cov, center = c(0, 0), scale = 1, npoints = 100) {
		theta <- (0:npoints) * 2 * pi/npoints
		Circle <- cbind(cos(theta), sin(theta))
		(center + scale * t(Circle %*% chol(cov)))
	}

	ellipses <- list()
	for( g in levels(NMDS$Aggregate.Fraction)) {
		newEllipse <- cbind(as.data.frame(t(with(NMDS[NMDS$Aggregate.Fraction==g,],
		veganCovEllipse(cov.wt(cbind(MDS1,MDS2),wt=rep(1/length(MDS1),length(MDS1)))$cov,center=c(mean(MDS1),mean(MDS2))))))
		,group=g)

		ellipses[[length(ellipses) + 1]] <- newEllipse
	}

	X1<-ggplot(data = NMDS, aes(MDS1, MDS2))
	X1 <- X1 + geom_point(aes(color = Aggregate.Fraction),size=5) +
		scale_color_manual(values=COLORS) + 
		theme(axis.text.x=element_text(size=20),axis.text.y=element_text(size=20),axis.title.x=element_text(size=20),axis.title.y=element_text(size=20)) + 
		theme(legend.title=element_text(size=15),legend.text=element_text(size=15))+theme_bw()+theme(aspect.ratio=1)
	for (ellipse in 1:length(ellipses)) {	
		X1 <- X1 + geom_path(data=ellipses[[ellipse]], aes(x=MDS1, y=MDS2, colour=group), size=0, linetype=5)
	}
	X1
}

#
#	Step 5: Check treatment effect on dissimilarity
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# 'normalized_order_counts.csv' is generated by 'normalize.R'
normalizedData <- read.csv("data/normalized_order_counts.csv", row.names=1, stringsAsFactors=T)

# 'decostand' in the package 'vegan' normalized community data (counts of orders per sample)
# The method 'total' performs this normalization be dividing each order's count by the sample total
# i.e. each order's datum is now the relative abundance of that order in the sample
library(vegan)
decostandData <- decostand(normalizedData[, -c(1:6)], "total")

# for the sake of replication, our -random- ANOVA must not actually be random.
set.seed(2)

# 'adonis' in the package 'vegan' performas a permutational analysis of variance
# In short, it shuffles community data between treatments 1000 times
#  For each shuffle, it quantifies, proportionally, how much of the variance in dissimilarity is explained by various treatment combinations
#  Finally, it gives the proportion of random shuffles for which the amount of variance explained by the treatments was equal to or greater than the amount of variance explained when the data were not shuffled
# The proprition of random shuffles with equally -effective- treatments is the p-value, the probability of obtaining our observed results under the null hypothesis that treatment doesn't affect dissimilarity
anova <- adonis(decostandData~normalizedData$Plants*normalizedData$Water*normalizedData$Season, permutations=9999)
write.csv(as.matrix(anova$aov.tab), "data/tmp_bc_anova.csv")

# 'adonis' will show with more rigor which factors affect dissimilarity
# 'mds' is, in this case, simply an effective way to visualize that effect
# Here we include every combination of treatments, to contrast an effective vs. an ineffective treatment combination
mds <- metaMDS(decostandData, autotransform=FALSE, k=3)
ggplot.NMDS(mds, normalizedData$Plants, rainbow(2))
ggplot.NMDS(mds, normalizedData$Water, rainbow(2))
ggplot.NMDS(mds, normalizedData$Season, rainbow(2))
ggplot.NMDS(mds, paste(normalizedData$Plants, normalizedData$Water), rainbow(4))
ggplot.NMDS(mds, paste(normalizedData$Plants, normalizedData$Season), rainbow(4))
ggplot.NMDS(mds, paste(normalizedData$Water, normalizedData$Season), rainbow(4))
ggplot.NMDS(mds, paste(normalizedData$Plants, normalizedData$Water, normalizedData$Season), rainbow(8))

q(save="no")