ANÁLISIS DE COMPONENTES PRINCIPALES (PCA) y ESCALAMIENTO MULTIDIMENSIONAL NO MÉTRICO (NMDS)
junio 24, 2022
Autor: Álvaro Alonso Fernández
Departamento de Ciencias de la Vida
Universidad de Alcalá (España)
1 Nuestros datos
Queremos realizar un ANÁLISIS DE COMPONENTES PRINCIPALES (PCA) para los datos ambientales recogidos en doce estaciones de muestreo (en cada columna se introduce una variable ambiental y en cada fila una estación de muestreo). Con un PCA reducimos las diferentes variables a unas pocas (componentes principales), permitiendo interpretar de forma más sencilla nuestros resultados. Es importante destacar que normalmente los primeros componentes principales recogen gran parte de la varianza y el resto (normalmente a partir del tercero) ya tienen poco peso a la hora de explicar la ordenación.
Entre los resultados más destacados de este análisis:
- Permite ordenar en el espacio nuestras estaciones o puntos de muestreo en función del conjunto de variables explicativas introducidas (por ejemplo las propiedades físico-químicas del agua o la composición de la comunidad).
- Los “eigenvectores” nos indican la importancia de la variable correspondiente (longitud del vector)
- Los “eigenvalues” nos indican la importancia de cada componente (porcentaje de la varianza explicada)
- Nuestros datos deben cumplir la normalidad (normalidad multivariada)
Veamos como son nuestros datos:
#Fijamos el directorio de trabajo
setwd(dir = "F:/R/MARKDOWN/pca")
#leemos el fichero de datos indicando que la primera columna no es una variable row.names
datospca<-read.csv("pca.csv", sep=";",row.names = 1)
str(datospca)## 'data.frame': 12 obs. of 6 variables:
## $ pH : num 7.2 7.6 7.8 7.1 7 8.1 8.2 8.3 8.5 8.4 ...
## $ distance : int 10 12 16 19 21 22 26 28 29 32 ...
## $ temperature: int 10 12 14 12 11 10 9 15 18 16 ...
## $ altitude : int 1000 1110 1121 1212 1001 1001 988 956 945 985 ...
## $ growth : int 10 11 12 13 15 14 16 19 18 19 ...
## $ diversity : int 10 12 15 15 16 10 21 22 24 25 ...head(datospca)## pH distance temperature altitude growth diversity
## Alto 7.2 10 10 1000 10 10
## Bajo 7.6 12 12 1110 11 12
## Intermedio 7.8 16 14 1121 12 15
## Muy alto 7.1 19 12 1212 13 15
## Medio alto 7.0 21 11 1001 15 16
## Alto bajo 8.1 22 10 1001 14 102 Normalidad de nuestros datos
Los datos deben cumplir la normalidad multiparamétrica, para ello existe el paquete MVN. También han desarrollado una web para testarlo on line.
Vamos a optar en primer lugar por la normalidad variable a variable:
par(mfrow=c(2,3))
qqnorm(datospca[,1])
qqline(datospca[,1])
qqnorm(datospca[,2])
qqline(datospca[,2])
qqnorm(datospca[,3])
qqline(datospca[,3])
qqnorm(datospca[,4])
qqline(datospca[,4])
qqnorm(datospca[,5])
qqline(datospca[,5])
qqnorm(datospca[,6])
qqline(datospca[,6])
Realizamnos un test de normalidad univariante (Shapiro-Wilk):
shapiro.test(datospca[,1])##
## Shapiro-Wilk normality test
##
## data: datospca[, 1]
## W = 0.89506, p-value = 0.137shapiro.test(datospca[,2])##
## Shapiro-Wilk normality test
##
## data: datospca[, 2]
## W = 0.93861, p-value = 0.4802shapiro.test(datospca[,2])##
## Shapiro-Wilk normality test
##
## data: datospca[, 2]
## W = 0.93861, p-value = 0.4802shapiro.test(datospca[,3])##
## Shapiro-Wilk normality test
##
## data: datospca[, 3]
## W = 0.91301, p-value = 0.2331shapiro.test(datospca[,3])##
## Shapiro-Wilk normality test
##
## data: datospca[, 3]
## W = 0.91301, p-value = 0.2331shapiro.test(datospca[,4])##
## Shapiro-Wilk normality test
##
## data: datospca[, 4]
## W = 0.88528, p-value = 0.1025En este caso todas nuestras variables son normales.
Testamos la normalidad multivariante (para MVN requiere https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/). Hay varios test multivariantes como “mardia”, etc. revisar ventajas e inconveninetes de cada uno de ellos antes de aplicarlos:
library(MVN)
mnvdemisdatospca <- mvn(datospca, mvnTest = "mardia")
mnvdemisdatospca## $multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 59.6116454400866 0.345715909912825 YES
## 2 Mardia Kurtosis -1.2499443330085 0.21131988309138 YES
## 3 MVN <NA> <NA> YES
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling pH 0.5229 0.1447 YES
## 2 Anderson-Darling distance 0.2816 0.5726 YES
## 3 Anderson-Darling temperature 0.3843 0.3355 YES
## 4 Anderson-Darling altitude 0.6523 0.0664 YES
## 5 Anderson-Darling growth 0.3975 0.3104 YES
## 6 Anderson-Darling diversity 0.3924 0.3199 YES
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th
## pH 12 7.916667 0.5654175 8.15 7 8.6 7.50 8.325
## distance 12 23.166667 7.7440104 24.00 10 33.0 18.25 29.250
## temperature 12 13.750000 3.6212755 13.00 9 19.0 10.75 16.500
## altitude 12 1011.750000 91.1553370 994.00 901 1212.0 953.25 1028.250
## growth 12 15.500000 3.5032452 15.50 10 20.0 12.75 19.000
## diversity 12 18.083333 5.7911898 18.50 10 26.0 14.25 22.500
## Skew Kurtosis
## pH -0.45081721 -1.4924827
## distance -0.32907432 -1.4254322
## temperature 0.23755956 -1.6168728
## altitude 0.84377740 -0.4746792
## growth -0.16281200 -1.6617798
## diversity -0.08805807 -1.67160013 Realizamos el PCA
Veamos a realizar el PCA, para ello cargamos los paquetes correspondientes (en este caso factoextra) y ejecutamos el comando prcom:
library(factoextra)## Loading required package: ggplot2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBapruebapca <- prcomp(datospca, scale = TRUE)3.1 Representación Eigenvalues
Representamos por columnas los eigenvalues para cada una de las dimensiones:
fviz_eig(pruebapca)
Podemos ver como la primera dimensión recoge casi el 80% de la varianza explicada
3.2 Representación de las dos primeras dimensiones
Representamos por columnas los eigenvalues para cada una de las dimensiones:
fviz_pca_ind(pruebapca,
col.ind = "cos2", # el color va en funci?n de la calidad de la representaci?n
gradient.cols = c("#800000", "#FFA500", "#008000"),
repel = TRUE # con esto se evita el solapado del texto
)
3.3 Representación de las dos primeras dimensiones y peso variables
Representamos por columnas los eigenvalues para cada una de las dimensiones con vectores para las variables:
fviz_pca_var(pruebapca,
col.var = "contrib", # el color va en funci?n de la contribuci?n al PC
gradient.cols = c("#800000", "#FFA500", "#008000"),
repel = TRUE # con esto se evita el solapado del texto
)
3.4 Representación de las dos primeras dimensiones, peso variables y estaciones
Representamos por columnas los eigenvalues para cada una de las dimensiones con vectores para las variables:
fviz_pca_biplot(pruebapca, repel = TRUE,
col.var = "#D2691E", # color para las variables
col.ind = "#000000" # colores de las estaciones
)
3.5 Representación gráfica de la importancia de cada variable
library("corrplot")## corrplot 0.92 loadedvar <- get_pca_var(pruebapca)
corrplot(var$cos2, is.corr = FALSE)
3.6 Datos numéricos del PCA
Para ello:
# Los Eigenvalues
eigvalues <- get_eigenvalue(pruebapca)
eigvalues## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 4.65938232 77.6563719 77.65637
## Dim.2 0.55662687 9.2771144 86.93349
## Dim.3 0.34691750 5.7819584 92.71544
## Dim.4 0.29646558 4.9410930 97.65654
## Dim.5 0.12429862 2.0716437 99.72818
## Dim.6 0.01630911 0.2718185 100.00000# Obtenemos los resultados de las variables utilizadas
resultsvariables <- get_pca_var(pruebapca)
resultsvariables$coord # Las coordenadas de las variables## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## pH 0.8613955 0.12299813 -0.1318622 0.47478978 0.002652093
## distance 0.9441953 -0.02226187 0.2527031 -0.02933266 0.195353927
## temperature 0.8034514 -0.41138152 -0.4065948 -0.13165566 0.050469355
## altitude -0.7666448 -0.57623010 0.1853844 0.20753104 0.050777368
## growth 0.9710016 -0.01322403 0.1812505 -0.09548320 0.068312723
## diversity 0.9218204 -0.19887679 0.1820235 -0.02569934 -0.276290257
## Dim.6
## pH -0.007016895
## distance 0.071535256
## temperature 0.005599322
## altitude -0.014135054
## growth -0.101712909
## diversity 0.023788936resultsvariables$contrib # Las contribuciones al PCA## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
## pH 15.92491 2.71789593 5.012035 76.0376072 0.005658629 0.3018976
## distance 19.13354 0.08903463 18.407511 0.2902208 30.702799842 31.3768959
## temperature 13.85450 30.40362642 47.653792 5.8466191 2.049222908 0.1922386
## altitude 12.61421 59.65237268 9.906494 14.5275320 2.074311934 1.2250805
## growth 20.23539 0.03141691 9.469615 3.0752443 3.754368366 63.4339650
## diversity 18.23746 7.10565343 9.550553 0.2227766 61.413638321 3.4699223resultsvariables$cos2 # Calidad de la representaci?n ## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## pH 0.7420022 0.0151285389 0.01738763 0.2254253349 7.033598e-06
## distance 0.8915047 0.0004955907 0.06385888 0.0008604047 3.816316e-02
## temperature 0.6455342 0.1692347530 0.16531934 0.0173332132 2.547156e-03
## altitude 0.5877442 0.3320411327 0.03436736 0.0430691324 2.578341e-03
## growth 0.9428442 0.0001748750 0.03285175 0.0091170410 4.666628e-03
## diversity 0.8497528 0.0395519760 0.03313254 0.0006604559 7.633631e-02
## Dim.6
## pH 4.923682e-05
## distance 5.117293e-03
## temperature 3.135241e-05
## altitude 1.997997e-04
## growth 1.034552e-02
## diversity 5.659135e-04# Resultados para cada una de las estaciones
resultsstations <- get_pca_ind(pruebapca)
resultsstations$coord # Bis de los anteriores## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Alto -2.8915052 0.91222143 -0.68590352 -0.47119937 -0.32880095
## Bajo -2.4433505 -0.31252282 -0.54035159 0.39181831 -0.14284698
## Intermedio -1.5639105 -0.81050071 -0.49327306 0.51037847 -0.10818946
## Muy alto -2.3202275 -1.49772327 0.73340522 -0.12605611 0.21718171
## Medio alto -1.2279571 0.34931076 0.57492723 -1.21794398 -0.03865547
## Alto bajo -1.0687580 1.09968843 -0.02260747 0.63726607 0.76461636
## Regular 0.2435975 0.85940498 1.06783280 0.58630839 -0.38582532
## Perdido 1.6275012 0.17718852 0.20166889 0.03386286 -0.02415986
## Encontrado 2.1954366 -0.21894332 -0.41309698 0.12279750 -0.17468594
## Casi 2.1350386 -0.34513384 0.45328105 0.19029430 -0.05723272
## Nulo 2.3131184 -0.13771531 -0.41844880 -0.56042374 0.57041415
## Acierto 3.0010165 -0.07527485 -0.45743376 -0.09710269 -0.29181551
## Dim.6
## Alto 0.07652281
## Bajo -0.09004704
## Intermedio 0.05985942
## Muy alto -0.01717443
## Medio alto -0.04122504
## Alto bajo -0.05353939
## Regular 0.12888077
## Perdido -0.27454619
## Encontrado 0.11971020
## Casi 0.07846180
## Nulo 0.15560141
## Acierto -0.14250431resultsstations$contrib # Bis de los anteriores## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Alto 14.9533449 12.45819518 11.30104515 6.24099567 7.24800662
## Bajo 10.6773123 1.46223797 7.01367125 4.31532908 1.36802673
## Intermedio 4.3743567 9.83470614 5.84476884 7.32199050 0.78473371
## Muy alto 9.6283427 33.58286502 12.92053617 0.44665508 3.16227003
## Medio alto 2.6968501 1.82674750 7.93995982 41.69645142 0.10017848
## Alto bajo 2.0429096 18.10480697 0.01227712 11.41526688 39.19580723
## Regular 0.1061295 11.05733129 27.39044297 9.66265974 9.98007711
## Perdido 4.7373233 0.47003017 0.97694554 0.03223232 0.03913283
## Encontrado 8.6204947 0.71765875 4.09917616 0.42386106 2.04582440
## Casi 8.1527073 1.78332166 4.93546028 1.01787874 0.21960447
## Nulo 9.5694313 0.28393507 4.20607682 8.82830885 21.81388482
## Acierto 16.1074643 0.08483096 5.02630655 0.26503732 5.70912023
## Dim.6
## Alto 2.9920601
## Bajo 4.1431189
## Intermedio 1.8308532
## Muy alto 0.1507138
## Medio alto 0.8683817
## Alto bajo 1.4646551
## Regular 8.4872049
## Perdido 38.5140615
## Encontrado 7.3223548
## Casi 3.1456105
## Nulo 12.3713056
## Acierto 10.3763466resultsstations$cos2 # Bis de los anteriores ## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Alto 0.83612968 0.0832197174 0.0470491482 0.0222041987 0.0108116465
## Bajo 0.91260830 0.0149305685 0.0446339930 0.0234683374 0.0031192847
## Intermedio 0.67530027 0.1813760459 0.0671812267 0.0719213483 0.0032317890
## Muy alto 0.65429542 0.2726314189 0.0653733506 0.0019312592 0.0057327006
## Medio alto 0.43744422 0.0353981219 0.0958919283 0.4303392042 0.0004334892
## Alto bajo 0.34140830 0.3614553387 0.0001527634 0.1213827386 0.1747440939
## Regular 0.02424592 0.3017786791 0.4659071077 0.1404575328 0.0608239018
## Perdido 0.94668487 0.0112210621 0.0145358566 0.0004098359 0.0002086179
## Encontrado 0.94537357 0.0094021044 0.0334707710 0.0029576087 0.0059851851
## Casi 0.92488257 0.0241685306 0.0416879244 0.0073472822 0.0006646054
## Nulo 0.86184108 0.0030548924 0.0282043523 0.0505899811 0.0524097448
## Acierto 0.96467352 0.0006069366 0.0224130014 0.0010099646 0.0091213755
## Dim.6
## Alto 5.856087e-04
## Bajo 1.239515e-03
## Intermedio 9.893232e-04
## Muy alto 3.584901e-05
## Medio alto 4.930358e-04
## Alto bajo 8.567660e-04
## Regular 6.786862e-03
## Perdido 2.693976e-02
## Encontrado 2.810761e-03
## Casi 1.249084e-03
## Nulo 3.899946e-03
## Acierto 2.175202e-033.7 Añadir más datos al PCA
Se puede añadir otra tanda de datos al PCA, en este caso 6 nuevas estaciones con las mismas variables ambientales:
pca_news <- read.csv("pca_news.csv", sep=";")
pca_news## pH distance temperature altitude growth diversity
## 1 6.9 8 21.0 425 22 12
## 2 7.1 7 19.0 456 21 15
## 3 7.0 6 18.0 564 21 16
## 4 6.8 5 17.5 690 20 19
## 5 6.9 9 17.9 578 18 20
## 6 6.9 10 18.9 594 19 19pruebapca_news <- predict(pruebapca, newdata = pca_news)
pruebapca_news## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] 1.8431420 3.876678 -3.599409 -4.675688 -0.5564863 -1.8714663
## [2,] 1.6942769 3.847688 -3.173456 -4.071518 -1.1167308 -1.7012468
## [3,] 1.1173740 3.013547 -2.572364 -3.708579 -1.1938257 -1.8749412
## [4,] 0.4702312 1.834606 -1.945953 -3.424225 -1.5499967 -1.7630414
## [5,] 1.0611739 2.700498 -2.349225 -3.700882 -1.6695437 -0.8557279
## [6,] 1.2127686 2.449783 -2.394728 -3.749617 -1.3425509 -1.0502277Ahora hacemos el ploteo de las 12 iniciales más las 6 nuevas:
# Ploteo de las estaciones las 12 iniciales
ploteo <- fviz_pca_ind(pruebapca, repel = TRUE)
# A lo anterior se a?aden las 6 nuevas estaciones
fviz_add(ploteo, pruebapca_news, color ="red")
3.8 PCA con estaciones por grupos
Ahora en nuestra matriz tenemos una columna categórica (A, B y C):
pca_groups <- read.csv("pca_groups.csv", sep=";")
pca_groups## pH distance temperature altitude growth diversity groups
## 1 7.2 10 10.0 1000 10 10 A
## 2 7.6 12 12.0 1110 11 12 A
## 3 7.8 16 14.0 1121 12 15 A
## 4 7.1 19 12.0 1212 13 15 A
## 5 7.0 21 11.0 1001 15 16 A
## 6 8.1 22 10.0 1001 14 10 B
## 7 8.2 26 9.0 988 16 21 B
## 8 8.3 28 15.0 956 19 22 B
## 9 8.5 29 18.0 945 18 24 B
## 10 8.4 32 16.0 985 19 25 B
## 11 8.2 33 19.0 921 19 21 B
## 12 8.6 30 19.0 901 20 26 B
## 13 6.9 8 21.0 425 22 12 C
## 14 7.1 7 19.0 456 21 15 C
## 15 7.0 6 18.0 564 21 16 C
## 16 6.8 5 17.5 690 20 19 C
## 17 6.9 9 17.9 578 18 20 C
## 18 6.9 10 18.9 594 19 19 CAhora seleccionamos las variables para el PCA (18) de la matriz en función de las tres estaciones:
groups <- as.factor(pca_groups$groups[1:18])
groups## [1] A A A A A B B B B B B B C C C C C C
## Levels: A B CSeleccionamos las estaciones y variables activas (la última columna es categórica)
pca_groups <- pca_groups[1:18, 1:6]
pca_groups## pH distance temperature altitude growth diversity
## 1 7.2 10 10.0 1000 10 10
## 2 7.6 12 12.0 1110 11 12
## 3 7.8 16 14.0 1121 12 15
## 4 7.1 19 12.0 1212 13 15
## 5 7.0 21 11.0 1001 15 16
## 6 8.1 22 10.0 1001 14 10
## 7 8.2 26 9.0 988 16 21
## 8 8.3 28 15.0 956 19 22
## 9 8.5 29 18.0 945 18 24
## 10 8.4 32 16.0 985 19 25
## 11 8.2 33 19.0 921 19 21
## 12 8.6 30 19.0 901 20 26
## 13 6.9 8 21.0 425 22 12
## 14 7.1 7 19.0 456 21 15
## 15 7.0 6 18.0 564 21 16
## 16 6.8 5 17.5 690 20 19
## 17 6.9 9 17.9 578 18 20
## 18 6.9 10 18.9 594 19 19pca_groups_graph <- prcomp(pca_groups, scale = TRUE)
pca_groups_graph## Standard deviations (1, .., p=6):
## [1] 1.7045225 1.5677318 0.5302477 0.4473621 0.3381012 0.2030071
##
## Rotation (n x k) = (6 x 6):
## PC1 PC2 PC3 PC4 PC5
## pH 0.3871287 -0.42178235 -0.5592271 0.02754150 -0.5785263
## distance 0.4024694 -0.43497817 -0.1705504 -0.12705540 0.6196472
## temperature -0.4426884 -0.34555648 -0.1550099 0.80075517 0.1311322
## altitude 0.5631818 0.06188552 0.2308650 0.37123449 0.2958464
## growth -0.4164644 -0.41643115 -0.1365669 -0.45005163 0.2931689
## diversity 0.0395024 -0.57986174 0.7498002 -0.03908995 -0.3011243
## PC6
## pH 0.15506402
## distance -0.46863539
## temperature -0.04674542
## altitude 0.63273825
## growth 0.58821552
## diversity -0.08829320Representamos la figura con las 18 estaciones agrupadas en las tres categorías:
fviz_pca_ind(pca_groups_graph,
col.ind = groups, # Cada grupo con un color
palette = c("#00AFBB", "#000000", "#FC4E07"),
addEllipses = TRUE, # elpises para cada grupo
ellipse.type = "confidence",
legend.title = "Groups",
repel = TRUE
)
Ahora añadimos los vectores con las variables en la figura anterior:
fviz_pca_biplot(pca_groups_graph,
col.ind = groups, # Cada grupo con un color
palette = c("#00AFBB", "#000000", "#FC4E07"),
addEllipses = TRUE, # elpises para cada grupo
ellipse.type = "confidence",
legend.title = "Groups",
repel = TRUE
)
4 Alternativa no paramétrica: NMDS (Escalamiento Multidimensional No Métrico)
No requiere de normalidad en los datos. Es un método multivariante adecuado para ordenar espacialmente estaciones en función de datos de comunidades o datos ambientales.
Vamos a utilizar los mismos datos de nuestro ejemplo anterior:
#Fijamos el directorio de trabajo
setwd(dir = "F:/R/MARKDOWN/pca")
#leemos el fichero de datos
datosnmds<-read.csv("pca.csv", sep=";",row.names = 1)
str(datosnmds)## 'data.frame': 12 obs. of 6 variables:
## $ pH : num 7.2 7.6 7.8 7.1 7 8.1 8.2 8.3 8.5 8.4 ...
## $ distance : int 10 12 16 19 21 22 26 28 29 32 ...
## $ temperature: int 10 12 14 12 11 10 9 15 18 16 ...
## $ altitude : int 1000 1110 1121 1212 1001 1001 988 956 945 985 ...
## $ growth : int 10 11 12 13 15 14 16 19 18 19 ...
## $ diversity : int 10 12 15 15 16 10 21 22 24 25 ...Cargamos el paquete necesario:
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-2Hacemos la ordenación espacial:
analizamosnmds <- metaMDS(datosnmds, k=2, trymax=100) #con k seleccionamos el n de dimensiones y con trymax incrementamos el n de iteraciones## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.02480033
## Run 1 stress 0.02480033
## ... Procrustes: rmse 6.033153e-06 max resid 1.136791e-05
## ... Similar to previous best
## Run 2 stress 0.02480041
## ... Procrustes: rmse 0.0002167966 max resid 0.0004851016
## ... Similar to previous best
## Run 3 stress 0.02626526
## Run 4 stress 0.0757909
## Run 5 stress 0.02626526
## Run 6 stress 0.02480034
## ... Procrustes: rmse 4.781706e-05 max resid 9.199704e-05
## ... Similar to previous best
## Run 7 stress 0.07579089
## Run 8 stress 0.02626527
## Run 9 stress 0.02626526
## Run 10 stress 0.07579089
## Run 11 stress 0.2518593
## Run 12 stress 0.02626526
## Run 13 stress 0.07579089
## Run 14 stress 0.02480032
## ... New best solution
## ... Procrustes: rmse 0.0001038054 max resid 0.0002248403
## ... Similar to previous best
## Run 15 stress 0.02480037
## ... Procrustes: rmse 8.604587e-05 max resid 0.0001790521
## ... Similar to previous best
## Run 16 stress 0.02480041
## ... Procrustes: rmse 0.0001202861 max resid 0.0002711023
## ... Similar to previous best
## Run 17 stress 0.02480034
## ... Procrustes: rmse 9.494851e-05 max resid 0.0002090252
## ... Similar to previous best
## Run 18 stress 0.02626526
## Run 19 stress 0.02626526
## Run 20 stress 0.02480034
## ... Procrustes: rmse 5.098908e-05 max resid 0.0001162622
## ... Similar to previous best
## *** Solution reachedstressplot(analizamosnmds) #Dibuja un Shepard diagram que es un plot de las distancias de la ordenación y el ajuste frente a las disimilitudes originales. Si los puntos se separan mucho de la línea no es bueno
analizamosnmds##
## Call:
## metaMDS(comm = datosnmds, k = 2, trymax = 100)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(datosnmds))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.02480032
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(datosnmds))'Por defecto se aplica una transformación de raíz cuadrada y se calculan las distancias de Bray-Curtis (adecuado para datos ecológicos).
Una buena regla para el stress: stress < 0.05 muy buena representación, con < 0.1 buena, < 0.2 es adecuada, y < 0.3 la representación es un poco pobre. A mayor stress peor representación.
Y la representamos:
plot(analizamosnmds)
Se pueden apreciar las 12 estaciones representadas por los puntos negros, y las seis variables ambientales representadas por las cruces rojas para los dos primeros ejes.
Se puede mejorar la representación con el nombre de las variables:
plot(analizamosnmds, type="t")
Podemos complicar la representación:
ordiplot(analizamosnmds,type="n")
orditorp(analizamosnmds,display="species",col="blue",air=0.01)
orditorp(analizamosnmds,display="sites",cex=1.2,air=0.01)
Ahora por ejemplo nuestras 6 primeras estaciones corresponden a un río y las otras 6 a otro y queremos que esto aparezca en la representación:
treat=c(rep("River1",6),rep("River2",6))
ordiplot(analizamosnmds,type="n")
ordihull(analizamosnmds,groups=treat,draw="polygon",col="grey90",label=F)
orditorp(analizamosnmds,display="species",col="red",air=0.05)
orditorp(analizamosnmds,display="sites",col=c(rep("red",6),rep("black",6)),
air=0.05,cex=1.2)
Lo mismo con gráfico de arañas:
ordiplot(analizamosnmds,type="n")
ordispider(analizamosnmds,groups=treat,col="black",label=F)
orditorp(analizamosnmds,display="species",col="red",air=0.05)
orditorp(analizamosnmds,display="sites",col=c(rep("red",6),rep("black",6)),
air=0.05,cex=1.2)
Con gráfico de elipses:
ordiplot(analizamosnmds,type="n")
ordiellipse(analizamosnmds,groups=treat,col="black",label=F)
orditorp(analizamosnmds,display="species",col="red",air=0.05)
orditorp(analizamosnmds,display="sites",col=c(rep("red",6),rep("black",6)),
air=0.05,cex=1.2)
5 CRÉDITOS
Departamento de Ciencias de la Vida
Universidad de Alcalá (España)
Truchas en un tramo de cabecera de un río