Category: Spanish

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

Yo tengo un PhViD

Luego de instalar y conocer los rudimentos de R y de instalar los componentes necesarios, podemos utilizar el paquete “PhViD” de R para detectar señales.Es conveniente primero, si no lo han hecho ya, revisar las nociones básicas de tablas de 2 x 2 y …

Yo tengo un PhViD

Luego de instalar y conocer los rudimentos de R y de instalar los componentes necesarios, podemos utilizar el paquete “PhViD” de R para detectar señales.
Es conveniente primero, si no lo han hecho ya, revisar las nociones básicas de tablas de 2 x 2 y el concepto de independencia en las mismas.
Bueno, voy a desarrollar el uso básico de los algoritmos para detección de señales en bases de datos de Eventos Adversos.

Lo primero que necesitamos es tener una tabla con datos para analizar. Las tablas que se utilizan habitualmente son sumamente extensas y esto hace que sea difícil manipular los datos de una forma entendible en primera instancia. Por eso, utilizaré una tabla pequeña con datos ficticios. La tabla puede descargarse del siguiente link: ListitaEventosDroga.
Descargamos y guardamos la tabla en cualquier carpeta, recordando dónde y vamos a importarla a R.
La importación y exportación de datos a y desde R es una de las tareas menos intuitivas en este software, pero con algo de práctica se le toma la mano.
Para hacerlo, usaremos la instrucción “read.csv”, la cual permite leer un archivo del tipo “csv” que consta de un archivo de texto común y que tiene los datos de las columnas separados por algún carácter (generalmente una coma, por eso “csv” quiere decir “comma separated variable”) y una línea para cada fila de la tabla.
El método read.csv recibe algunos parámetros. El más importante es el que indica la ruta del archivo a importar. Aquí, lo más importante es entrecomillar la ruta del archivo y que, a diferencia de la forma de indicar la ruta en Windows, se usa el carácter “/” en lugar de “\” para separar los nombres de las carpetas.
Los otros parámetros que usaremos serán “header = TRUE” que indica que nuestro archivo tiene los títulos de las columnas, y “sep = “,””, indicando que el separador de columnas es la coma. Finalmente, hay que saber que convertiremos la tabla en un objeto que sea manipulable por el nombre del mismo, asignando la instrucción para importar la tabla a una variable.
La estructura del comando debería seguir este patrón:
variableObjeto = read.csv(“rutaArchivo”, header=TRUE, sep = “,”)
donde “variableObjeto” será reemplazado por el nombre del objeto a crear (el mismo es totalmente arbitrario y pueden elegir el que quieran, siempre y cuando lo recuerden para luego poderlo invocar) y “rutaArchivo” por la ruta para llegar al archivo alojado en el disco. En mi caso la instrucción será la siguiente:
x1=read.csv(“C:/Users/Francisco/Documents/Repositorios/00-PhVid/listitaeventodroga.csv”, header = TRUE, sep = “;”)
Ahora podemos visualizar la tabla con solo escribir el nombre del objeto, como en el siguiente ejemplo:
> x1
        Droga    Evento  n
1 lamotrigina   prurito 10
2   omeprazol   prurito  1
3 lamotrigina       tos  1
4   omeprazol   diarrea  1
5    aspirina sialorrea  1
6  ibuprofeno depresion  1
Vemos que R muestra una tabla con 4 columnas. La primera de ellas es solamente una etiqueta con el número de línea, colocado por R, pero no tiene ningún valor durante el análisis a realizar. De tal modo, la primer columna útil es la que indica el nombre de la droga, la siguiente la que indica el evento ocurrido, y la última indicando el número de reportes para esa combinación Droga-Evento.
Todas las tablas que tengamos que analizar con el paquete “PhViD” deberán tener esta estructura, respetando esas 3 columnas. El número de asociaciones Droga-Evento no tiene limitaciones (a excepción de la memoria del sistema que estemos usando).
Como se puede ver fácilmente, la relación entre los pares Droga-Evento se halla intencionalmente exagerada, resaltando la ocurrencia de uno de ellos (lamotrigina-prurito). Es conveniente repetir que los datos son totalmente inventados a los fines de este tutorial.
Bueno, como expliqué previamente, al iniciar R se cargan solamente los paquetes básicos. Para dotar a R de las capacidades que necesitamos hay que cargar este paquete. Escribimos la siguiente instrucción:
> library(PhViD)
Así el paquete se carga en el espacio de trabajo actual y podemos utilizar sus métodos para manipular los objetos creados.

Si les llega a dar algún error avisando que no puede instalar “LBE”, tipean lo siguiente:

source("http://bioconductor.org/biocLite.R")
biocLite("LBE")

Y luego vuelven a tipear “library(PhViD)”. Con eso debería solucionarse.

