Sobredispersión

Están aquí.

Dos aclaraciones:

• El tipo de letra que uso es Windsor en homenaje a Woody Allen.
• Las presentaciones las construyo con una versión tuneada y ad hoc de revealjs.

[Quod si sal evanuerit in quo sallietur ad nihilum valet ultra nisi ut mittatur foras et conculcetur ab hominibus.]

Vuelvo con mi monotema de los últimos días: cómo hacer GLMs de Poisson robustos. Encuentro la tesis Robust Inference for Generalized Linear Models: Binary and Poisson Regression y pienso: ajá, será cuestión de copipegar.

Nada más lejos de la realidad. El método propuesto en la tesis está basado en asignaciones de pesos a las observaciones usando kernels con centros y anchuras basadas respectivamente en

[Este es un extracto, una píldora atómica, de mi charla del otro día sobre el modelo de Poisson y al sobredispersión.]

Aunque me guste expresar el modelo lineal de la forma

$$ y_i \sim N(a_0 + \sum_j a_j x_{ij}, \sigma_i)$$

hoy, para lo que sigue, es más conveniente la representación tradicional

$$ y_i = a_0 + \sum_j a_j x_{ij} + \epsilon_i$$

donde si no sabes lo que es cada cosa, más vale que no sigas leyendo.

Le he dado muchas vueltas en estos últimos tiempos al asunto de la sobredispersión, particularmente en dos tipos de modelos: Poisson y logístico. Así que, aunque solo sea por proximidad semántica, se me quedan pegados ejemplos y casos de ese fenómeno mucho menos frecuente que es el de la infradispersión.

Un ejemplo ilustrativo del fenómeno que se me ocurrió el otro día era

pero hace nada, ese señor lleno de paz y amor que es Putin, nos ha regalado otro:

[Esta entrada continúa el ciclo al que he dedicado esta y esta otra entradas durante los últimos días.]

Las dos entradas anteriores de la serie se resumen en que:

• el modelo de Poisson no recoge todas las fuentes de error que pueden existir en los datos y que
• las soluciones al uso (como, p.e., usar modelos quasi-Poisson) son puros remiendos.

Si el error en el modelo de Poisson entra (también) en el término lineal, podemos modelar ese error explícitamente. Podría haber implementado la solución INLA o Stan del problema, pero me conformaré con la lme4. Primero, generaré los datos (igual que en las entradas anteriores) y añadiré una variable categórica que identifique cada registro:

[Esta entrada abunda sobre la de ayer y sin la cual no se entiende.]

Generemos unos datos, las x:

n <- 1000
sigma <- .5
x <- rep(-2:2, each = n)
x_real <- -1 + .5 * x + rnorm(length(x), 0, sigma)

En el bloque anterior hemos creado una/la variable observada, x, el término lineal que operará en el modelo de Poisson, -1 + .5 * x, y el real, -1 + .5 * x + rnorm(length(x), 0, sigma), que agrega al anterior el impacto de otras variables no tenidas en cuenta a través de un error normal al uso.

El modelo de Poisson viene a decir que si y es una variable con valores 0, 1,… y x1,…, xn son variables explicativas tiene cierto sentido en algunos casos plantear un modelo de la forma

$$ y | x_i \sim \text{Pois}(\exp(a_0 + \sum_i a_i x_i) ),$$

Es decir , para cada combinación de las xi, el modelo proporciona el parámetro de una distribución de Poisson de la que y es una realización. Hay una incertidumbre (o un error irreductible) que reside en que de y solo conocemos la distribución.

Esta entrada muestra cómo afrontar (con Stan) un problema que encontré el otro día en un lugar que no puedo mencionar pero en el que sé que me leen (y los destinatarios sabrán que va por ellos).

El contexto es el siguiente: se hace un test A/B donde la variable de interés son unos conteos. Hay varios grupos (aquí los reduciré a dos) y los datos siguen aproximadamente (aquí omitiré la parte de la inflación de ceros) una distribución de Poisson. Pero solo aproximadamente: existe sobredispersión, es decir, la varianza de los datos excede su media.

La distribución de Poisson se utiliza de oficio cuando se quiere modelar datos relativos a conteos. Sin embargo, tiene un problema serio: la varianza está fijada a la media: ambas son $latex \lambda$, el parámetro de la distribución.

Muy frecuentemente se observan datos con sobredispersión. Si $latex \lambda$ es 1000, el número esperado de eventos está contenido en un intervalo demasiado estrecho,

qpois(c(0.025, 0.975), 1000)
#[1]  938 1062

como para ser realista en muchas aplicaciones.