###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))'
#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))'
#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 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
#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")