miércoles, 19 de octubre de 2011

Análisis de Correspondencia Canónica en R

El análisis de correspondencia canónica es una técnica multivariante que  maximiza la relación entre una serie de variables dependientes y una serie de variables independientes. Esta relación es hecha en base de regresión múltiple.
(mayor información).

Es muy común que se utilicen dos matrices; una de variables dependientes (por ejemplo, abundancia de especies en X's sitios) y otra de variables independientes o explicativas (por ejemplo, variables climáticas de los X's sitios).
Para nuestro ejemplo utilizaremos este tipo de datos.

Para nuestra primer matriz tenemos los datos de 10 sitios, con la abundancia de 5 especies para estos sitios, lo guardaremos como un archivo formato .CSV en la carpeta sobre la que estemos trabajando, con el nombre "especies.csv",asi:


sitios,sp1,sp2,sp3,sp4,sp5
1,1,0,0,1,0
2,0,0,0,0,0
3,0,0,0,0,0
4,1,8,1,5,1
5,1,0,3,1,0
6,3,5,0,1,0
7,117,288,141,100,14
8,0,0,0,0,1
9,0,0,0,0,0
10,1,2,26,3,2


Para la segunda matriz, tenemos los datos de los mismos 10 sitios con sus respectivos valores de 5 variables climáticas, en una archivo formato .CSV en la misma carpeta anterior, con el nombre "clima.csv", asi:

sitios,v1,v2,v3,v4,v5
1,28.4,80.1,30.7,54.3,32.4
2,26.8,85.8,19.4,23.7,51.2
3,25.5,85.2,56.2,66.2,60.0
4,27.3,90.4,43.3,52.1,53.0
5,24.6,83.8,12.5,22.6,69.0
6,24.5,83.8,12.5,25.0,66.3
7,26.0,81.2,32.0,77.0,65.0
8,25.1,86.6,11.3,45.2,83.1
9,27.5,84.5,23.7,41.2,55.7
10,26.0,94.1,20.0,80.0,60.0

Teniendo nuestras matrices de datos listas, nos dispondremos a trabajar sobre R:

Lo primero que debemos hacer, es descargar los paquetes que necesitemos para hacer nuestro analisis de correspondencia canónica, estos paquetes son "fields" y "CCA"; para esto ejecutamos los siguientes comandos:


>install.packages ("fields")
>install.packages ("CCA")

# Descargamos nuestro paquete del espejo deseado (yo, utilizo "USA CA2"), a nuestro computador.


>library (fields)
>library (CCA)

#Cargamos a nuestra area de trabajo en R, los paquetes recien descargados.

Una vez teniendo los paquetes cargados, introducimos en R, como objetos, las dos matrices que creamos que estan en nuestra carpeta de trabajo:


>clima <- read.csv("clima.csv",header=T)
# llamamos nuestra matriz de datos de clima


>especies <- read.csv("especies.csv",header=T)
#Llamamos nuestra matriz de datos de especies por sitios, deben recordar que el argumento "header=T", quiere decir que la primer fila de las matrices corresponde a los nombres de las columnas.

Como en nuestras dos matrices, tenemos enumerados los sitios del 1 al 10, es importante que esta columna, no entre en el análisis de correlación de las variables, para esto, para esto modificamos los dos objetos que acabamos de crear en R ("especies" y "clima"), eliminando la primera columna de cada matriz, utilizando los siguientes comandos:


>especies<-especies[2:6]
# Esta orden le dice a R, que tome desde la columna 2 hasta la 6 [2:6] , de la matriz original de especies.


>clima<-clima[2:6]
# Esta orden le dice a R, que tome desde la columna 2 hasta la 6 [2:6] , de la matriz original de clima.


> ñero<-matcor(clima,especies)
# Realizamos varias matrices de correlacion de las matrices originales, es decir; realizamos matrices de correlación entre matrices y dentro de cada matriz y nombramos estos resultados con el nombre de "ñero".

Si queremos ver los resultados de este análisis de correlación simplemente tenemos que darle a R el nombre del objeto, osea "ñero".


>ñero

$Xcor
            v1          v2          v3         v4         v5
v1  1.00000000 -0.07921803  0.36040283  0.2125409 -0.8404823
v2 -0.07921803  1.00000000  0.00531708  0.2648932  0.1915028
v3  0.36040283  0.00531708  1.00000000  0.5353831 -0.3988760
v4  0.21254086  0.26489317  0.53538314  1.0000000 -0.0983227
v5 -0.84048225  0.19150278 -0.39887598 -0.0983227  1.0000000

$Ycor
          sp1       sp2       sp3       sp4       sp5
