QK Box 16.2

Worked example of NMDS, ANOSIM, PERMANOVA, MV-ABUND, BIO-ENV: invertebrates in artificial ponds

We continue with the pond invertebrate example from Box 16.1. We will now use these data to produce an NMDS ordination plot showing the relationship among ponds based on invertebrate family abundances. We will compare this analysis to one where the abundances are converted to proportions of pond totals and where abundances are converted to presence-absence data. We will also compare the four management groups using ANOSIM, PERMANOVA, and MV-ABUND and examine relationships between the patterns among ponds based on invertebrate families and patterns based on environmental variables.

Preliminaries

Load graphics tweaks and core set of packages

Load additional packages: vegan, mvabund

Load lemminvert data and remove labels

lemminvert <- read_csv("../data/lemminvert.csv")

lemminvert1 <- lemminvert[,-(1:3)]

MDS on original invert abundances

lemminvert.bc <- vegdist(lemminvert1,'bray')    #Create Bray-Curtis matrix
lemminvert.mds <- metaMDS(lemminvert.bc,k=2,autotransform=FALSE)

Run 0 stress 0.1156702016 
Run 1 stress 0.1153228917 
... New best solution
... Procrustes: rmse 0.01897880188  max resid 0.07580301109 
Run 2 stress 0.1625649627 
Run 3 stress 0.1156700351 
... Procrustes: rmse 0.01894864552  max resid 0.07543301313 
Run 4 stress 0.1153230838 
... Procrustes: rmse 0.0001747205327  max resid 0.0006785231796 
... Similar to previous best
Run 5 stress 0.115670024 
... Procrustes: rmse 0.01894433962  max resid 0.07540411245 
Run 6 stress 0.1172330644 
Run 7 stress 0.1156700384 
... Procrustes: rmse 0.01893607684  max resid 0.07532678562 
Run 8 stress 0.117232991 
Run 9 stress 0.1156700522 
... Procrustes: rmse 0.01895515353  max resid 0.07548437274 
Run 10 stress 0.1156700458 
... Procrustes: rmse 0.01893148736  max resid 0.07527102811 
Run 11 stress 0.1172329882 
Run 12 stress 0.1156700954 
... Procrustes: rmse 0.01891740638  max resid 0.07517139284 
Run 13 stress 0.1156700384 
... Procrustes: rmse 0.01895070766  max resid 0.0754500985 
Run 14 stress 0.1153228735 
... New best solution
... Procrustes: rmse 0.0001002065992  max resid 0.0003940480001 
... Similar to previous best
Run 15 stress 0.1172331314 
Run 16 stress 0.1153230881 
... Procrustes: rmse 0.0002666670402  max resid 0.001041512682 
... Similar to previous best
Run 17 stress 0.1153228836 
... Procrustes: rmse 7.996045759e-05  max resid 0.000312427442 
... Similar to previous best
Run 18 stress 0.3826070117 
Run 19 stress 0.1156700262 
... Procrustes: rmse 0.01892727566  max resid 0.07538609972 
Run 20 stress 0.1172331017 
*** Best solution repeated 3 times

stressplot(lemminvert.mds, main="Shepard plot")

lemminvert.mds

Call:
metaMDS(comm = lemminvert.bc, k = 2, autotransform = FALSE) 

global Multidimensional Scaling using monoMDS

Data:     lemminvert.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.1153228735 
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 14 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

plot(lemminvert.mds$points,type="n")
text(lemminvert.mds$points,lab=lemminvert$manag,col="black")

points(lemminvert.mds$points,pch=16,col=as.numeric(lemminvert$managsymb))

Standardize abundance to proportional abundance for each pond

lemminvert1s <- decostand(lemminvert1, method = "total")
lemminvert1s.bc <- vegdist(lemminvert1s,'bray')
lemminvert1s.mds <- metaMDS(lemminvert1s.bc,k=2,try=20,trymax=40,maxit=200,autotransform=FALSE)

