library(ape)
library(ggplot2)
## Warning: le package 'ggplot2' a été compilé avec la version R 4.4.3
library(corrplot)
## corrplot 0.94 loaded
library(vegan)
## Le chargement a nécessité le package : permute
## Le chargement a nécessité le package : lattice
## This is vegan 2.6-8
library(readr)
## Warning: le package 'readr' a été compilé avec la version R 4.4.3

MULTIDIMENSIONAL SCALING

Exercice 1

Aim: Create an index summarizing flower morphology to compare species.
=> Focus on variables : PCA

Aim: Determine whether individuals within a species are more similar in terms of flower morphology than individuals of different species. Hypothesis: Iris species differ in petal and sepal dimensions.
=> Focus on individuals : MDS

  • Distances between points => Classic Multidimensional Scaling = PCoA

  • Rank order => Non-Metric Multidimensional Scaling (NMDS)

data(iris)
View(iris)
# • Quantitative continuous variables
# • Same measurement scales
# • No missing data
# • Natural 0 (non existant dans le tableau)
## Chioce of metric : Tu compares tes sites (ou tes échantillons) un par un. Tu veux savoir si le Site A ressemble au Site B en fonction de leur pH, de leur taux de calcaire, etc. Chaque point sur ton graphique final représente un seul site.  => Euclidian distance 

Build a similarity/distance matrix based on the chosen distance coefficient

# Distance matrix between individuals
data(iris)
iris_num <- iris[, 1:4]
iris_num <- scale(iris_num) # Les variances peuvent différer entre les pétales et les sépales. La standardisation (Z-score) permet de donner le même poids à chaque variable dans le calcul de la distance.
dist_eucl <- dist(iris_num, method = "euclidean")

PCoA

# Transformation de matrice de distances en coordonnées spatiales
pcoa_iris <- cmdscale(dist_eucl, 
                      k = 4,      # 4 variables initiales
                      eig = TRUE) # valeurs propres (eigenvalues)
# eig and points values 

Analyser l’Inertie (Variance). Les deux premiers axes suffisent à résumer les données ?

# Récupérer les valeurs propres (eig) = Inertie // Variance Totale
pcoa_eig <- pcoa_iris$eig[1:4]
percent_inertie <- (pcoa_eig / sum(pcoa_eig)) * 100
cumul_inertie <- cumsum(percent_inertie)

inertia_table <- data.frame(Axe = 1:4,
                            Inertie = pcoa_eig,
                            Pourcentage_Inertie = percent_inertie,
                            Cumululative_Inertie = cumul_inertie)
inertia_table
##   Axe    Inertie Pourcentage_Inertie Cumululative_Inertie
## 1   1 434.856175          72.9624454             72.96245
## 2   2 136.190540          22.8507618             95.81321
## 3   3  21.866774           3.6689219             99.48213
## 4   4   3.086511           0.5178709            100.00000
scree_plot_data <- data.frame(Axe = factor(1:4),         # Factor 
                              Inertie = percent_inertie,
                              Cumululative_Inertie = cumul_inertie)

ggplot(scree_plot_data, aes(x = Axe)) +
  geom_bar(aes(y = Inertie), stat = "identity", fill = "lightblue") +  geom_line(aes(y = Cumululative_Inertie, group = 1), color = "blue") +  geom_point(aes(y = Cumululative_Inertie), color = "darkblue") +
  labs(title = "Scree Plot de la PCoA (Données Iris Standardisées)",
       x = "Dimension Principale (Axe)", 
       y = "Pourcentage d'Inertie Expliquée (%)")

# point value 
pcoa_points <- as.data.frame(pcoa_iris$points)
colnames(pcoa_points) <- c("Axe1", "Axe2", "Axe3", "Axe4")
pcoa_points$Species <- iris$Species