sp1 1.0000000 0.9996739 0.9834022 0.9987807 0.9870899
sp2 0.9996739 1.0000000 0.9830768 0.9994826 0.9885262
sp3 0.9834022 0.9830768 1.0000000 0.9858946 0.9938675
sp4 0.9987807 0.9994826 0.9858946 1.0000000 0.9911531
sp5 0.9870899 0.9885262 0.9938675 0.9911531 1.0000000

$XYcor
             v1          v2          v3         v4         v5         sp1
v1   1.00000000 -0.07921803  0.36040283  0.2125409 -0.8404823 -0.05392549
v2  -0.07921803  1.00000000  0.00531708  0.2648932  0.1915028 -0.37144081
v3   0.36040283  0.00531708  1.00000000  0.5353831 -0.3988760  0.13278602
v4   0.21254086  0.26489317  0.53538314  1.0000000 -0.0983227  0.46217977
v5  -0.84048225  0.19150278 -0.39887598 -0.0983227  1.0000000  0.14388437
sp1 -0.05392549 -0.37144081  0.13278602  0.4621798  0.1438844  1.00000000
sp2 -0.04569713 -0.35776621  0.14543384  0.4681405  0.1429093  0.99967391
sp3 -0.06149136 -0.23766354  0.10945076  0.5577507  0.1510460  0.98340218
sp4 -0.03504070 -0.33865725  0.15226692  0.4835602  0.1339936  0.99878072
sp5 -0.05218253 -0.23469160  0.12520778  0.5517128  0.1810415  0.98708989
            sp2         sp3        sp4         sp5
v1  -0.04569713 -0.06149136 -0.0350407 -0.05218253
v2  -0.35776621 -0.23766354 -0.3386572 -0.23469160
v3   0.14543384  0.10945076  0.1522669  0.12520778
v4   0.46814048  0.55775070  0.4835602  0.55171276
v5   0.14290927  0.15104604  0.1339936  0.18104151
sp1  0.99967391  0.98340218  0.9987807  0.98708989
sp2  1.00000000  0.98307679  0.9994826  0.98852620
sp3  0.98307679  1.00000000  0.9858946  0.99386751
sp4  0.99948262  0.98589456  1.0000000  0.99115306
sp5  0.98852620  0.99386751  0.9911531  1.00000000


Después procedemos a hacer un análisis de correlación canónica con el siguiente comando:


>gala <- cc(clima,especies)
#con este comenado hacemos el análisis de correlación canónica y llamamos a este análisis "gala".

Si queremos ver los resultados de este análisis, solo tenemos que nombrarlo, es decir, escribir "gala".


>gala

$cor
[1] 1.0000000 0.9873988 0.8237346 0.5648062 0.2742222

$names
$names$Xnames
[1] "v1" "v2" "v3" "v4" "v5"

$names$Ynames
[1] "sp1" "sp2" "sp3" "sp4" "sp5"

$names$ind.names
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"


$xcoef
           [,1]         [,2]         [,3]        [,4]         [,5]
v1  0.336685062 -0.230574895 -1.009579966 -0.04601321  0.985822606
v2  0.230731161 -0.024787992  0.085323153  0.07411478 -0.025693901
v3  0.017698536 -0.061736475 -0.033217231 -0.02906274 -0.043634878
v4 -0.006286382  0.053495860  0.007531747 -0.02614144  0.001067509
v5  0.021817542 -0.008533328 -0.144485604  0.01789087  0.028632651

$ycoef
           [,1]        [,2]       [,3]         [,4]       [,5]
sp1 -0.98163323  0.64368084  0.2355047  0.441119087  1.2186888
sp2  0.42111846 -0.29567457 -0.2818513 -0.009289717 -0.9252762
sp3  0.08332968 -0.01161349  0.1318968 -0.061027074 -0.2518835
sp4 -0.19737790 -0.07755685  0.5919466 -0.634888351  1.3265043
sp5  0.02968457  1.51438054 -1.8122917  1.499392109  1.9448609

$scores
$scores$xscores
            [,1]        [,2]       [,3]        [,4]        [,5]
 [1,] -1.0541235 -0.12954861  1.1004450 -1.27018272  1.36831073
 [2,]  0.1248876 -1.00169802  0.6306688  0.69057810  0.64324152
 [3,]  0.1248876 -0.76049941 -0.2818393 -1.31715393 -1.93133861
 [4,]  1.6383268 -1.20258957 -0.3216990 -0.39631331  0.05694330
 [5,] -0.8041345 -0.22961427  0.3301686  1.19132365 -0.66461283
 [6,] -0.9117977 -0.05512673  0.8393139  1.08488015 -0.83794123
 [7,] -1.0168043  1.25247646 -0.9651501 -1.12617449 -0.12499530
 [8,]  0.1545721  0.74846218 -1.7628854  1.07217854  0.23656347
 [9,]  0.1248876 -0.49856533 -0.8481710  0.06008361  1.32661776
