jueves, 15 de diciembre de 2011

Entrada y manejo de datos en R (guía básica)

Hace unos días, un muy buen amigo mio (Jeffrey Vega), me pregunto acerca de la entrada de datos a R, específicamente acerca de archivos .CSV.
En ese momento me dio la luz para escribir este post (no sabia sobre que escribir, con tantas cosas que se pueden hacer!), ya que la entrada de datos en R, es sumamente importante y teniendo claro esto, muchas de las cosas que uno hará en R serán "breves*" (como decimos los "ñeros" jejeje).


*BREVE= en el lenguaje "ñeristico", significa, fácil, sin complique, suave, etc!

Para empezar, aclarare donde deben estar nuestros archivos de datos en nuestro computador y como decirle a R, donde estan o como decirle cuando cambiar de lugar para trabajar.


Entonces, como bien sabemos en nuestro computador podriamos tener 3.548.125 carpetas en lugares como "C", "D", el Escritorio u otros dentro de estos o afuera de estos, por ejemplo:

Podría tener un archivo en formato ".txt" (mas adelante explicare lo de los formatos!) , llamado "arc_1.txt" en mi escritorio:




O podria tener otro archivo llamado "arc_2.txt" en una carpeta personal, en este caso (una carpeta mia llamada "morphometria"):




Entonces, como puedo tener archivos en muchas carpetas (no solo en estas dos), debo indicarle a R donde se encuentran los archivos que voy a utilizar en esta oportunidad; esto lo podemos hacer via linea de comandos, o lo podemos hacer graficamente.
En windows es muy sencillo hacerlo, puesto que solo se debe dar click, en "cambiar el directorio".



En linux es aun mas sencillo porque solo se debe utilizar un comando con la ruta que estemos trabajando, este comando es "setwd".
Por ejemplo, si quiero decirle que vamos a trabajar sobre el escritorio en donde esta mi archivo "arc_1.txt", le damos:


> setwd ("/home/ambrosio/Escritorio")

Y si por ejemplo, queremos trabajar en mi carpeta "morphometria", donde tengo mi archivo "arc_2.txt", solo le tenemos que indicar la ruta para acceder a la carpeta:


> setwd ("/home/ambrosio/morphometria")




Y listo!!!, de este modo le indicamos a R, donde están nuestros archivos!!
Muchas veces, a muchas personas se les hace difícil, saber la ruta, porque es una ruta muy larga o porque no la conocen, para esto solo deben pararse sobre el archivo que tenga sus datos, le dan click derecho, después le damos click a propiedades y aparecerá la ruta del archivo:




Una vez estemos en el directorio sobre el que queramos trabajar, tenemos que incluir los archivos de datos que queramos utilizar dentro de el directorio (aunque en el mismo R podemos escribir los datos con los que vamos a trabajar, como vectores, tablas, etc etc etc. pero este no es el objeto de este post, ya que estamos trabajando acerca incluir datos a R).

