Poisson - regresion
Alfredo Aro Terleira
2024-09-04
Subimos la data
library(rio)
data = import("data_final.xlsx")Revisamos la información de la data
str(data)## 'data.frame': 1870 obs. of 11 variables:
## $ Ubigeo : chr "010101" "010102" "010103" "010104" ...
## $ Departamento : chr "AMAZONAS" "AMAZONAS" "AMAZONAS" "AMAZONAS" ...
## $ Provincia : chr "CHACHAPOYAS" "CHACHAPOYAS" "CHACHAPOYAS" "CHACHAPOYAS" ...
## $ Distrito : chr "CHACHAPOYAS" "ASUNCION" "BALSAS" "CHETO" ...
## $ PEA_desocupada : num 601 12 10 3 0 19 3 64 33 40 ...
## $ poblacion : num 32589 262 1136 642 585 ...
## $ educacion : num 11 7.5 6.5 6.8 5.6 7.4 6 7.9 6.5 7.6 ...
## $ bibliotecas : chr "1" "0" "0" "0" ...
## $ instituciones : num 57 4 19 5 9 24 7 4 40 28 ...
## $ hogares_dependencia: num 1.9 0 8.5 1.4 0 6.5 3.4 2.4 4.7 7.3 ...
## $ plan_familiar : chr "1" "1" "0" "0" ...Explciación Poisson
Cuando nos referimos a crear una Regresión Poisson asumimos que la variable dependiente de nuestra regresión resulta ser de carácter numérica discreta, es decir, una variable de conteos. Ya sean por unidad de tiempo o espacio, por ejemplo: peso, altura, número de fallecidos, presupuestos o cantidad de votos a favor. (Esto también incluye a las tasas)
Asimismo, otros requisitos de la regresión Poisson son:
Independencia: Nuetras filas no deben tener relación alguna entre ellas
Equidispersión: Para que la Regresión Poisson funcione se necesita que la media y la varianza resulten ser iguales. Si este supuesto no se cumplen podemos encontrarnos frente a una Sobredispersión o Subdispersión
Linealidad: El logaritmo de la media de los datos, log(λ), debe ser una función lineal de los datos
Empecemos con el ejemplo
1. Regresión
Nuestra hipótesis es la siguiente: Nosotros creemos que la cantidad de la población económica desempleada (PEA_desocupada) se encuentra influenciada por otros factores como: el promedio de años de estudios de la población entre 15 años y más (educacion), si la municipalidad respectiva posee una biblioteca (bibliotecas), la cantidad de instituciones educativas públicas (instituciones), el porcentaje de hogares con alta dependencia económica (hogares_dependencia) y por si la municipalidad implementó o no un programa de control de planificación familiar (plan_familiar)
hp = formula(PEA_desocupada ~ educacion + bibliotecas + instituciones + hogares_dependencia + plan_familiar)
reg_p=glm(hp, data = data, #censo es nuestra base de datos
offset=log(poblacion), #no es obligatorio el offset
family = poisson(link = "log"))
summary(reg_p)##
## Call:
## glm(formula = hp, family = poisson(link = "log"), data = data,
## offset = log(poblacion))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.790e+00 1.143e-02 -331.577 < 2e-16 ***
## educacion 1.600e-02 1.028e-03 15.567 < 2e-16 ***
## bibliotecas1 4.599e-02 2.965e-03 15.515 < 2e-16 ***
## instituciones 3.365e-05 8.107e-06 4.151 3.31e-05 ***
## hogares_dependencia -1.087e-02 4.776e-04 -22.759 < 2e-16 ***
## plan_familiar1 -5.194e-03 2.491e-03 -2.085 0.0371 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 106098 on 1869 degrees of freedom
## Residual deviance: 100372 on 1864 degrees of freedom
## AIC: 112156
##
## Number of Fisher Scoring iterations: 5El offset no resulta obligatorio, pero es importante ya que representa nuestra variable de control. En el presente caso sería la población total de cada distrito. Colocamos el offset ya que, por ejemplo, no es lo mismo sostener que el impacto del covid fue igual de fuerte en dos ciudades solo porque ambas tuvieron alrededor de 2000 fallecidos. El impacto tiene que estar en relación con la cantidad total de población, por lo que si comprendemos que la primera ciudad contaba con alrededor de 4000 habitantes, y la segunda ciudad con 23000 no podríamos señalar que el impacto fue igual.
Si queremos una mejor visualización
library(kableExtra)
library(modelsummary)## `modelsummary` 2.0.0 now uses `tinytable` as its default table-drawing
## backend. Learn more at: https://vincentarelbundock.github.io/tinytable/
##
## Revert to `kableExtra` for one session:
##
## options(modelsummary_factory_default = 'kableExtra')
## options(modelsummary_factory_latex = 'kableExtra')
## options(modelsummary_factory_html = 'kableExtra')
##
## Silence this message forever:
##
## config_modelsummary(startup_message = FALSE)modelP=list("POISSON asegurados"=reg_p)
modelsummary(modelP, title = "Regresión Poisson",
stars = TRUE,
output = "kableExtra")| POISSON asegurados | |
|---|---|
| (Intercept) | -3.790*** |
| (0.011) | |
| educacion | 0.016*** |
| (0.001) | |
| bibliotecas1 | 0.046*** |
| (0.003) | |
| instituciones | 0.000*** |
| (0.000) | |
| hogares_dependencia | -0.011*** |
| (0.000) | |
| plan_familiar1 | -0.005* |
| (0.002) | |
| Num.Obs. | 1870 |
| AIC | 112155.6 |
| BIC | 112188.8 |
| Log.Lik. | -56071.796 |
| F | 1096.432 |
| RMSE | 205.14 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
¡PERO!: a diferencia de la anterior regresión, en este caso debemos exponenciar los coeficientes para calcular el efecto real sobre nuestra variable dependiente
modelsummary(modelP, title = "Regresión Poisson Exponenciada",
stars = TRUE,
statistic = 'conf.int',
exponentiate = T,
output = "kableExtra")| POISSON asegurados | |
|---|---|
| (Intercept) | 0.023*** |
| [0.022, 0.023] | |
| educacion | 1.016*** |
| [1.014, 1.018] | |
| bibliotecas1 | 1.047*** |
| [1.041, 1.053] | |
| instituciones | 1.000*** |
| [1.000, 1.000] | |
| hogares_dependencia | 0.989*** |
| [0.988, 0.990] | |
| plan_familiar1 | 0.995* |
| [0.990, 1.000] | |
| Num.Obs. | 1870 |
| AIC | 112155.6 |
| BIC | 112188.8 |
| Log.Lik. | -56071.796 |
| F | 1096.432 |
| RMSE | 205.14 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
¿Qué sigue?: Prueba de Equidispersión
Se supone que deberíamos proceder con el análisis de impacto de cada variable; sin embargo, si recuerdas lo que dijimos al inicio, debemos confirmar si nuestra media y nuestra varianza son iguales. Para ello empleamos la prueba de Equidispersión
library(magrittr)
library(AER)## Cargando paquete requerido: car## Cargando paquete requerido: carData## Cargando paquete requerido: lmtest## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Cargando paquete requerido: sandwich## Cargando paquete requerido: survivaloverdispersion=AER::dispersiontest(reg_p,alternative='greater')$ p.value<0.05
underdispersion=AER::dispersiontest(reg_p,alternative='less')$ p.value<0.05
# tabla
testResult=as.data.frame(rbind(overdispersion,underdispersion))
names(testResult)='¿Es probable?'
testResult%>%kable(caption = "Test de Equidispersión")%>%kableExtra::kable_styling()| ¿Es probable? | |
|---|---|
| overdispersion | TRUE |
| underdispersion | FALSE |
En este caso notamos que nuestra regresión posee sobredispersión, por lo que tendremos que emplear la regresión Quasipoisson y la Binomial Negativa
RECUERDA: La Quasipoisson se usa cuando existe Sobredispersión (overdispersion) o Subdispersión (underdispersion). Mientras que la Binomial Negativa solo se puede emplear cuando existe Sobredispersión.
Quasipoisson
reg_qp = glm(hp, data = data,
offset=log(poblacion),
family = quasipoisson(link = "log"))
summary(reg_qp)##
## Call:
## glm(formula = hp, family = quasipoisson(link = "log"), data = data,
## offset = log(poblacion))
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.790e+00 9.602e-02 -39.465 < 2e-16 ***
## educacion 1.600e-02 8.635e-03 1.853 0.06407 .
## bibliotecas1 4.599e-02 2.491e-02 1.847 0.06497 .
## instituciones 3.365e-05 6.811e-05 0.494 0.62132
## hogares_dependencia -1.087e-02 4.013e-03 -2.709 0.00682 **
## plan_familiar1 -5.194e-03 2.093e-02 -0.248 0.80405
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 70.59156)
##
## Null deviance: 106098 on 1869 degrees of freedom
## Residual deviance: 100372 on 1864 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5Binomial Negativa
library(MASS)
hp2=formula(PEA_desocupada ~ educacion + bibliotecas + instituciones + hogares_dependencia + plan_familiar + offset(log(poblacion))) #nueva fórmula necesaria!!!
reg_bn=glm.nb(hp2,data=data)
summary(reg_bn)##
## Call:
## glm.nb(formula = hp2, data = data, init.theta = 2.425123544,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.0903003 0.1102158 -37.112 < 2e-16 ***
## educacion 0.0549957 0.0122235 4.499 6.82e-06 ***
## bibliotecas1 0.0935495 0.0372028 2.515 0.0119 *
## instituciones -0.0009381 0.0002896 -3.239 0.0012 **
## hogares_dependencia 0.0026304 0.0037317 0.705 0.4809
## plan_familiar1 0.0046177 0.0364908 0.127 0.8993
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(2.4251) family taken to be 1)
##
## Null deviance: 2082.8 on 1869 degrees of freedom
## Residual deviance: 2034.2 on 1864 degrees of freedom
## AIC: 20775
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 2.4251
## Std. Err.: 0.0805
##
## 2 x log-likelihood: -20761.1780Si queremos ver nuestras 3 regresiones exponenciadas:
formatoNum <- function(x) format(x, digits = 4, scientific = FALSE)
models_total=list('Poisson asegurados'=reg_p,
'QuasiPoisson asegurados'=reg_qp,
'Binomial Negativa asegurados'=reg_bn)
modelsummary(models_total,fmt=formatoNum,
exponentiate = T,
statistic = 'conf.int',
title = "EXP() de la Regresiones Poisson, Quasi Poisson y Binomial Negativa",
stars = TRUE,
output = "kableExtra")| Poisson asegurados | QuasiPoisson asegurados | Binomial Negativa asegurados | |
|---|---|---|---|
| (Intercept) | 0.02261*** | 0.02261*** | 0.01673*** |
| [0.0221, 0.02312] | [0.01873, 0.02729] | [0.01345, 0.02082] | |
| educacion | 1.01613*** | 1.01613+ | 1.05654*** |
| [1.0141, 1.01818] | [0.99906, 1.03346] | [1.03106, 1.08275] | |
| bibliotecas1 | 1.04707*** | 1.04707+ | 1.09806* |
| [1.0410, 1.05317] | [0.99729, 1.09958] | [1.02107, 1.18149] | |
| instituciones | 1.00003*** | 1.00003 | 0.99906** |
| [1.0000, 1.00005] | [0.99990, 1.00017] | [0.99853, 0.99962] | |
| hogares_dependencia | 0.98919*** | 0.98919** | 1.00263 |
| [0.9883, 0.99011] | [0.98138, 0.99694] | [0.99563, 1.00977] | |
| plan_familiar1 | 0.99482* | 0.99482 | 1.00463 |
| [0.9900, 0.99969] | [0.95474, 1.03638] | [0.93509, 1.08021] | |
| Num.Obs. | 1870 | 1870 | 1870 |
| AIC | 112155.6 | 20775.2 | |
| BIC | 112188.8 | 20813.9 | |
| Log.Lik. | -56071.796 | -10380.589 | |
| F | 1096.432 | 15.532 | 10.696 |
| RMSE | 205.14 | 205.14 | 374.19 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
¿Qué nos toca hacer ahora?: Debemos confirmar si es que se eliminó la Sobredispersión para saber que regresión escoger
Debemos recordar que si bien los AIC/BIC son importantes, en este caso buscaremos a la regresión que pueda eliminar la sobredispersión o la subdispersión si es el caso. Asimismo, recordar que la tabla anova si se emplea para comparar debe ser con cuidado, ya que no se tratan de modelos anidados
#poisson case
performance::check_overdispersion(reg_p)## # Overdispersion test
##
## dispersion ratio = 70.592
## Pearson's Chi-Squared = 131582.675
## p-value = < 0.001## Overdispersion detected.#quasipoisson case
performance::check_overdispersion(reg_qp)## # Overdispersion test
##
## dispersion ratio = 70.592
## Pearson's Chi-Squared = 131582.675
## p-value = < 0.001## Overdispersion detected.#negative binomial case
performance::check_overdispersion(reg_bn)## # Overdispersion test
##
## dispersion ratio = 0.186
## p-value = < 0.001## Underdispersion detected.Lo ideal sería que nos salga un mensaje así en una de las tres pruebas (ya que así sabremos cuál regresión escoger): ## No overdispersion detected.
Sin embargo, en ocasiones la vida se complica y debemos continuar :´(
En mi caso optaré por utilizar la regresión QuasiPoisson:
Interpretación de la Regresión QuasiPoisson
modelQP=list("QUASIPOISSON asegurados"=reg_qp)
modelsummary(modelQP, title = "Regresión Poisson Exponenciada",
stars = TRUE,
statistic = 'conf.int',
exponentiate = T,
output = "kableExtra")| QUASIPOISSON asegurados | |
|---|---|
| (Intercept) | 0.023*** |
| [0.019, 0.027] | |
| educacion | 1.016+ |
| [0.999, 1.033] | |
| bibliotecas1 | 1.047+ |
| [0.997, 1.100] | |
| instituciones | 1.000 |
| [1.000, 1.000] | |
| hogares_dependencia | 0.989** |
| [0.981, 0.997] | |
| plan_familiar1 | 0.995 |
| [0.955, 1.036] | |
| Num.Obs. | 1870 |
| F | 15.532 |
| RMSE | 205.14 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Como observamos, nuestra tabla nos muestra algo nuevo: los intervalos de confianza Estos resultan importantes ya que nos van a confimar cuáles son realmente las variables que resultan significativas. En nuestro caso debemos tener claro que para que nuestra variable independiente resulta significativa, su intervalo de confianza no debe contener al 1.
Si no lo ves muy bien (como yo jeje) podemos usar una tabla
library(ggplot2)
library(dotwhisker)
dotwhisker::dwplot(list(Poisson=reg_p,CuasiPoisso=reg_qp,BinomialNegativa=reg_bn),exp=T) + scale_y_discrete(labels=c("plan_familiar","hogares_dependencia","instituciones","bibliotecas","educacion")) + scale_color_discrete(name="Modelos para:\nCantidad de Asegurados") + geom_vline(
xintercept = 1,
colour = "grey60",
linetype = 2
)
Esta gráfica nos demuestrra que nuestra única variable significante
resulta ser: hogares_dependencia
Interpretación
Debemos tener claro que nuestra interpretación siempre se encontrará en función del número 1, por lo que debemos comprender el resultado de la siguiente forma:
Hogares dependencia:
exp(coef(reg_qp)[['hogares_dependencia']])## [1] 0.9891885((1-0.9891885)*-1)*100## [1] -1.08115–> Encontramos que por cada unidad que aumenta nuestra variable HOGARES_DEPENDENCIA, nuestra variable dependiente disminuye en 1.08%.
Todo esto no es necesario si es que notamos que la interpretación se encuentra en función del 1
Pero si sale positivo? QUÉ HAGO??
Aquí por fines educativos vamos a intepretar otra variable
exp(coef(reg_qp)[['educacion']])## [1] 1.016126((1-1.016126)*-1)*100## [1] 1.6126Señalaremos, pues, que por cada unidad en la que aumenta nuestra variable EDUCACION, la variable dependiente aumentará en 1.61%