Category: Spanish

Si es Bayes, es bueno…

Los métodos de inferencia vistos hasta el momento forman parte de lo que se da en llamar “inferencia frecuentista”. El modelo subyacente es que, verificando ciertos supuestos, es posible ajustar una distribución (normal u otra) a mis datos para inferir conclusiones acerca de la población.
La inferencia bayesiana tiene otra visión. En general parte de que la distribución “real” de la población es algo imposible de conocer. En lugar de tomar una distribución como la verdadera, se toma lo que conoce “a priori” el investigador para elegir una distribución que se ajuste a ello y trabajar inicialmente bajo ese supuesto. Luego se toma una muestra de la población y con los datos obtenidos de ella se realiza un “ajuste” de la distribución empírica, llevando a una distribución que podría ser diferente a la anterior. Es lo que se conoce como distribución “a posteriori”.
O sea que, proponemos un modelo de acuerdo a lo que conocemos, tomamos evidencia para contrastarla contra el modelo y corregimos el modelo de acuerdo a lo aprendido.
Para tener un ejemplo intentaré estimar cuál es la proporción de estudiantes de una universidad que duerme al menos 8 horas.
Este ejemplo está tomado por completo del libro “Bayesian computation with R” de Jim Albert.
Supongamos que un investigador está interesado en conocer esa proporción. El verdadero valor de la proporción le es desconocido.
De acuerdo a la información de otras publicaciones y experiencia previa, piensa 2 cosas:
·         La mejor predicción que puede hacer es de 30%
·         Es muy probable que se halle por debajo del 50%
Construimos un rango de valores posibles para la proporción. Llamaremos “p” a la proporción. En R escribimos:
> p=seq(0.05,0.95,by=0.1)
> p
 [1] 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95
Esta instrucción construye un rango entre 0.05 (5%) y 0.95 (95%) que aumenta de a 0.1 (10%). Es un rango de valores posibles de p que nos servirá para construir una distribución a priori.
Ahora a cada uno de estos valores le otorgamos un puntaje entre 0 y 10 de acuerdo a lo que pensamos (esta escala es arbitraria). Estos puntajes son llamados “pesos”.
> pesos = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
> pesos
[1] 1.0 5.2 8.0 7.2 4.6 2.1 0.7 0.1 0.0 0.0
Y convertimos esos pesos a un formato de probabilidad. Para hacerlo, dividimos cada peso por la suma total de los mismos. A estas probabilidades se las llama “priors” o probabilidades “a priori”.
> prior = pesos/sum(pesos)
> prior
Esto nos permite crear una distribución que indica una probabilidad para cada valor de p. Para visualizarla podemos usar un gráfico:
> plot(p, prior, type = "h")
El gráfico obtenido nos muestra nuestra distribución “a priori”. Ahora la tarea será obtener una muestra de estudiantes. Lo hacemos y obtenemos datos de 27 estudiantes. Luego de revisar las encuestas notamos que hay 11 estudiantes que duermen más de 8 horas. Esto nos da un número de “éxitos” (11) y “fracasos” (16). Colocamos esos valores dentro de un objeto:
> data=c(11,16)
> data
[1] 11 16
El uso de la instrucción “c” es para concatenar varios valores y en este caso, asignarlos a un solo objeto.
Ahora hay que cargar un paquete de R llamado “LearnBayes”. Por supuesto, tiene que estar previamente instalado.
> library("LearnBayes")
Si no estuviese instalado nos daría un error diciendo que no encuentra el paquete. Procedemos como en el post que explica sobre instalación de paquetes de R y lo instalamos. Luego volvemos a tipear la instrucción anterior.
Luego de cargado el paquete escribimos el siguiente comando:
> post = pdisc(p, prior, data)
> post
Esto nos da las probabilidades a posteriori. Es bastante claro que son diferentes a las probabilidades a priori.
Para visualizarlas en forma comparativa escribimos una serie de instrucciones que no explicaré en detalle, pero nos permiten ver una gráfica sobre la otra. En caso de que a alguien le interese la explicación detallada me lo puede pedir. No obstante, se puede obtener ayuda de R escribiendo “?” seguido de la instrucción, o help(“instrucción”).
> library(lattice)
> PRIOR=data.frame("prior",p,prior)
> POST=data.frame("posterior",p,post)
> names(PRIOR)=c("Type","P","Probability")
> names(POST)=c("Type","P","Probability")
> data=rbind(PRIOR,POST)
>xyplot(Probability~P|Type,data=data,layout=c(1,2),type="h",lwd=3,col="black")
Esto nos muestra que las 2 distribuciones, la a priori construida con las suposiciones del investigador y la a posteriori, calculada usando la a priori y corrigiendo de acuerdo a las observaciones, son diferentes.
Existen funciones matemáticas para calcular las distribuciones a posteriori a partir de las a priori. Estas funciones son diferentes en cada situación y mi objetivo en este post no es ser extenso en éste tema.
Baste decir, que con esta metodología es posible estimar distribuciones pero también parámetros, como medias, sd y otros.
Con estos conocimientos podemos ir, ahora sí, a continuar con el aprendizaje del paquete PhViD para detección de señales.