Para esto lo primero que tenemos que tener en cuenta son los formatos de los archivos, personalmente los archivos que yo utilizo, estan en formato .csv o .txt, es decir archivos separados por comas "," o archivos separados por espacios " ".
son muy faciles de utilizar porque pueden ser modificados en excel (para los que aman excel!!) o en editores de texto como Geany (para los que amamos geany!!) y es sencillo manejar grandes cantidades de datos con estos formatos!!!
Para aprender a crear archivos .csv  y .txt he generado unos datos !IMAGINARIOS¡ acerca de algunas variables "medidas" en las 7 etapas del barrio Zapamanga de Floridablanca/Santander. Los datos fueron incluidos dentro de una tabla utilizando LibreOffice  3.3  (Abajo excel!), de allí vamos a guardarlos en formato .csv y en geany los transformaremos en formato .txt (para observar como es de fácil transformar en geany formatos sin tener que ir de nuevo a hojas de calculo y guardar en formato .txt .

Los datos originales en una tabla son estos:



Entonces guardamos la tabla en formato .csv yendo a archivo y dando clic en "guardar como..", allí escojemos la opción de guardar en formato .csv, como se muestra en la imagen:





De este modo ya tenemos guardado nuestros datos en formato .csv, pero hay algo muy importante para tener en cuenta y es que ni openoffice, ni excel, ni libreoffice utiliza puntos para los decimales (Y R SI!!!), por lo que en este sentido ninguno es eficiente porque R leera las comas de los decimales como si fueran dos numeros de dos columnas distintas y no de la misma; es decir, supongamos que tenemos el decimal 24,5 en una hoja de calculo, al guardarlo en formato .csv se guardara como tal (24,5), pero al leerlo en R no se leera como un solo numero de una columna sino como dos numero apartes de dos columnas contiguas 24 y 5.

Al abrirlo en geany se vería así y no se distinguirán los decimales:




Sí es un archivo con pocos datos (como en este caso) se pueden cambiar las comas de los decimales por puntos manualmente uno por uno, pero esto no seria muy practico al momento de tener archivos con miles de datos, así que tomaremos la vía de edición de formatos con editor de texto (valga la redundancia).

Entonces lo único que tenemos que hacer es seleccionar y copiar los datos de la hoja de calculo y pegarlos en una hoja en blanco de geany (o con el editor de texto que estemos trabajando), nos vamos a la solapa que dice "Buscar" damos clic ahi, y despues damos clic en "Reemplazar", finalmente le decimos que reemplace las comas por puntos en el documento y ya tenemos nuestros decimales definidos por puntos (ahora lo que tenemos que hacer es separar nuestras columnas por comas o por espacios, segun lo que queramos hacer (o .txt o .csv).


Para separar nuestras columnas entonces seleccionamos un espacio en blanco entre dos números: 



Vamos de nuevo a "Buscar/Reemplazar" y automáticamente en el espacio de "Buscar por:" aparece la selección que hemos hecho, lo único que tenemos que hacer es ponerle una coma "," en el espacio de "Reemplazar con:" y le damos reemplazar en el documento: 



Finalmente lo que obtnemos es nuestro archivo en formato .csv (es decir las columnas separadas por comas) listo para ser llamado a R!!!
Si hubiésemos querido el archivo en formato .txt lo unico que teniamos que hacer era reemplazar la seleccion hecha en el paso anterior por un espacio en blanco y ya!!!.
Obviamente una vez terminada nuestra edición debemos guardar el archivo en la carpeta que estemos trabajando, entonces, nos vamos en geany a "Archivo", le damos "Guardar como" y nombramos el archivo como "luna.csv" (en honor a la nueva mascota de Juancho) y sí la tenemos en formato .txt la nombramos "luna.txt". Este seria el aspecto final de los archivos en ambos formatos:


 Archivo en formato .txt (luna.txt)



Archivo en formato .csv (luna.csv)

Una vez teniendo listos nuestros datos en la carpeta que estemos utilizando, los incluiremos en R, para esto necesitamos simples comandos de lectura en R, dependiendo si el formato que vamos a  leer es .csv o .txt.
empezaremos llamando nuestro formato .csv

> chorizo <- read.csv("luna.csv",header=T)
#Este comando indica que crearemos un objeto con el nombre "chorizo" en R, a partir de una matriz de datos externa (nuestro "luna.csv"), para esto se utiliza la función "read.csv" y lo que indica que el nuevo objeto en R se llamara chorizo es la flecha "<-", el argumento header=T o header=TRUE, significa que la primer linea o fila de la matriz, son los nombres de las variables, de nos ser así, se tendría que poner, header=F o header=FALSE.

Para ver el nuevo objeto que hemos creado en R, solo tenemos que darle el nombre del objeto es decir "chorizo"

> chorizo
#
Y veremos algo así en la ventana de R:

          Etapa Temperatura Humedad Calentura Muertos.año Poblacion Estrato
1   Zapamanga_I        24.5      82         7           3      5300       2
2  Zapamanga_II        24.3      81         5           1      4300       2
3 Zapamanga_III        26.0      84         7           3      6200       2
4  Zapamanga_IV        27.8      83         8           6      8260       2
5   Zapamanga_V        24.0      86         6           2      4260       2
6  Zapamanga_VI        25.3      84         6           1      4500       2
7 Zapamanga_VII        25.0      82         6           2      5230       2
  No_canchas No_escuelas
1          2           1
2          0           2
3          2           1
4          2           3
5          1           1
6          0           1
7          1           1


Después leeremos esta misma matriz de datos pero en formato .txt y utilizaremos el siguiente comando:

> galletas <- read.table("luna.txt",header=T)
#Con este comando indicamos que queremos crear un objeto en R con el nombre de "galletas", para eso utilizamos la función read.table (que es para leer archivos en formato.txt, el argumento "header=T", se utiliza de la misma forma como con la función read.csv, lo mismo que el uso de la flecha.
para ver el objeto hacemos lo mismo que hicimos antes, pero en vez de escribir en R "chorizo", escribimos "galletas" y veremos nuestros datos.

> galletas
#
asi veremos los mismos datos que mostré anteriormente.


Una vez teniendo nuestros objetos creados en R, vamos a ver un poco acerca de como manejarlos y un tanto "jugar con ellos".
Primero que todo, debemos saber (si queremos!) que objetos hemos creado en R, para esto solo le damos el comando:

>objects ()
#esto nos muestra los objetos que hemos creado en R:

[1] "chorizo"  "galletas"

En ocasiones necesitamos hacer analisis con solo algunas variables de las matrices, o solo necesitamos tomar algunos datos de esta; para eso podemos utilizar submatrices, que son muy facil de crear, o se pueden tomar los datos directamente an el analisis, utilizando los mismos argumentos de creacion de submatrices.
Por ejemplo si queremos tomar solo las variables temperatura y humedad, le indicamos con un comando que esas variables corresponden a las columnas 2 y 3 y el damos el comando:




> choricito<-(chorizo[,2:3])
#con esto le decimos que tome de la matriz original "chorizo" las columnas 2 y 3 y cree una submatriz llamada "choricito", asi que al darle el siguiente comando, nos aparecera la nueva mariz:


>choricito

  Temperatura Humedad
1        24.5      82
2        24.3      81
3        26.0      84
4        27.8      83
5        24.0      86
6        25.3      84
7        25.0      82


Y si queremos decirle que tome las 4 primeras variables solo tenemos que darle el sigueinte comando:

>a<-(chorizo[,1:4])

>a


          Etapa Temperatura Humedad Calentura
1   Zapamanga_I        24.5      82         7
2  Zapamanga_II        24.3      81         5
3 Zapamanga_III        26.0      84         7
4  Zapamanga_IV        27.8      83         8
5   Zapamanga_V        24.0      86         6
6  Zapamanga_VI        25.3      84         6
7 Zapamanga_VII        25.0      82         6

Ahora bien, esto lo realizamos para tomar el numero de columnas deseado o las variables que queramos incluir en el análisis, pero si por ejemplo queremos mas bien escoger el numero de filas o el numero de lugares a evaluar lo único que tenemos que hacer es cambiar la posición de la coma en el comando y ponerla al final, es decir:



>alto<-(chorizo[1:4,])
>alto

          Etapa Temperatura Humedad Calentura Muertos.año Poblacion Estrato
1   Zapamanga_I        24.5      82         7           3      5300       2                                                                                                                                                      
2  Zapamanga_II        24.3      81         5           1      4300       2                                                                                                                                                      
3 Zapamanga_III        26.0      84         7           3      6200       2                                                                                                                                                      
4  Zapamanga_IV        27.8      83         8           6      8260       2                                                                                                                                                      
  No_canchas No_escuelas                                                                                                                                                                                                         
1          2           1                                                                                                                                                                                                         
2          0           2                                                                                                                                                                                                         
3          2           1                                                                                                                                                                                                         
4          2           3    

> galle<-(luna[2:4,])
>galle


          Etapa Temperatura Humedad Calentura Muertos.año Poblacion Estrato
2  Zapamanga_II        24.3      81         5           1      4300       2
3 Zapamanga_III        26.0      84         7           3      6200       2
4  Zapamanga_IV        27.8      83         8           6      8260       2
  No_canchas No_escuelas
2          0           2
3          2           1
4          2           3



Con esta pequeña introduccion al manejo de datos en R, es muy facil empezar a hacer nuestros analisis, otra herramienta muy util en el manejo de datos son los archivos multiples tipo  "Array", que explique en uno de mis post anteriores:
Array's en R

Espero que este post, haya sido de gran utilidad, especialmente para las personas que estan empezando a utilizar R.











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

miércoles, 3 de agosto de 2011

Graficas de Filomorfoespacios en R

En diciembre de 2010 Liam Revell escribió la version 1 de la función para trazar gráficos de filomorfoespacios siguiendo el método de Sidlaukas (2008); Recientemente (agosto, 2011) el incorporo la version 3 de la función a la versión beta del paquete "phytools", desarrollado por Revell.


Al principio tuve varios inconvenientes para probar la función, pero ahora que ha sido mejorada y que la he cacharreado bastante, espero explicar como llevar a cabo el análisis y obtener la gráfica.


Para empezar, definiré las características que deben tener las entradas ("imputs") de nuestros datos, ya que sin duda alguna la gran mayoría de los errores que podemos tener están en los formatos y en los argumentos que utilizamos.


Lo primero que debemos tener es un árbol en formato .phylo en el que no deben ir las etiquetas de los nodos ya que esto puede hacer que la función confunda los valores de la matriz que utilicemos de estados de caracteres ancestrales (sí no utilizamos las mismas etiquetas para la matriz). Es decir:


Consideremos el siguiente fragmento de un árbol:


(Struthioniformes:21.8,Tinamiformes:21.8)78:4.1
en este fragmento podemos ver que la etiqueta del nodo es "78", esto quiere decir que en la matriz de estado de caracteres ancestrales este nodo debe estar etiquetado con el mismo numero ="78".
Pero para evitarnos ese problema podemos dejar el árbol sin etiquetas de los nodos, así:


(Struthioniformes:21.8,Tinamiformes:21.8):4.1


Para nuestro ejemplo, utilizaremos un árbol que viene con el paquete "ape" = data(bird.orders).
y esta construido de la siguiente forma:


(((Struthioniformes:21.8,Tinamiformes:21.8):4.1,((Craciformes:21.6,Galliformes:21.6):1.3,Anseriformes:22.9):3):2.1,(Turniciformes:27,(Piciformes:26.3,((Galbuliformes:24.4,((Bucerotiformes:20.8,Upupiformes:20.8):2.6,(Trogoniformes:22.1,Coraciiformes:22.1):1.3):1):0.6,(Coliiformes:24.5,(Cuculiformes:23.7,(Psittaciformes:23.1,(((Apodiformes:21.3,Trochiliformes:21.3):0.6,(Musophagiformes:20.4,Strigiformes:20.4):1.5):0.6,((Columbiformes:20.8,(Gruiformes:20.1,Ciconiiformes:20.1):0.7):0.8,Passeriformes:21.6):0.9):0.6):0.6):0.8):0.5):1.3):0.7):1);


debido a que la función de filomorfoespacios, trabaja con objetos .phylo lo único que tenemos que hacer es copiar el árbol anterior (en paréntesis) en un editor de texto (geany, worpad, gedit, etc) y guardarlo como "prueba.phylo".

Para cargar el árbol a R solo necesitamos darle el comando, (obviamente debemos tener el archivo en el directorio que estemos trabajando!!!):


>prueba<-read.tree (file="prueba.phylo")

y debe tener el siguiente aspecto, sí le damos:


>plot (prueba)



Después de esto, necesitaremos cargar nuestros valores para los estados  de caracteres de las terminales, estos caracteres pueden ser discretos o continuos y pueden ser medidas (como peso, altura, longitud rostro-cloaca, etc) o pueden ser valores de análisis de reducción de variables (PCA, PCO, relative warps (RW), etc).
Para este ejemplo utilizaremos valores que se pueden generar en R aleatoriamente con la función:


>x <- rnorm(23)
#genera 23 valores al azar y los llama X.


suponiendo que generamos dos veces (para tener valores de dos variables) valores, pondremos en una matriz los valores con el nombre de las etiquetas de las terminales, el archivo debe ser guardado en formato .csv (aunque puede estar en otro formato como .txt) con el nombre de "tips.csv" y debe lucir asi:


"V1","V2"
"Struthioniformes",-0.292614262,-1.1963108
"Tinamiformes",-0.057614850,1.0345608
"Craciformes",-0.006225218,0.4230613
"Galliformes",-0.475368798,-0.1927969
"Anseriformes",0.668049937,0.2761842
"Turniciformes",0.769136135,-0.1066864
"Piciformes",1.376599742,0.1124578
"Galbuliformes",-0.063333022,-0.5463322
"Bucerotiformes",0.195140568,-0.1388508
"Upupiformes",1.406734039,0.7346157
"Trogoniformes",1.509326716,0.8676560
"Coraciiformes",0.230836319,0.3352026
"Coliiformes",-0.049803844,1.6605007
"Cuculiformes",-0.433148495,0.1290520
"Psittaciformes",0.544576328,1.0001100
"Apodiformes",0.348601645,-0.9107667
"Trochiliformes",-0.701743108,2.4540648
"Musophagiformes",-0.496092895,0.1223289
"Strigiformes",-0.417999500,1.3746733
"Columbiformes",-0.429004577,-0.2578424
"Gruiformes",0.642364170,0.5630810
"Ciconiiformes",1.599520049,1.1040882
"Passeriformes",1.010477618,-0.1672313

cargaremos estos datos utilizando el comando:


> terminales <-read.csv ("tips.csv", header=TRUE)
# en donde header=TRUE, significa que le estamos dando en la matriz los nombres de las filas y columnas, nombramos el objeto en R como "terminales".

al ver nuestros datos con el siguiente comando se deben ver asi:


>terminales


                          V1         V2
Struthioniformes -0.292614262 -1.1963108
Tinamiformes     -0.057614850  1.0345608
Craciformes      -0.006225218  0.4230613
Galliformes      -0.475368798 -0.1927969
Anseriformes      0.668049937  0.2761842
Turniciformes     0.769136135 -0.1066864
Piciformes        1.376599742  0.1124578
Galbuliformes    -0.063333022 -0.5463322
Bucerotiformes    0.195140568 -0.1388508
Upupiformes       1.406734039  0.7346157
Trogoniformes     1.509326716  0.8676560
Coraciiformes     0.230836319  0.3352026
Coliiformes      -0.049803844  1.6605007
Cuculiformes     -0.433148495  0.1290520
Psittaciformes    0.544576328  1.0001100
Apodiformes       0.348601645 -0.9107667
Trochiliformes   -0.701743108  2.4540648
Musophagiformes  -0.496092895  0.1223289
Strigiformes     -0.417999500  1.3746733
Columbiformes    -0.429004577 -0.2578424
Gruiformes        0.642364170  0.5630810
Ciconiiformes     1.599520049  1.1040882
Passeriformes     1.010477618 -0.1672313

Adicionalmente debemos incluir los valores de los estados ancestrales de los caracteres (nodos) que estemos trabajando, esto no es necesario ya que la función escrita por Revell, calcula los estados de los caracteres ancestarles utilizando la función "ace"  (Ancestral Character Estimation) del paquete "ape".
Pero supongamos que tenemos los valores de los estados de los caracteres ancestrales en una matriz diferente, en un archivo .txt Tendremos los datos de esta forma, el archivo sera llamado "nodos.txt":


"V1" "V2"
"24" 0.3086038 0.2754068
"25" 0.2238189 0.2228665
"26" 0.1147769 0.1398433
"27" 0.1824843 0.2085591
"28" 0.1370078 0.198521
"29" 0.3489775 0.3004267
"30" 0.3663461 0.3284963
"31" 0.3486656 0.3913058
"32" 0.3755762 0.3750678
"33" 0.4384152 0.3857682
"34" 0.5109196 0.3681905
"35" 0.4838538 0.4084688
"36" 0.31944 0.428996
"37" 0.2847358 0.4490893
"38" 0.276882 0.4722624
"39" 0.262075 0.4817261
"40" 0.2075108 0.5081077
"41" 0.1870265 0.5221627
"42" 0.1223112 0.5389269
"43" 0.3217109 0.4563515
"44" 0.3492107 0.4568929
"45" 0.399463 0.4814211

Y para cargarlos en R solo necesitamos utilizar el comando:
>nodos<-read.table ("nodos.txt", header=TRUE)
# De la misma forma que cargamos la matriz de valores de estado de caracteres de las terminales.


De este modo, tendremos cargados todos nuestros datos necesarios en R, para utilizar la función de gráfica de filomorfoespacios.


>phylomorphospace(prueba,terminales,nodos,label=TRUE)
#el argumento label=TRUE, indica que queremos que las terminales en la grafica lleven las etiquetas.


Es posible tener una lista en forma de vector con los colores que queremos que utilicen las terminales en el gráfico, si tenemos la lista solo hay que agregar los argumentos





>phylomorphospace(prueba,terminales,nodos,label=TRUE,control=list(X))


Sí no tenemos la información del estado ancestral de los caracteres en los nodos, es decir, la matriz "nodos.txt", la función de filomorfoespacios puede calcular por si misma esta matriz, tomando solo la información del arbol y de la matriz de los terminales ("tips.csv"), para hacer esto solo necesitamos agregar al comando lo siguiente:



>phylomorphospace(prueba,terminales,nodos=NULL,label=TRUE) #con el argumento "nodos=NULL" le estamos diciendo que no poseemos esta matriz y que debe calcularla.


Debido a que nosotros tenemos ambas matrices, utilizaremos el primer comando que describimos y tendremos esl siguiente resultado:

>phylomorphospace(prueba,terminales,nodos,label=TRUE)





Gracias a Liam Revell por su ayuda en la utilización de la función!!
Thanks Liam Revell!


Espero que les guste y les sea util!


martes, 19 de julio de 2011

Test de Mantel en R

En su sentido mas básico el Test de Mantel es un análisis estadístico de correlación entre dos matrices del mismo rango, es usado usualmente para comparar dos matrices de distancia (por ejemplo una matriz de distancia genética entre individuos y una matriz de distancia geográfica entre estos mismos individuos).
Lógicamente en R es muy sencillo realizar este Test, en principio bajo cualquier coeficiente de correlación, pero es comúnmente mayor utilizado el coeficiente de correlación producto-momento de Pearson.


Nathan Mantel (un gordito con cara de buena onda)


Aquí utilizaremos los datos de cantidad de germenes en 4 heladerías!!! (en honor a Juliette fanatica de los helados); Se midió la distancia geográfica entre las heladerías para obtener una matriz y la otra matriz de distancia se obtuvo de los valores de especies compartidas entre pares de heladerías.


los datos los tenemos en dos archivos en formato .txt (aunque podrían estar en cualquier otro formato como .csv) de la siguiente manera:

germen.dist.txt

"hel1" "hel2" "hel3" "hel4" 
"hel1" 0 16.987 107.913 139.049
"hel2" 16.987 0 104.433 129.082 
"hel3" 107.913 104.433 0 62.435
"hel4" 139.049 129.082 62.435 0 


geo.dit.txt



"hel1" "hel2" "hel3" "hel4" 
"hel1" 0 2.54 10.735 14.054
"hel2" 2.54 0 8.254 11.624 
"hel3" 10.735 8.254 0 4.325
"hel4" 14.054 11.624 4.325 0 


Teniendo estos datos podemos empezar a calcular nuestro Test de Mantel.



> gerd <- read.table("germen.dist.txt",header=T)
# llamo mis datos de la matriz de distancia de germenes presentes. Header=T o header=TRUE quiere decir que mis datos poseen los nombres de las heladerias en la primer fila y la primer columna del archivo.

> geod <- read.table("geo.dist.txt",header=T)
#llamo mis datos de la matriz de distancia geográfica entre heladerías.
Header=T o header=TRUE quiere decir que mis datos poseen los nombres de las heladerias en la primer fila y la primer columna del archivo.

> gerd


#visualizo mis datos



        hel1    hel2    hel3    hel4
hel1   0.000  16.987 107.913 139.049
hel2  16.987   0.000 104.433 129.082
hel3 107.913 104.433   0.000  62.435
hel4 139.049 129.082  62.435   0.000

> geod
#visualizo mis datos

       hel1   hel2   hel3   hel4
hel1  0.000  2.540 10.735 14.054
hel2  2.540  0.000  8.254 11.624
hel3 10.735  8.254  0.000  4.325
hel4 14.054 11.624  4.325  0.000

> install.packages ("vegan")
#instalo el paquete vegan.


> library (vegan)
#cargo el paquete vegan a R, para utilizar la función de test de mantel.

>mantel.germenes<-mantel(geod, gerd)
# Realizo el test de mantel para examinar la correlación entre las dos matrices de datos. Sí no se especifica el método o el coeficiente de correlación a utilizar, se utilizara por defecto el coeficiente de correlación producto-momento de Pearson.

>mantel.germenes
#visualizo mis resultados.

Mantel statistic based on Pearson's product-moment correlation 
Call:                                                                                                                                                                                                                            
mantel(xdis = geod, ydis = gerd)                                                                                                                                                                          
      Mantel statistic r: 0.9634                                                                                                                                                                                                       
      Significance: 0.036                                                                                                                 
Empirical upper confidence limits of r:                                                                                                                                                                                          
  90%   95% 97.5%   99%                                                                                                                                                                                                          
0.865 0.931 0.963 0.963                                                                                                                  
Based on 999 permutations

También se pueden utilizar dos coeficientes de correlacion mas (Kendall y Spearman)
mas información. 


Adjunto tambien otra forma de hallar matrices de distancia y hacer un test de mantel, a partir de datos que tengamos:




distclim.csv



estacion,SAL,TEM,O,NH4,N3,N2,P4,S04, CLA,COL,FIT,TUR,ZOT
1,5.0,22.0,2.1,1.6,1.2,2.3,2.2,16.0,1.4,100,121,9,60
2,9.0,20.0,2.2,2.3,1.4,2.5,2.3,12.0,1.4,99,125,6,75
3,11.0,18.0,2.6,0.1,2.0,0.6,1.7,13.0,1.4,89,228,9,62
4,15.0,16.0,3.2,0.0,1.3,0.5,0.7,9.0,0.8,97,245,14,121
5,20.0,18.0,3.1,0.1,0.1,0.0,0.3,6.0,0.8,97,258,14,159
6,25.0,15.0,3.5,0.1,0.5,0.2,0.4,4.0,0.6,75,287,15,186
7,24.0,18.0,4.2,0.1,0.3,0.6,0.4,3.0,1.2,72,175,16,98
8,6.0,19.0,2.7,0.1,0.6,2.5,1.4,13.0,1.6,70,103,22,421
9,10.0,19.0,3.1,0.1,1.1,1.8,1.6,9.0,1.8,85,99,21,369
10,15.0,21.0,4.9,0.1,1.3,1.6,0.6,7.0,1.0,72,145,19,251
11,18.0,18.0,4.6,0.1,2.0,1.5,0.3,5.0,0.6,56,158,16,187
12,22.0,17.0,5.2,0.1,2.0,1.2,0.3,3.0,0.4,55,167,12,111
13,24.0,16.0,4.8,0.1,2.1,1.1,0.2,1.0,0.4,42,262,10,105
14,31.0,16.0,4.2,0.1,1.7,0.2,0.2,6.0,0.6,49,112,8,192
15,30.0,16.0,4.7,0.1,1.5,0.2,0.2,8.0,0.8,23,156,8,214
16,32.0,15.0,4.5,0.5,1.3,0.2,0.3,10.0,0.6,26,128,10,227
17,33.0,16.0,5.1,0.4,1.0,0.4,0.4,6.0,0.2,15,137,12,269
18,35.0,16.0,5.2,0.3,0.7,0.8,0.3,1.0,0.2,15,139,15,301
19,35.0,15.0,4.9,0.1,0.6,0.7,0.2,3.0,0.2,13,241,18,441
20,35.0,16.0,5.1,0.1,0.9,0.3,0.2,2.0,0.3,19,258,19,379

distesp.csv


estacion,AC,OI,SP,EA,NA,PA,TY,FR,OA,CD,CO,CA
1,15,0,0,5,12,0,0,0,0,0,0,0
2,22,3,0,3,14,0,0,0,0,0,0,0
3,37,5,0,8,17,0,0,0,0,0,0,0
4,21,3,0,1,9,0,0,0,0,0,1,0
5,14,12,1,8,16,0,2,3,0,0,0,0
6,21,6,0,14,6,6,3,5,0,2,2,2
7,19,8,0,12,11,5,1,8,0,3,0,3
8,9,3,5,3,9,11,4,0,0,0,0,5
9,15,0,9,12,7,15,10,0,0,0,0,4
10,11,0,8,9,8,18,15,0,0,0,0,7
11,8,0,14,6,5,17,6,0,0,0,0,8
12,13,0,11,5,6,9,14,0,0,0,0,3
13,24,1,9,7,4,3,5,10,0,0,0,9
14,8,2,13,4,2,5,6,12,0,0,0,1
15,6,2,18,12,7,1,9,14,0,26,0,9
16,0,8,16,13,18,2,3,13,0,23,10,14
17,0,9,14,15,14,0,0,20,18,19,2,10
18,3,10,35,18,20,0,0,18,15,17,15,15
19,0,9,22,11,14,0,0,22,19,13,13,18
20,0,7,24,15,11,0,0,17,23,18,9,17



                                             
>  distspe <- read.csv("distesp.csv",header=T) 
# (llamo mis datos de las especies en las estaciones)

>  distclima <- read.csv("distclim.csv",header=T) 
# (llamo mis datos de las observaciones climaticas por estaciones)

>  mdclima <- dist(cbind(distclima$SAL, distclima$TEM, distclima$O, distclima$NH4, distclima$N3, distclima$N2, distclima$P4, distclima$SO4, distclima$CLA, distclima$CO2, distclima$FIT, distclima$TUR, distclima$ZOT)) 
# (creo mi matriz de distancias a partir de mis datos climaticos, esta formula es la via larga en donde se escojen las variables que van a ser tomadas de los datos iniciales con el comando "cbind", lo unico que se necesita hacer es escojer las variables que deso incluir para crear mi matriz)

> mdspe <- dist(cbind(distspe$AC,distspe$OI, distspe$SP, distspe$EA, distspe$NA., distspe$PA, distspe$TY, distspe$FR, distspe$OA, distspe$CD, distspe$CO, distspe$CA )) 
# (realizo otra matriz de distamcia de mis datos de especies por estaciones por la via larga escojiendo mis variables con "cbind")

>mdclima <- dist(distclima)
 # (si tengo pereza y ademas quiero incluir todas las variables de mis datos iniciales de clima, solo utilizo la via rapida para crear mi matriz de distancias)

> mdspe <- dist(distspe) 
# (lo mismo para mis datos de especie)
 
> mantelhilton<-mantel.rtest(mdspe, mdclima, nrepet = 9999) 
# (por ultimo realizo el test de mantel con 9999 replicas simuladas) 


Hasta la vista!!!