Si bien antes dije que la tabla debe tener una estructura determinada, PhViD necesita cambiar el formato que le entregamos para realizar el análisis que queremos. Para esto, tiene el método “as.PhViD”. Como siempre, asignamos el resultado del método a un nuevo objeto, con la siguiente instrucción:
> l1=as.PhViD(x1)
Ahora el objeto “l1” contiene mis datos en un formato que es analizable mediante los métodos que aporta PhViD.
Voy a mostrar los pasos usando la metodología usada por el Uppsala Monitoring Centre, con el algoritmo llamado “Bayes Confidence Propagation Neural Network” (BCPNN). El método a usar es justamente el “BCPNN”. Los mismos pasos pueden usarse para aplicar cualquiera de los otros algoritmos disponibles en el paquete.
Crearemos un objeto usando el método BCPNN. Podemos nombrar al objeto de cualquier manera, yo lo llamaré “bl1.1”. Para ello simplemente pasamos el objeto “l1” previamente creado, por el método BCPNN asignándolo a la variable nueva.
> bl1.1=BCPNN(l1)
El nuevo objeto creado tiene varios miembros capaces de entregarnos los resultados obtenidos. Si queremos ver los miembros de bl1.1 usamos el método “names” de la siguiente manera.
> names(bl1.1)
[1] "INPUT.PARAM" "ALLSIGNALS"  "SIGNALS"     "NB.SIGNALS"
Los miembros que nos interesan son “SIGNALS” y “ALLSIGNALS”, que nos entregan la lista de pares detectados como señales positivas ordenados por ranking y la lista ordenada por ranking de los valores obtenidos para todos los pares, respectivamente.
> bl1.1$SIGNALS
Nos muestra la tabla con las señales positivas:

par
drug code
event effect
count
expected count
post.H0
n11/E
drug margin
event margin
FDR
FNR
Se
Sp
lamotrigina prurito
lamotrigina
prurito
10
8.066666667
0.180741235
1.239669421
11
11
0.180741235
0.735913946
0.222650697
1

Que en nuestro caso consta solamente del par “lamotrigina prurito”. 

A su vez: 
> bl1.1$ALLSIGNALS

par
drug code
event effect
count
expected count
post.H0
n11/E
drug margin
event margin
FDR
FNR
Se
Sp
lamotrigina prurito
lamotrigina
prurito
10
8.066666667
0.180741235
1.239669421
11
11
0.180741235
0.735913946
0.222650697
1
aspirina sialorrea
aspirina
sialorrea
1
0.066666667
0.295550161
15
1
1
0.238145698
0.715077741
0.414099668
0.922108741
ibuprofeno depresion
ibuprofeno
depresion
1
0.066666667
0.295550161
15
1
1
0.257280519
0.718620375
0.605548639
0.794740052
omeprazol diarrea
omeprazol
diarrea
1
0.133333333
0.330254358
7.5
2
1
0.275523979
0.725705643
0.787566019
0.667371363
lamotrigina tos
lamotrigina
tos
1
0.733333333
0.541284471
1.363636364
11
1
0.328676077
0.781665645
0.912231554
0.525046743
omeprazol prurito
omeprazol
prurito
1
1.466666667
0.677049884
0.681818182
2
11
0.386738379
Inf
1
0.29177773

Nos da la tabla con todos los pares ordenados por intensidad de señal. 

Obviamente todos estos números son difíciles de leer, pero analizaremos alguna de las columnas y les mostraré, que en el corazón del método residen mecanismos que ya conocemos. 
En primer lugar, la primera columna contiene las etiquetas con los pares combinados. La segunda y tercera columnas son los nombres de las drogas y eventos, respectivamente. La tercera columna es solamente el número de veces en las que aparece esa combinación. Hasta aquí, salvo por el orden, no es más que lo que se hallaba contenido en la tabla original. 
La columna “expected count” corresponde a lo que conocemos como “valor esperado” en las tablas de 2 x 2 desarrollado previamente. Para verlo más fácil, construiremos la tabla para el par “lamotrigina prurito”. Necesitamos las marginales y el total general. Sabemos que el total de reportes es de 15. El número de reportes para ese par es de 10. Continuamos colocando el número de reportes que contienen “prurito” y el número de reportes que  contienen “lamotrigina”. La tabla quedaría así:


lamotrigina
si
no
prurito
si
10
11
no
11
15

Podemos ver que estos dos últimos valores son los que están presentes en las columnas “drug margin” y “event margin”. Es interesante que con sólo conocer éstos valores, podemos calcular el resto sin necesidad de volver a los datos y es por esto que BCPNN solamente guarda estos valores. La tabla completa sería:

lamotrigina
si
no
prurito
si
10
1
11
no
1
3
4
11
4
15

Sabemos que estos son los “valores observados” y que podemos calcular a través de ellos los “valores esperados”. Esta vez lo haremos sin redondear los números. La tabla es:

lamotrigina
si
no
prurito
si
8.06666667
2.93333333
11
no
2.93333333
1.06666667
4
11
4
15

Chan! Resulta que el valor esperado para la celda correspondiente al par “lamotrigina prurito” es el mismo que el de la columna expected count”. Exactamente lo mismo obtendremos calculando las tablas correspondientes a cada par. 
Agregaré que la columna “n11/E” no es nada más que el valor obtenido dividiendo “count” sobre “expected.count”. El resto de las columnas serán analizadas en entregas posteriores. 

Con esto ya podemos usar el paquete PhViD para generación de señales. Claro que otro tema es que hacer luego de tenerlas, ¿no? 
Bueno, los que se animen, a probar. Espero sus comentarios.