Si es Bayes, es bueno…

Los métodos de inferencia vistos hasta el momento formanparte de lo que se da en llamar “inferencia frecuentista”. El modelo subyacentees que, verificando ciertos supuestos, es posible ajustar una distribución(normal u otra) a mis datos para infe…

Yo infiero, tu infieres

Luego de repasar algunos rudimentos de distribuciones estadísticas, es conveniente ver algo de inferencia estadística. Es la parte de la estadística que permite deducir propiedades de una población, a partir de una muestra de la misma.
Acá hay algunos conceptos que habitualmente conducen a error. Básicamente, como decíamos antes, cuando estudiamos una muestra una de las cosas que hay que realizar es una serie de pruebas para establecer si la misma proviene de una población donde la variable tiene una distribución normal. Esto implica que existe una población teórica (e imposible de aprehender) desde la cual se toma la muestra. En esta muestra verificaremos propiedades similares a la de la población, pero no idénticas. Por ejemplo, si calculamos la media en la muestra de glucemias con la que trabajamos anteriormente veremos que:
> mean(glucemias)
[1] 90.03842
Su promedio no es exactamente 90, que es lo esperado conociendo la media real, pero está muy cerca. Esto nos puede llevar a una pregunta. Si tomamos muchas muestras de la misma población y calculamos sus promedios, ¿cuánto se alejarán éstos de la media poblacional? Es aquí que surge la idea conceptualmente diferente de la que hablaba antes. Existe una distribución que se puede ajustar a los valores de nuestra muestra para describirla, pero lo que nos interesa en realidad, es una distribución teórica que muestra cómo se distribuirían los promedios si tomáramos muchas muestras respecto de la media poblacional. Inicialmente, es fácil imaginar que la mayoría de los promedios estarán cerca de la media y si bien obtendremos por azar valores alejados, éstos se irán haciendo menos y menos probables a medida que nos alejamos de ella. Esto es muy importante, existe una distribución teórica de medias muestrales alrededor de la media poblacional. Esta distribución tiene a su vez una media (la media poblacional o “verdadera” media) y un desvío estándar, que en general se denotan usando las letras griegas µ y σ para diferenciarlas y resaltar que en general son valores teóricos desconocidos y a estimar.
Éste es el razonamiento básico subyacente en lo que solemos ver como Intervalo de Confianza (ojo que su acrónimo en español “IC” puede confundirse con el de Information Component que nos interesa, pero no es lo mismo).
Si tuviéramos información acerca de la media poblacional y su desvío estándar, podríamos calcular que tan probable es obtener un promedio determinado para una muestra tomada. Como no lo tenemos, lo inferimos. Quiero decir, vamos al revés. En lugar de ir de la información de la población para evaluar la muestra, vamos de la información que hay en la muestra para evaluar cómo será la población.
Si tuviéramos los datos de la población, calcularíamos la media y el desvío estándar de ella, mediríamos a qué distancia en términos de sd se encuentra la media de nuestra muestra y buscaríamos la probabilidad de encontrar un valor así de alejado de la media poblacional. Como no conocemos la media y el sd poblacional, usamos la media muestral para estimar la poblacional y para estimar el sd poblacional hacemos una corrección del sd muestral. Esta corrección hará que el sd sea más estrecho, como sería de esperar si tuviéramos toda la población y obviamente, está relacionada al tamaño de la muestra (n).
 Este nuevo valor del sd es llamado Error Standard (se) y se trata del valor que usamos para estimar el sd poblacional. Entonces con esto podemos calcular un intervalo de confianza para una muestra. SI luego del análisis preliminar concluimos que nuestra muestra proviene de una población con distribución de valores normal, podemos usar aproximadamente 2 errores estándar (sabemos por lo visto en el post anterior que el valor exacto es de 1.959964 sd) hacia la derecha y 2 errores estándar hacia la izquierda del promedio para establecer un intervalo de confianza de 95%. Entre las características que podemos usar para verificar si se cumple el principio de normalidad, uno muy importante es el tamaño de la muestra. Cuanto mayor el tamaño muestral, menos dependemos del resto de los supuestos.