ggplot(pcoa_points, aes(x = Axe1, y = Axe2, color = Species)) +
  geom_point() +
  stat_ellipse(level = 0.95, linetype = 2) + # Ellipses de confiance à 95%
  scale_color_manual(values = c("setosa" = "#E41A1C", 
                                "versicolor" = "#377EB8", 
                                "virginica" = "#4DAF4A")) +
  labs(title = "Score Plot de la PCoA (Mode Q)",
       x = paste0("Axe 1 (", round(percent_inertie[1], 1), "%)"),
       y = paste0("Axe 2 (", round(percent_inertie[2], 1), "%)")) +
  theme_minimal()

## Information on:
# - Intraspecific variation
# - Interspecific variation
# - How different individuals and species are 

## Species differ in petal dimensions and sepal length
# Intraspecific variation mainly in sepal width
# Overlap between versicolor and virginica species -> difficulites to identify the species based on flower morphology

Calculer la matrice de corrélation

# On compare iris_num (variables originales) et pcoa_scores (nouveaux axes)
pcoa_points <- pcoa_iris$points
cor <- cor(iris_num, pcoa_points)
cor
##                    [,1]       [,2]        [,3]        [,4]
## Sepal.Length  0.8901688 0.36082989  0.27565767  0.03760602
## Sepal.Width  -0.4601427 0.88271627 -0.09361987 -0.01777631
## Petal.Length  0.9915552 0.02341519 -0.05444699 -0.11534978
## Petal.Width   0.9649790 0.06399985 -0.24298265  0.07535950
col_palette <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF",
                                  "#77AADD", "#4477AA"))
corrplot(cor, 
         method = "circle",
         col = col_palette(200), addCoef.col = "black", tl.col = "black")

# Species differ in petal dimensions and sepal length
# Intraspecific variation mainly in sepal width
# Overlap between versicolor and virginica species -> difficulites to identify the species based on flower morphology

Exercice 2

Context: Investigating feeding specificity of marine herbivores on different macroalgal habitats in Sydney Harbour.
Aim: Determine whether microalgae habitats are similar based on the composition of marine herbivore feeding on them. Hypothesis: microalgal habitats differ in herbivore species composition
.

  • Five species of microalgae, 20 replicates

  • Seven species of herbivorous crustacean

  • Abundance of custaceans recorded

  • 2 sampling times: Day and Night

herb <- read_delim("Herbivore_specialisation.csv", 
                   delim = ";", escape_double = FALSE, trim_ws = TRUE)
