Schefferville - transitions; phylobetadiversity

nmds ordinations

Date: February 7, 2015

R version 3.1.0


nmds calculations for four different analyses


###nmds for all species based on abundance data

#combine all plots so that they can be in a nmds
all.plots<-rbind(abd.sp.allsp, SMF.ab.allsp, TS.ab.allsp.nozeros)
dim(all.plots)
## [1] 378 125
#Plot an nmds ordination to check how plots fall out
trans.ord<-metaMDS(all.plots)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1026377 
## Run 1 stress 0.1026762 
## ... procrustes: rmse 0.002830603  max resid 0.03441154 
## Run 2 stress 0.1026941 
## ... procrustes: rmse 0.00253594  max resid 0.02328229 
## Run 3 stress 0.1026581 
## ... procrustes: rmse 0.001134633  max resid 0.01488761 
## Run 4 stress 0.1026542 
## ... procrustes: rmse 0.001154571  max resid 0.01953637 
## Run 5 stress 0.1027068 
## ... procrustes: rmse 0.002434354  max resid 0.0221711 
## Run 6 stress 0.1026517 
## ... procrustes: rmse 0.001342974  max resid 0.0158842 
## Run 7 stress 0.1026801 
## ... procrustes: rmse 0.00149698  max resid 0.02197543 
## Run 8 stress 0.1026782 
## ... procrustes: rmse 0.002346658  max resid 0.01587941 
## Run 9 stress 0.102658 
## ... procrustes: rmse 0.001221854  max resid 0.01857605 
## Run 10 stress 0.1027368 
## ... procrustes: rmse 0.002410674  max resid 0.02356667 
## Run 11 stress 0.1026844 
## ... procrustes: rmse 0.002701931  max resid 0.02913972 
## Run 12 stress 0.1026647 
## ... procrustes: rmse 0.0007436732  max resid 0.006343857 
## *** Solution reached
trans.ord
## 
## Call:
## metaMDS(comm = all.plots) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(all.plots)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1026377 
## Stress type 1, weak ties
## Two convergent solutions found after 12 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(all.plots))'
#Remove TSI9.4a outlier. This plot only has one species with an abundance of 0.25, also remove TSF6.3a, TSJ10.2a for the same reason
TS.ab.allsp.nolows<-TS.ab.allsp.nozeros[c(1:48, 50:76,78:81, 83:92),]

#again, combine all plots so that they can be in a nmds
all.plots<-rbind(abd.sp.allsp, SMF.ab.allsp, TS.ab.allsp.nolows)
dim(all.plots)
## [1] 375 125
#Plot an nmds ordination to check how plots fall out
trans.ord<-metaMDS(all.plots)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2140895 
## Run 1 stress 0.2165046 
## Run 2 stress 0.2173625 
## Run 3 stress 0.2239977 
## Run 4 stress 0.2193066 
## Run 5 stress 0.214149 
## ... procrustes: rmse 0.008717089  max resid 0.1194009 
## Run 6 stress 0.2229946 
## Run 7 stress 0.2174829 
## Run 8 stress 0.2160215 
## Run 9 stress 0.2221839 
## Run 10 stress 0.2150742 
## Run 11 stress 0.215303 
## Run 12 stress 0.2244825 
## Run 13 stress 0.2237545 
## Run 14 stress 0.219337 
## Run 15 stress 0.2155598 
## Run 16 stress 0.2159295 
## Run 17 stress 0.2241671 
## Run 18 stress 0.218478 
## Run 19 stress 0.2158725 
## Run 20 stress 0.2183839
trans.ord
## 
## Call:
## metaMDS(comm = all.plots) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(all.plots)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.2140895 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(all.plots))'

nmds based on angiosperm only data