Este intervalo de confianza tiene el siguiente significado: si tomáramos muchas muestras con el mismo método que tomamos la primera, el 95% de las medias se encontraría en ese rango y además, lo que es más importante, en ese rango se encuentra la media poblacional con un 95% de probabilidad.
Veamos esto con las glucemias. Usaremos la función “mean” y “length” ya vistas, la función “sd” que calcula el desvío estándar y la función “sqrt” que nos da la raíz cuadrada.
> nG=length(glucemias)
> meanG=mean(glucemias)
> sdG=sd(glucemias)
> seG=sdG/sqrt(nG)
> sdG
[1] 9.519673
> seG
[1] 0.9519673
Vemos que el sd es de 9.52, cercano a 10 que es el valor poblacional, pero el se es mucho menor, cercano a 0.92. Para obtener el intervalo de confianza acerca de nuestra media hacemos lo siguiente:
> LimInf=meanG-1.959964*seG
> LimSup=meanG+1.959964*seG
> LimInf
[1] 88.1726
> LimSup
[1] 91.90424
El límite inferior es de 88.17 y el superior de 91.90. Como dijimos antes, esto quiere decir entre otras cosas que la media poblacional “real” se encuentra en ese rango con un 95% de probabilidad. Sabemos que esto es lo que ocurre en la realidad, pues la media poblacional es de 90 mg/dl, según lo predeterminado por la simulación con la que generamos los datos muestrales.
También, nos dice otra cosa muy importante: Si en otro momento tomamos otra muestra y al calcular el promedio encontramos que se encuentra fuera de ese rango, podemos tener la “confianza” de que es improbable que esos valores provengan de la misma población. Mejor dicho, la probabilidad de que provengan de la misma población es inferior al 5%.
Bueno, hasta aquí los conceptos de estadística “clásica” que son importantes para entender lo que viene de generación de señales con el paquete PhViD.
No obstante, para entenderlo más acabadamente hay que tener nociones de inferencia Bayesiana.

Yo infiero, tu infieres

Luego de repasar algunos rudimentos de distribuciones estadísticas, esconveniente ver algo de inferencia estadística. Es la parte de la estadísticaque permite deducir propiedades de una población, a partir de una muestra de lamisma.Acá hay algunos…

Se me caen las medias!

Voya hablar un poco acerca de distribuciones de probabilidad. Este es un temaextenso así que solamente trataré los puntos que me parecen necesarios paraentender los demás.Comovimos en el post sobre nociones de independencia, en las tablas de 2 x 2 e…

Se me caen las medias!