## Rows: 100 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr (2): Habitat, DayNight
## dbl (7): Peramphithoe_parmerong, Ampithoe_caddi, Ampithoe_kava, Ampithoe_nga...
## 
## ℹ 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.
# => Ce concentre sur des rangs => NMDS car comptage 
# • Quantitative discrete variables (exemple :1,2,3etc) 
# pas du presence/abscence mais abscene et plusieur presence donc continue
# • Same measurement scales
# • Missing data -> 0 = True absence of a species car pas observé
## Na meaning : Missing data or Abscence of observation ?
# Meaning of 0 : Missed individual or True abscence ?
# petite parcelle donc on sait qu'elle n'etait pas présente dans cette parcelle car elle etait petite (not a missed individual). signifie une abscence d'observation et donc remplace les NAs par 0.
# • 0 in 2 sites : 
# - Similar sites -> symmetrical coefficient (X)
# - Not similar -> asymmetrica (V)
# on ne se concentre sur 7 especes uniquement donc 0;0 ne seront pas considerer comme similaire
# • Toute nos sp sont sur le meme point d'egalité et ignore les double 0= bray-curtis
herb[is.na(herb)] <- 0
summary(herb)
##    Habitat            DayNight         Peramphithoe_parmerong Ampithoe_caddi  
##  Length:100         Length:100         Min.   : 0.00          Min.   :  0.00  
##  Class :character   Class :character   1st Qu.: 0.00          1st Qu.:  0.00  
##  Mode  :character   Mode  :character   Median : 0.00          Median :  5.00  
##                                        Mean   :10.90          Mean   : 17.37  
##                                        3rd Qu.:16.25          3rd Qu.: 19.25  
##                                        Max.   :96.00          Max.   :179.00  
##  Ampithoe_kava   Ampithoe_ngana  Cymadusa_munnu Exampithoe_kutti
##  Min.   :  0.0   Min.   :  0.0   Min.   :0.00   Min.   :  0.00  
##  1st Qu.:  0.0   1st Qu.:  1.0   1st Qu.:0.00   1st Qu.:  1.00  
##  Median :  4.0   Median : 10.0   Median :0.00   Median :  4.00  
##  Mean   : 16.5   Mean   : 26.3   Mean   :0.25   Mean   : 14.67  
##  3rd Qu.: 22.0   3rd Qu.: 40.5   3rd Qu.:0.00   3rd Qu.: 12.00  
##  Max.   :159.0   Max.   :136.0   Max.   :6.00   Max.   :167.00  
##  Plumithoe_quadrimanus
##  Min.   : 0.00        
##  1st Qu.: 0.00        
##  Median : 0.00        
##  Mean   : 3.28        
##  3rd Qu.: 0.00        
##  Max.   :50.00
herb <- transform(herb,  Habitat=factor(Habitat), DayNight=factor(DayNight))
herb_num <- herb[,3:9]
# herb_num <- sqrt(herb_num) # Regler le probleme d'asymetrie pour éviter qu'une espèce ultra-dominante n'écrase les autres.

NMDS: Distance matrix between individuals if needed

set.seed(42) # Pour que les résultats soient reproductibles (méthode itérative)
nmds_herb <- metaMDS(comm = herb_num, distance = "bray", k = 2, 
                     try = 40, trymax = 50,
                     trace = FALSE, autotransform = T)
# solution obtenue par un algorithme itératif → il faut plusieurs essais (`try`, `trymax`) pour vérifier la convergence vers une solution stable.

MDS: Score plot with habitats

