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.115322946 
... New best solution
... Procrustes: rmse 0.0189448746  max resid 0.07577272532 
Run 2 stress 0.1153228722 
... New best solution
... Procrustes: rmse 0.0001332765475  max resid 0.0005149485731 
... Similar to previous best
Run 3 stress 0.1172333437 
Run 4 stress 0.1156700569 
... Procrustes: rmse 0.01891880149  max resid 0.07525363734 
Run 5 stress 0.1156700851 
... Procrustes: rmse 0.01895538217  max resid 0.0755501726 
Run 6 stress 0.1153230967 
... Procrustes: rmse 0.0002702610702  max resid 0.001057340867 
... Similar to previous best
Run 7 stress 0.1153229703 
... Procrustes: rmse 0.0001739346127  max resid 0.0006742127962 
... Similar to previous best
Run 8 stress 0.1156700194 
... Procrustes: rmse 0.01892698264  max resid 0.07536101802 
Run 9 stress 0.1153229521 
... Procrustes: rmse 0.0001457488044  max resid 0.0005644363622 
... Similar to previous best
Run 10 stress 0.1172331317 
Run 11 stress 0.1156701841 
... Procrustes: rmse 0.01890121632  max resid 0.07510464765 
Run 12 stress 0.1156701456 
... Procrustes: rmse 0.01890185452  max resid 0.07512527764 
Run 13 stress 0.1153229505 
... Procrustes: rmse 0.0001607374236  max resid 0.0006251011904 
... Similar to previous best
Run 14 stress 0.1153230711 
... Procrustes: rmse 0.0002520185237  max resid 0.0009852698337 
... Similar to previous best
Run 15 stress 0.1153228854 
... Procrustes: rmse 6.752969035e-05  max resid 0.0002566588905 
... Similar to previous best
Run 16 stress 0.1153230349 
... Procrustes: rmse 0.0002138194223  max resid 0.0008367394134 
... Similar to previous best
Run 17 stress 0.1156700584 
... Procrustes: rmse 0.01894639128  max resid 0.07549075399 
Run 18 stress 0.1153229199 
... Procrustes: rmse 0.0001164512099  max resid 0.0004509812175 
... Similar to previous best
Run 19 stress 0.1156700244 
... Procrustes: rmse 0.01893168798  max resid 0.07538220799 
Run 20 stress 0.1156701597 
... Procrustes: rmse 0.0188997287  max resid 0.07511089486 
*** Best solution repeated 9 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.1153228722 
Stress type 1, weak ties
Best solution was repeated 9 times in 20 tries
The best solution was from try 2 (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.1649065617 
... New best solution
... Procrustes: rmse 0.01494844861  max resid 0.05093303219 
Run 2 stress 0.1763011963 
Run 3 stress 0.1969958637 
Run 4 stress 0.1648801187 
... New best solution
... Procrustes: rmse 0.008275086618  max resid 0.03127678786 
Run 5 stress 0.1709487893 
Run 6 stress 0.2054917683 
Run 7 stress 0.170948787 
Run 8 stress 0.1783397316 
Run 9 stress 0.2008724331 
Run 10 stress 0.2197541405 
Run 11 stress 0.1783397314 
Run 12 stress 0.1790250486 
Run 13 stress 0.1763011964 
Run 14 stress 0.1793740312 
Run 15 stress 0.1736187217 
Run 16 stress 0.2005321408 
Run 17 stress 0.2031756856 
Run 18 stress 0.1847458372 
Run 19 stress 0.1742414652 
Run 20 stress 0.1856856402 
Run 21 stress 0.1649421052 
... Procrustes: rmse 0.01127697205  max resid 0.04881201499 
Run 22 stress 0.1832527941 
Run 23 stress 0.1707635696 
Run 24 stress 0.1736187271 
Run 25 stress 0.1779554366 
Run 26 stress 0.1832527928 
Run 27 stress 0.1835084232 
Run 28 stress 0.1742453837 
Run 29 stress 0.1648801121 
... New best solution
... Procrustes: rmse 1.485100402e-05  max resid 2.888931322e-05 
... Similar to previous best
*** 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.1648801121 
Stress type 1, weak ties
Best solution was repeated 1 time in 29 tries
The best solution was from try 29 (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.2107007937 
Run 2 stress 0.1996370419 
Run 3 stress 0.1972852048 
Run 4 stress 0.1676823255 
... New best solution
... Procrustes: rmse 0.0002381066046  max resid 0.0007976072167 
... Similar to previous best
Run 5 stress 0.1676824713 
... Procrustes: rmse 0.0002957253877  max resid 0.000990869083 
... Similar to previous best
Run 6 stress 0.167682628 
... Procrustes: rmse 0.0003271587264  max resid 0.001111299241 
... Similar to previous best
Run 7 stress 0.2080583428 
Run 8 stress 0.1676826526 
... Procrustes: rmse 0.0003639721309  max resid 0.001238599788 
... Similar to previous best
Run 9 stress 0.2057095586 
Run 10 stress 0.2153016562 
Run 11 stress 0.1676824279 
... Procrustes: rmse 0.0001916372554  max resid 0.0006558477933 
... Similar to previous best
Run 12 stress 0.167682576 
... Procrustes: rmse 0.0003141864367  max resid 0.001071523201 
... Similar to previous best
Run 13 stress 0.1676824319 
... Procrustes: rmse 0.0001718822342  max resid 0.0005840181324 
... Similar to previous best
Run 14 stress 0.2096575865 
Run 15 stress 0.1676827318 
... Procrustes: rmse 0.0004701871926  max resid 0.001586789546 
... Similar to previous best
Run 16 stress 0.2045142119 
Run 17 stress 0.1926151338 
Run 18 stress 0.167682351 
... Procrustes: rmse 6.396017731e-05  max resid 0.0002169851584 
... Similar to previous best
Run 19 stress 0.1676823616 
... Procrustes: rmse 0.0001270003117  max resid 0.000416667006 
... Similar to previous best
Run 20 stress 0.2016022116 
*** Best solution repeated 10 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.1676823255 
Stress type 1, weak ties
Best solution was repeated 10 times in 20 tries
The best solution was from try 4 (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, ...
ℹ 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.04739254 -0.99887634 0.17444  0.147  
silt     0.86955032  0.49384434 0.12896  0.265  
ph       0.29026870 -0.95694518 0.12647  0.253  
temp    -0.16698561 -0.98595933 0.00161  0.989  
surface  0.01866995 -0.99982570 0.21048  0.106  
cond    -0.19976409 -0.97984402 0.08806  0.432  
chla     0.59023073 -0.80723459 0.08147  0.411  
tn      -0.49309516 -0.86997538 0.19841  0.117  
tp      -0.36949074 -0.92923441 0.34855  0.021 *
reed    -0.67315093  0.73950512 0.21729  0.091 .
emerg   -0.62436360  0.78113385 0.12489  0.281  
lchla    0.15591532 -0.98777043 0.15802  0.182  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999

plot(fit1)

Error in strwidth(labels, cex = cex, ...) : 
  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.97880776 -0.20478126 0.15296  0.204  
silt    -0.10985270 -0.99394788 0.02816  0.718  
ph       0.94740494  0.32003732 0.08464  0.419  
temp     0.20181041 -0.97942461 0.02851  0.787  
surface  0.86028214 -0.50981824 0.13431  0.214  
cond     0.46289144 -0.88641498 0.08084  0.419  
chla     0.64368593  0.76528976 0.02808  0.763  
tn       0.33297274  0.94293645 0.15822  0.179  
tp       0.57504149  0.81812425 0.14517  0.179  
reed    -0.70398692  0.71021293 0.26946  0.058 .
emerg   -0.79831105  0.60224536 0.15104  0.191  
lchla    0.98631265  0.16488587 0.11666  0.314  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999

plot(fit2)

Error in strwidth(labels, cex = cex, ...) : 
  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.67050446 -0.74190550 0.02424  0.782
silt    -0.81424703 -0.58051853 0.00658  0.929
ph       0.05657893 -0.99839813 0.04594  0.609
temp     0.56727107 -0.82353114 0.03545  0.688
surface -0.01791782 -0.99983946 0.06936  0.473
cond    -0.07819187 -0.99693833 0.05738  0.496
chla     0.11377452 -0.99350660 0.01580  0.824
tn       0.10835065 -0.99411274 0.04330  0.668
tp       0.24001764 -0.97076853 0.10239  0.336
reed     0.36901297 -0.92942425 0.06176  0.531
emerg    0.89198410  0.45206677 0.19599  0.103
lchla    0.43934019 -0.89832076 0.01534  0.851
Permutation: free
Number of permutations: 999

plot(fit3)

Error in strwidth(labels, cex = cex, ...) : 
  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.02 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0939 0.1264 0.1606 0.1922 

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
Permutation: free
Number of permutations: 999

adonis2(formula = lemminvert.bc ~ manag, data = lemminvert, permutations = 999)
         Df  SumOfSqs         R2       F Pr(>F)   
Model     3 1.3565598 0.25111301 2.12366  0.002 **
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)

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

Model: lemminvertmv ~ lemmenv$manag

Multivariate test:
              Res.Df Df.diff      Dev Pr(>Dev)    
(Intercept)       22                              
lemmenv$manag     19       3 232.0596    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.

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.07 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0958 0.1363 0.1644 0.2052 

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
Permutation: free
Number of permutations: 999

adonis2(formula = lemminvert1s.bc ~ manag, data = lemminvert, permutations = 999)
         Df  SumOfSqs         R2       F Pr(>F)  
Model     3 0.8296405 0.23689983 1.96614   0.03 *
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 'data2' not found

plot(lemminvert1smv.mv)

Error: object 'lemminvert1smv.mv' not found

anova(lemminvert1smv.mv)

Error: object 'lemminvert1smv.mv' not found

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.077 

Permutation: free
Number of permutations: 999

Upper quantiles of permutations (null model):
   90%    95%  97.5%    99% 
0.0936 0.1144 0.1389 0.1722 

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
Permutation: free
Number of permutations: 999

adonis2(formula = lemminvert.jac ~ manag, data = lemminvert, permutations = 999)
         Df   SumOfSqs         R2       F Pr(>F)  
Model     3 0.54394094 0.21570186 1.74183  0.013 *
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 0 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.025 *
---
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.622 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.187 0.242 0.284 0.346 
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.515 

Upper quantiles of permutations (null model):
  90%   95% 97.5%   99% 
0.233 0.297 0.364 0.404 
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
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)
Model    11 3.0380838 0.56238019 1.28509  0.107
Residual 11 2.3641047 0.43761981               
Total    22 5.4021885 1.00000000               

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.156  
silt      1 0.1810559 0.03351528 0.84244  0.565  
ph        1 0.2562631 0.04743690 1.19237  0.299  
temp      1 0.2256969 0.04177879 1.05015  0.376  
surface   1 0.2482610 0.04595563 1.15514  0.329  
cond      1 0.0572557 0.01059861 0.26641  0.989  
lchla     1 0.3096201 0.05731384 1.44064  0.156  
tn        1 0.1236755 0.02289360 0.57545  0.795  
tp        1 0.3417611 0.06326345 1.59019  0.133  
reed      1 0.2324913 0.04303650 1.08176  0.396  
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.010 **
lemmenv1$depth      6.18052     0.389   
lemmenv1$silt       6.81522     0.234   
lemmenv1$ph         6.41864     0.352   
lemmenv1$temp       6.45409     0.382   
lemmenv1$surface    7.05475     0.239   
lemmenv1$cond       4.69355     0.769   
lemmenv1$lchla      7.41865     0.221   
lemmenv1$tn         5.29815     0.642   
lemmenv1$tp         6.96154     0.318   
lemmenv1$reed       6.07868     0.513   
lemmenv1$emerg      5.45911     0.592   
--- 
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Test statistic:  28.80876, p-value: 0.105 
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