NMDS analysis
V. Gårdman, B. Aulin, M. Molander & M. Jonsell
2025-02-23
This is the code for running the NMDS and Permanova models and from the article The potential of semiochemical blends to monitor saproxylic beetle communities published in Journal of Insect Conservation
Note: All packages used in the Rmarkdown script are automatically downloaded if needed, and unpacked, in the first (hidden) code chunk. This includes a package hosted on Github.
Step 1. Load and prepare data. In this example, the file is expected to be inside a folder named “Data” within the project directory. If your copy of the repository does not contain this folder, or if the file is stored elsewhere, you may need to adjust the path accordingly.
SpeciesMatrix <- read.csv(here("Data", "SpeciesMatrix.csv"),
sep=";", header=TRUE, as.is=TRUE,
na.strings=c(NA, "na", " NA"))
SpeciesMatrix <- SpeciesMatrix %>%
select(-Trap_Count, -Month) %>%
mutate(across(c("Harbor", "Blend", "Year"), as.factor)) %>%
group_by(Harbor, Blend, Year) %>%
summarise(across(where(is.numeric), sum), .groups = "drop")
#Standardize data
SpeciesMatrix <- SpeciesMatrix %>%
mutate(row_sum = rowSums(across(where(is.numeric))), # Dynamically select numeric columns
across(where(is.numeric), ~ . / row_sum)) %>% # Standardize
select(-row_sum) # Remove the temporary sum column
#Subset the dataframe on which to base the ordination (dataframe 1)
species_data <- SpeciesMatrix[,4:147]
#Identify the columns that contains the descriptive/environmental data (dataframe 2)
envi_data <- SpeciesMatrix[,1:3]Step 2: Fit multivariate model (NMDS). The standard two dimensions is used, with 1000 iterations.
NMDS_mod <- metaMDS(species_data, distance = "bray", k = 2, trymax = 1000)## Run 0 stress 0.1907318
## Run 1 stress 0.1853752
## ... New best solution
## ... Procrustes: rmse 0.04306565 max resid 0.1838044
## Run 2 stress 0.1941006
## Run 3 stress 0.1850025
## ... New best solution
## ... Procrustes: rmse 0.01321224 max resid 0.09390426
## Run 4 stress 0.1957739
## Run 5 stress 0.2089876
## Run 6 stress 0.201888
## Run 7 stress 0.1850025
## ... New best solution
## ... Procrustes: rmse 3.482412e-05 max resid 0.0001843095
## ... Similar to previous best
## Run 8 stress 0.1937152
## Run 9 stress 0.1915718
## Run 10 stress 0.1853752
## ... Procrustes: rmse 0.01319673 max resid 0.09373753
## Run 11 stress 0.2251218
## Run 12 stress 0.2003089
## Run 13 stress 0.2037894
## Run 14 stress 0.1983763
## Run 15 stress 0.2201509
## Run 16 stress 0.2239627
## Run 17 stress 0.2181841
## Run 18 stress 0.2005883
## Run 19 stress 0.1983921
## Run 20 stress 0.2045145
## *** Best solution repeated 1 times# Run twice to avoid local optimum.
NMDS_mod <- metaMDS(species_data, distance = "bray", k = 2, trymax = 1000,
previous.best = NMDS_mod)## Starting from 2-dimensional configuration
## Run 0 stress 0.1850025
## Run 1 stress 0.1957253
## Run 2 stress 0.2066428
## Run 3 stress 0.1941079
## Run 4 stress 0.1966414
## Run 5 stress 0.2272336
## Run 6 stress 0.1962833
## Run 7 stress 0.2048511
## Run 8 stress 0.1908475
## Run 9 stress 0.2170085
## Run 10 stress 0.2140594
## Run 11 stress 0.2014603
## Run 12 stress 0.2004468
## Run 13 stress 0.2051993
## Run 14 stress 0.2037236
## Run 15 stress 0.2122028
## Run 16 stress 0.1998858
## Run 17 stress 0.1938884
## Run 18 stress 0.221409
## Run 19 stress 0.1915587
## Run 20 stress 0.2132116
## *** Best solution repeated 1 timesStep 3: Plot the NMDS
datascores = as.data.frame(scores(NMDS_mod)$sites) #extract the site scores
datascores$Harbor = SpeciesMatrix$Harbor
datascores$Blend = SpeciesMatrix$Blend
datascores$Year = SpeciesMatrix$Year
#plot by blend and harbor
NMDS_plot <- ggplot(datascores, aes(x = NMDS1, y = NMDS2, colour = Blend)) +
geom_point(size = 2, aes(shape = Harbor, colour = Blend)) + # Points by Blend
coord_fixed() +
theme_bw() +
stat_ellipse(aes(color = Blend, linetype = Blend), level=0.9) +
theme(legend.position="right",legend.text=element_text(size=12),
legend.title=element_text(size=14),
legend.direction='vertical',
panel.grid.minor = element_blank(),
panel.grid.major = element_blank()) +
scale_color_brewer(palette = "Dark2") +
theme(legend.text = element_text(size = 12),
legend.title = element_text(size = 14))
Step 4: Fit Permanova model.
#Treat year as a random factor
permutations <- with(envi_data, how(nperm=999, blocks = Year))
#Calculate distance matrix
dist_matrix <- vegdist(species_data, method = "bray")
###Fit permanova model
Permanova_mod <- adonis2(dist_matrix ~ Blend * Harbor, data=envi_data,
permutations=permutations, method="bray")## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: Year
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist_matrix ~ Blend * Harbor, data = envi_data, permutations = permutations, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Blend 2 2.0508 0.13588 5.8164 0.001 ***
## Harbor 2 3.8455 0.25480 10.9065 0.001 ***
## Blend:Harbor 4 1.2627 0.08366 1.7906 0.009 **
## Residual 45 7.9333 0.52565
## Total 53 15.0923 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Check assumption of homogeneity of multivariate dispersion
BlendEffect <- anova(betadisper(dist_matrix, envi_data$Blend))
HarborEffect <- anova(betadisper(dist_matrix, envi_data$Harbor))## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.04747 0.023737 1.1264 0.3321
## Residuals 51 1.07472 0.021073## Analysis of Variance Table
##
## Response: Distances
## Df Sum Sq Mean Sq F value Pr(>F)
## Groups 2 0.08757 0.043783 3.0618 0.0555 .
## Residuals 51 0.72930 0.014300
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Step 5: Run pairwise comparisons between blends and harbors
#Extract the axes scores
datascores = as.data.frame(scores(NMDS_mod)$sites) #extract the site scores
# Combine the NMDS scores with the factors in 'envi_data'
datascores$Harbor = envi_data$Harbor
datascores$Blend = envi_data$Blend
# Create a new interaction factor for the pairwise comparisons
# Add the group_factor to envi_data
envi_data$group_factor <- interaction(envi_data$Harbor, envi_data$Blend)
# Fit the pairwise comparisons
Pairwise_mod <- pairwise.adonis2(dist_matrix ~ group_factor, data=envi_data,
permutations=999, method="bray", p.adjust.m="none")## [[1]]
## $parent_call
## [1] "dist_matrix ~ group_factor , strata = Null , permutations 999"
##
## $Gothenburg.Ips_vs_Gothenburg.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.92612 0.32812 4.8837 0.005 **
## Residual 10 1.89636 0.67188
## Total 11 2.82248 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Ips_vs_Gothenburg.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.91182 0.34149 5.1859 0.005 **
## Residual 10 1.75829 0.65851
## Total 11 2.67011 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Ips_vs_Monsteras.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.4614 0.53679 11.588 0.002 **
## Residual 10 1.2611 0.46321
## Total 11 2.7225 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Ips_vs_Monsteras.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.1264 0.40187 6.7189 0.002 **
## Residual 10 1.6764 0.59813
## Total 11 2.8028 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Ips_vs_Monsteras.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.1972 0.43176 7.5981 0.002 **
## Residual 10 1.5757 0.56824
## Total 11 2.7729 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Ips_vs_Norrkoping.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.5057 0.54642 12.047 0.003 **
## Residual 10 1.2499 0.45358
## Total 11 2.7556 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Ips_vs_Norrkoping.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.1732 0.42059 7.2589 0.003 **
## Residual 10 1.6162 0.57941
## Total 11 2.7895 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Ips_vs_Norrkoping.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.1095 0.36972 5.8658 0.002 **
## Residual 10 1.8914 0.63028
## Total 11 3.0008 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Monochamus_vs_Gothenburg.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.54247 0.19578 2.4344 0.026 *
## Residual 10 2.22834 0.80422
## Total 11 2.77081 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Monochamus_vs_Monsteras.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.7465 0.50222 10.089 0.002 **
## Residual 10 1.7311 0.49778
## Total 11 3.4777 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Monochamus_vs_Monsteras.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.1612 0.35107 5.4099 0.003 **
## Residual 10 2.1465 0.64893
## Total 11 3.3077 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Monochamus_vs_Monsteras.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.0585 0.34099 5.1742 0.002 **
## Residual 10 2.0457 0.65901
## Total 11 3.1042 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Monochamus_vs_Norrkoping.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.8227 0.51451 10.598 0.003 **
## Residual 10 1.7199 0.48549
## Total 11 3.5427 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Monochamus_vs_Norrkoping.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.1397 0.35328 5.4626 0.006 **
## Residual 10 2.0863 0.64672
## Total 11 3.2260 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Monochamus_vs_Norrkoping.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.2601 0.34794 5.3361 0.003 **
## Residual 10 2.3614 0.65206
## Total 11 3.6215 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Pissodes_vs_Monsteras.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.6935 0.51528 10.63 0.002 **
## Residual 10 1.5931 0.48472
## Total 11 3.2866 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Pissodes_vs_Monsteras.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.2452 0.38271 6.1998 0.004 **
## Residual 10 2.0084 0.61729
## Total 11 3.2536 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Pissodes_vs_Monsteras.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.74403 0.28059 3.9002 0.004 **
## Residual 10 1.90765 0.71941
## Total 11 2.65169 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Pissodes_vs_Norrkoping.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.7192 0.5208 10.868 0.006 **
## Residual 10 1.5819 0.4792
## Total 11 3.3011 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Pissodes_vs_Norrkoping.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.3979 0.41776 7.1752 0.002 **
## Residual 10 1.9482 0.58224
## Total 11 3.3461 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Gothenburg.Pissodes_vs_Norrkoping.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 1.1379 0.33853 5.1178 0.003 **
## Residual 10 2.2234 0.66147
## Total 11 3.3612 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Ips_vs_Monsteras.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.8248 0.35308 5.4579 0.007 **
## Residual 10 1.5112 0.64692
## Total 11 2.3360 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Ips_vs_Monsteras.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.41565 0.22762 2.947 0.048 *
## Residual 10 1.41044 0.77238
## Total 11 1.82609 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Ips_vs_Norrkoping.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.06142 0.05359 0.5663 0.724
## Residual 10 1.08466 0.94641
## Total 11 1.14608 1.00000
##
## $Monsteras.Ips_vs_Norrkoping.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.58564 0.28755 4.0361 0.02 *
## Residual 10 1.45101 0.71245
## Total 11 2.03666 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Ips_vs_Norrkoping.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.20958 0.10827 1.2141 0.261
## Residual 10 1.72616 0.89173
## Total 11 1.93573 1.00000
##
## $Monsteras.Monochamus_vs_Monsteras.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.44734 0.19679 2.4501 0.031 *
## Residual 10 1.82580 0.80321
## Total 11 2.27315 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Monochamus_vs_Norrkoping.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.82409 0.35458 5.4938 0.007 **
## Residual 10 1.50002 0.64542
## Total 11 2.32411 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Monochamus_vs_Norrkoping.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.17435 0.08543 0.9342 0.46
## Residual 10 1.86638 0.91457
## Total 11 2.04073 1.00000
##
## $Monsteras.Monochamus_vs_Norrkoping.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.41768 0.16321 1.9504 0.089 .
## Residual 10 2.14152 0.83679
## Total 11 2.55920 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Pissodes_vs_Norrkoping.Ips
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.4948 0.26124 3.5361 0.024 *
## Residual 10 1.3993 0.73876
## Total 11 1.8941 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Pissodes_vs_Norrkoping.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.50138 0.22117 2.8397 0.044 *
## Residual 10 1.76561 0.77883
## Total 11 2.26700 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Monsteras.Pissodes_vs_Norrkoping.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.27666 0.11938 1.3557 0.212
## Residual 10 2.04076 0.88062
## Total 11 2.31741 1.00000
##
## $Norrkoping.Ips_vs_Norrkoping.Monochamus
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.50687 0.26037 3.5203 0.039 *
## Residual 10 1.43983 0.73963
## Total 11 1.94670 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Norrkoping.Ips_vs_Norrkoping.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.13762 0.07428 0.8024 0.567
## Residual 10 1.71498 0.92572
## Total 11 1.85259 1.00000
##
## $Norrkoping.Monochamus_vs_Norrkoping.Pissodes
## Df SumOfSqs R2 F Pr(>F)
## group_factor 1 0.25756 0.11012 1.2375 0.255
## Residual 10 2.08133 0.88988
## Total 11 2.33889 1.00000
##
## attr(,"class")
## [1] "pwadstrata" "list"References:
## @Manual{,
## title = {vegan: Community Ecology Package},
## author = {Jari Oksanen and Gavin L. Simpson and F. Guillaume Blanchet and Roeland Kindt and Pierre Legendre and Peter R. Minchin and R.B. O'Hara and Peter Solymos and M. Henry H. Stevens and Eduard Szoecs and Helene Wagner and Matt Barbour and Michael Bedward and Ben Bolker and Daniel Borcard and Gustavo Carvalho and Michael Chirico and Miquel {De Caceres} and Sebastien Durand and Heloisa Beatriz Antoniazi Evangelista and Rich FitzJohn and Michael Friendly and Brendan Furneaux and Geoffrey Hannigan and Mark O. Hill and Leo Lahti and Dan McGlinn and Marie-Helene Ouellette and Eduardo {Ribeiro Cunha} and Tyler Smith and Adrian Stier and Cajo J.F. {Ter Braak} and James Weedon},
## year = {2024},
## note = {R package version 2.6-6.1},
## url = {https://CRAN.R-project.org/package=vegan},
## }## @Manual{,
## title = {here: A Simpler Way to Find Your Files},
## author = {Kirill Müller},
## year = {2020},
## note = {R package version 1.0.1},
## url = {https://CRAN.R-project.org/package=here},
## }## @Article{,
## title = {Welcome to the {tidyverse}},
## author = {Hadley Wickham and Mara Averick and Jennifer Bryan and Winston Chang and Lucy D'Agostino McGowan and Romain François and Garrett Grolemund and Alex Hayes and Lionel Henry and Jim Hester and Max Kuhn and Thomas Lin Pedersen and Evan Miller and Stephan Milton Bache and Kirill Müller and Jeroen Ooms and David Robinson and Dana Paige Seidel and Vitalie Spinu and Kohske Takahashi and Davis Vaughan and Claus Wilke and Kara Woo and Hiroaki Yutani},
## year = {2019},
## journal = {Journal of Open Source Software},
## volume = {4},
## number = {43},
## pages = {1686},
## doi = {10.21105/joss.01686},
## }## @Manual{,
## title = {pairwiseAdonis: Pairwise Multilevel Comparison using Adonis},
## author = {Pedro {Martinez Arbizu}},
## year = {2017},
## note = {R package version 0.4.1, commit cb190f7668a0c82c0b0853927db239e7b9ec3e83},
## url = {https://github.com/pmartinezarbizu/pairwiseAdonis},
## }## @Manual{,
## title = {RColorBrewer: ColorBrewer Palettes},
## author = {Erich Neuwirth},
## year = {2022},
## note = {R package version 1.1-3},
## url = {https://CRAN.R-project.org/package=RColorBrewer},
## }