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.