domingo, 9 de agosto de 2015

Análisis de Componentes Principales en R [Principal Component Analysis using R]

En esta entrada explicare de la forma mas clara y simple posible que es un Análisis de Componentes Principales y como se hace en R. De  igual manera presento un forma bonita de graficar los resultados.


Cuando se va a campo, por lo general se toman muchos datos, pero son solo unos pocos las que contienen la mayor cantidad de información o la mas importante. El PCA es una técnica multivariada que se encarga de reducir el numero de variables, donde a partir de una matriz que contiene las variables originales genera una nueva con la misma cantidad de variables (generalmente) y casos, las nuevas variables recibirán el nombre de componentes (scores). El análisis genera 4 matrices en total: la primera corresponde a los Eingevalues de cada uno de ellos, la segunda a la cantidad de varianza de cada componente, una tercera con las correlaciones (loadings) entre cada una de las variables originales y los componentes, y finalmente la ultima matriz son las nuevas variables (componentes) creadas por el PCA para cada uno de los casos. Los componentes del PCA son los que se usaran posteriormente en otras técnicas como Análisis de Factores o regresiones múltiples. Para elegir los componentes principales a se usa  la cantidad de varianza acumulada (> 80% en general) o  componentes con Eingenvalues mayores a 1. Una vez determinamos los componentes principales, en la matriz de correlación las variables con valores mas altos serán las que estarán representadas en cada componente. Esta forma de reducir las variables con la matriz de correlación "a ojito" es algo subjetiva, por lo que existen otros métodos para eliminar esta subjetividad (pero... no hablare de ellos en este post).



Ahora veamos un ejemplo en R.

Para realizar un PCA en R existen varios paquetes: 


- prcomp () [stats]

- princomp () [stats]
- PCA () [FactoMineR]
- dudi.pca () [ade4]
- acp () [amap]

Si desea saber mas sobre cada una de las funciones puede visitar: 5 functions to do Principal Components Analysis in R


Tomare los datos morfológicos de el trabajo de Aaron M. Ellison & Elizabeth J. Fransworth.


Cargamos la librerías necesarias:


>library(FactoMineR)

>library(ggplot2)

Leo los datos y extraigo solo las variables morfológicas:


>Darli <- read.table("Darlingtonia.csv", header=T,sep=",")

>morfo <- Darli[,2:14]
>morfo <- as.data.frame(mapply(as.numeric,morfo))
>morfo <- as.data.frame(mapply(log,morfo)) ## Transformo las variables.
>head(morfo,5L) ##Miramos que tiene morfo.


En esta matriz tenemos 14 variables y 87 casos, es una matriz algo grande, por lo que haremos un PCA para reducirla. Usaremos la función prcomp del paquete stats.
 
  
>pca <- prcomp(morfo,scale=T)
>summary(pca)
 

El componente que acumula la mayor cantidad de varianza con casi el 98% es el PC1, seguido del PC2 y PC3 que juntos acumulan el 99%, ahora veamos los loadings para ver que variables son las mejor representadas en cada componente.
 
 
>pca$rotation

 
La selección de las variables es algo subjetiva, pero podemos ver que height es la variable mas importante para el PC1 con una correlación -0.49, otras variables, aunque con correlaciones no tan altas son: Wing1_length, Wing2_length y Wingsprea. Por otro lado para los componentes PC2 y PC3 son wingmass_g y tubemass_g respectivamente las variables mejor representadas, esto también lo podemos ver de forma gráfica:
 
 
>circle <- function(center=c(0,0),npoints=100){
    r <- 1
    tt <- seq(0,2*pi,length=npoints)
    xx <- center[1] + r * cos(tt)
    yy <- center[1] + r * sin(tt)
    return(data.frame(x=xx,y=yy))
  }

>corcir <- circle(c(0,0),npoints=100)

>corr <- as.data.frame(cor(morfo,pca$x))

>flechas <- data.frame(x1=rep(0,length(corr$PC1)),y1=rep(0,length(corr$PC2)  ),x2=corr$PC1,y2=corr$PC2)
  
>g <- ggplot() + geom_path(data=corcir,aes(x=x,y=y),colour="gray65")

>g <- g +  geom_segment(data=flechas,aes(x=x1,y=y1,xend=x2,yend=y2),col  our="gray65")

>g <- g + geom_text(data=corr,aes(x=PC1,y=PC2,label=rownames(corr)),size   =10)

>g <- g + geom_hline(yintercept=0,colour="gray65")+geom_vline(xintercep=0  ,colour="gray65")

>g <- g + labs(title="Loadings de los PC1 y PC2")+xlab("PC1")+ylab("PC2")
>g


Con esta gráfica los que se pretende es que las variables con una alta correlación en los loadings aparezcan mas cerca al eje del componente, pero entes caso no ocurre, pues height que es la variables mas importante para el PC1 aparece muy lejos de su eje, pero no hay que alarmarse una forma de poder solucionar esto es aplicar una rotación a los componentes con la función varimax( ) [stats] y listo !!.  Con un PCA también podemos ver como se distribuyen los casos (sitios en este ejemplo) dentro del morfoespacio usando los scores de los componentes.
 
 
> m <- ggplot()

> m <- m + geom_point(aes(x=pca$x[,1],y=pca$x[,2],shape=Darli$Sitios,colour=Darli$Sitios,group=Darli$Sitios),size=10)

> m <- m + xlab("PC1")+ylab("PC2")+labs(title="PCA para datos morfológicos")

> m <- m + theme(plot.title=element_text(size=30,face="bold"),
           axis.title.x=element_text(size=20),axis.title.y=element_text(size=20))

> m <- m+ geom_hline(aes(yintercept=0),linetype=2)

> m <-m+ geom_vline(aes(yintercept=0),linetype=2)

> m <- m + scale_colour_manual(name="Sitios",values=c(DG="red",HD="blue",LEH="green",TJH="orange"))

> m <- m + scale_shape_manual(name="Sitios",values=c(DG=7,HD=8,LEH=9,TJH=10))

> m



En esta grafica podemos ver como se agrupan los sitios dentro del morfoespacio del PCA, por ejemplo vemos que HD es el que mejor se diferencia y se agrupa en su mayoria hacia los valores altos del eje X (PC1), por lo que se puede decir que los sitios HD presentan valores altos del PC1, mientras que LEH y DG en su mayoria presentan valores bajos del PC1. Estas tendencias tambien pueden verse  usando el componente 2 y 3. 
 
Para profundizar mas en PCA, en espcieal en su matemática recomiendo estos dos libros: Experimental Design and Data Analysis for Biologist y Applied Multivariate Statistical Analysis