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
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
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
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 :
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