#Set up nmds ordination without ferns and gymnosperms; only with angiosperms; all.plots have zero values removed so can continue analysis as is
names(all.plots)
##   [1] "AchMil" "AgrMer" "AlnVir" "AmeBar" "AndPol" "AnePar" "AntAlp"
##   [8] "AntMon" "ArcAlp" "ArnAng" "BarAlp" "BetGla" "BisViv" "CalCan"
##  [15] "CarAqu" "CarBig" "CarBru" "CarCap" "CarCat" "CarDef" "CarDis"
##  [22] "CarGyn" "CarLep" "CarLim" "CarMag" "CarSci" "CarTri" "CarUtr"
##  [29] "CarVag" "CasSep" "CerAlp" "ChaAng" "CopTri" "CorCan" "CysMon"
##  [36] "DasFru" "AveFle" "DiaLap" "DipCom" "DryExp" "DryInt" "ElyTra"
##  [43] "EmpNig" "EpiHor" "EquArv" "EquSci" "EquSyl" "EquVar" "EriVir"
##  [50] "EurRad" "FraVir" "GauHis" "GeoLiv" "HupApr" "JunCom" "JunTri"
##  [57] "KalPol" "LinBor" "LisCor" "LonVil" "LuzPar" "LycAno" "MaiTri"
##  [64] "MinBif" "MitNud" "MoeMac" "MonUni" "MyrGal" "OrtSec" "PacAur"
##  [71] "ParKot" "PetFri" "PhyCae" "PoaArc" "PyrAsa" "PyrGra" "RhiMin"
##  [78] "RhoGro" "RhoLap" "RibGla" "RubArc" "RubCha" "RubIda" "SalArc"
##  [85] "SalArg" "SalGla" "SalHum" "SalPed" "SalPla" "SalUva" "SalVes"
##  [92] "SchPur" "SelSel" "SolMac" "SolMul" "SteBor" "SteLon" "TofPus"
##  [99] "TriAlp" "TriBor" "TriCes" "TriSpi" "VacCes" "VacMyr" "VacOxy"
## [106] "VacUli" "VacVit" "VibEdu" "VioAdu" "VioRen" "AbiBal" "LarLar"
## [113] "PicGla" "PicMar" "AreHum" "BetMin" "BetPum" "CarBel" "CarGla"
## [120] "CopLap" "CorTri" "LuzSpi" "PoaAlp" "SalHer" "SalMyr"
angios.trans<-all.plots[,c(1:34,36:38, 41:44,49:53, 56:61, 63:92, 94:110, 115:125)]
names(angios.trans)
##   [1] "AchMil" "AgrMer" "AlnVir" "AmeBar" "AndPol" "AnePar" "AntAlp"
##   [8] "AntMon" "ArcAlp" "ArnAng" "BarAlp" "BetGla" "BisViv" "CalCan"
##  [15] "CarAqu" "CarBig" "CarBru" "CarCap" "CarCat" "CarDef" "CarDis"
##  [22] "CarGyn" "CarLep" "CarLim" "CarMag" "CarSci" "CarTri" "CarUtr"
##  [29] "CarVag" "CasSep" "CerAlp" "ChaAng" "CopTri" "CorCan" "DasFru"
##  [36] "AveFle" "DiaLap" "DryInt" "ElyTra" "EmpNig" "EpiHor" "EriVir"
##  [43] "EurRad" "FraVir" "GauHis" "GeoLiv" "JunTri" "KalPol" "LinBor"
##  [50] "LisCor" "LonVil" "LuzPar" "MaiTri" "MinBif" "MitNud" "MoeMac"
##  [57] "MonUni" "MyrGal" "OrtSec" "PacAur" "ParKot" "PetFri" "PhyCae"
##  [64] "PoaArc" "PyrAsa" "PyrGra" "RhiMin" "RhoGro" "RhoLap" "RibGla"
##  [71] "RubArc" "RubCha" "RubIda" "SalArc" "SalArg" "SalGla" "SalHum"
##  [78] "SalPed" "SalPla" "SalUva" "SalVes" "SchPur" "SolMac" "SolMul"
##  [85] "SteBor" "SteLon" "TofPus" "TriAlp" "TriBor" "TriCes" "TriSpi"
##  [92] "VacCes" "VacMyr" "VacOxy" "VacUli" "VacVit" "VibEdu" "VioAdu"
##  [99] "VioRen" "AreHum" "BetMin" "BetPum" "CarBel" "CarGla" "CopLap"
## [106] "CorTri" "LuzSpi" "PoaAlp" "SalHer" "SalMyr"
dim(angios.trans)
## [1] 375 110
#Plot an nmds ordination to check how plots fall out
angios.trans.ord<-metaMDS(angios.trans)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2120982 
## Run 1 stress 0.2153257 
## Run 2 stress 0.2198586 
## Run 3 stress 0.2137071 
## Run 4 stress 0.2129187 
## Run 5 stress 0.2217382 
## Run 6 stress 0.2130166 
## Run 7 stress 0.2146896 
## Run 8 stress 0.2194485 
## Run 9 stress 0.2131429 
## Run 10 stress 0.2152398 
## Run 11 stress 0.2144363 
## Run 12 stress 0.2205705 
## Run 13 stress 0.2194085 
## Run 14 stress 0.2148005 
## Run 15 stress 0.2119797 
## ... New best solution
## ... procrustes: rmse 0.01178555  max resid 0.1550561 
## Run 16 stress 0.2130615 
## Run 17 stress 0.2190742 
## Run 18 stress 0.2120996 
## ... procrustes: rmse 0.01190457  max resid 0.1548617 
## Run 19 stress 0.2142114 
## Run 20 stress 0.2146334
angios.trans.ord
## 
## Call:
## metaMDS(comm = angios.trans) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(angios.trans)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.2119797 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(angios.trans))'