nmds_herb$points
##            MDS1         MDS2
## 1   -1.29283831  0.134526635
## 2   -1.15478777 -0.162559913
## 3   -1.22528257 -0.923218578
## 4   -1.10281157 -1.160884210
## 5   -1.63993724 -0.545861932
## 6   -1.01509413 -0.758422137
## 7   -1.35783303  0.008862767
## 8   -1.26207145 -0.313466879
## 9   -1.84175224 -0.156551593
## 10  -0.94959439  0.054171670
## 11  -0.71308094 -0.150273649
## 12  -1.29167821  0.150256875
## 13  -0.86118045  1.332021316
## 14  -1.33478560 -0.705717212
## 15  -1.22534021 -0.225918345
## 16  -1.16393058 -0.624534194
## 17  -0.82098866 -1.023846701
## 18  -1.40990529 -0.064783969
## 19  -0.98517061 -0.107021981
## 20  -1.25778530 -0.224739696
## 21   0.19705869  0.248179572
## 22   0.12960516  0.288338797
## 23   0.15070226  0.199340648
## 24   0.08050478  0.266318974
## 25   0.17751161  0.270261628
## 26   0.10740119  0.193953140
## 27   0.19364230  0.059409144
## 28   0.25068178  0.180259625
## 29  -0.05396935  0.398234005
## 30   0.05785329  0.272319563
## 31  -0.12117428  0.435627757
## 32   0.15966497  0.362020312
## 33   0.05187010  0.245256368
## 34  -0.05929267  0.374109328
## 35   0.04172454  0.279690113
## 36  -0.51090844  0.967256048
## 37  -0.25103070  0.594505386
## 38  -0.09917941  0.426595332
## 39   0.16595938  0.249574749
## 40   0.04929399  0.274225486
## 41   0.54502391 -0.093041952
## 42   0.62355475  0.122757481
## 43   0.32159756 -0.087230039
## 44   0.33890522  0.015307599
## 45   0.48738107 -0.142335598
## 46   0.63993564 -0.074222783
## 47   0.50483292 -0.002656278
## 48   0.28516422 -0.232401675
## 49   0.38589343  0.031514042
## 50   0.54487296  0.201265058
## 51   0.39467131  0.069203142
## 52   0.46507310  0.232690079
## 53   0.51103601 -0.073172200
## 54   0.33082402  0.005404667
## 55   0.48729396 -0.037101593
## 56   0.50198394  0.096746658
## 57   0.51111593  0.045056548
## 58   0.50059279 -0.217534533
## 59   0.59177971  0.054159918
## 60   0.42199748 -0.132345307
## 61   0.50301778 -0.458558731
## 62   0.53194196 -0.571388625
## 63   0.65797479 -0.367359480
## 64   0.40373896 -0.413852754
## 65   0.66806670 -0.480575479
## 66   0.79866252 -0.391985614
## 67   0.65140787 -0.241588826
## 68   0.40131935 -0.121975515
## 69   0.39075042 -0.443267392
## 70   0.79347233 -0.493857746
## 71   0.47757716 -0.371878120
## 72   0.87706670 -0.959207594
## 73   0.80636807 -0.843606174
## 74   0.78616799 -0.822874885
## 75   0.71727275 -0.732486492
## 76   0.63765675 -0.290427399
## 77   0.44371938 -0.235394817
## 78   0.77787355 -1.051756087
## 79   0.44703626 -0.152235323
## 80   1.07427255 -1.328740752
## 81   0.16124854  0.420584199
## 82  -0.24756611 -0.051613926
## 83   0.59861241  0.494758873
## 84  -0.05667038  0.917346314
## 85   0.15327790  1.008992284
## 86   0.29573548  0.623715691
## 87   0.45479012  0.576923591
## 88  -0.86029655 -0.137971549
## 89   0.26802632  0.285494635
## 90  -0.16392466  0.864860536
## 91   0.06787251  0.633920570
## 92  -0.21778162  0.046481063
## 93   0.10922969  0.980136221
## 94  -0.16013919 -0.035574372
## 95  -0.44209401  0.340411887
## 96   0.37156393  0.539018119
## 97   0.09614012  0.643978962
## 98   0.07915572  0.622903847
## 99   0.34007796  0.368999513
## 100  0.10177741  0.730073861
## attr(,"centre")
## [1] TRUE
## attr(,"pc")
## [1] TRUE
## attr(,"halfchange")
## [1] TRUE
## attr(,"internalscaling")
## [1] 1.183399
# Extraction des coordonnées (sites)
nmds_points <- as.data.frame(scores(nmds_herb)$sites)
nmds_points$Habitat <- herb$Habitat
nmds_points$DayNight <- herb$DayNight

# Graphique par HABITAT
ggplot(nmds_points, aes(x = NMDS1, y = NMDS2, color = Habitat)) +
  geom_point() +
  stat_ellipse(level = 0.95) +
  theme_minimal() +
  labs(title = "NMDS : Similarité des Habitats",
       subtitle = paste("Stress =", round(nmds_herb$stress, 3)))

# • Uninterpretable axes: on a un comptage d'sp pour chaque habitat donc peut pas interpreter les axes. There are no eigenvalues or percentage of variation explained by axis because the axes are not combinations of variables.​
# • Samples are more likely to be similar within a habitat than between habitats
# • Species composition of herbivores differs across microalgal habitats

MDS: Score plot with sampling time

# Graphique par Jour/Nuit
ggplot(nmds_points, aes(x = NMDS1, y = NMDS2, color = DayNight)) +
  geom_point() +
  stat_ellipse(level = 0.95) +
  theme_minimal() +
  labs(title = "NMDS : Différence Jour vs Nuit",
       subtitle = paste("Stress =", round(nmds_herb$stress, 3)))

