miércoles, 5 de octubre de 2011

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.

No hay comentarios:

Publicar un comentario