nmds ordinations using phylogenetic betadiversity for all plot data



#first, upload maximum clade credibility phylogeny
phy<-read.nexus("trans.one.tree.nex")

x<-colnames(all.plots)
y<-phy$tip.label
setdiff(x,y)
## character(0)
setdiff(y,x)
## [1] "ComUmb" "CorAlt" "CorSto" "HydArb" "TheHum" "XimAme"
#drop unneeded tips from phylogeny
phy.trans<-drop.tip(phy, c("ComUmb", "CorAlt", "CorSto", "HydArb", "TheHum", "XimAme"))

dis.trans<-phylosor(all.plots,phy.trans)

#nmds for phySor data
phydist.nmds=metaMDS(dis.trans)
## Run 0 stress 0.3926988 
## Run 1 stress 0.3925827 
## ... New best solution
## ... procrustes: rmse 0.04286876  max resid 0.1055854 
## Run 2 stress 0.3951701 
## Run 3 stress 0.3922048 
## ... New best solution
## ... procrustes: rmse 0.04284769  max resid 0.1068832 
## Run 4 stress 0.3929441 
## Run 5 stress 0.3933658 
## Run 6 stress 0.3938603 
## Run 7 stress 0.3894258 
## ... New best solution
## ... procrustes: rmse 0.03773136  max resid 0.1097528 
## Run 8 stress 0.3934369 
## Run 9 stress 0.3908544 
## Run 10 stress 0.3955425 
## Run 11 stress 0.3932786 
## Run 12 stress 0.3953333 
## Run 13 stress 0.3897825 
## ... procrustes: rmse 0.03014595  max resid 0.1257717 
## Run 14 stress 0.3914681 
## Run 15 stress 0.3912269 
## Run 16 stress 0.3903247 
## Run 17 stress 0.3898219 
## ... procrustes: rmse 0.03249144  max resid 0.1192162 
## Run 18 stress 0.391699 
## Run 19 stress 0.3934622 
## Run 20 stress 0.3929618
phydist.nmds
## 
## Call:
## metaMDS(comm = dis.trans) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     dis.trans 
## Distance: user supplied 
## 
## Dimensions: 2 
## Stress:     0.3894258 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation 
## Species: scores missing

nmds phylsor with angiosperms only



#nmds with phylosor values among plots (angiosperms only)
#drop unneeded tips from phylogeny
phy.trans.angios<-drop.tip(phy.trans, c("CysMon", "DryExp", "EquArv", "EquSyl", "EquVar", "EquSci", "LycAno", "SelSel", "DipCom", "HupApr", "PicGla", "PicMar", "AbiBal",
"LarLar", "JunCom"))

#Calculate phylogenetic distances with PhyloSor index
dis.trans.angios<-phylosor(angios.trans,phy.trans.angios)

#nmds for phySor data
phydist.angios.nmds=metaMDS(dis.trans.angios)
## Run 0 stress 0.3974195 
## Run 1 stress 0.3969613 
## ... New best solution
## ... procrustes: rmse 0.047248  max resid 0.09608611 
## Run 2 stress 0.3957487 
## ... New best solution
## ... procrustes: rmse 0.04437489  max resid 0.1034775 
## Run 3 stress 0.3956532 
## ... New best solution
## ... procrustes: rmse 0.03919195  max resid 0.1134524 
## Run 4 stress 0.3955308 
## ... New best solution
## ... procrustes: rmse 0.03898343  max resid 0.1118756 
## Run 5 stress 0.3972252 
## Run 6 stress 0.3952598 
## ... New best solution
## ... procrustes: rmse 0.04299911  max resid 0.1067321 
## Run 7 stress 0.397012 
## Run 8 stress 0.3955803 
## ... procrustes: rmse 0.04248626  max resid 0.1078232 
## Run 9 stress 0.3964919 
## Run 10 stress 0.3959786 
## Run 11 stress 0.395442 
## ... procrustes: rmse 0.04261351  max resid 0.1085131 
## Run 12 stress 0.3955846 
## ... procrustes: rmse 0.04153832  max resid 0.1037393 
## Run 13 stress 0.398003 
## Run 14 stress 0.3952999 
## ... procrustes: rmse 0.04084545  max resid 0.1118642 
## Run 15 stress 0.3959211 
## Run 16 stress 0.3963693 
## Run 17 stress 0.3953508 
## ... procrustes: rmse 0.04135994  max resid 0.1114875 
## Run 18 stress 0.3978303 
## Run 19 stress 0.3940138 
## ... New best solution
## ... procrustes: rmse 0.04062714  max resid 0.1095988 
## Run 20 stress 0.3946411
phydist.angios.nmds
## 
## Call:
## metaMDS(comm = dis.trans.angios) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     dis.trans.angios 
## Distance: user supplied 
## 
## Dimensions: 2 
## Stress:     0.3940138 
## Stress type 1, weak ties
## No convergent solutions - best solution after 20 tries
## Scaling: centring, PC rotation 
## Species: scores missing