Run 0 stress 0.1649421061 
Run 1 stress 0.1968737718 
Run 2 stress 0.1878371053 
Run 3 stress 0.1790320285 
Run 4 stress 0.2033229782 
Run 5 stress 0.1787485915 
Run 6 stress 0.1790320285 
Run 7 stress 0.1649065586 
... New best solution
... Procrustes: rmse 0.0149550513  max resid 0.05093557104 
Run 8 stress 0.2169582729 
Run 9 stress 0.1649421055 
... Procrustes: rmse 0.01495701528  max resid 0.05108815652 
Run 10 stress 0.1648801132 
... New best solution
... Procrustes: rmse 0.00829849936  max resid 0.03137433076 
Run 11 stress 0.170948786 
Run 12 stress 0.1861307271 
Run 13 stress 0.1820631099 
Run 14 stress 0.1649065583 
... Procrustes: rmse 0.008276015597  max resid 0.03133371479 
Run 15 stress 0.1648801129 
... New best solution
... Procrustes: rmse 2.258851733e-05  max resid 6.742058596e-05 
... Similar to previous best
Run 16 stress 0.2042870927 
Run 17 stress 0.1790250489 
Run 18 stress 0.1707635842 
Run 19 stress 0.1894802214 
Run 20 stress 0.1856856478 
*** Best solution repeated 1 times

stressplot(lemminvert1s.mds, main="Shepard plot")

lemminvert1s.mds

Call:
metaMDS(comm = lemminvert1s.bc, k = 2, try = 20, trymax = 40,      autotransform = FALSE, maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     lemminvert1s.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.1648801129 
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 15 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

plot(lemminvert1s.mds$points,type="n")
text(lemminvert1s.mds$points,lab=lemminvert$manag,col="black")

points(lemminvert1s.mds$points,pch=16,col=as.numeric(lemminvert$managsymb))

MDS on presence-absence using Jaccard transformed data

lemminvert.jac <- vegdist(lemminvert1,binary=TRUE,method='jaccard')
lemminvert.jac.mds <- metaMDS(lemminvert.jac,k=2,autotransform=FALSE)

Run 0 stress 0.1676824162 
Run 1 stress 0.1676823749 
... New best solution
... Procrustes: rmse 0.0003612348514  max resid 0.001233603919 
... Similar to previous best
Run 2 stress 0.1775569077 
Run 3 stress 0.1978442817 
Run 4 stress 0.167682573 
... Procrustes: rmse 0.000499013239  max resid 0.001703304889 
... Similar to previous best
Run 5 stress 0.1676825406 
... Procrustes: rmse 0.0004497567292  max resid 0.001533184971 
... Similar to previous best
Run 6 stress 0.1776146578 
Run 7 stress 0.200584873 
Run 8 stress 0.2159557664 
Run 9 stress 0.167682381 
... Procrustes: rmse 1.962723544e-05  max resid 5.857404159e-05 
... Similar to previous best
Run 10 stress 0.1818470917 
Run 11 stress 0.1676825091 
... Procrustes: rmse 0.0004518002795  max resid 0.001531750207 
... Similar to previous best
Run 12 stress 0.2122791599 
Run 13 stress 0.1978441145 
Run 14 stress 0.1999713348 
Run 15 stress 0.1775570557 
Run 16 stress 0.1676824875 
... Procrustes: rmse 0.0001202097147  max resid 0.0003915126096 
... Similar to previous best
Run 17 stress 0.1938968838 
Run 18 stress 0.1676823497 
... New best solution
... Procrustes: rmse 4.464386131e-05  max resid 0.0001540719556 
... Similar to previous best
Run 19 stress 0.1676825948 
... Procrustes: rmse 0.0004712925821  max resid 0.001598077215 
... Similar to previous best
Run 20 stress 0.1676824373 
... Procrustes: rmse 0.000110061048  max resid 0.0003686477231 
... Similar to previous best
*** Best solution repeated 3 times

stressplot(lemminvert.jac.mds, main="Shepard plot")

lemminvert.jac.mds

Call:
metaMDS(comm = lemminvert.jac, k = 2, autotransform = FALSE) 

global Multidimensional Scaling using monoMDS

Data:     lemminvert.jac 
Distance: binary jaccard 

Dimensions: 2 
Stress:     0.1676823497 
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 18 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

plot(lemminvert.jac.mds$points,type="n")
text(lemminvert.jac.mds$points,lab=lemminvert$manag,col="black")

points(lemminvert.jac.mds$points,pch=16,col=as.numeric(lemminvert$managsymb))

Plot env variables on ordinations

Select 11 variables (same as ch15)

lemmenv <- read_csv("../data/lemmenv.csv")

Rows: 23 Columns: 17── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): site, manag
dbl (15): depth, silt, ph, o2, temp, surface, trans, cond, chla, tn, tp, susps, subm, reed, emerg
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