Voy a hablar un poco acerca de distribuciones de probabilidad. Este es un tema extenso así que solamente trataré los puntos que me parecen necesarios para entender los demás.
Como vimos en el post sobre nociones de independencia, en las tablas de 2 x 2 existe un valor observado para cada celda interna, que es el que anotamos a partir de las observaciones, y un valor esperado que podemos calcular y es a su vez el valor más probable bajo la hipótesis de que las variables son independientes.
Entonces, queda configurada una “distancia” entre el valor esperado y el observado. Intuitivamente puede verse que a mayor diferencia del valor observado respecto del esperado, menos probable se hace el mismo si las variables fueran independientes y por lo tanto se hace cada vez más difícil aceptar esta hipótesis. Cuando la distancia es demasiada (y para ello se define un punto de corte), rechazamos la hipótesis de independencia y aceptamos que hay asociación entre las variables, pues de otro modo sería sumamente improbable obtener esa configuración de celdas.
Ahora bien: ¿Cómo se mide esa distancia y cómo se hace para saber el valor de probabilidad de obtener esa diferencia? Aquí entran las distribuciones de probabilidad. Usándolas, podemos traducir entre distancia medida en términos comunes y probabilidad asociada a la misma.
Para tablas de 2 x 2, en general se usa la distribución de chi cuadrado, porque las propiedades que tienen estas tablas establecen que este método sea adecuado la mayor parte de las veces. Pero otras veces, lo adecuado es elegir otras distribuciones para medir la distancia observada y eso depende de la situación a la que nos enfrentemos.
Para dar una noción más intuitiva del uso de las distribuciones de probabilidad, voy a ejemplificar con la distribución llamada Normal que es posiblemente la más escuchada.
Es de esperar que la mayoría de nosotros haya visto alguna vez una de esas reglas que sirven para medir en diferentes escalas. Por ejemplo, una regla que de un lado tiene una escala en centímetros y del otro en pulgadas.
Una regla así, nos podría servir como método para traducir entre una escala y otra. Es fácil ver en la imagen que si medimos algo que mide 1 pulgada, medirá aproximadamente 2,5 centímetros. Algo así pasa con las distribuciones de probabilidad, vale decir que podríamos medir una distancia en centímetros y traducir a probabilidades. Para eso, vamos a usar una herramienta que se consigue en cualquier ferretería estadística. 😉
Esta regla tendrá, al menos inicialmente, una escala en centímetros. Pero posee algunas características particulares. En primer lugar, el valor 0 no está colocado en un extremo sino en el centro de la misma. En segundo lugar, los centímetros son positivos hacia la derecha del cero pero se hacen cada vez más negativos hacia la izquierda del mismo.
Como esta no es una regla cualquiera, tiene además la particularidad de que se le pueden acoplar diferentes escalas en el borde opuesto al de los números, de forma de poder traducir de escala en centímetros a cualquier otra.
En esta oportunidad, le acoplaremos una curva de probabilidad Normal:
Como ven, esta nueva escala es bastante rara. Básicamente, a diferencia de lo que ocurre con la escala en centímetros, su grosor es máximo en el centro o valor 0 y éste se va adelgazando a medida que nos alejamos del mismo hasta hacerse mínimo más allá de los 3 centímetros de distancia respecto del centro. Obviamente, el área cubierta por la figura en cada segmento también disminuye. También podemos ver que la figura es simétrica hacia ambos lados.
Lo interesante (como veremos) es que si asumimos que el área total de la figura es 100%, su forma es tal que es posible saber exactamente qué porcentaje del área queda encerrada en cada segmento medido en centímetros. En primer lugar, es fácil darse cuenta que cualquiera de las 2 mitades contiene un área de 50%. Entre 0 y 1 cm queda determinada un área de 34% del total y dada su simetría lo mismo ocurre entre 0 y -1 cm. Por lo tanto, entre -1 y 1 cm queda un área de 68%.
Varios rangos que podríamos verificar serían:
Rango en cm
Porcentaje
del área
0 a 1
34%
-1 a 1
68%
0 a 2
47,7%
-2 a 2
95,4%
0 a 3
49,9%
-3 a 3
99,7%
-1 a 3
84%
Todo el rango menor o igual a cero
50 %
Todo el rango menor o igual a 2
97,7%
La verdad es que si bien estuvimos hablando en términos de porcentajes, también estuvimos hablando en términos de probabilidades. El principio fundamental es que no debemos nunca usar un destornillador para sacar un clavo. Vale decir, cada herramienta tiene una situación de aplicación para lo cual fue diseñada.
Así, lo más importante antes de hacer una medición con la herramienta que mostramos, es saber si nuestros datos pueden ser medidos con la misma o debe usarse otra. Sabemos por lo que leímos en muchas publicaciones, que muchos valores en medicina tendrían una gráfica de distribución normal y por lo tanto simétrica alrededor de la media. (La palabra “tendría” es condicional aquí, pues hay situaciones en las que esto no se cumple).
Vamos a usar R para simular datos con distribución normal y reforzar lo que venimos hablando.  Pretendamos que estamos midiendo valores de glucemia en una población sana en mg/dl y que el rango normal va de 70 a 110 mg/dl. Usaremos el valor intermedio del rango como media (90 mg/dl) y un desvío estándar de 10 mg/dl.
La función “rnorm” sirve para generar números aleatorios según la distribución normal y toma como argumentos la media y el desvío estándar (sd), además del número de observaciones o “casos” que queremos. Pediremos 100 casos:
> n=100
> media=90
> sd=10
> glucemias=rnorm(n,media,sd)
Luego, graficaremos la densidad de distribución de nuestros datos de la siguiente manera:
plot(density(glucemias))
Obtendremos una gráfica más o menos así:

Posiblemente la gráfica sea medio torcida, pero vemos que es bastante parecida a la normal (si en lugar de un n de 100 usáramos 1000 o 10000, veríamos que cada vez es más parecida).
Como suponemos que es aceptable usar la distribución normal para estudiar los datos (aquí este supuesto es obviamente correcto pues los datos fueron generados por simulación) mediremos la probabilidad de obtener valores dentro de determinados rangos y los compararemos con la tabla de probabilidad normal aportada más arriba.
Para obtener los valores dentro de un determinado rango, tomaremos nuestro objeto “glucemias” que contiene nuestros datos y lo “cortaremos” según necesidad. Para tomar segmentos de una lista de valores contenida en un objeto, se usa el nombre del objeto seguido de corchetes (“[ ]”), y dentro de los corchetes irá la instrucción que indica dónde cortaremos.
> glucemias[glucemias<=90]
 
Esto nos dará la lista de valores de glucemia que tienen valor menor o igual a nuestra media, que es el valor de 90 mg/dl. Como tenemos 100 casos totales, el número de casos en un segmento será también el porcentaje del total de casos contenido en ese segmento. Para saber el número de casos que hay en la lista de valores sin tener que contarlos a mano usamos la instrucción “length”, envolviendo a la instrucción anterior de la siguiente manera:
 
> length(glucemias[glucemias<=90])
 
Rango en sd
Porcentaje
del área esperada
Cantidad obtenida en el rango
Instrucción
0 a 1
34%
31
length(glucemias[glucemias<=100])-
length(glucemias[glucemias<=90])
-1 a 1
68%
64
length(glucemias[glucemias<=100])-
length(glucemias[glucemias<=80])
0 a 2
47,7%
49
length(glucemias[glucemias<=110])-
length(glucemias[glucemias<=90])
-2 a 2
95,4%
99
length(glucemias[glucemias<=110])-
length(glucemias[glucemias<=70])
0 a 3
49,9%
49
length(glucemias[glucemias<=120])-
length(glucemias[glucemias<=90])
-3 a 3
99,7%
100
length(glucemias[glucemias<=120])-
length(glucemias[glucemias<=60])
-1 a 3
84%
82
length(glucemias[glucemias<=120])-
length(glucemias[glucemias<=80])
Todo el rango menor o igual a la media
50 %
51
length(glucemias[glucemias<=90])
Todo el rango menor o igual a 2 sd
97,7%
100
length(glucemias[glucemias<=110])
Es interesante verificar que al obtener una muestra de una población, a pesar de tener la variable medida una distribución expresamente normal en esa población, su distribución en la muestra será aproximada pero no exactamente normal. Habitualmente el procedimiento será inverso, se tomará una muestra de una población y deberán realizarse algunas evaluaciones para determinar si se cumplen los supuestos que indican que la distribución normal es pasible de ser aplicada sin incurrir en violaciones de los mismos que, en general, conducirían a resultados y conclusiones erróneas.

