Category: Spanish

Acerca de… si es que el R sirve para todo.

En el blog sobre ecología y evolución EBB Flow. Nos han felicitado Halloween con una tarjeta hecha en R. Os pego el código: Ale, a destriparlo y a aprender cosas sobre gráficos en R. No os perdáis el segundo comentario. jaume Filed under: General, open source, R Tagged: Halloween, R, R cran

Acerca de… RStudio

Hola Dicen que rectificar es de sabios. Ya voy por el tercer intermediario (GUI) para R. Todo empezó con el JGR, necesitaba hacer el uso de R un poco más fácil en linux (vedlo aquí y aquí) La verdad es que no iba demasiado bien y era un lío trabajar de dos formas diferentes en […]

Web scraping rudimentario

La gracia de los métodos computacionales de análisis de datos reside en parte en la enorme cantidad de datos con la que disponemos en estos tiempos digitales, y la relativa facilidad para acceder a ellos.  Pero en general los datos no vienen ser…

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 lasalida del paquete PhViD para detección …

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.

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…

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.