lemmenv1 <- lemmenv[,c(3:5,7:8,10:13,16:17)]
# standardize predictors
lemmenv1$lchla <- log10(lemmenv1$chla)
lemmenv1$depth <- scale(lemmenv1$depth)
lemmenv1$silt <- scale(lemmenv1$silt)
lemmenv1$ph <- scale(lemmenv1$ph)
lemmenv1$temp <- scale(lemmenv1$temp)
lemmenv1$surface <- scale(lemmenv1$surface)
lemmenv1$cond <- scale(lemmenv1$cond)
lemmenv1$lchla <- scale(lemmenv1$lchla)
lemmenv1$tn <- scale(lemmenv1$tn)
lemmenv1$tp <- scale(lemmenv1$tp)
lemmenv1$reed <- scale(lemmenv1$reed)
lemmenv1$emerg <- scale(lemmenv1$emerg)
# plot using original abundances
fit1 <- envfit(lemminvert.mds, lemmenv1)
fit1

***VECTORS

              NMDS1       NMDS2      r2 Pr(>r)  
depth    0.04741568 -0.99887524 0.17443  0.165  
silt     0.86935396  0.49418993 0.12897  0.253  
ph       0.29031669 -0.95693062 0.12645  0.254  
temp    -0.16690777 -0.98597251 0.00161  0.989  
surface  0.01868586 -0.99982540 0.21044  0.093 .
cond    -0.19972348 -0.97985230 0.08807  0.395  
chla     0.59035813 -0.80714143 0.08145  0.441  
tn      -0.49295639 -0.87005402 0.19845  0.123  
tp      -0.36943509 -0.92925654 0.34857  0.021 *
reed    -0.67318080  0.73947793 0.21728  0.074 .
emerg   -0.62429867  0.78118574 0.12491  0.247  
lchla    0.15594189 -0.98776623 0.15798  0.185  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999

plot(fit1)

Error in strwidth(labels, ...) : plot.new has not been called yet

# plot using proportional abundances
fit2 <- envfit(lemminvert1s.mds, lemmenv1)
fit2

***VECTORS

              NMDS1       NMDS2      r2 Pr(>r)  
depth    0.97880991 -0.20477100 0.15296  0.215  
silt    -0.10984606 -0.99394861 0.02816  0.755  
ph       0.94740555  0.32003550 0.08464  0.365  
temp     0.20178633 -0.97942957 0.02851  0.757  
surface  0.86028323 -0.50981641 0.13431  0.235  
cond     0.46289026 -0.88641559 0.08084  0.420  
chla     0.64367295  0.76530068 0.02808  0.738  
tn       0.33297260  0.94293650 0.15821  0.170  
tp       0.57503048  0.81813198 0.14516  0.217  
reed    -0.70398668  0.71021318 0.26945  0.047 *
emerg   -0.79831881  0.60223507 0.15104  0.217  
lchla    0.98631286  0.16488465 0.11665  0.293  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999

plot(fit2)

Error in strwidth(labels, ...) : plot.new has not been called yet

# plot using pa
fit3 <- envfit(lemminvert.jac.mds, lemmenv1)
fit3

***VECTORS

              NMDS1       NMDS2      r2 Pr(>r)
depth   -0.67142734 -0.74107039 0.02419  0.793
silt    -0.81534590 -0.57897415 0.00657  0.924
ph       0.05623164 -0.99841775 0.04591  0.620
temp     0.56637607 -0.82414692 0.03547  0.702
surface -0.01819445 -0.99983447 0.06929  0.498
cond    -0.07848917 -0.99691497 0.05734  0.567
chla     0.11334636 -0.99355554 0.01579  0.832
tn       0.10793391 -0.99415807 0.04332  0.645
tp       0.23971575 -0.97084312 0.10234  0.333
reed     0.36892078 -0.92946084 0.06170  0.550
emerg    0.89227883  0.45148475 0.19596  0.110
lchla    0.43900893 -0.89848270 0.01533  0.849
Permutation: free
Number of permutations: 999

