Veni, ViDi, Vici

Bueno, luego de realizar una pequeña introducción a las distribuciones de probabilidad, inferencia frecuentista e inferencia Bayesiana, necesarias para comprender lo que sigue, continuaré con el análisis de la salida del paquete PhViD para detección de señales en farmacovigilancia.

Teníamos nuestro espacio de trabajo cargado con el paquete PhViD y los datos de la tablita descargada con datos inventados a los fines de poder trabajar. Si no está configurado así, repetir los pasos según sugerí en el post anterior.

Ahora bien, había explicado el significado de las columnas “count”, “expected count”, “drug margin”, “event margin” y “n11/E”. Ahora nos adentraremos en la utilidad de las demás.

En primer lugar, hablaré un poco del Information Component o “IC”. Como vimos antes, la columna “n11/E” nos entrega la razón entre el valor observado y el esperado. El IC se expresa como el logaritmo natural del siguiente cociente: la probabilidad de encontrar el evento y la droga asociados, dividida por la multiplicación de la probabilidad de hallar la droga sola por la probabilidad de hallar el evento solo. Habitualmente se encuentra notada de la siguiente manera:

 Pero puede encontrarse escrito de otras maneras que tienen el mismo significado. El mismo puede calcularse también obteniendo el logaritmo natural del valor de n11/E, que es equivalente a la expresión entre paréntesis en la ecuación anterior.

En este cociente, si la ocurrencia observada en forma conjunta es igual a la esperada, el resultado será igual a 1 y su logaritmo igual a 0. Con este resultado encontraríamos que la asociación entre Droga – Evento no ocurre más frecuentemente que lo esperable y por lo tanto no se generaría una señal. Por el contrario, si es mayor a 1, tendríamos una posible señal.

El algoritmo BCPNN va más allá. Utilizando métodos de inferencia Bayesiana estima un nuevo valor para el IC, tomando como evidencia para reformularlo, los valores encontrados en la tabla. Este valor es el valor esperado del IC luego de tomar evidencia que permite estimar su valor “a posteriori”. También estima el sd “a posteriori”. Con estos valores asume (en la función por defecto, pero puede indicarse que lo haga de otra manera) que la distribución de probabilidad del IC es normal y calcula la probabilidad de obtener un valor igual o menor que 0, que sería el valor obtenido para el IC si no se genera señal alguna. Si esta probabilidad es baja se acepta la hipótesis de que hay asociación entre Droga y Evento. Esta probabilidad esta expresada en la columna postH0.

Si tenemos ganas de ver el corazón del algoritmo, podemos escribir en R:

> BCPNN

Luego de presionar “Enter” aparecerá una ristra de instrucciones que son las que se ejecutan cuando llamamos a la función. La verdad es que la mayoría es bastante críptica para alguien que no se halla acostumbrado a ver este tipo de salidas, pero hay un par de cosas que pueden notarse.

Puede verse que la función calcula dentro de su operatoria un objeto “EICb” que es la media esperada para el IC y supongo que su acrónimo significará algo como “la esperanza del IC luego de aplicar Bayes”. También, a continuación de la anterior se observa un objeto denominado “VICb” que es la varianza del IC. La raíz cuadrada de la varianza es el sd. Aquí se usa sd como equivalente del se porque estamos estimando una distribución teórica que sería equivalente a la “poblacional”, y en ese caso el sd es usado como el estimador de σ.

También, en ese listado se usa el símbolo “<-” en lugar de “=” para asignar el contenido a un objeto. Hay que saber que el primero es el nativo de R. Yo utilizo “=” por comodidad y para mayor entendimiento, pero puede pasar infrecuentemente que alguna función no arroje nada o nos de error si no utilizamos la otra forma.

Las otras columnas importantes son las denominadas “FDR” y “FNR” que son acrónimos de False Discovery Rate y False Negative Rate. Estos son conceptos introducidos por los autores del paquete y pueden verse desarrollados en los papers publicados por ellos^[1].

Básicamente, los métodos para detección de señales son muy sensibles y tienden a producir un gran número de falsos positivos, entonces se genera una variable que mide la probabilidad de obtener un falso positivo (FDR) y la de obtener un falso negativo (FNR).