# Si les ellipses Jour et Nuit sont l'une sur l'autre, le moment de l'échantillonnage n'influence pas la communauté trouvée sur les algues.

MDS: Eigenvalues / % of inertia

# Piège de l'exercice : En NMDS, il n'y a pas de valeurs propres (Eigenvalues) car ce n'est pas une décomposition spectrale comme la PCA ou la PCoA.
# En NMDS, on mesure la qualité par le Stress. Si ton professeur demande quand même l'inertie, c'est peut-être un test pour voir si tu as compris la différence, ou alors il veut que tu fasses une PCoA en parallèle. Pour le NMDS, on regarde le Shepard Plot.

nmds_herb$stress # Si stress < 0.2, représentation des points acceptable.
## [1] 0.1260513
stressplot(nmds_herb) # Shepard Plot (Qualité de l'ajustement) = mesure l’écart entre les dissimilarités réelles et celles de la configuration réduite.

MDS: Correlations between axes and original variables

Exercice 3

Context: Study of Boloria eunomia abundance. Compare sites based on the plant species they contain.

  • 4 sampled sites divided into plots in which a quadrat is created to record the abundance of each species

  • A quadrat is composed of 25 squares

  • Species abundance = number of squares in which the spcies is present

Sites :

  • Pisserote: bog (very humid)
    Bihain: humid meadow
    Grande Fange: strong gradient from humid meadow to bog
    Petite Langlire: one invaded by « eagle fern », one bog

Aim: Determine whether sites are similar based on the plant composition

vegetation <- read_delim("RelevesVegetation_Clean.csv", 
    delim = ";", escape_double = FALSE, trim_ws = TRUE)