[10,]  1.6192984  1.87670328  1.2791486  0.01078041 -0.07278881

$scores$yscores
            [,1]       [,2]        [,3]        [,4]       [,5]
 [1,] -1.0541235 -0.1229924  0.88336995  0.00975034  1.5515810
 [2,]  0.1248876 -0.6891164  0.05591873  0.20351960 -0.9936121
 [3,]  0.1248876 -0.6891164  0.05591873  0.20351960 -0.9936121
 [4,]  1.6383268 -1.2958493 -0.68404897 -1.16575576  1.1483660
 [5,] -0.8041345 -0.1578328  1.27906044 -0.17333088  0.7959305
 [6,] -0.9117977 -0.3140035 -0.05487719  0.84553993 -0.6374223
 [7,] -1.0168043  1.2754077 -1.14318226 -1.96314865 -0.5236562
 [8,]  0.1545721  0.8252642 -1.75637292  1.70291171  0.9512489
 [9,]  0.1248876 -0.6891164  0.05591873  0.20351960 -0.9936121
[10,]  1.6192984  1.8573552  1.30829476  0.13347450 -0.3052117

$scores$corr.X.xscores
         [,1]       [,2]        [,3]       [,4]       [,5]
v1 0.18557161 -0.2821097  0.12445671 -0.5554653  0.7496307
v2 0.94032126  0.1960841  0.12974856  0.2071041 -0.1327409
v3 0.23588079 -0.3601799 -0.10955360 -0.8382693 -0.3161224
v4 0.32283979  0.5725302 -0.09871644 -0.7387230 -0.1119242
v5 0.01251680  0.3692552 -0.56373702  0.5710781 -0.4685769

$scores$corr.Y.xscores
          [,1]      [,2]       [,3]       [,4]        [,5]
sp1 -0.3635856 0.4433349 -0.3264553 -0.3897888 -0.04971763
sp2 -0.3447263 0.4348985 -0.3358075 -0.3950264 -0.04891287
sp3 -0.2560776 0.5608868 -0.2563731 -0.3898132 -0.05383898
sp4 -0.3246888 0.4420926 -0.3268741 -0.4034079 -0.04464397
sp5 -0.2355762 0.5332581 -0.3337258 -0.3857010 -0.04100433

$scores$corr.X.yscores
         [,1]       [,2]        [,3]       [,4]        [,5]
v1 0.18557161 -0.2785548  0.10251931 -0.3137303  0.20556540
v2 0.94032126  0.1936132  0.10687838  0.1169737 -0.03640051
v3 0.23588079 -0.3556412 -0.09024310 -0.4734597 -0.08668778
v4 0.32283979  0.5653156 -0.08131615 -0.4172354 -0.03069211
v5 0.01251680  0.3646021 -0.46436971  0.3225485 -0.12849421

$scores$corr.Y.yscores
          [,1]      [,2]       [,3]       [,4]       [,5]
sp1 -0.3635856 0.4489927 -0.3963112 -0.6901284 -0.1813042
sp2 -0.3447263 0.4404487 -0.4076646 -0.6994016 -0.1783694
sp3 -0.2560776 0.5680449 -0.3112326 -0.6901715 -0.1963334
sp4 -0.3246888 0.4477346 -0.3968196 -0.7142413 -0.1628021
sp5 -0.2355762 0.5400635 -0.4051375 -0.6828908 -0.1495296




Después de esto solo tenemos que proceder a realizar los resultados gráficos del análisis de correspondencia canónica.


>img.matcor(ñero, type = 2)
#Obtenemos el resultado de la función de correlación entre y dentro de las matrices, por ejemplo en la siguiente gráfica podemos observar una alta correlación para la matriz de especies, pero en la matriz de variables climáticas vemos una correlación menos marcada, asi mismo, podemos ver el diagrama de correlación cruzada (es decir, de ambas matrices).
















>barplot(gala$cor, xlab = "Dimension",ylab = "Canonical correlations",
names.arg = 1:5, ylim = c(0,1))

#Realizamos el gráfico de explicación de cada dimensión a las que fueron reducidas las variables utilizadas, ya que sobre este es que se decidirá con cuales dimensiones o "componentes" se trabajara.


>plt.cc(gala)
#Finalmente, utilizando las dos dimensiones que mas "poder explicativo tienen", hacemos nuestro diagrama de correspondencia el cual nos servirá para agrupar los sitios, las variables y las especies correlacionadas (por ejemplo, en la siguiente gráfica podemos observar la alta similaridad entre los sitios 9,3 y 2, y como las especies se encuentra muy cercanas (muy correlacionadas))







