---
title: "NMDS analysis"
author: "V. Gårdman, B. Aulin, M. Molander & M. Jonsell"
date: "2025-02-23"
output: html_document
---

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

```{r, echo = FALSE, warning=FALSE, message=FALSE}
# ---- Setup: install and load packages ----

# Install/load CRAN packages
load_or_install <- function(packages) {
  for (pkg in packages) {
    if (!requireNamespace(pkg, quietly = TRUE)) {
      install.packages(pkg)
    }
    library(pkg, character.only = TRUE)
  }
}

# CRAN packages
cran_packages <- c(
  "tidyverse", "here", "vegan",
  "RColorBrewer"
)
load_or_install(cran_packages)

# GitHub package
if (!requireNamespace("devtools", quietly = TRUE))
  install.packages("devtools")

if (!requireNamespace("pairwiseAdonis", quietly = TRUE)) {
  devtools::install_github("pmartinezarbizu/pairwiseAdonis")
}
library(pairwiseAdonis)
```



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

```{r, warning=FALSE, message=FALSE}
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.

```{r, warning=FALSE, message=FALSE}
NMDS_mod <- metaMDS(species_data, distance = "bray", k = 2, trymax = 1000)
# Run twice to avoid local optimum.
NMDS_mod <- metaMDS(species_data, distance = "bray", k = 2, trymax = 1000, 
                      previous.best = NMDS_mod)
```

### Step 3: Plot the NMDS

```{r, warning=FALSE, message=FALSE}
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))
```
```{r, echo = FALSE, warning=FALSE, message=FALSE, fig.width=12, fig.height=8}
NMDS_plot
```


### Step 4: Fit Permanova model.


```{r, warning=FALSE, message=FALSE}
#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")
```
```{r, echo = FALSE, warning=FALSE, message=FALSE}
Permanova_mod
```
```{r, warning=FALSE, message=FALSE}
#Check assumption of homogeneity of multivariate dispersion
BlendEffect <- anova(betadisper(dist_matrix, envi_data$Blend))

HarborEffect <- anova(betadisper(dist_matrix, envi_data$Harbor))
```
```{r, echo = FALSE, warning=FALSE, message=FALSE}
BlendEffect
HarborEffect
```
### Step 5: Run pairwise comparisons between blends and harbors

```{r, warning=FALSE, message=FALSE}
#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")
```
```{r, echo = FALSE, warning=FALSE, message=FALSE}
list(Pairwise_mod)
```

## References:
```{r, echo = FALSE, warning=FALSE, message=FALSE}
toBibtex(citation("vegan"))
toBibtex(citation("here"))
toBibtex(citation("tidyverse"))
toBibtex(citation("pairwiseAdonis"))
toBibtex(citation("RColorBrewer"))
```