# hypothesis testing and various other things you can do with
#   species composition data

# illustration uses Park Grass data set

# Mantel test
# need to construct environmental distance matrix

# use Eucl distance for continuous variables
park.DpH <- vegdist(park.env$pH,'eucl')

# could extend to multiple env var by specifying a matrix with
#   the relevant variables as the first argument

# use 0 = within, 1 = between for factors
# construct the matrix for all pairs, then convert to distance object
#  the factor variable is repeated twice, the '!=' should not be changed

temp <- outer(park.env$N,park.env$N,'!=')
park.DN <- as.dist(temp)

park.mantel.pH <- mantel(park.bray, park.DpH)
park.mantel.N <- mantel(park.bray, park.DN)

# methods only for factors
# anosim (Clarke's R)
#  can specify a distance matrix or a data matrix and distance= (default bray)
# can increase number of permutations by specifying perm = # (e.g. perm=9999)

park.anosim.N <- anosim(park.bray, park.env$N)

# has print, summary and plot methods
# plot gives boxplots of the ranks for each group

plot(park.anosim.N)

# mrpp

park.mrpp.N <- mrpp(park.bray, park.env$N)
# print this result to see details of the computation and the p-value

# associated with mrpp is the meandist() function
#  produces matrix of mean distances between all pairs of groups
temp <- meandist(park.bray, park.env$N)
temp

# ploting the meandist result gives a clustering of groups (average link)
plot(temp)

# upward links (e.g. 48) are groups that are more variable within than
#  between 

# any linear model, potentially with blocks
# adonis 
# formula interface, Y is distance object or data matrix
#  variables are 'looked up' within the specified data frame

park.adonis <- adonis(park.bray ~ pH + N, data=park.env)
# only a print method
park.adonis

# if have blocks, permute within blocks.  specify by strata = groups
# ie to permute within N groups (because no blocks in these data)

park.adonis2 <- adonis(park.bray ~ pH, data=park.env, strata=park.env$N)

# remember, these are sequential tests so order of terms matters
#   unless data are balanced.

# alternative ordination methods

#  all have a print (short output table), 
#             summary (long output table),
#    and      plot (species and site locations)

#  can use any of the 'ordiXXX' functions, e.g. ordirgl(), 
#    to plot or decorate a plot

# Correspondence analysis = reciprocal averaging
#   cca() function with only a composition matrix

park.ca <- cca(park.m)

# one of many ways to do principal components analysis

park.pca <- rda(park.m)

# add constraints by providing a formula 

# canonical correspondence analysis
park.cca <- cca(park.m ~ pH + N, data=park.env)

# redundancy analysis
park.rda <- rda(park.m ~ pH + N, data=park.env)


# procrustes rotation
#  below recreates my example: two different standardizations of the data

park.mds <- metaMDS(decostand(park.m,'total'), 'bray', autotransform=F, expand=T)
park.mdsw <- metaMDS(wisconsin(park.m), 'bray', autotransform=F, expand=T)

temp <- procrustes(park.mds, park.mdsw)

# print() gives small output, summary() longer, 
#  plot() plots superposition and arrows indicating errors


# comparison of beta diversity by examining dissimilarity within groups
park.bd <- betadisper(park.bray, park.env$N)

# has print(), anova(), plot() methods
# boxplot() plots within group dissimilarities
boxplot(park.bd)