Eso fue todo por hoy, para una mejor interpretación de sus resultados puede ir a este link.
Recuerden que este es un blog de implementación de métodos estadísticos en R, y por ello la explicación de los resultados es muy superficial, eso corre por cuenta de uds. y de cada uno de sus trabajos de investigación y datos.!!!
Nos pillamos, mano..... 

14 comentarios:

  1. Hola, una pregunta, tengo problemas con R para descargar los paquetes siempre me sale error descargue la última versión para windows 7 de 64 bits. Sabes que puedo hacer? o como los descargo? Gracias

    ResponderEliminar
  2. Es muy extraño que tenga errores descargando los paquetes, generalmente es super sencillo!, y sobre todo en windows; pueden ser muchas cosas las que esten pasando, pero basicamente trate de cambiar de espejo (de CRAN mirror), yo personalmente utilizo "CA2", que es de USA, tambien puede ser que tenga errores el R instalado en su compu, podria eliminarlo y volver a instalarlo de nuevo!;
    lo mas sencillo para descargar paquetes es definitivamente

    >install.packages ("nombre del paquete")
    # el nombre del paquete debe ir entre comillas y dentro de parentesis.

    ejemplo:
    >install.packages ("phytools")

    despues solo tiene que cargarlo a su area de trabajo:

    >library (phytools)

    cualquier duda, aqui estare!

    ResponderEliminar
    Respuestas
    1. Que tal gracias por hacer este blog ¿como soluciono el siguiente error?

      Error: could not find function "matcor"

      Eliminar
  3. Análisis de Correspondencia Canónica en R.......
    creo que debemos modificar el nombre de la entrada, es ANÁLISIS DE CORRELACIÓN CANÓNICA EN R....
    PD: muy bueno tu blog, me encanta....
    :)

    ResponderEliminar
  4. Felicitaciones por esta iniciativa!!! Aunque los blogs y recursos para R online abundan, siempre es bueno encontrarse con uno que es por y para chibchombianos.

    ResponderEliminar
  5. Creo que lo mas conveniente para la representacion grafica de un CCA (principal producto de un CCA) es un triplot de puntos para las especies y , vectores flecha para las variables ambientales (ver ter Braak 1995),ademas de la prueba de significancia de Monte Carlo, sin embargo te felicito, excelente tutorial

    ResponderEliminar
  6. Que tal; Felicidades y gracias por compartir el blog. He estado tratando de usar el script para correspondencias canónicas pero me emite el siguiente mensaje ¿cómo puedo solucionarlo?

    Error: could not find function "matcor"

    ResponderEliminar
  7. Alguien podria auxiliarme con este error:
    Error in chol.default(Bmat) :
    the leading minor of order 63 is not positive definite

    ResponderEliminar
  8. Hola, que pena no poder haber contestado antes...
    la función "matcor" esta dentro del paquete CCA. lo más posible es que si tu R no encuentra la función es porque el paquete no esta instalado o se te pifio y se quebro, intenta inatalarlo o reinstalarlo.

    install.packages ("CCA") ## instalar el paquete
    library(CCA) ## cargar el paquete a tu area de trabajo
    ?matcor ## con este comando puedes ver si matcor esta presente y como funciona.

    Gracias por leer el blog. espero te sirva de ayuda.

    pd: si tienes problemas para instalar el paqquete busca en google, ya que este paquete puede tener algunas dependencias o tener requerimientos particulares como por ejemplo la versión de R que estes usando.

    ResponderEliminar
    Respuestas
    1. Ambrosio:
      Pues gracias por la respuesta, ya había podido solucionar este error al que haces referencia. Pero lo que no logro decifrar es lo siguiente: Alguien podria auxiliarme con este error:
      Error in chol.default(Bmat) :
      the leading minor of order 63 is not positive definite

      Este mensaje aparece cuando quiero aplicar la función: xxcc <- cc(clima,especies)
      ¿sabes como puedo rebasar este error?
      Saludos y gracias nuevamente.

      Eliminar
  9. Hola. Me puedes explicar por qué al poner la selección de las columnas, el análisis no me toma en cuenta la última columna, ni de las especies ni de las variables? Muchas gracias

    ResponderEliminar
  10. Hola, quisiera saber como le agregas a la grafica de ordenacion del CCA otros datos, por ejemplo los meses de muestreo, he tratado de varias formas pero no he podido. tengo las var ambientales y las de los bichos y los puntos de los sitios pero necesito agregar a la grafica la columna de datos de los meses. Te agradezco si me puedes colaborar.

    ResponderEliminar
  11. Hola me gustariia poner etiquetas a los puntosen los graficos de puntos y triangulos

    ResponderEliminar
  12. Hola si me sale este error alguien me puede ayudar!
    Error en chol.default(Bmat) :
    la submatriz de orden 4 no es definida positiva
    > gala
    Error: objeto 'gala' no encontrado

    ResponderEliminar