Separación perfecta en el modelo de Poisson

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.

Voy a ilustrarlo con el segundo ejemplo más sencillo que se me ocurre. Supongamos que

$$Y \sim \text{Pois}(X)$$

donde $latex X$ es una variable aleatoria que toma los valores $latex a$ y $latex b$. Supongamos que tenemos una muestra de tamaño $latex 2N$ donde a cada nivel de $latex X$ le corresponden $latex N$ casos. Los estimadores por máxima verosimilitud de los coeficientes correspondientes a esos valores son $latex \log n_a /N$ y $latex \log n_b /N$ respectivamente.

Pero, ¿qué pasa si $latex n_a = 0$? El estimador es $latex -\infty$; aunque, en realidad, acabo de ver que R se come la tostada:

set.seed(1)
N <- 100
x <- rep(c("a", "b"), each = N)
y <- c(rep(0, N), rpois(N, 1))
modelo <- glm(y ~ -1 + x, family = poisson)
summary(modelo)
# Call:
#   glm(formula = y ~ -1 + x, family = poisson)
#
# Deviance Residuals:
#   Min        1Q    Median        3Q       Max
# -1.42127  -0.00997  -0.00006  -0.00006   2.24293
#
# Coefficients:
#   Estimate Std. Error z value Pr(>|z|)
# xa  -20.30259 1554.18637  -0.013     0.99
# xb    0.00995    0.09950   0.100     0.92
#
# (Dispersion parameter for poisson family taken to be 1)
#
# Null deviance: 292.635  on 200  degrees of freedom
# Residual deviance:  92.625  on 198  degrees of freedom
# AIC: 252.98
#
# Number of Fisher Scoring iterations: 18

El coeficiente xb es, efectivamente

log(sum(y) / N)
# [1] 0.009950331

aunque glm nos engaña y da por convergida una regresión que no lo está. De todos modos, el coeficiente xa tiene un valor de -20 y un error estándar de 1554, nada menos, lo que debería hacer saltar alarmas donde hubiere luces.

Para saber más, y para que quede constancia de de dónde he sacado todo lo anterior, Bias Reduction as a Remedy to the Consequences of Infinite Estimates in Poisson and Tobit Regression del, entre otros, genial y nunca suficientemente apreciado A. Zeileis.

Coda: Apenas acabo lo anterior, me doy cuenta de que ya había hablado del tema de pasada hace un par de años largos.