## Rows: 165 Columns: 55
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ";"
## chr  (1): Site
## dbl (53): Sample ID, Bambou, Achillea millefolium, Achillea ptarmica, Ajuga ...
## lgl  (1): Digitalis purpurea
## 
## ℹ 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.
View(vegetation)
summary(vegetation)
##    Sample ID          Site               Bambou       Achillea millefolium
##  Min.   :  1.00   Length:165         Min.   :  1.00   Min.   :6           
##  1st Qu.: 45.00   Class :character   1st Qu.: 52.00   1st Qu.:6           
##  Median : 86.00   Mode  :character   Median :100.00   Median :6           
##  Mean   : 85.84                      Mean   : 99.81   Mean   :6           
##  3rd Qu.:127.00                      3rd Qu.:148.00   3rd Qu.:6           
##  Max.   :168.00                      Max.   :199.00   Max.   :6           
##                                                       NA's   :164         
##  Achillea ptarmica Ajuga reptans Anemonoides nemorosa Angelica sylvestris
##  Min.   :2         Min.   :2     Min.   : 1.000       Min.   : 1.000     
##  1st Qu.:2         1st Qu.:2     1st Qu.: 2.000       1st Qu.: 4.000     
##  Median :2         Median :2     Median : 4.000       Median : 7.000     
##  Mean   :2         Mean   :2     Mean   : 5.526       Mean   : 8.487     
##  3rd Qu.:2         3rd Qu.:2     3rd Qu.: 7.500       3rd Qu.:10.500     
##  Max.   :2         Max.   :2     Max.   :20.000       Max.   :24.000     
##  NA's   :164       NA's   :164   NA's   :146          NA's   :126        
##  Athyrium filix-femina Betula pendula Betula pubescens Bistorta officinalis
##  Min.   :12            Min.   :5      Min.   :1.00     Min.   : 0.00       
##  1st Qu.:12            1st Qu.:5      1st Qu.:1.25     1st Qu.:13.00       
##  Median :12            Median :5      Median :2.50     Median :22.00       
##  Mean   :12            Mean   :5      Mean   :2.90     Mean   :18.51       
##  3rd Qu.:12            3rd Qu.:5      3rd Qu.:3.75     3rd Qu.:25.00       
##  Max.   :12            Max.   :5      Max.   :7.00     Max.   :25.00       
##  NA's   :164           NA's   :164    NA's   :155      NA's   :3           
##  Caltha palustris Carex rostrata Chamaenerion angustifolium Cirsium palustre
##  Min.   : 2.000   Min.   :25     Min.   : 1.000             Min.   : 1.00   
##  1st Qu.: 3.500   1st Qu.:25     1st Qu.: 2.250             1st Qu.: 5.00   
##  Median : 5.000   Median :25     Median : 5.500             Median : 9.00   
##  Mean   : 8.333   Mean   :25     Mean   : 8.667             Mean   :10.73   
##  3rd Qu.:11.500   3rd Qu.:25     3rd Qu.:11.750             3rd Qu.:17.50   
##  Max.   :18.000   Max.   :25     Max.   :25.000             Max.   :25.00   
##  NA's   :162      NA's   :164    NA's   :159                NA's   :98      
##  Comarum palustre Deschampsia cespitosa Digitalis purpurea
##  Min.   : 1.00    Min.   : 2.00         Mode:logical      
##  1st Qu.: 5.00    1st Qu.: 6.00         NA's:165          
##  Median :15.50    Median :11.00                           
##  Mean   :14.27    Mean   :11.71                           
##  3rd Qu.:22.00    3rd Qu.:18.00                           
##  Max.   :25.00    Max.   :25.00                           
##  NA's   :139      NA's   :110                             
##  Dryopteris carthusiana Epilobium palustre Equisetum fluviatile Frangula alnus 
##  Min.   : 1.000         Min.   :3.00       Min.   :4            Min.   :1.000  
##  1st Qu.: 2.500         1st Qu.:4.25       1st Qu.:4            1st Qu.:1.500  
##  Median : 4.000         Median :5.50       Median :4            Median :2.000  
##  Mean   : 5.815         Mean   :5.50       Mean   :4            Mean   :1.667  
##  3rd Qu.: 8.000         3rd Qu.:6.75       3rd Qu.:4            3rd Qu.:2.000  
##  Max.   :18.000         Max.   :8.00       Max.   :4            Max.   :2.000  
##  NA's   :138            NA's   :163        NA's   :164          NA's   :162    
##  Galeopsis tetrahit Galium aparine   Galium palustre Holcus lanatus 
##  Min.   :1.000      Min.   : 2.000   Min.   :8       Min.   : 2.00  
##  1st Qu.:1.000      1st Qu.: 3.000   1st Qu.:8       1st Qu.: 6.50  
##  Median :2.000      Median : 7.000   Median :8       Median : 8.00  
##  Mean   :2.143      Mean   : 9.947   Mean   :8       Mean   :10.36  
##  3rd Qu.:2.500      3rd Qu.:15.500   3rd Qu.:8       3rd Qu.:13.50  
##  Max.   :5.000      Max.   :25.000   Max.   :8       Max.   :25.00  
##  NA's   :158        NA's   :146      NA's   :164     NA's   :154    
##  Juncus acutiflorus Juncus effusus Lotus pedunculatus Lysimachia europaea
##  Min.   : 2.00      Min.   : 1.0   Min.   : 1.00      Min.   : 1.0       
##  1st Qu.:22.00      1st Qu.: 4.0   1st Qu.: 4.75      1st Qu.: 3.0       
##  Median :25.00      Median :11.0   Median :11.50      Median : 5.0       
##  Mean   :22.19      Mean   :13.2   Mean   :11.36      Mean   : 6.2       
##  3rd Qu.:25.00      3rd Qu.:25.0   3rd Qu.:16.25      3rd Qu.: 9.0       
##  Max.   :25.00      Max.   :25.0   Max.   :25.00      Max.   :15.0       
##  NA's   :30         NA's   :130    NA's   :137        NA's   :140        
##  Lysimachia vulgaris Molinia caerulea Potentilla erecta Pteridium aquilinum
##  Min.   : 2.00       Min.   : 1.00    Min.   : 1.000    Min.   :25         
##  1st Qu.: 7.00       1st Qu.:10.00    1st Qu.: 3.500    1st Qu.:25         
##  Median :11.00       Median :21.00    Median : 5.000    Median :25         
##  Mean   :12.45       Mean   :17.75    Mean   : 7.533    Mean   :25         
##  3rd Qu.:20.50       3rd Qu.:25.00    3rd Qu.: 9.000    3rd Qu.:25         
##  Max.   :24.00       Max.   :25.00    Max.   :23.000    Max.   :25         
##  NA's   :154         NA's   :109      NA's   :150       NA's   :161        
##  Quercus robur  Ranunculus acris Ranunculus repens Rhinanthus minor
##  Min.   :1.00   Min.   :2.00     Min.   :3.0       Min.   :1.00    
##  1st Qu.:1.25   1st Qu.:2.25     1st Qu.:3.5       1st Qu.:1.75    
##  Median :1.50   Median :2.50     Median :4.0       Median :2.50    
##  Mean   :1.50   Mean   :2.50     Mean   :4.0       Mean   :2.50    
##  3rd Qu.:1.75   3rd Qu.:2.75     3rd Qu.:4.5       3rd Qu.:3.25    
##  Max.   :2.00   Max.   :3.00     Max.   :5.0       Max.   :4.00    
##  NA's   :163    NA's   :163      NA's   :162       NA's   :161     
##  Rumex acetosa    Salix aurita    Senecio ovatus Silene flos-cuculi
##  Min.   : 1.00   Min.   : 1.000   Min.   :7      Min.   : 1.00     
##  1st Qu.: 6.25   1st Qu.: 3.500   1st Qu.:7      1st Qu.: 4.75     
##  Median :12.00   Median : 6.500   Median :7      Median : 6.50     
##  Mean   :10.97   Mean   : 8.083   Mean   :7      Mean   : 6.00     
##  3rd Qu.:15.00   3rd Qu.: 9.750   3rd Qu.:7      3rd Qu.: 7.75     
##  Max.   :25.00   Max.   :25.000   Max.   :7      Max.   :10.00     
##  NA's   :67      NA's   :153      NA's   :164    NA's   :161       
##  Stachys officinalis Stellaria alsine Stellaria graminea Stellaria holostea
##  Min.   :16.00       Min.   :1.00     Min.   :2.00       Min.   :7         
##  1st Qu.:18.25       1st Qu.:1.75     1st Qu.:3.50       1st Qu.:7         
##  Median :20.50       Median :2.50     Median :4.00       Median :7         
##  Mean   :20.50       Mean   :2.50     Mean   :4.25       Mean   :7         
##  3rd Qu.:22.75       3rd Qu.:3.25     3rd Qu.:4.75       3rd Qu.:7         
##  Max.   :25.00       Max.   :4.00     Max.   :7.00       Max.   :7         
##  NA's   :163         NA's   :163      NA's   :161        NA's   :164       
##  Succisa pratensis Urtica dioica    Vaccinium myrtillus Vaccinium uliginosum
##  Min.   : 2.00     Min.   : 1.000   Min.   : 7.00       Min.   :6           
##  1st Qu.: 7.50     1st Qu.: 1.000   1st Qu.: 9.75       1st Qu.:6           
##  Median : 9.50     Median : 3.000   Median :12.50       Median :6           
##  Mean   :11.00     Mean   : 4.727   Mean   :12.50       Mean   :6           
##  3rd Qu.:13.75     3rd Qu.: 5.000   3rd Qu.:15.25       3rd Qu.:6           
##  Max.   :23.00     Max.   :21.000   Max.   :18.00       Max.   :6           
##  NA's   :159       NA's   :154      NA's   :163         NA's   :164         
##  Valeriana officinalis Veronica chamaedrys  Vicia sepium Viola palustris 
##  Min.   : 2.00         Min.   : 2.000      Min.   :2     Min.   : 1.000  
##  1st Qu.: 7.00         1st Qu.: 6.500      1st Qu.:2     1st Qu.: 4.750  
##  Median :16.00         Median :11.000      Median :2     Median : 7.000  
##  Mean   :14.56         Mean   : 8.667      Mean   :2     Mean   : 8.396  
##  3rd Qu.:20.00         3rd Qu.:12.000      3rd Qu.:2     3rd Qu.:10.250  
##  Max.   :25.00         Max.   :13.000      Max.   :2     Max.   :24.000  
##  NA's   :140           NA's   :162         NA's   :164   NA's   :117
vegetation <- transform(vegetation, Site=factor(Site))
vegetation[is.na(vegetation)] <- 0 # N'ont pas été observé
vegetation_num <- vegetation[4:55]

