NMDS of 2011 NWCA Coastal Plains Estuarine Wetland Sites
# Site Data Start with subset of 2011 NWCA data. The dataset includes estuarine sites in the Coastal Plains ecoregion that had soil data available. It was neccessary to select only species that occured on greater than or equal to 5% of the total sites in order to reach a solution when running the NMDS (reducing species helpt ot reduce 0’s in the dataset). When possible, species that did not meet the 5% cutoff were grouped into genus groups.
# CPL_5 is species cover data for CPL estuarine sites that only includes species that occured on 5% (16) or more of the plots.
cpl5 <- read.csv("cpl_v1_soilplots_soilonly.csv")
cpl5 <- cpl5[,-c(1:3)]The function metaMDS()in the R package Vegan was used to perform NMDS on the dataset. 2 dimensions were selected with a maximum number of 1000 tries # Running NMDS
library("vegan")## Warning: package 'vegan' was built under R version 3.4.4## Loading required package: permute## Warning: package 'permute' was built under R version 3.4.4## Loading required package: lattice## This is vegan 2.5-1# Trying to get convergence with NMDS
cpl5.ord <- metaMDS(cpl5, trymax = 1000, dim = 2)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1487978
## Run 1 stress 0.1623657
## Run 2 stress 0.1705917
## Run 3 stress 0.1636854
## Run 4 stress 0.1555586
## Run 5 stress 0.1556092
## Run 6 stress 0.1601869
## Run 7 stress 0.1630541
## Run 8 stress 0.1611683
## Run 9 stress 0.1482454
## ... New best solution
## ... Procrustes: rmse 0.01390773 max resid 0.0548812
## Run 10 stress 0.1479615
## ... New best solution
## ... Procrustes: rmse 0.006493519 max resid 0.04766277
## Run 11 stress 0.1475582
## ... New best solution
## ... Procrustes: rmse 0.01105201 max resid 0.06282129
## Run 12 stress 0.1640063
## Run 13 stress 0.1475028
## ... New best solution
## ... Procrustes: rmse 0.003072261 max resid 0.02651801
## Run 14 stress 0.1563938
## Run 15 stress 0.1597209
## Run 16 stress 0.1596905
## Run 17 stress 0.1640657
## Run 18 stress 0.1475897
## ... Procrustes: rmse 0.002395905 max resid 0.02519296
## Run 19 stress 0.1479753
## ... Procrustes: rmse 0.009703945 max resid 0.06394385
## Run 20 stress 0.1546606
## Run 21 stress 0.1543823
## Run 22 stress 0.1487477
## Run 23 stress 0.1486235
## Run 24 stress 0.1487695
## Run 25 stress 0.1640017
## Run 26 stress 0.1540043
## Run 27 stress 0.1575688
## Run 28 stress 0.1506196
## Run 29 stress 0.158597
## Run 30 stress 0.1589604
## Run 31 stress 0.1475675
## ... Procrustes: rmse 0.002798101 max resid 0.02632182
## Run 32 stress 0.1491882
## Run 33 stress 0.1604627
## Run 34 stress 0.1482257
## Run 35 stress 0.1474597
## ... New best solution
## ... Procrustes: rmse 0.004866118 max resid 0.04264182
## Run 36 stress 0.1656298
## Run 37 stress 0.1487021
## Run 38 stress 0.1493223
## Run 39 stress 0.1476154
## ... Procrustes: rmse 0.005927215 max resid 0.04289392
## Run 40 stress 0.1686082
## Run 41 stress 0.1564759
## Run 42 stress 0.168702
## Run 43 stress 0.1542984
## Run 44 stress 0.1506788
## Run 45 stress 0.1667068
## Run 46 stress 0.1694386
## Run 47 stress 0.1487652
## Run 48 stress 0.147463
## ... Procrustes: rmse 0.0007057016 max resid 0.003534712
## ... Similar to previous best
## *** Solution reachedThe ordination resulted in a stress level of 0.1473803, convergent solutions were reached after 20 tries.
cpl5.ord##
## Call:
## metaMDS(comm = cpl5, trymax = 1000, dim = 2)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(cpl5))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.1474597
## Stress type 1, weak ties
## Two convergent solutions found after 48 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(cpl5))'# Plot NMDS Simple plot of the NMDS ordination. Because there are only 2 dimensions, we only have to plot NMDS1 vs NMDS2.
plot(cpl5.ord, type = "n")
points(cpl5.ord, display = "sites", cex = 0.8, pch=21, col="red", bg="yellow")
text(cpl5.ord, display = "spec", cex=0.7, col="black")
# Adding Environmental Vectors ## Soil Vectors
#read in csv
soil.env<- read.csv("soilplot_env_ARKH.csv")
soil.env <- soil.env[,-1]
#Scale dataset
soil.env <- as.data.frame(scale(soil.env))
# use envfit() to fit soil environmental variables to the ordination
soil.env.fit <- envfit(cpl5.ord, soil.env)Output of vector fit values and significance
soil.env.fit <- envfit(cpl5.ord, soil.env, permu = 999)
soil.env.fit##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## SAND_10 0.53356 -0.84576 0.3650 0.001 ***
## SILT_10 -0.28387 0.95886 0.2179 0.001 ***
## CLAY_mean_10 -0.66435 0.74742 0.3610 0.001 ***
## pH_H2O_max_10 0.34981 -0.93682 0.0876 0.004 **
## EC_mean_10 -0.97378 0.22751 0.3418 0.001 ***
## InorgC_gkg_10 -0.35313 -0.93557 0.1271 0.001 ***
## OrgC_gkg_mean_10 0.72509 0.68866 0.0231 0.159
## TOTN_gkg_10 0.26228 0.96499 0.0569 0.006 **
## TOTP_mgkg_10 0.08156 0.99667 0.0906 0.001 ***
## MEHLICH_P_mgkg_10 -0.73897 0.67374 0.0420 0.035 *
## Pox_mgkg_10 0.01804 0.99984 0.1398 0.001 ***
## S_gkg_10 -0.91788 0.39686 0.3168 0.001 ***
## CEC_mean_10 -0.14422 0.98955 0.1070 0.001 ***
## BaseSat_10 -0.98657 0.16331 0.1731 0.001 ***
## Ca_cmolkg_10 -0.84734 -0.53105 0.0471 0.020 *
## K_cmolkg_10 -0.89762 0.44077 0.4683 0.001 ***
## Mg_cmolkg_10 -0.93280 0.36039 0.2435 0.001 ***
## Na_cmolkg_10 -0.97900 0.20385 0.3518 0.001 ***
## BD_mean_10 0.63031 -0.77634 0.1125 0.001 ***
## TOT_CN_10 -0.06116 -0.99813 0.1654 0.001 ***
## TOTC_TOTP_10 0.11220 -0.99369 0.1146 0.002 **
## DPS_10 -0.29882 -0.95431 0.0550 0.009 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999Plot of all soil vectors
plot(cpl5.ord, type = "n")
points(cpl5.ord, display = "sites", cex = 0.8, pch=21, col="red", bg="yellow")
text(cpl5.ord, display = "spec", cex=0.7, col="black")
plot(soil.env.fit, col = "blue")
Select most important soil parameters
soil.select <- read.csv("soil_select.csv")
soil.select <- as.data.frame(scale(soil.select))
soil.select.fit <- envfit(cpl5.ord, soil.select)Plot
plot(cpl5.ord, type = "n")
points(cpl5.ord, display = "sites", cex = 0.8, pch=21, col="red", bg="yellow")
text(cpl5.ord, display = "spec", cex=0.7, col="black")
plot(soil.select.fit, col = "blue")
##Landscape Vectors
land.env <- read.csv("cpl_v1_soilplots_landenv.csv")
land.env <- land.env[,-1]
land.env <- as.data.frame(scale(land.env))
land.env.fit <- envfit(cpl5.ord, land.env)land.env.fit <- envfit(cpl5.ord, land.env, permu = 999)
land.env.fit##
## ***VECTORS
##
## NMDS1 NMDS2 r2 Pr(>r)
## DEVELOPED_PCT -0.05396 -0.99854 0.0726 0.003 **
## AGRICULTURE_PCT 0.61145 0.79128 0.0609 0.006 **
## FOREST_PCT 0.24139 -0.97043 0.0132 0.366
## IMPERVIOUS_PCT -0.17587 -0.98441 0.0903 0.005 **
## POPDEN -0.34271 -0.93944 0.1265 0.001 ***
## ROADDEN 0.03306 -0.99945 0.0586 0.007 **
## ELEVMAX 0.99822 0.05963 0.0944 0.004 **
## ELEVMIN 0.95181 0.30668 0.0410 0.043 *
## PCVPY_PT -0.22641 -0.97403 0.0655 0.008 **
## PIP_PT 0.92039 0.39100 0.0028 0.811
## PMAX_PT -0.26080 -0.96539 0.0950 0.001 ***
## Water_Cover -0.97287 0.23136 0.1615 0.001 ***
## SWater_Depth -0.67618 0.73674 0.1273 0.001 ***
## SP_Water_Depth -0.96034 0.27882 0.0717 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999plot(cpl5.ord, type = "n")
points(cpl5.ord, display = "sites", cex = 0.8, pch=21, col="red", bg="yellow")
text(cpl5.ord, display = "spec", cex=0.7, col="black")
plot(land.env.fit, col = "red")
Select most important Landscape Vectors
land.select <- read.csv("landscape_select.csv")
land.select <- as.data.frame(scale(land.select))
land.select.fit <- envfit(cpl5.ord, land.select)plot
plot(cpl5.ord, type = "n")
points(cpl5.ord, display = "sites", cex = 0.8, pch=21, col="red", bg="yellow")
text(cpl5.ord, display = "spec", cex=0.7, col="black")
plot(land.select.fit, col = "red")
# Combine Soil and Landscape Vectors
plot(cpl5.ord, type = "n")
points(cpl5.ord, display = "sites", cex = 0.8, pch=21, col="red", bg="yellow")
text(cpl5.ord, display = "spec", cex=0.7, col="black")
plot(soil.select.fit, col = "blue")
plot(land.select.fit, col = "red")
# Site Factors
siteinfo <- read.csv("cpl_v1_soilonly_siteinfo.csv")## States
colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.vec[siteinfo$STATE],
bg = colors.vec[siteinfo$STATE])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$STATE), bty = "n",
col = colors.vec, pch = 21, pt.bg = colors.vec, title = "STATE")
colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.vec[siteinfo$STATE],
bg = colors.vec[siteinfo$STATE])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$STATE), bty = "n",
col = colors.vec, pch = 21, pt.bg = colors.vec, title = "STATE")
plot(soil.select.fit, col = "blue")
colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.vec[siteinfo$STATE],
bg = colors.vec[siteinfo$STATE])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$STATE), bty = "n",
col = colors.vec, pch = 21, pt.bg = colors.vec, title = "STATE")
plot(land.select.fit, col = "red")
## Nonnative Stress
colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$STRESS_NONNATIVE],
bg = colors.scale[siteinfo$STRESS_NONNATIVE])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$STRESS_NONNATIVE), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "STRESS_NONNATIVE")
#plot(soil.select.fit, col = "blue")
#plot(land.select.fit, col = "red")colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$STRESS_NONNATIVE],
bg = colors.scale[siteinfo$STRESS_NONNATIVE])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$STRESS_NONNATIVE), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "STRESS_NONNATIVE")
plot(soil.select.fit, col = "blue")
#plot(land.select.fit, col = "red")colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$STRESS_NONNATIVE],
bg = colors.scale[siteinfo$STRESS_NONNATIVE])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$STRESS_NONNATIVE), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "STRESS_NONNATIVE")
#plot(soil.select.fit, col = "blue")
plot(land.select.fit, col = "red")
## Stress from Surface Hardening
colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$STRESS_HARD],
bg = colors.scale[siteinfo$STRESS_HARD])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$STRESS_HARD), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "STRESS_HARD")
#plot(soil.select.fit, col = "blue")
#plot(land.select.fit, col = "red")colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$STRESS_HARD],
bg = colors.scale[siteinfo$STRESS_HARD])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$VEGCOND), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "STRESS_HARD")
plot(soil.select.fit, col = "blue")
#plot(land.select.fit, col = "red")colors.vec <- rainbow(13)
colors.scale <- c("red", "green", "orange", 'dark red')
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$STRESS_HARD],
bg = colors.scale[siteinfo$STRESS_HARD])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$STRESS_HARD), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "STRESS_HARD")
#plot(soil.select.fit, col = "blue")
plot(land.select.fit, col = "red")
## Veg Condition
colors.vec <- rainbow(13)
colors.scale <- c("orange", "green", "red")
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$VEGCOND],
bg = colors.scale[siteinfo$VEGCOND])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$VEGCOND), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "VEGCOND")
#plot(soil.select.fit, col = "blue")
#plot(land.select.fit, col = "red")colors.vec <- rainbow(13)
colors.scale <- c("orange", "green", "red")
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$VEGCOND],
bg = colors.scale[siteinfo$VEGCOND])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$VEGCOND), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "VEGCOND")
plot(soil.select.fit, col = "blue")
#plot(land.select.fit, col = "red")colors.vec <- rainbow(13)
colors.scale <- c("orange", "green", "red")
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$VEGCOND],
bg = colors.scale[siteinfo$VEGCOND])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$VEGCOND), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "VEGCOND")
#plot(soil.select.fit, col = "blue")
plot(land.select.fit, col = "red")
## Wetland Class
colors.vec <- rainbow(13)
colors.scale <- c("green", "purple")
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$CLASS_FIELD_FWSST],
bg = colors.scale[siteinfo$CLASS_FIELD_FWSST])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$CLASS_FIELD_FWSST), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "CLASS_FIELD_FWSST")
#plot(soil.select.fit, col = "blue")
#plot(land.select.fit, col = "red")colors.vec <- rainbow(13)
colors.scale <- c("green", "purple")
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$CLASS_FIELD_FWSST],
bg = colors.scale[siteinfo$CLASS_FIELD_FWSST])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$CLASS_FIELD_FWSST), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "CLASS_FIELD_FWSST")
plot(soil.select.fit, col = "blue")
#plot(land.select.fit, col = "red")colors.vec <- rainbow(13)
colors.scale <- c("green", "purple")
plot(cpl5.ord, display = "sites", type = "n")
points(cpl5.ord, display = "sites", cex=0.7, pch = 21,
col = colors.scale[siteinfo$CLASS_FIELD_FWSST],
bg = colors.scale[siteinfo$CLASS_FIELD_FWSST])
text(cpl5.ord, display = "spec", cex=0.7, col="black")
legend("topright", legend = levels(siteinfo$CLASS_FIELD_FWSST), bty = "n",
col = colors.scale, pch = 21, pt.bg = colors.scale, title = "CLASS_FIELD_FWSST")
#plot(soil.select.fit, col = "blue")
plot(land.select.fit, col = "red")