Si no indicamos nada en contrario el ránking es ordenado según postH0, de menor a mayor. Esto puede cambiarse mediante el argumento “RANKSTAT”. Por defecto toma el valor igual a 1 y ordena por postH0. Acepta también el valor 2 y en este caso, produce un doble efecto: arroja una columna adicional llamada “Q_0.025(log(IC))” que indica el extremo inferior de un intervalo de confianza del 95% para el IC a posteriori, lo segundo, es que ordena el ránking por esa columna.

No obstante el orden, la regla para decidir cuáles son señales y cuáles no, es por defecto el FDR, con un punto de corte de 0.05.

Haremos ahora unas pruebas con el set de datos que viene en el paquete PhViD.

El set de datos es “PhViDdata.frame”. Lo asignaremos a un objeto que tenga un nombre un poco más manejable, por ejemplo:

> lista=PhViDdata.frame

En segundo lugar, convertiremos el objeto “lista” a un objeto en el formato que maneja el paquete:

> data=as.PhViD(lista)

Ahora usamos “BCPNN” para buscar señales, asignando la salida al objeto “Bdata”:

> Bdata=BCPNN(data)

Por último, asignaremos las señales a el objeto “signals”:

> signals=Bdata$SIGNALS

Ahora, dado que la salida será enorme el objeto signals es muy difícil de visualizar en R, y mostraré algunas funciones para poder enterarnos de lo que contiene. En primer lugar la instrucción “head”:

> head(signals)

Nos muestra el contenido de las 6 primeras filas de la tabla.

La instrucción “summary”:

> summary(signals)

Nos muestra resúmenes de las columnas, cuando son numéricos nos da estadísticos descriptivos y cuando son categóricos nos muestra las frecuencias de ocurrencias de los primeros valores.

La instrucción “str”:

> str(signals)

Nos muestra también los primeros valores pero puestos como líneas. Adicionalmente nos dice si una variable es numérica marcándola como “num” o si es categórica, marcándola como “factor”. Otra información que nos da es la dimensión de la tabla, que en este caso es de 19130 líneas x 12 columnas. Esto es bastante útil, porque el número de líneas de “signals” será igual al número de señales producidas. Si hiciéramos lo mismo con el objeto “lista” que contiene la lista de pares Droga – Evento:

> str(lista)

Obtenemos que el número total es de 102483. Vale decir, que obtuvimos:

> 19130/102483*100

[1] 18.66651

Un 18% de pares con señal positiva. Esto es lo que expresaba antes cuando decía que estos métodos daban gran cantidad de señales y que la forma de manejarlos es tomar las más intensas. Por ello, es importante el ránking ordenado de más intensa a menos intensa.

Podemos abrir una ventana con la tabla para navegarla con:

> View(signals)

Posiblemente, la forma más práctica sea exportar la tabla

> write.table(signals,"C:/PhVid/listaExp.csv",quote=F,row.names=F,sep=";")

A continuación, muestro como se genera una lista de señales dándole a “RANKSTAT” el valor de 2, para mostrar cómo se llama la función con este argumento:

> Bdata2=BCPNN(data,RANKSTAT=2)

> signals2=Bdata2$SIGNALS

Si hacemos un “str” de “signals2”, se verá que el número de señales es distinto al anterior. Esto ocurre porque el orden de la lista interviene en la creación de la variable “FDR”, que es la que determina el punto de corte.

Hay otro argumento que es importante y es “MC”. Por defecto su valor es “FALSE”. Si seteamos su valor a “MC=TRUE”, la distribución de probabilidad del IC es calculada por simulación por un método llamado Monte Carlo. Adicionalmente, el argumento “NB.MC” indica el número de simulaciones a realizar. Por defecto es 10000.

Si nuestra base es pequeña (y la mayoría los son) conviene usar este argumento para generar los parámetros calculados con menor error. Si bien el método es más lento, prioriza la precisión a la velocidad.

Finalmente, quiero decirles que no soy estadístico ni matemático ni informático. Estoy aprendiendo sobre esto casi al mismo tiempo que ustedes.

Les digo esto porque me gustaría que vuelquen sus opiniones o inquietudes en los comentarios al pie de cada post; todos serán bienvenidos. Contestaré todas las preguntas cuya respuesta sepa, y las que no, buscaremos la forma de aprenderlas.

Muchas gracias y disculpen la demora en postear. Nos vemos.

1.    Ahmed I, Thiessard F, Miremont-Salamé G, Bégaud B, Tubert-Bitter P. Pharmacovigilance Data Mining With Methods Based on False Discovery Rates: A Comparative Simulation Study. Clin Pharmacol Ther 2010 Sep;88(4):492-498.