¿Pero, para qué utilicé la analogía de los centímetros? Y ¿Por qué en una tabla puse los rangos en centímetros y en otra en sd (desvíos estándar)? Existe una distribución normal especial, denominada “Normal(0,1)”, en donde la media es igual a cero y el desvío estándar igual a 1. Si bien cualquier probabilidad puede ser calculada computacionalmente sin recurrir al uso de tablas, es útil saber que las tablas de probabilidad normal que a veces se ven publicadas y algunos tuvimos que usar alguna vez, están expresadas según la Normal(0,1). Aunque tengamos una media y sd diferentes, cualquier valor puede ser fácilmente convertido a esa escala mediante una cuenta en donde restamos el valor de la media y luego lo dividimos por el desvío estándar:
Así, sabemos a cuántos desvíos estándar está nuestro valor, respecto de la media y podemos buscarlo en las tablas. En R podemos usar el comando “pnorm” con el resultado de ese valor, para obtener la probabilidad de obtenerlo. El resultado será evaluado como si se tratara de uno obtenido de una Normal(0,1) (salvo que indiquemos otros argumentos). Por ejemplo, si la cuenta anterior diera un resultado de 2 (sería el caso si con las glucemias usáramos el valor de 110 mg/dl), sabríamos que este valor está a 2 sd respecto de la media. Usamos “pnorm” para establecer la probabilidad de obtener un valor así, tomado de una distribución normal:
 
> pnorm(2)
[1] 0.9772499
Tal vez hubiéramos esperado que el valor estuviera cerca del 95% (a esta altura la mayoría se imaginará que cuando hablo de 95% y 0.95 me refiero a lo mismo. En realidad las probabilidades se expresan en valores de 0 a 1. Para transformarlas en porcentaje se multiplican por 100), de acuerdo a los conocimientos que suelen transmitirnos. En realidad, la función da la probabilidad de obtener un valor igual o menor a 2, no entre ±2sd o lo que es lo mismo, entre 2 y -2. Si queremos obtener esto último, debemos restarle la probabilidad de obtener un valor menor o igual a -2, de la siguiente manera:
> pnorm(2)-pnorm(-2)
[1] 0.9544997
Como ven, está cercana a 95%, aunque no exacta. Si quisiéramos obtener el valor para ±sd que nos da exactamente el 95% de probabilidad, usamos “qnorm” indicándole el valor 0.975 (o sea 97,5%, que es el valor que deja 2,5% hacia la derecha, pues queremos que nos quede un 2,5% a la izquierda).
> qnorm(0.975)
[1] 1.959964
Ahora volvemos a utilizar “pnorm” con ese valor:
> pnorm(1.959964)-pnorm(-1.959964)
[1] 0.95
Exactamente 95%.

Bueno, acá quiero que respiren un poco, se estiren y despejen un poco la cabeza. Vamos a pasar a un punto importante.

Cuando hacemos cualquier estudio, la mayoría de las veces queremos generalizar el resultado. Vale decir, puede ser que las conclusiones nos interesen solamente para ser aplicadas para esos pacientes y ese momento particular, pero habitualmente queremos que sean útiles para usarlas en otros pacientes. Esto implica que las conclusiones serán utilizables para aplicar a otros pacientes similares en otras regiones, en bases de datos históricas, y en pacientes que aún no han llegado, pero llegarán en el futuro. Esto quiere decir que cuando tomo una muestra, no solamente lo hago porque es impracticable estudiar a toda la población de interés por una cuestión de recursos sino que esa población es una población teórica que la mayoría de las veces no existe al momento de realizar el estudio.
Hasta el momento, con lo que tenemos, solamente podemos describir nuestra muestra en términos de media y desvío estándar. Ahora, si queremos evaluar cuán aplicable es mi resultado a otros pacientes no incluidos en el estudio, debemos usar otro razonamiento y realizar un salto conceptual para entrar en el terreno de la inferencia estadística.
 

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