NMDS: Distance matrix

set.seed(123)
bray_dist <- vegdist(vegetation_num, method = "bray")

nmds_veg <- metaMDS(bray_dist, k = 2, trymax = 100)
## Run 0 stress 0.1963054 
## Run 1 stress 0.2255869 
## Run 2 stress 0.2247842 
## Run 3 stress 0.1943476 
## ... New best solution
## ... Procrustes: rmse 0.01673462  max resid 0.2065978 
## Run 4 stress 0.196522 
## Run 5 stress 0.194374 
## ... Procrustes: rmse 0.002115744  max resid 0.02330085 
## Run 6 stress 0.2321771 
## Run 7 stress 0.1943537 
## ... Procrustes: rmse 0.0007199693  max resid 0.006693355 
## ... Similar to previous best
## Run 8 stress 0.197081 
## Run 9 stress 0.1944078 
## ... Procrustes: rmse 0.00468915  max resid 0.05835511 
## Run 10 stress 0.1943739 
## ... Procrustes: rmse 0.002131619  max resid 0.02331612 
## Run 11 stress 0.2270577 
## Run 12 stress 0.1944323 
## ... Procrustes: rmse 0.002656075  max resid 0.02890492 
## Run 13 stress 0.2233056 
## Run 14 stress 0.1944324 
## ... Procrustes: rmse 0.005093689  max resid 0.05838784 
## Run 15 stress 0.197229 
## Run 16 stress 0.194408 
## ... Procrustes: rmse 0.004684189  max resid 0.0583173 
## Run 17 stress 0.2330493 
## Run 18 stress 0.1963081 
## Run 19 stress 0.2228328 
## Run 20 stress 0.1963064 
## *** Best solution repeated 1 times
nmds_veg
## 
## Call:
## metaMDS(comm = bray_dist, k = 2, trymax = 100) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     bray_dist 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1943476 
## Stress type 1, weak ties
## Best solution was repeated 1 time in 20 tries
## The best solution was from try 3 (random start)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: scores missing
# k = 2 signifie que l'on reduit les donnees a 2 dimensions (pour une visualisation en 2D)
# trymax = 100 indique le nombre maximal d'essais pour trouver une solution stable

Score plot with habitats

nmds_point_veg <- as.data.frame(nmds_veg$points)
nmds_point_veg$Site <- vegetation$Site

# Visualiser avec ggplot2
ggplot(nmds_point_veg, aes(x = MDS1, y = MDS2, color = Site, fill = Site)) +
  geom_point(size = 1) +
  stat_ellipse(geom = "polygon", alpha = 0.3) +  
  theme_minimal() +
  labs(title = "Analyse nMDS avec indice de Bray-Curtis",
       x = "MDS1",
       y = "MDS2")

nmds_veg$stress # 0.1943476 ~ 0.2
## [1] 0.1943476
stressplot(nmds_veg) #stabilité de la nMDS