nmds Plots


#Plot four plots for vegetation data side-by-side
#dev.new(width=11.8, height=8)
par(mfrow=c(2,2))
plot(trans.ord, type="n", main="nmds ordination of Schefferville transition \nvegetation plots", cex.main=1, xlim=c(-2,2))
#points(trans.ord, display="sites", cex=0.2, pch=21, col="black", bg="black")
points(trans.ord$points[c(1:176),c(1:2)], cex=0.8, col="black",bg="black", pch=21)
points(trans.ord$points[c(177:286),c(1:2)], col="grey70", pch=16, cex=0.8)
points(trans.ord$points[c(287:375),c(1:2)], col="black",bg="white", pch=21, cex=0.8)
legend("bottomleft", c("Large grid", "Spruce-moss forest - fen transition", "Tunda-shrub transition"), col = c("black","grey70", "white" ), cex=0.5,
      lty = c(1, 1, 1), pch = c(NA, NA, NA), lwd=3 , bg = "gray99", bty="n")

plot(angios.trans.ord, type="n", main="nmds ordination of Schefferville transition \nvegetation plots (angiosperms only)", cex.main=1, xlim=c(-2,2))
#points(trans.ord, display="sites", cex=0.2, pch=21, col="black", bg="black")
points(angios.trans.ord$points[c(1:176),c(1:2)], cex=0.8, col="black",bg="black", pch=21)
points(angios.trans.ord$points[c(177:286),c(1:2)], col="grey70", pch=16, cex=0.8)
points(angios.trans.ord$points[c(287:375),c(1:2)], col="black",bg="white", pch=21, cex=0.8)
legend("bottomleft", c("Large grid", "Spruce-moss forest - fen transition", "Tunda-shrub transition"), col = c("black","grey70", "white" ), cex=0.5,
      lty = c(1, 1, 1), pch = c(NA, NA, NA), lwd=3 , bg = "gray99", bty="n")

      #Plot nmds ordination; color depending on which grid they are taken from
plot(phydist.nmds, type="n", main="nmds ordination of PhySor values \nfor Schefferville transition vegetation plots", cex.main=1, xlim=c(-2,2))
## Warning in ordiplot(x, choices = choices, type = type, display = display,
## : Species scores not available
#points(trans.ord, display="sites", cex=0.2, pch=21, col="black", bg="black")
points(phydist.nmds$points[c(1:176),c(1:2)], cex=0.8, col="black",bg="black", pch=21)
points(phydist.nmds$points[c(177:286),c(1:2)], col="grey70", pch=16, cex=0.8)
points(phydist.nmds$points[c(287:375),c(1:2)], col="black", bg="white", pch=21, cex=0.8)
legend("bottomleft", c("Large grid", "Spruce-moss forest - fen transition", "Tunda-shrub transition"), col = c("black","grey70", "white" ), cex=0.5,
      lty = c(1, 1, 1), pch = c(NA, NA, NA), lwd=3 , bg = "white", bty="n")

plot(phydist.angios.nmds, type="n", main="nmds ordination of PhySor values for Schefferville \ntransition vegetation plots (angiosperms only)", cex.main=1, xlim=c(-2,2))
## Warning in ordiplot(x, choices = choices, type = type, display = display,
## : Species scores not available
#points(trans.ord, display="sites", cex=0.2, pch=21, col="black", bg="black")
points(phydist.angios.nmds$points[c(1:176),c(1:2)], cex=0.8, col="black",bg="black", pch=21)
points(phydist.angios.nmds$points[c(177:286),c(1:2)], col="grey70", pch=16, cex=0.8)
points(phydist.angios.nmds$points[c(287:375),c(1:2)], col="black", bg="white", pch=21, cex=0.8)
legend("bottomleft", c("Large grid", "Spruce-moss forest - fen transition", "Tunda-shrub transition"), col = c("black","grey70", "white" ), cex=0.5,
      lty = c(1, 1, 1), pch = c(NA, NA, NA), lwd=3 , bg = "white", bty="n")