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.

Análisis de componentes principales (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)

Escalamiento multidimensional no métrico (NMDS)

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