plot(fit3)

Error in strwidth(labels, ...) : plot.new has not been called yet

Model/test management types

Check homogeneity of dispersions on original abundances; if OK, for comparison run anosim, permanova and mvabund

lemminvert.disp <- betadisper(lemminvert.bc,lemminvert$manag)
anova(lemminvert.disp)

Analysis of Variance Table

Response: Distances
          Df     Sum Sq     Mean Sq F value Pr(>F)
Groups     3 0.05364661 0.017882203 1.06078  0.389
Residuals 19 0.32029544 0.016857655               

# anosim on original abundances
lemminvert.ano <- anosim(lemminvert.bc, lemminvert$manag,permutations=999)
summary(lemminvert.ano)

Call:
anosim(x = lemminvert.bc, grouping = lemminvert$manag, permutations = 999) 
Dissimilarity: bray 

ANOSIM statistic R: 0.1729514 
      Significance: 0.017 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0953 0.1275 0.1502 0.1828 

Dissimilarity ranks between and within classes:
        0%   25%   50%    75% 100%   N
Between  1 68.00 133.0 200.00  253 197
li      10 15.75  43.5 105.50  181  10
nf      34 93.00 127.0 151.75  193  10
nm      25 64.00 132.0 177.00  236  21
yf       6 67.50  85.0 163.00  195  15

# permanova on original abundances
lemminvert.ado <- adonis2(lemminvert.bc~manag,data=lemminvert,permutations=999)
print(lemminvert.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = lemminvert.bc ~ manag, data = lemminvert, permutations = 999)
         Df  SumOfSqs         R2       F Pr(>F)   
manag     3 1.3565598 0.25111301 2.12366  0.007 **
Residual 19 4.0456287 0.74888699                  
Total    22 5.4021885 1.00000000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# mv-abund on original abundances
lemminvertmv <- mvabund(lemminvert1)
lemminvertmv.mv <- manyglm(lemminvertmv~lemmenv$manag,family="negative_binomial")
plot(lemminvertmv.mv)
anova(lemminvertmv.mv)

Repeat process for proportional abundance data

# check homog of dispersions on prop abundances
lemminvert1s.disp <- betadisper(lemminvert1s.bc,lemminvert$manag)
anova(lemminvert1s.disp)

Analysis of Variance Table

Response: Distances
          Df     Sum Sq      Mean Sq F value  Pr(>F)
Groups     3 0.02873894 0.0095796473 0.51922 0.67413
Residuals 19 0.35055280 0.0184501474                

# anosim on proportional abundances
lemminvert1s.ano <- anosim(lemminvert1s.bc, lemminvert$manag,permutations=999)
summary(lemminvert1s.ano)

Call:
anosim(x = lemminvert1s.bc, grouping = lemminvert$manag, permutations = 999) 
Dissimilarity: bray 

ANOSIM statistic R: 0.117839 
      Significance: 0.062 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0906 0.1278 0.1474 0.1847 

Dissimilarity ranks between and within classes:
        0%   25%   50%    75% 100%   N
Between  1  63.0 132.0 200.00  253 197
li      15  46.0  65.5  81.75  159  10
nf      42  88.0 102.5 186.75  191  10
nm       7  83.0 149.0 173.00  240  21
yf       2 104.5 123.0 161.00  236  15

# permanova on proportional abundances
lemminvert1s.ado <- adonis2(lemminvert1s.bc~manag,data=lemminvert,permutations=999)
print(lemminvert1s.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = lemminvert1s.bc ~ manag, data = lemminvert, permutations = 999)
         Df  SumOfSqs         R2       F Pr(>F)  
manag     3 0.8296405 0.23689983 1.96614  0.035 *
Residual 19 2.6724326 0.76310017                 
Total    22 3.5020732 1.00000000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mv-abund on proportional abundances using composition=TRUE (this either doesnt work or is very slow!!!) Note: when composition=TRUE, gives error message; runs with composition=FALSE, but that isn’t relative abundance

lemminvert1smv.mv <- manyglm(lemminvertmv~lemmenv$manag,composition=TRUE,family="negative_binomial")

Error in manyglm(lemminvertmv ~ lemmenv$manag, composition = TRUE, family = "negative_binomial") : 
  object '"negative_binomial"' not found

