# PCoA and nMDS from a distance matrix

library(vegan)
library(MASS)

# linear scaling (PCoA), k specifies number of dimensions
#  cmdscale is part of vanilla R, returns nsite x k matrix with coordinates
park.pcoa <- cmdscale(park.bray, k=2)

plot(park.pcoa)
# labeled plot, for non-vegan functions
plot(park.pcoa, type='n', xlab='MDS axis 1', ylab='MDS axis 2')
text(park.pcoa, dimnames(park.pcoa)[[1]])


# nMDS 'by hand', uses isoMDS from MASS library, 
#  returns list with $points = coordinates and $stress = stress value
# generate a random starting configuration of points 
start <- initMDS(park.bray)
park.mds <- isoMDS(park.bray, start)
plot(park.mds$points)


plot(park.mds$points, type='n')
text(park.mds$points, dimnames(park.mds$points)[[1]] )


# semi-automatic 'black box' function
#  metaMDS in vegan package
#  does various things to try to generate a useful analysis without
#   thought.  
#  Always good things include 20 random starts,
#  almost always good things include step-across distances (W lecture)
#  not so always good things include automatic transformation of data
#   and expanding species scores.

# can give it a data matrix and the name of a distance measure
#  note vegdist uses method=, but metaMDS uses distance=
park.mds2 <- metaMDS(park.m,distance='bray',autotransform=F, expand=F)
plot(park.mds2)
plot(park.mds2,type='t')
plot(park.mds2, disp='sites', type='t')  # or species

# or a distance matrix, but then does not compute species scores
park.mds3 <- metaMDS(park.bray, expand=F)
plot(park.mds3)

# finding solutions with more (or fewer) dimensions
park.mds1 <- metaMDS(park.m,k=1,distance='bray',autotransform=F)
park.mds3 <- metaMDS(park.m,k=3,distance='bray',autotransform=F)

# plotting stress vs number of dimensions
plot(1:3,c(park.mds1$stress,park.mds2$stress,park.mds3$stress), type='b')
plot(1:3,
     c(park.mds1$stress,park.mds2$stress,park.mds3$stress),
     type='b', ylim=c(0,park.mds1$stress), xlab='# dimensions', 
     ylab='Stress')

# interactive 3D plot, requires rgl package
ordirgl(park.mds3, type='t')