Más sobre la presunta sobredispersión en el modelo de Poisson

Estadística

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

Generamos las y:

y <- sapply(x_real, function(lambda) rpois(1, exp(lambda)))

Realidad, y según la teoría bajo la que cabe aplicar el modelo de Poisson, los vectores

tapply(y, x, var)
#        -2        -1         0         1         2
# 0.1750711 0.2701892 0.4400561 0.8555716 1.5842683

y

exp(-1 + .5 * (-2:2))
# [1] 0.1353353 0.2231302 0.3678794 0.6065307 1.0000000

deberían ser aproximadamente iguales. De hecho, casi lo son si forzamos $latex \sigma = 0$:

y <- sapply(-1 + .5 * x, function(lambda) rpois(1, exp(lambda)))
tapply(y, x, var)
exp(-1 + .5 * (-2:2))

Si hacemos

datos <- data.frame(
    x = x,
    y = y
)

modelo_glm <- glm(y ~ x, data = datos, family = poisson)
summary(modelo_glm)

Obtenemos

Call:
glm(formula = y ~ x, family = poisson, data = datos)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.4731  -0.9109  -0.5633   0.3742   4.4913

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.87974    0.02415  -36.44   <2e-16 ***
x            0.48065    0.01598   30.07   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 6184.4  on 4999  degrees of freedom
Residual deviance: 5145.1  on 4998  degrees of freedom
AIC: 9051.2

Number of Fisher Scoring iterations: 6

Mientras que si nos decantamos por el modelo quasi-Poisson,

modelo_quasi <- glm(y ~ x, data = datos, family = quasipoisson)
summary(modelo_quasi)

obtenemos algo parecido,

Call:
glm(formula = y ~ x, family = quasipoisson, data = datos)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.4731  -0.9109  -0.5633   0.3742   4.4913

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.87974    0.02646  -33.24   <2e-16 ***
x            0.48065    0.01752   27.44   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 1.201217)

    Null deviance: 6184.4  on 4999  degrees of freedom
Residual deviance: 5145.1  on 4998  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

con la salvedad de que el modelo quasi nos advirte que existe una sobredispersión del 20%, i.e., que la varianza estimada de los datos es un 20% superior a la esperada por el modelo de Poisson.

De hecho,

tapply(y, x, var) / predict(modelo_glm, data.frame(x = -2:2), type = "response")
#        -2        -1         0         1         2
# 1.1426982 0.9928202 1.1523015 1.3737663 1.3469143

que es un vector cuya media coincide prácticamente con el parámetro de dispersión, 1.2017. Pero es que con este planteamiento en particular, la sobredispersión ni siquiera es proporcional al valor de $latex \hat{\mu}$, como demuestra la siguiente simulación,

x <- runif(100000, -2, 5)

y <- exp(x + rnorm(length(x), 0, .2))
y <- sapply(y, function(lambda) rpois(1, lambda))

#plot(x, y)

rejilla_x <- seq(min(x), max(x), length.out = 100)

get_local_variance <- function(x, y, centro, d = .1){
    tmp <- y[abs(x - centro) < d]
    var(tmp) / exp(centro)
}

local_variance <- sapply(
    rejilla_x,
    function(centro) get_local_variance(x, y, centro))

plot(rejilla_x, local_variance, type = "l",
        xlab = "", ylab = "dispersión local",
        main = "dispersión local en función del\nvalor de la expresión lineal")

que produce

[Nota aclaratoria: la gráfica anterior debería ser aproximadamente constante para que pudiésemos creernos el ajuste que la quasi-Poisson hace al modelo de Poisson en el caso en cuestión.]

Y el próximo día trato de ajustar el modelo sin las restricciones que impone la tecnología de los setenta.