Regresión de Poisson
Maria Durban
Introducción
Si la variable respues ta es de conteo:
Acotadad⇒ Binomial
Grandes⇒ Normal
Otro caso⇒ Poisson
Introducción
Si la variable respues ta es de conteo:
Acotadad⇒ Binomial
Grandes⇒ Normal
Otro caso⇒ Poisson
Si el conteo es pequeño:
solo predictores categóricos (contingency table) ⇒ modelos log-lineales
predictores categóricos/continuos ⇒ regresión de Poisson
Introducción
La regresión de Poisson también se usa para modelizar tasas
Una tasa es: un conteo de sucesos que ocurren en una unidad de observación particular
dividido por alguna medida de esa unidad de observación exposición
el número de sucesos se puede ver afectado por la población en riesgo de sufrir ese
evento (demografía, epidemiología, seguros, etc.)
Se incluye la exposición (exp) en el denominador y se calcula ta tasa:
tasa = Y /exp ⇒ Y = tasa × exp
Introducción
Los modelos de regresión lineal (varianza constante y errores normales) no son apropiados para datos
de conteo por 4 razones fundamentales:
Los modelos de regresión lineal pueden dar lugar a predicciones de conteo negativas
En el caso de datos de conteo la varianza de la variable respuesta suele depender de la media
Los errores no se distribuyen como una Normal
Es difícil usar transformaciones cuando hay zeros en la variable respuesta
La distribución de Poisson
Dada una variable Y que se distribuye como una Poisson con parámetro μ > 0 , para
k = 0, 1, 2, … su función de probabilidd viene dada por:
−μ k
e μ
P r(Y = k) =
k!
Además, E[Y] = V ar[Y] = μ .
La distribución de Poisson
Estimación de parámetros e inferencia
En este caso, tampoco modelizamos directamente la media, μ,sino una transformación de la
misma: log() (así aseguramos valores positivos)
log(μ) = β 0 + β 1 X1 + … + β k Xk = Xβ
Estimación de parámetros e inferencia
En este caso, tampoco modelizamos directamente la media, μ,sino una transformación de la
misma: log() (así aseguramos valores positivos)
log(μ) = β 0 + β 1 X1 + … + β k Xk = Xβ
La verosimilitud:
n −μ y
e μ i
L(y1 , … , yn ) = ∏
yi !
i=1
Estimación de parámetros e inferencia
En este caso, tampoco modelizamos directamente la media, μ,sino una transformación de la
misma: log() (así aseguramos valores positivos)
log(μ) = β 0 + β 1 X1 + … + β k Xk = Xβ
La verosimilitud:
n −μ y
e μ i
L(y1 , … , yn ) = ∏
yi !
i=1
La log-verosimilitud:
n n
′ ′
log L(y1 , … , yn ) = ∑ −μ + yi log(μ) − log(yi !) = ∑ exp(X β) + yi log(X β) − log(yi !)
i i
i=1 i=1
Estimación de parámetros e inferencia
En este caso, tampoco modelizamos directamente la media, μ,sino una transformación de la
misma: log() (así aseguramos valores positivos)
log(μ) = β 0 + β 1 X1 + … + β k Xk = Xβ
La verosimilitud:
n −μ y
e μ i
L(y1 , … , yn ) = ∏
yi !
i=1
La log-verosimilitud:
n
∂log L
′
= ∑(yi − exp(X β))Xi
i
∂β
i=1
Resolvemos por Newton-Rapson
Estimación de parámetros e inferencia
es una función lineal de los predictores ⇒ μ es una función multiplicativa de los
log(μ)
mismos:
Estimación de parámetros e inferencia
es una función lineal de los predictores ⇒ μ es una función multiplicativa de los
log(μ)
mismos:
β β X1 β Xk Xβ
μ = eo × e 1
×… ×e k
= e
e
βj
representa el efecto multiplicativo del j-ésimo predictor sobre la media
un incremento en 1 unidad de x tiene un efecto multiplicativo e
j
βj
sobre μ
Estimación de parámetros e inferencia
Usamos el deviance para comparar modelos:
verosimilitud del modelo ajustado
D = −2 log[ ]
lverosimilitud del modelo perfecto
Para muestras grandes D ∼ χ 2
n−k+1
Estimación de parámetros e inferencia
Bondad de ajuste:
Usamos el estadístico chi-cuadrado de Pearson:
n 2
∑i=1 (yi − μ
^ )
i
2
χ =
^
μi
comparamos este estadístico con una χ 2
k+1
.
Este test es similar al test del Deviance que también podemos usar.
Estimación de parámetros e inferencia
Goodness of fit:
Usamos el estadístico chi-cuadrado de Pearson:
n 2
∑i=1 (yi − μ
^ )
i
2
χ =
^
μi
comparamos este estadístico con una χ 2
k+1
.
Este test es similar al test del Deviance que también podemos usar.
Gráficos de residuos:
Usamos los residuos de Pearson:
^
yi − μi
ri = ri ≈ N (0, 1)
−−
^
√μ
i
Caso de estudio: riqueza de especies
Experimeto agrícola que consta de 90 parcelas de 25m x
25m acda una.
Se mide la biomasa, el pH del suelo y la riqueza de
especies (número de especies diferentes en cada parcela).
La hipótesis de partida es que el número de especies
desciende con la biomasa.
¿Depende la pendiente de la relación del tipo de pH del
suelo?
Caso de estudio: riqueza de especies
pH Biomass Species
2 0.4692972 30
2 1.7308704 39
2 2.0897785 44
2 3.9257871 35
2 4.3667927 25
2 5.4819747 29
2 6.6846859 23
2 7.5116506 18
2 8.1322025 19
2 9.5721286 12
Caso de estudio: riqueza de especies
Primero: ajustamos un modelo de regresión lineal para estos datos:
El modelo daría predicciones negativas del número de especies si la biomasa es 7.5 en suelos
con pH bajo
Caso de estudio: riqueza de especies
Empezamos por ajustar cada predictor por separado y comprobamos si son significativos:
Biomasa
m0=glm(Species~Biomass,family=poisson,data=species)
m0
##
## Call: glm(formula = Species ~ Biomass, family = poisson, data = species)
##
## Coefficients:
## (Intercept) Biomass
## 3.18409 -0.06444
##
## Degrees of Freedom: 89 Total (i.e. Null); 88 Residual
## Null Deviance: 452.3
## Residual Deviance: 407.7 AIC: 830.9
anova(m0,test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: Species
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 89 452.35
## Biomass 1 44.673 88 407.67 2.328e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Caso de estudio: riqueza de especies
Empezamos por ajustar cada predictor por separado y comprobamos si son significativos:
pH
species$pH=factor(species$pH)
m1=glm(Species~pH,family=poisson,data=species)
m1
##
## Call: glm(formula = Species ~ pH, family = poisson, data = species)
##
## Coefficients:
## (Intercept) pH1 pH2
## 2.4481 0.5476 0.8403
##
## Degrees of Freedom: 89 Total (i.e. Null); 87 Residual
## Null Deviance: 452.3
## Residual Deviance: 265.1 AIC: 690.3
anova(m1,test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: Species
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 89 452.35
## pH 2 187.23 87 265.12 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Caso de estudio: riqueza de especies
Veamos si el efecto de la biomasa en el número de especies depende del pH del suelo y si la
fuerza de dicha relación cambia con el pH
m2 <- glm(Species~Biomass+pH,family=poisson,data=species)
m3 <- glm(Species~Biomass*pH,family=poisson,data=species)
anova(m2,m3,test="Chisq")
## Analysis of Deviance Table
##
## Model 1: Species ~ Biomass + pH
## Model 2: Species ~ Biomass * pH
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 86 99.242
## 2 84 83.201 2 16.04 0.0003288 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Caso de estudio: riqueza de especies
Caso de estudio: riqueza de especies
El modelo ajustado:
log(Especies) = 2.95 − 0.261 × biomass + 0.484pHmedium + 0.815 × pHhigh
+0.123 × pHmedium × biomass + 0.155 × pHmedium × biomass
Un aumento en una unidad en la biomasa resultaría en un número de especies exp(−0.262) = 0.769
veces menor en suelos de pH bajo
Un aumento en una unidad en la biomasa resultaría en un número de especies
exp(−0.262 + 0.123) = 0.87 veces menor en suelos de pH medio
Un aumento en una unidad en la biomasa resultaría en un número de especies
exp(−0.262 + 0.155) = 0.898 veces menor en suelos de pH alto
Caso de estudio: riqueza de especies
Usamos los gráficos de residuos para comprobar las hipótesi del modelo (recuerda usar los
residuos de Pearson)
Los gráficos son razonablemente buenos
Caso de estudio: riqueza de especies
Comprobamos la bomdad de ajuste del modelo:
H0 : Model fits well the data
H1 : Model does not fit well the data
#Pearson's chi-squared statistic
pchisq(m3$deviance,m3$[Link],[Link]=FALSE)
## [1] 0.5041192
Regresión de Poisson para tasas de incidencia
El número de eventos observados puede depender del tamaño de una variable que
determina el número de oportunidades para que ese evento ocurra.
Ejemplos:
Número de robos denunciados en dioferentes ciudades depende del número de viviendas
en cada ciudad
Número de inidentes violentos cometidos en urgencias de un hospital psiquiátrico en un
periodo de 6 meses depende de la oportunidad de cometerlos (el número de días dentro
de ese periodo en que un paciente está en el área cercana al hospital)
Regresión de Poisson para tasas de incidencia
Y = conteo (e.g., número de ataques con violencia)
t = índice de espacio, tiempo o población (por ejemplo, número de días)
La tasa muestral de incidencia es Y /t.
El valor esperado sería:
E[Y /t] = E[Y ]/t = μ/t.
Regresión de Poisson para tasas de incidencia
Y = conteo (e.g., número de ataques con violencia)
t = índice de espacio, tiempo o población (por ejemplo, número de días)
La tasa muestral de incidencia es Y /t.
El valor esperado sería:
E[Y /t] = E[Y ]/t = μ/t
No podemo modelizar la tasa directamente, lo hacemos a través de Y
Regresión de Poisson para tasas de incidencia
El modelo de regresión de Poisson para la tasa esperada de ocurrencia de un evento es:
log(μ/t) = β 0 + β 1 x
log(μ) − log(t) = β 0 + β 1 x ⇒ log(μ) = β 0 + β 1 x + log(t)
Regresión de Poisson para tasas de incidencia
El modelo de regresión de Poisson para la tasa esperada de ocurrencia de un evento es:
log(μ/t) = β 0 + β 1 x
log(μ) − log(t) = β 0 + β 1 x ⇒ log(μ) = β 0 + β 1 x + log(t)
log(t) se llama offset
Es conocido y no hay que estimarlo, tampoco es una covariable por lo que no le
acompaña ningún parámetro
Modifica el valor del conteo estimado:
^ ^
β0 β 1x
^ = t×e
μ e
Regresión de Poisson para tasas de incidencia
Interpretación de los parámetros:
Riesgo relativo (RR)
Sea λ 0 = μ0 /t la tasa de incidencia cuando x = x y λ cuando x = x
0 1 0 +1 ,
λ1 λ1
β
log ( ) = log(λ1 ) − log(λ0 ) = β ⇒ RR = = e
λ0 λ0
e
β
es el riesgo relativo entre individuos con x = x y x = x . Es por cuánto se multiplica
0 1
el riesgo de que ocurra el evento para individuoa con x = x frente a individuos con
1
x = x . 0
Caso de estudio: Diabetes en la Comunidad de Madrid
Los datos corresponden al registro de incidencias de
diabetes mellitus en niños de 0 a 14 años de la
Comunidad de Madrid entre 1997 y 2005
Caso de estudio: Diabetes en la Comunidad de Madrid
sexo: 1=Varón, 2=Mujer
edad: 1=0-4,2=5-9, 3=10-15
mes: 1- 12
periodo: 1997-2003
casos: número de casos por estrato
poblacion: número de personas en el estrato
estación: 1=fría (Octubre a Marzo), 2= caliente (Abril a Septiembre)
Caso de estudio: Diabetes en la Comunidad de Madrid
Leemos los datos y definimos las variables categóricas como factores
diabetes=[Link]("./datos/[Link]",header=TRUE)
names(diabetes)
## [1] "edad" "sexo" "mes" "periodo" "casos" "poblacion"
## [7] "estacion"
diabetes$periodo=factor(diabetes$periodo)
diabetes$estacion=factor(diabetes$estacion)
diabetes$sexo=factor(diabetes$sexo)
diabetes$mes=factor(diabetes$mes)
Por ejemplo, para niños de edad 0 a 4, en Diciembre 1997 , y = 2 y E = 110324 , Por lo tanto:
2
λ = = 0.0000181 ó 1.81 por cada 100000 niños
110324
Caso de estudio: Diabetes en la Comunidad de Madrid
# Calculamos el offset
logpobla=log(diabetes$poblacion)
diabetes1=glm(casos~periodo+offset(logpobla),family=poisson,data=diabetes)
summary(diabetes1)
##
## Call:
## glm(formula = casos ~ periodo + offset(logpobla), family = poisson,
## data = diabetes)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0439 -0.7709 -0.3837 0.7293 2.7097
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.21951 0.08944 -125.438 <2e-16 ***
## periodo1998 0.02386 0.12649 0.189 0.850
## periodo1999 0.04955 0.12599 0.393 0.694
## periodo2000 0.13387 0.12391 1.080 0.280
## periodo2001 -0.15710 0.13238 -1.187 0.235
## periodo2002 0.08106 0.12369 0.655 0.512
## periodo2003 0.03106 0.12527 0.248 0.804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 686.27 on 503 degrees of freedom
## Residual deviance: 680.44 on 497 degrees of freedom
## AIC: 1716.9
##
## Number of Fisher Scoring iterations: 5
El año no parece significativo, lo comprobamos:
# Calculamos el offset
anova(diabetes1,test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: casos
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 503 686.27
## periodo 6 5.8389 497 680.44 0.4415
Qué estaríamos calculando si hacemos: e (0.04955−0.02386)
Caso de estudio: Diabetes en la Comunidad de Madrid
Calculamos intervalos de confianza para el riesgo relativo:
exp(confint(diabetes1))
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 1.119437e-05 0.0000159
## periodo1998 7.990119e-01 1.3127114
## periodo1999 8.206769e-01 1.3456605
## periodo2000 8.968513e-01 1.4585479
## periodo2001 6.584309e-01 1.1071044
## periodo2002 8.511035e-01 1.3829685
## periodo2003 8.068865e-01 1.3192888
Para calcular las tasas estimadas:
[Link]=predict(diabetes1,type="response")
[Link]=[Link]/diabetes$poblacion
Caso de estudio: Diabetes en la Comunidad de Madrid
Ejercicio
Comprueba si hay cambios en la incidencia de diabetes por edad, mes (o estación) y sexo.
¿Qué ocurre si metemos estación junto con mes?
A qué esdad es mayor la incidencia?, es mayor en niños o niñas?
