Poisson

Estamos acostumbrados a la caracterización habitual de la distribución binomial negativa como el aburrido número de fracasos en una serie de ensayos de Bernoulli hasta lograr $r$ éxitos. Esto, junto con un poco de matemáticas de primero de BUP —todo aquello de combinaciones, etc.— lleva a la expresión conocida de su función de probabilidad,

$$\binom{n + x - 1}{x} p^r (1 - p)^x.$$

Pero esta caracterización, muy útil para resolver problemas de probabilidad construidos artificialmente para demostrar que los alumnos han estudiado la lección con aprovechamiento, se queda muy corta a la hora de proporcionar intuiciones sobre cómo, cuándo y por qué utilizarla en el ámbito en el que es más útil: el análisis de los procesos puntuales.

El asunto de la separación perfecta en el modelo logístico es sobradamente conocido. Solo quiero añadir al respecto dos cosas que no se suelen decir:

• Es un dolor que solo duele a los frecuentistas que no usan regularización (y van quedando cada vez menos de esos).
• Que no es malo sino bueno: ¿qué cosa mejor que tus datos puedan responder categóricamente las preguntas que les planteas (supuesto, claro, está, un N suficientemente grande).

Lo que es menos conocido es que el problema de la separación perfecta también puede afectar a la regresión de Poisson.

Una de las consencuencias del coronavirus es que vamos a tener que replantearnos lo que significa ajustar series temporales. Es decir, comenzar a ajustar series temporales y no repetir la consabida teoría que subyace a los modelos ARIMA simplemente porque es guay.

También tendremos que replantearnos qué hacer con los outliers que la pandemia va dejando tras de sí. Y tratar de hacerlo más elegantemente que cierta gente, por supuesto. En particular, habrá que ver cuál y cómo es el efecto de los outliers en determinados modelos. En particular, en esos en los que yo más trabajo últimamente, que son los de Poisson.

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.

Este es un anuncio de una charla que daré este viernes (2020-09-18) dentro del congreso virtual EncuentRos en la fase R. Ni que decir tiene que los detalles logísticos pueden consultarse en el enlace anterior.

Hablaré de cuestiones relativas al modelo de Possion (gran parte de las cuales pueden trasladarse también al logístico) de las que se habla poco y sobre las que la teoría que uno tropieza por ahí no es del todo clara pero que se manifiestan claramente en datos como los de la monitorización de la mortalidad, que será discutida también de pasada.

[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.

Fueron mis modelos favoritos un tiempo, cuando modelaba visitas y revisitas de usuarios a cierto malhadado portal.

Si las visitas fuesen aleatorias (en cierto sentido), tendrían un aspecto no muy distinto del que se obtiene haciendo

library(IHSEP)

suppressWarnings(set.seed(exp(pi * complex(imaginary = 1))))

tms <- simPois(int = function(x) .1, cens = 1000)
hist(tms, breaks = 100, main = "Proceso homogéneo de Poisson",
      xlab = "", ylab = "frecuencia")

Es decir,

o bien una distribución uniforme en el tiempo. Pero bien puede ocurrir que una visita incremente la probabilidad de otra inmediatamente después, por lo que las visitas tenderían a arracimarse en determinados momentos. Con el paquete [IHSEP](https://cran.r-project.org/package=IHSEP) de R pueden simularse (y ajustarse) este tipo de modelos. Por ejemplo,

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.

Escribo esta entrada con cierta prevención porque soy consciente de que dan pábulo a determinadas teorías conspiranoicas de las que soy declarado enemigo. Pero es que los números de muertos en carretera por accidente en España en los últimos años,

(extraídos de aquí) dan que pensar: la varianza de las observaciones correspondientes a los años 2013, 2014 y 2015 es muy baja, demasiado baja. Al menos, si se da como bueno un modelo de Poisson para modelar esos conteos.