plot(lemminvert1smv.mv)

anova(lemminvert1smv.mv)

Time elapsed: 0 hr 0 min 1 sec
Analysis of Deviance Table

Model: lemminvert1smv ~ lemmenv$manag

Multivariate test:
              Res.Df Df.diff    Dev Pr(>Dev)    
(Intercept)       22                            
lemmenv$manag     19       3 8.3276    0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Arguments:
 Test statistics calculated assuming uncorrelated response (for faster computation) 
 P-value calculated using 999 iterations via PIT-trap resampling.

try mv-abund on standardized abunds from earlier

lemminvert1smv <- mvabund(lemminvert1s)
lemminvert1smv.mv <- manyglm(lemminvert1smv~lemmenv$manag,family="negative_binomial")

Warning: Non-integer data are fitted to the negative.binomial model.

plot(lemminvert1smv.mv)

anova(lemminvert1smv.mv)

Time elapsed: 0 hr 0 min 1 sec
Analysis of Deviance Table

Model: lemminvert1smv ~ lemmenv$manag

Multivariate test:
              Res.Df Df.diff    Dev Pr(>Dev)    
(Intercept)       22                            
lemmenv$manag     19       3 8.3276    0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Arguments:
 Test statistics calculated assuming uncorrelated response (for faster computation) 
 P-value calculated using 999 iterations via PIT-trap resampling.

Finally, run analysis on presence-absence data

# check homog of dispersions on pa
lemminvert.jac.disp <- betadisper(lemminvert.jac,lemminvert$manag)
anova(lemminvert.jac.disp)

Analysis of Variance Table

Response: Distances
          Df      Sum Sq      Mean Sq F value  Pr(>F)
Groups     3 0.058952285 0.0196507616  2.1726 0.12472
Residuals 19 0.171851500 0.0090448158                

# anosim on pa
lemminvert.jac.ano <- anosim(lemminvert.jac, lemminvert$manag,permutations=999)
summary(lemminvert.jac.ano)

Call:
anosim(x = lemminvert.jac, grouping = lemminvert$manag, permutations = 999) 
Dissimilarity: binary jaccard 

ANOSIM statistic R: 0.1024293 
      Significance: 0.074 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0905 0.1205 0.1445 0.1833 

Dissimilarity ranks between and within classes:
          0%     25%    50%     75%  100%   N
Between  1.0  72.000 131.50 186.500 253.0 197
li      21.5  53.250  72.75 113.125 204.5  10
nf       5.5  34.625  66.00 100.500 219.5  10
nm      41.0 117.000 195.50 230.500 249.0  21
yf      21.5  41.000  72.00 128.250 192.5  15

# permanova on pa
lemminvert.jac.ado <- adonis2(lemminvert.jac~manag,data=lemminvert,permutations=999)
print(lemminvert.jac.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = lemminvert.jac ~ manag, data = lemminvert, permutations = 999)
         Df   SumOfSqs         R2       F Pr(>F)   
manag     3 0.54394094 0.21570186 1.74183  0.008 **
Residual 19 1.97778484 0.78429814                  
Total    22 2.52172579 1.00000000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mv-abund on pa

lemminvert1pa <- decostand(lemminvert1, method = "pa")
lemminvert1pamv <- mvabund(lemminvert1pa)
lemminvert1pamv.mv <- manyglm(lemminvert1pamv~lemmenv$manag,family="binomial")
plot(lemminvert1pamv.mv)

anova(lemminvert1pamv.mv)

Time elapsed: 0 hr 0 min 1 sec
Analysis of Deviance Table

Model: lemminvert1pamv ~ lemmenv$manag

Multivariate test:
              Res.Df Df.diff     Dev Pr(>Dev)  
(Intercept)       22                           
lemmenv$manag     19       3 127.529    0.034 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Arguments:
 Test statistics calculated assuming uncorrelated response (for faster computation) 
 P-value calculated using 999 iterations via PIT-trap resampling.

Use bio-env to explore relationships with environmental variables

# bio-env raw data
bioenv(lemminvert.bc, lemmenv1, method="spearman", metric="euclidean")

4095 possible subsets (this may take time...)

Call:
bioenv(comm = lemminvert.bc, env = lemmenv1, method = "spearman",      metric = "euclidean") 

