Estadística Bayesiana

Imagina que tienes que generar (reitero: generar) datos compatibles con el siguiente modelo:

• Tienes n sujetos a los que se proporciona un remedio para dormir en distintas dosis (conocidas) en distintos días.
• El número adicional de horas que duerme cada sujeto es lineal con una pendiente que depende de la dosis (una serie de dosis fijas).
• Esa recta tiene un término independiente (el número de horas que duerme el sujeto con una dosis igual a cero del remedio).

Argumento que para generar los términos independientes usarías algo así como una normal de media igual a 8 horas. Seguro que usarías alguna otra distribución razonable para las pendientes (p.e., que prohibiese que con dosis pequeñas se durmiese, p.e., 80 horas).

A primeros de julio (de 2018) impartiré un curso de 15 horas de estadística bayesiana aplicada con Stan en la UPC (Barcelona). La información relevante está aquí y aquí.

El proyecto y su definición es un tanto contradictorio en sus propios términos, lo reconozco. Es muy difícil hacer algo aplicado y, a la vez, bayesiano. Y más, con Stan. Además, podrían acusarme de hipócrita: ¿cuándo fue la última vez que facturé (recuérdese: facturable es el grado máximo de aplicado) por algo hecho con Stan? Porque la idea, en el fondo, es otra: esencialmente, cómo replantear modelos y estrategias de modelización, aunque se implenten con herramientas métodos de índole frecuentista, para enriquecerlos con la visión bayesiana.

ABC significa, entre otras cosas, approximate bayesian computation. Por lo que parece, consiste en calcular $latex P(\theta ,|, \text{datos})$ por el tradicional y directo método del rechazo. Es decir:

• Planteas un modelo generativo, con sus prioris y todo.
• Simulas casos, casos y casos.
• Te quedas con los que cumplen un criterio de aceptación.

La distribución empírica de los parámetros en el subconjunto de los casos aceptados representa, en los libros está escrito, la distribución a posteriori. Sin MCMC ni historias.

Se ve que hay arqueólogos bayesianos. Un problema con el que se encuentran es que tropiezan con cacharros antiguos y quieren estimar su antigüedad.

Así que prueban distintos métodos (¿químicos?), cada uno de los cuales con su precisión, y acaban recopilando una serie de estimaciones y errores. Obviamente, tienen que combinarlas de alguna manera.

El modelo más simple es

$$ M_i \sim N(\mu, \sigma_i)$$

donde $latex \mu$ es la antigüedad (desconocida) del artefacto y los $latex \sigma_i$ son las varianzas distintas de los distintos métodos de medida, que arrojan las estimaciones $latex M_i$.

Allá por el 2007 publicó The Independent una portada en que se retractaba. El diario había sido un histórico defensor de la legalización de la marihuana. Ese día hizo público su cambio de postura. Al parecer, motivada por las evidencias sobre los efectos sobre la salud mental.

Este fin de semana he asistido a una serie de conferencias. En una de ellas participaba el representante de una organización que:

• Adoptaba de partida una posición militante, de parte, en cierto asunto de interés público.
• Se definía como data driven, evidence driven, etc.

La pregunta obvia y que no tuve ocasión de plantear (por eso la traigo aquí) es la siguiente: si los datos y la evidencia se obstinaran en subrayar la bondad de la posición contraria a la que actualmente mantienen, ¿cuál de sus dos principios abandonarían primero?

El contexto

Cada día $latex i$ ocurren eventos de cierta naturaleza (transacciones, fallecimientos, infartos, etc.) que interesa contar.

El problema

El número de eventos $latex n_i$ que ocurren el día $latex i$ no se conoce el día $latex i$ sino que va siendo conocido progresivamente los días $latex i+1, \dots$. Pero hace falta una estimación de $latex n_i$ antes del fin del mundo.

Los datos

• La distribución de los $latex n_i$ (basados en el histórico).
• La proporción (probabilidad) $latex p_\Delta$ de eventos del día $latex i$ que se conocen el día $latex i+\Delta$.

La solución prebayesiana

• La estadística tiene que estar totalmente integrada en la investigación: el diseño experimental es fundamental.
• Que no te asuste utilizar métodos modernos
• Preserva toda la información disponible en los datos: evita categorizar los predictores continuos y los valores predichos
• No asumas que algo opera linealmente
• Ten en cuenta la incerditumbre sobre el (no del) modelo y trata de minimizarlo usando conocimiento previo sobre la materia
• Usa remuestreos
• Considera el tamaño muestral una variable aleatoria cuando sea posible
• Usa estadística bayesiana siempre que sea posible
• Usa buenos gráficos frecuentemente
• Para que sea creíble, la investigación tiene que ser reproducible
• Toda la manipulación de datos y el análisis estadístico tiene que ser reproducible

Así se titula un libro que no he leído y que, pese a lo cual, como los malos críticos, voy a comentar. Los libros suelen estar plagados de hojarasca, tal vez porque de otra manera no se puede hacer crecer un par de ideas más o menos originales a las cientoypico páginas como mínimo que uno espera encontrar entre dos tapas. El relato corto no da caché. Y yo ando corto de tiempo.

Esto que cuento hoy puede ser muy útil: cómo mejorar los GLMs mediante la introducción de prioris (casi) a voluntad sobre los coeficientes. Usando el paquete arm de R, claro.

De momento y porque aún tengo sucios los datos sobre los que me gustaría aplicar el modelo, extraeré un ejemplo de la ayuda de la función principal del paquete, bayesglm.

Primero, preparo unos datos:

n <- 100
x1 <- rnorm (n)
x2 <- rbinom (n, 1, .5)
b0 <- 1
b1 <- 1.5
b2 <- 2
y <- rbinom (n, 1, invlogit(b0+b1*x1+b2*x2))

Comenzamos con un glm de toda la vida.

Tengo la sensación de que un lenguaje funcional (como Scala) está particularmente bien adaptado al tipo de operaciones que exige MCMC.

Juzguen Vds.

Primero, genero datos en R:

datos <- rnorm(500, 0.7, 1)
writeLines(as.character(datos), "/tmp/datos.txt")

Son de una normal con media 0.7. En el modelo que vamos a crear, suponemos conocida (e igual a 1) la varianza de la normal y trataremos de estimar la media suponiéndole una distribución a priori normal estándar. Y con Scala, así:

Porque ahí puedes tomarte una foto tal que

o

y luego tuitear cosas como

After #StrataHadoop: @jayusor and me in front of Bayes grave with my new @OReillyMedia book (emulating @gilbellosta) pic.twitter.com/SlguJIeLw0

— Antonio Sánchez Chinchón (@aschinchon) June 4, 2016

Para mayor referencia (y por tenerlo a mano cuando vuelva),

Las diapositivas de mi charla Datos, modelos y parámetros en el grupo Machine Learning Spain pueden verse/bajarse de aquí.

Porque voy a dar una charla en él. Es este jueves, por la tarde, en el Campus de Google de Madrid (los detalles).

Se tratará de una introducción a y justificación de aproximaciones más bayesianas de lo habitual a problemas reales del análisis de datos. Que comenzará con una explicación sobre cuándo 100% no significa 100% para terminar con lo que viene siéndome habitual últimamente: un ejemplo en rstan con su discusión.