Introducción
Las técnicas de ordenación más utilizadas son las de tipo métrico basadas en el método de los Eigenvectores, dentro de estas el análisis de componentes principales (PCA) y el análisis de correspondencias (CA) son los más utilizados.
El PCA está basado en la distancia euclídea, mientras que el CA utiliza la como distancia el ?? cuadrado, el análisis de coordenadas principales puede usar cualquier tipo de medida de distancia, pero cuando utiliza la euclídea su resultado es idéntico al del PCA.
Esta técnica rota los datos en el espacio vectorial, de tal modo que la mayor parte de las variaciones se acumula en los primeros ejes.
Con frecuencia, el análisis de componentes principales revela relaciones de las que no se sospechaba inicialmente, y por tanto este análisis permite interpretaciones de los datos que no podrían ser derivadas directamente de las variables originales.
Aunque se necesitan las p componentes principales para reproducir toda la variabilidad del sistema, generalmente la mayor parte de esa variabilidad es explicada por un número pequeño k de componentes principales. En estos casos las k primeras componentes principales reemplazan las p variables originales, logrando una reducción del sistema original.
El desarrollo del procedimiento de componentes principales no requiere del supuesto de la normalidad multivariada.
Que debemos saber?
data(iris)
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
str(iris)
## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
pairs(iris[,1:4])
Análisis exploratorios
tapply(iris$Sepal.Length, iris$Species, mean)
## setosa versicolor virginica
## 5.006 5.936 6.588
tapply(iris$Sepal.Length, iris$Species, sd)
## setosa versicolor virginica
## 0.3524897 0.5161711 0.6358796
O mediante una selección precisa.
by(iris[,1:4] , iris$Species , colMeans)
## iris$Species: setosa
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.006 3.428 1.462 0.246
## --------------------------------------------------------
## iris$Species: versicolor
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.936 2.770 4.260 1.326
## --------------------------------------------------------
## iris$Species: virginica
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 6.588 2.974 5.552 2.026
Visualización de diagrama de cajas
boxplot(Sepal.Length ~ Species, data=iris)
o mostrarlo de una forma mas elegante con “ggplot2”
library(ggplot2)
ggplot(iris, aes(x = Species, y = Sepal.Length)) +
geom_boxplot(fill = "grey80", colour = "gray20") +
scale_x_discrete() + xlab("Tratamiento") +
ylab("Sepal.Length (cm)")
pca1 <- prcomp(iris[, 1:4], scale = T)
summary(pca1)
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion 0.7296 0.9581 0.99482 1.00000
En otras palabras, cada eje explica una cierta proporción de varianza. El primer eje explica 72.96%, el segundo 22. 85%, y así sucesivamente segun indica la “Proportion of Variance”.
pca1$rotation[, 1:2]
## PC1 PC2
## Sepal.Length 0.5210659 -0.37741762
## Sepal.Width -0.2693474 -0.92329566
## Petal.Length 0.5804131 -0.02449161
## Petal.Width 0.5648565 -0.06694199
#Esto es la interpretacion de los componentes. Notamos que el Petal.Width,y Sepal.Lengt, está correlacionado con el PC1. Mientras que el segundo componente está relacionado de forma negativa Sepal.Width.
Es conveniente dibujar los componentes del PCA
biplot(pca1, scale=0)
reg1Petal.Width<-lm(iris$Petal.Width ~ pca1$x[, 1:2])
summary(reg1Petal.Width)
##
## Call:
## lm(formula = iris$Petal.Width ~ pca1$x[, 1:2])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.46380 -0.12686 -0.00901 0.11421 0.59337
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.199333 0.015940 75.24 < 2e-16 ***
## pca1$x[, 1:2]PC1 0.430555 0.009362 45.99 < 2e-16 ***
## pca1$x[, 1:2]PC2 -0.051026 0.016729 -3.05 0.00271 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1952 on 147 degrees of freedom
## Multiple R-squared: 0.9353, Adjusted R-squared: 0.9344
## F-statistic: 1062 on 2 and 147 DF, p-value: < 2.2e-16
Vemos que el modelo explica un 93.44% de la variabilidad. Y ambos ejes son significativos, pero la primera variable (PC1) es positivo.
Es importante revisar los supuestos del modelo.
par(mfrow = c(2, 2))
plot(reg1Petal.Width)
Normalidad
shapiro.test(residuals(reg1Petal.Width))
##
## Shapiro-Wilk normality test
##
## data: residuals(reg1Petal.Width)
## W = 0.98946, p-value = 0.3215
Heterocedasticidad
library(car)
ncvTest(reg1Petal.Width)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 39.23071 Df = 1 p = 3.765705e-10
#Vemos que falla el supuesto
Aplicacmos un modelo glm (no requiere de linealidad)
glm1Petal.Width<- glm(iris$Petal.Width ~ pca1$x[, 1:2], family = poisson)
summary(glm1Petal.Width)
##
## Call:
## glm(formula = iris$Petal.Width ~ pca1$x[, 1:2], family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.96302 -0.23760 0.01564 0.14754 0.50619
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.07314 0.09636 -0.759 0.448
## pca1$x[, 1:2]PC1 0.44898 0.05949 7.547 4.45e-14 ***
## pca1$x[, 1:2]PC2 0.06118 0.08392 0.729 0.466
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 85.907 on 149 degrees of freedom
## Residual deviance: 10.350 on 147 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 4
Ahora solo una variable es signi???cativas.
Vamos a ver si esta vez los residuos son adecuados
par(mfrow = c(2, 2))
plot(glm1Petal.Width)
shapiro.test(residuals(glm1Petal.Width))
##
## Shapiro-Wilk normality test
##
## data: residuals(glm1Petal.Width)
## W = 0.97972, p-value = 0.02581
library(devtools)
install_github("vqv/ggbiplot")
## Skipping install of 'ggbiplot' from a github remote, the SHA1 (7325e880) has not changed since last install.
## Use `force = TRUE` to force installation
ir.species <- iris[, 5]
library(ggbiplot)
## Loading required package: plyr
## Loading required package: scales
## Loading required package: grid
g <- ggbiplot(pca1, obs.scale = 1, var.scale = 1,
groups = ir.species, ellipse = TRUE,
circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal',
legend.position = 'top')
print(g)
plot(pca1)
Es una técnica multivariante de interdependencia que trata de representar en un espacio geométrico de pocas dimensiones a las proximidades existentes entre un conjunto de objetos.
*El NMDS es un m´etodo de ordenación adecuado para datos que no son normales o que están en una escala discontinua o arbitraria
Una ventaja del NMDS frente a otras técnicas de ordenación es que, al estar basada en rangos de distancias, tiende a linealizar la relación entre las distancias ambientales y las distancias biológicas. Una de las desventajas de esta técnica es la di???cultad para alcanzar una solución estable única.
La ventaja del NMDS es que nos permite, al igual que el PCA, reducir la dimensionalidad de nuestros datos originales. El resultado de la ordenación se puede visualizar en un grá???co de ordenación. Posteriormente podemos relacionar los ejes resultantes de dicha ordenación con distintas variables ambientales para determinar de manera indirecta el efecto de éstas sobre la matriz de sitios x especies.
Aunque en ecología se utiliza típicamente esta técnica para analizar datos de comunidades biológicas (matriz de sitios x especies) también se puede aplicar a otro tipo de datos, como por ejemplo múltiples variables físico-químicas medidas en distintos cuerpos de agua (ríos, embalses, pantanos).
bio<-read.csv("bio.csv", h=T, sep = ",")
head(bio)
## X Alnusacuminata Anonaceaesp. Arbutusxalapensis Buddleiacordata
## 1 Bazom1 0 0 0 0
## 2 Bazom2 0 0 8 0
## 3 Bazom3 0 0 0 0
## 4 Bazom4 0 0 1 0
## 5 Bazom5 3 0 0 0
## 6 Bazom6 0 0 0 0
## Buddleianitida Calyptranthespallens Clethramacrophylla Clethraoleoides
## 1 0 0 25 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## Clethrasuaveolens Cleyeratheaeoides Clusiarosea Cornusdisciflora
## 1 0 18 0 0
## 2 0 0 0 0
## 3 0 8 0 0
## 4 0 0 0 0
## 5 0 3 0 15
## 6 0 19 0 9
## Cornusexcelsa Cupaniadentata Cupressuslusitanica Cyatheafulva
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## Daphnopsisselerolum Dendropanaxarboreus Deppeagrandiflora
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## Drimysgranadensis Erythrinachiapasana Eugeniacapuli Eugeniacapulioides
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## Eupatoriumligustrinum Eupatoriumnubigenum Garryalaurifolia
## 1 0 10 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 4 0
## Holodiscusargenteus Ingasp. Juniperusgamboana Lippiasubstrigosa
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## Liquidambarstyraciflua Litseaglaucescens Magnoliasharpii Meliosmadives
## 1 0 0 15 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 10 0
## Miconiaglaberrima Microtropiscontracta Myricacerifera Myrtaceaesp1.
## 1 0 5 0 0
## 2 0 0 0 0
## 3 0 2 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## Nyssasilvatica Ocoteaeffusa Ocoteahelicterifolia Olmediellabetschleriana
## 1 0 0 0 2
## 2 0 0 0 0
## 3 0 0 0 2
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## Oreopanaxxalapensis Ostryavirginiana Parathesisbelizensis
## 1 5 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 2 0 0
## Perseaamericana Pinusayacahuite Pinusmontezumae Pinusoocarpa
## 1 14 0 0 0
## 2 0 0 19 0
## 3 2 0 0 0
## 4 0 0 12 0
## 5 2 0 0 0
## 6 26 0 0 0
## Pinuspseudostrobus Pinustecunumani Prunuslundelliana Prunusrhamnoides
## 1 0 0 0 0
## 2 37 13 0 0
## 3 5 5 0 2
## 4 1 30 0 0
## 5 0 32 0 0
## 6 0 0 0 0
## Prunusserotina Quercusacatenangensis Quercuscandicans Quercuscrassifolia
## 1 0 0 0 0
## 2 0 0 0 0
## 3 1 71 0 2
## 4 0 0 0 0
## 5 1 0 0 9
## 6 2 0 0 3
## Quercuscrispipilis Quercuslaurina Quercuspeduncularis Quercuspolymorpha
## 1 0 6 0 0
## 2 0 0 0 0
## 3 0 19 0 0
## 4 0 10 0 0
## 5 0 30 0 0
## 6 0 13 0 0
## Quercusrugosa Quercussapotaefolia Quercussegoviensis Quercusskutchii
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## Rapaneajuergensenii Rapaneamyricoides Rhamnussharpii
## 1 9 0 0
## 2 0 0 0
## 3 0 0 4
## 4 0 0 0
## 5 26 0 0
## 6 48 0 2
## Rondeletiastenosiphon Saurauialatipetala Saurauiascabrida
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## Sebastianiacruenta Solanumnudum Styraxargenteus Symplococarponpurpusii
## 1 0 0 4 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 1 0
## 6 0 0 0 0
## Symplocoslimoncillo Tecomastans Ternstroemialineata Ternstroemiaoocarpa
## 1 0 0 2 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 0 0 0 0
## 5 0 0 10 0
## 6 0 0 25 0
## Verbesinaperymenioides Vernonialeiocarpa Vernoniascorpioides
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 2 0 0
## 6 6 0 0
## Viburnumaff.lautum Viburnumjucundum Weinmanniapinnata
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 0 0
## Zanthoxylummelanostictum
## 1 5
## 2 0
## 3 0
## 4 0
## 5 0
## 6 2
env<-read.csv("env.csv")
head(env)
## X Forest.type Productivity Elevation
## 1 1 CF 0.642 2482
## 2 2 PF 0.621 2353
## 3 3 POF 0.899 2398
## 4 4 PF 0.584 2409
## 5 5 POF 0.655 2461
## 6 6 CF 0.616 2451
library(vegan)
## Loading required package: permute
##
## Attaching package: 'permute'
## The following object is masked from 'package:devtools':
##
## check
## Loading required package: lattice
## This is vegan 2.4-1
set.seed(0)#Para que los resultados no se brinden aleatorios
nmds1 <- metaMDS(bio[,2:86])
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.181971
## Run 1 stress 0.1903725
## Run 2 stress 0.1919633
## Run 3 stress 0.183261
## Run 4 stress 0.1894679
## Run 5 stress 0.1769346
## ... New best solution
## ... Procrustes: rmse 0.01052617 max resid 0.1464001
## Run 6 stress 0.1933279
## Run 7 stress 0.1975294
## Run 8 stress 0.1820899
## Run 9 stress 0.1769395
## ... Procrustes: rmse 0.002968075 max resid 0.03952345
## Run 10 stress 0.1769321
## ... New best solution
## ... Procrustes: rmse 0.0002113219 max resid 0.002012466
## ... Similar to previous best
## Run 11 stress 0.2007416
## Run 12 stress 0.187395
## Run 13 stress 0.1878302
## Run 14 stress 0.1817502
## Run 15 stress 0.1882357
## Run 16 stress 0.1812257
## Run 17 stress 0.1935675
## Run 18 stress 0.1769233
## ... New best solution
## ... Procrustes: rmse 0.001570212 max resid 0.02055457
## Run 19 stress 0.1848603
## Run 20 stress 0.1822167
## *** No convergence -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
plot(nmds1)
#Este plot no brinda mucha información
plot(nmds1, type = "n")
points(nmds1$points, pch = as.numeric(env$Forest.type), col = as.numeric(env$Forest.type), cex = 1.5)
legend(x = "topright", legend = c("Cloud forest", "Oak forest", "Pine forest", "Pine-oak forest", "Transitional forest"), pch = c(1:5), col = c(1:5))
ordiplot (nmds1, display = 'sites', type = 'n')
points (nmds1, col = as.numeric(env$Forest.type), pch = as.numeric(env$Forest.type))
for (i in unique (env$Forest.type)) ordihull (nmds1, groups = as.numeric(env$Forest.type, show.group = i, col = i), draw = 'polygon', label = T)
Vamos a ver un detalle importante de los PCA´s
library(vegan)
pca2 <-decorana(bio[,2:86])
library(vegan)
ordiplot (pca2, display = 'sites', type = 'n')
points (pca2, col = as.numeric(env$Forest.type), pch = as.numeric(env$Forest.type))
for (i in unique (env$Forest.type)) ordihull (pca2, groups = as.numeric(env$Forest.type, show.group = i, col = i), draw = 'polygon', label = T)
library (vegan3d)
ordirgl (pca2)
orglspider (pca2, groups = as.numeric(env$Forest.type))