Subset of environmental variables with best correlation to community data.

Correlations:    spearman 
Dissimilarities: bray 
Metric:          euclidean 

Best model has 3 parameters (max. 12 allowed):
depth tp emerg
with correlation  0.1053499931 

# bio-env prop data
bioenv(lemminvert1s.bc, lemmenv1,method="spearman", metric="euclidean")

4095 possible subsets (this may take time...)

Call:
bioenv(comm = lemminvert1s.bc, env = lemmenv1, method = "spearman",      metric = "euclidean") 

Subset of environmental variables with best correlation to community data.

Correlations:    spearman 
Dissimilarities: bray 
Metric:          euclidean 

Best model has 3 parameters (max. 12 allowed):
depth emerg lchla
with correlation  0.1969941504 

# bio-env on pa data
bioenv(lemminvert.jac, lemmenv1,method="spearman", metric="euclidean")

4095 possible subsets (this may take time...)

Call:
bioenv(comm = lemminvert.jac, env = lemmenv1, method = "spearman",      metric = "euclidean") 

Subset of environmental variables with best correlation to community data.

Correlations:    spearman 
Dissimilarities: binary jaccard 
Metric:          euclidean 

Best model has 4 parameters (max. 12 allowed):
surface tp reed emerg
with correlation  0.0649229655 

# mantel test
lemmenv.euc <- vegdist(lemmenv1,"euclidean")
mantel(lemminvert.bc,lemmenv.euc,method="spearman")

Mantel statistic based on Spearman's rank correlation rho 

Call:
mantel(xdis = lemminvert.bc, ydis = lemmenv.euc, method = "spearman") 

Mantel statistic r: -0.04405329 
      Significance: 0.562 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.169 0.234 0.278 0.328 
Permutation: free
Number of permutations: 999

mantel(lemminvert.bc,lemmenv.euc,method="pearson")

Mantel statistic based on Pearson's product-moment correlation 

Call:
mantel(xdis = lemminvert.bc, ydis = lemmenv.euc, method = "pearson") 

Mantel statistic r: -0.01516855 
      Significance: 0.49 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.212 0.274 0.324 0.392 
Permutation: free
Number of permutations: 999

model/test continuous variables with permanova and mv-abund

# raw abundances
lemminvert2.ado <- adonis2(lemminvert.bc~depth+silt+ph+temp+surface+cond+lchla+tn+tp+reed+emerg,data=lemmenv1,permutations=999)
print(lemminvert2.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = lemminvert.bc ~ depth + silt + ph + temp + surface + cond + lchla + tn + tp + reed + emerg, data = lemmenv1, permutations = 999)
         Df  SumOfSqs         R2       F Pr(>F)  
depth     1 0.3151774 0.05834254 1.46650  0.193  
silt      1 0.3977300 0.07362387 1.85061  0.071 .
ph        1 0.2925690 0.05415750 1.36130  0.192  
temp      1 0.2549017 0.04718489 1.18604  0.296  
surface   1 0.2704965 0.05007165 1.25860  0.277  
cond      1 0.1157972 0.02143524 0.53880  0.834  
lchla     1 0.2177099 0.04030031 1.01299  0.429  
tn        1 0.2588641 0.04791838 1.20448  0.298  
tp        1 0.2439905 0.04516512 1.13527  0.323  
reed      1 0.2886512 0.05343227 1.34307  0.230  
emerg     1 0.3821963 0.07074842 1.77833  0.090 .
Residual 11 2.3641047 0.43761981                 
Total    22 5.4021885 1.00000000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lemminvert3.ado <- adonis2(lemminvert.bc~depth+silt+ph+temp+surface+cond+lchla+tn+tp+reed+emerg,data=lemmenv1,by="margin",permutations=999)
print(lemminvert3.ado)

Permutation test for adonis under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 999

adonis2(formula = lemminvert.bc ~ depth + silt + ph + temp + surface + cond + lchla + tn + tp + reed + emerg, data = lemmenv1, permutations = 999, by = "margin")
         Df  SumOfSqs         R2       F Pr(>F)  
depth     1 0.3294949 0.06099286 1.53311  0.147  
silt      1 0.1810559 0.03351528 0.84244  0.569  
ph        1 0.2562631 0.04743690 1.19237  0.274  
temp      1 0.2256969 0.04177879 1.05015  0.366  
surface   1 0.2482610 0.04595563 1.15514  0.327  
cond      1 0.0572557 0.01059861 0.26641  0.991  
lchla     1 0.3096201 0.05731384 1.44064  0.166  
tn        1 0.1236755 0.02289360 0.57545  0.814  
tp        1 0.3417611 0.06326345 1.59019  0.134  
reed      1 0.2324913 0.04303650 1.08176  0.345  
emerg     1 0.3821963 0.07074842 1.77833  0.084 .
Residual 11 2.3641047 0.43761981                 
Total    22 5.4021885 1.00000000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lemminvertmv <- mvabund(lemminvert1)
lemminvertmv.mv <- manyglm(lemminvertmv~lemmenv1$depth+lemmenv1$silt+lemmenv1$ph+lemmenv1$temp+lemmenv1$surface+lemmenv1$cond+lemmenv1$lchla+lemmenv1$tn+lemmenv1$tp+lemmenv1$reed+lemmenv1$emerg,family="negative_binomial")
# use summary to get Wald tests (at least vaguely comparable to Type III Permanova results)
summary(lemminvertmv.mv)

Test statistics:
                 wald value Pr(>wald)   
(Intercept)        52.51585     0.009 **
lemmenv1$depth      6.18052     0.376   
lemmenv1$silt       6.81522     0.268   
lemmenv1$ph         6.41864     0.364   
lemmenv1$temp       6.45409     0.364   
lemmenv1$surface    7.05475     0.262   
lemmenv1$cond       4.69355     0.786   
lemmenv1$lchla      7.41865     0.219   
lemmenv1$tn         5.29815     0.635   
lemmenv1$tp         6.96154     0.311   
lemmenv1$reed       6.07868     0.518   
lemmenv1$emerg      5.45911     0.615   
--- 
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Test statistic:  28.80876, p-value: 0.111 
Arguments:
 Test statistics calculated assuming response assumed to be uncorrelated 
 P-value calculated using 999 resampling iterations via pit.trap resampling (to account for correlation in testing).

plot(lemminvertmv.mv)

Generate graphs

br=c("nm","li","nf","yf")
la=c("None", "Light", "No fish", "Young fish")
a<-as.data.frame(lemminvert.mds$points)
a<-cbind(lemminvert[c(1:3)],a)   #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=MDS1, y=MDS2, shape=manag, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1", title="Raw")+
  scale_shape_manual(values=sym4,
                     name="Management",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()
as<-as.data.frame(lemminvert1s.mds$points)
as<-cbind(lemminvert[c(1:3)],as)   #Add site names & symbols from original data file
p2<-ggplot(data=as, aes(x=MDS1, y=MDS2, shape=manag, ) )+
  geom_point()+
  labs(y=NULL, x="MDS1", title="Standardized")+
  scale_shape_manual(values=sym4,
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()
aj<-as.data.frame(lemminvert.jac.mds$points)
aj<-cbind(lemminvert[c(1:3)],aj)   #Add site names & symbols from original data file
p3<-ggplot(data=aj, aes(x=MDS1, y=MDS2, shape=manag, ) )+
  geom_point()+
  labs(y=NULL, x="MDS1", title="Jaccard")+
  scale_shape_manual(values=sym4,
                      breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()
p4<-p1+p2+p3+plot_layout(guides='collect')

p4

Colour version

p1a<-ggplot(data=a, aes(x=MDS1, y=MDS2, color=manag, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1", title="Raw")+
  scale_color_uchicago(
                      breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()
p2a<-ggplot(data=as, aes(x=MDS1, y=MDS2, color=manag, ) )+
  geom_point()+
  labs(y=NULL, x="MDS1", title="Standardized")+
  scale_color_uchicago(
                     name="Management",
                     breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()
p3a<-ggplot(data=aj, aes(x=MDS1, y=MDS2, color=manag, ) )+
  geom_point()+
  labs(y=NULL, x="MDS1", title="Jacard")+
  scale_color_uchicago(
                      breaks=br,
                     labels=la,
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )+
  theme_qk()
p4a<-p1a+p2a+p3a+plot_layout(guides='collect')&theme(axis.text=element_blank())

p4a