Sesión 3
Profesor:Dr. José Manuel MAGALLANES REYES, Ph.D
Profesor Principal del Departamento de Ciencias Sociales, Sección de Ciencia Política y Gobierno.
Oficina 223 - Edificio CISEPA / ECONOMIA / CCSS
Telefono: (51) 1 - 6262000 anexo 4302
Correo Electrónico: jmagallanes@pucp.edu.pe
Regresionando conteos
1 Data como conteos
Las variables que representan conteos son prominentes en las ciencias sociales. Podemos encontrar muchas en los indicadores que se crean en los censos. Por ejemplo, el Instituto Nacional de Estadística del Perú (INEI), organizó muy bien diversos indicadores para su XI Censo de Población y VI Censo de Vivienda (INEI 2007) en su portal para el 2007. Si queriamos descargar los indicadores disponibles para ese censo a nivel distrital, simplemente se accede a su página y se selecciona todos los distritos (es lo que queremos en este momento):
Descarguemos los siguientes indicadores a nivel de distrito:
- Conteos:
- Número de personas que tienen algún tipo de seguro.
- Total de habitantes (2007).
- Procentajes:
- Porcentaje de Trabajadores Independientes o por cuenta propia no profesionales.
- Porcentaje de personas analfabetas de 15 años y más.
Para ello, en el portal debes seleccionar esos indicadores de esta manera:
Al pedir “procesar” aparece este resultado en otra página:
Esa página permite que reorganices los datos. Entonces, busca darle este formato:
Luego, descargalo con CSV (añade .csv al nombre de archivo), ábrelo en Excel y verás algo así:
Nota que todo ha aperecido en una sola columna. Como se ve, cada columna tiene una coma que la separa de otra; sepáralas de manera adecuaday guardalo como .xlsx (aprovecha de eliminar las primeras dos columnas) .
El archivo guardado lo podemos guardar en un repositorio de GitHub, y abrirlo desde ahí.
rm(list = ls())
library(rio)
censo2007=import("https://github.com/Estadistica-AnalisisPolitico/Sesion3/raw/main/data2007peru.xlsx")
head(censo2007)## Pais Departamento Provincia Distrito Tema
## 1 Perú Amazonas Bagua Aramango Demográfico
## 2 Perú Amazonas Bagua Aramango Demográfico
## 3 Perú Amazonas Bagua Aramango Demográfico
## 4 Perú Amazonas Bagua Aramango Demográfico
## 5 Perú Amazonas Bagua Aramango Demográfico
## 6 Perú Amazonas Bagua Aramango Económico-Laboral
## Sub Tema
## 1 General
## 2 General
## 3 General
## 4 General
## 5 General
## 6 Población Económicamente Activa
## Descripcion
## 1 Total de habitantes del censo 2007
## 2 Total de habitantes del censo 2007
## 3 Total de habitantes del censo 2007
## 4 Total de habitantes del censo 2007
## 5 Total de habitantes del censo 2007
## 6 Porcentaje de Trabajadores Independientes o por cuenta propia no profesionales
## Clase Indicadores Valor
## 1 Total NA 11442
## 2 Área Urbana NA 2657
## 3 Área Rural NA 8785
## 4 Sexo - Hombre NA 6141
## 5 Sexo - Mujer NA 5301
## 6 Total NA 54.61Veamos qué tipos de datos tenemos, y los nombres de columnas:
str(censo2007)## 'data.frame': 36620 obs. of 10 variables:
## $ Pais : chr "Perú" "Perú" "Perú" "Perú" ...
## $ Departamento: chr "Amazonas" "Amazonas" "Amazonas" "Amazonas" ...
## $ Provincia : chr "Bagua" "Bagua" "Bagua" "Bagua" ...
## $ Distrito : chr "Aramango" "Aramango" "Aramango" "Aramango" ...
## $ Tema : chr "Demográfico" "Demográfico" "Demográfico" "Demográfico" ...
## $ Sub Tema : chr "General" "General" "General" "General" ...
## $ Descripcion : chr "Total de habitantes del censo 2007" "Total de habitantes del censo 2007" "Total de habitantes del censo 2007" "Total de habitantes del censo 2007" ...
## $ Clase : chr "Total" "Área Urbana" "Área Rural" "Sexo - Hombre" ...
## $ Indicadores : logi NA NA NA NA NA NA ...
## $ Valor : chr "11442" "2657" "8785" "6141" ...Lo más notorio es que la columna Valor ha sido interpretada como texto y no cómo tipo numérico. R hará eso si es que alguna celda tiene un carácter que no sea número o punto. Veamos:
# qué celda tiene algo que no sea número
censo2007$Valor[!grepl("[0-9]", censo2007$Valor)]## [1] "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-"
## [20] "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-"
## [39] "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-"
## [58] "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-" "-"
## [77] "-" "-"Sabiendo ello, podemos convertir la columna a numérico sin temor a perder valores:
censo2007$Valor=as.numeric(censo2007$Valor)## Warning: NAs introduced by coercionAhora, podemos quedarnos las columnas que utilizaremos en el ejemplo :
censo2007_sub=censo2007[,-c(1,5,6,9)] # -c() elimina esas columnas
# queda
str(censo2007_sub)## 'data.frame': 36620 obs. of 6 variables:
## $ Departamento: chr "Amazonas" "Amazonas" "Amazonas" "Amazonas" ...
## $ Provincia : chr "Bagua" "Bagua" "Bagua" "Bagua" ...
## $ Distrito : chr "Aramango" "Aramango" "Aramango" "Aramango" ...
## $ Descripcion : chr "Total de habitantes del censo 2007" "Total de habitantes del censo 2007" "Total de habitantes del censo 2007" "Total de habitantes del censo 2007" ...
## $ Clase : chr "Total" "Área Urbana" "Área Rural" "Sexo - Hombre" ...
## $ Valor : num 11442 2657 8785 6141 5301 ...Este último data frame está en formato long:
library(magrittr)
head(censo2007_sub,20)%>%
rmarkdown::paged_table(options = )Reconocemos el formato long cuando, sabiendo que la unidad de análisis es el distrito, vemos que el distrito se repite en muchas filas. Arriba vemos que todas las variables de Aramango aparecen al lado de este distrito. Nota además que los nombres de las variables están en las filas (no como titulos de las columnas) y que los valores de la variable aparecen en la última fila.
Pasémoslo a formato wide:
library(reshape2)
censo2007_wide=dcast(data=censo2007_sub,
formula=Departamento+Provincia+Distrito ~ Descripcion+Clase,
value.var="Valor")Los nombres de columnas son muy largos luego de esta conversión serán muy largos:
# columnas del formato wide
names(censo2007_wide)## [1] "Departamento"
## [2] "Provincia"
## [3] "Distrito"
## [4] "Personas que tienen algún tipo de seguro_Área Rural"
## [5] "Personas que tienen algún tipo de seguro_Área Urbana"
## [6] "Personas que tienen algún tipo de seguro_Sexo - Hombre"
## [7] "Personas que tienen algún tipo de seguro_Sexo - Mujer"
## [8] "Personas que tienen algún tipo de seguro_Total"
## [9] "Porcentaje de personas analfabetas de 15 años y más_Área Rural"
## [10] "Porcentaje de personas analfabetas de 15 años y más_Área Urbana"
## [11] "Porcentaje de personas analfabetas de 15 años y más_Sexo - Hombre"
## [12] "Porcentaje de personas analfabetas de 15 años y más_Sexo - Mujer"
## [13] "Porcentaje de personas analfabetas de 15 años y más_Total"
## [14] "Porcentaje de Trabajadores Independientes o por cuenta propia no profesionales_Área Rural"
## [15] "Porcentaje de Trabajadores Independientes o por cuenta propia no profesionales_Área Urbana"
## [16] "Porcentaje de Trabajadores Independientes o por cuenta propia no profesionales_Sexo - Hombre"
## [17] "Porcentaje de Trabajadores Independientes o por cuenta propia no profesionales_Sexo - Mujer"
## [18] "Porcentaje de Trabajadores Independientes o por cuenta propia no profesionales_Total"
## [19] "Total de habitantes del censo 2007_Área Rural"
## [20] "Total de habitantes del censo 2007_Área Urbana"
## [21] "Total de habitantes del censo 2007_Sexo - Hombre"
## [22] "Total de habitantes del censo 2007_Sexo - Mujer"
## [23] "Total de habitantes del censo 2007_Total"Pero podemos cambiarlos por nombres más cortos así:
# uso de gsub()
# donde diga "Area Rural" cambiar a "rural":
names(censo2007_wide)=gsub(pattern = "Área Rural",
replacement = "rural",
x = names(censo2007_wide))
names(censo2007_wide)=gsub(pattern = "Área Urbana",
replacement = "urban",
x = names(censo2007_wide))
names(censo2007_wide)=gsub(pattern = "Sexo - Hombre",
replacement = "Hombre",
x = names(censo2007_wide))
names(censo2007_wide)=gsub(pattern = "Sexo - Mujer",
replacement = "Mujer",
x = names(censo2007_wide))
names(censo2007_wide)=gsub(pattern = "Sexo - Mujer",
replacement = "Mujer",
x = names(censo2007_wide))
names(censo2007_wide)=gsub(pattern = "Personas que tienen algún tipo de seguro",
replacement = "conSeguro",
x = names(censo2007_wide))
names(censo2007_wide)=gsub(pattern = "Porcentaje de personas analfabetas de 15 años y más",
replacement = "analfa",
x = names(censo2007_wide))
names(censo2007_wide)=gsub(pattern = "Porcentaje de Trabajadores Independientes o por cuenta propia no profesionales",
replacement = "indep",
x = names(censo2007_wide))
names(censo2007_wide)=gsub(pattern = "Total de habitantes del censo 2007",
replacement = "poblacion",
x = names(censo2007_wide))
# nuevos nombres
names(censo2007_wide)## [1] "Departamento" "Provincia" "Distrito" "conSeguro_rural"
## [5] "conSeguro_urban" "conSeguro_Hombre" "conSeguro_Mujer" "conSeguro_Total"
## [9] "analfa_rural" "analfa_urban" "analfa_Hombre" "analfa_Mujer"
## [13] "analfa_Total" "indep_rural" "indep_urban" "indep_Hombre"
## [17] "indep_Mujer" "indep_Total" "poblacion_rural" "poblacion_urban"
## [21] "poblacion_Hombre" "poblacion_Mujer" "poblacion_Total"¿Cómo nos quedaríamos por ahora con las columnas de “Totales”?
# grep() devuelve las posiciones
grep(pattern = "Total",names(censo2007_wide))## [1] 8 13 18 23Procedamos:
colsTotal=grep(pattern = "Total",names(censo2007_wide))
censo2007_wide_Totales=censo2007_wide[,c(1,2,3,colsTotal)]
# ahora:
head(censo2007_wide_Totales)## Departamento Provincia Distrito conSeguro_Total analfa_Total indep_Total
## 1 Amazonas Bagua Aramango 5057 13.81 54.61
## 2 Amazonas Bagua Copallin 2593 13.46 35.05
## 3 Amazonas Bagua El Parco 524 12.17 20.91
## 4 Amazonas Bagua Imaza 11294 17.29 47.94
## 5 Amazonas Bagua La Peca 14310 7.62 36.51
## 6 Amazonas Bongará Chisquilla 206 3.29 42.95
## poblacion_Total
## 1 11442
## 2 6126
## 3 1274
## 4 21409
## 5 31506
## 6 346Planteemos la primera hipótesis:
A nivel distrital, la cantidad de personas con algún seguro de salud está afectada por el nivel analfabetismo.
De lo visto en las sesiones anteriores (MAGALLANES 2022a, 2022b), podríamos intentar una regresión lineal, que nos daría estos resultados:
library(knitr)
library(modelsummary)
h1=formula(conSeguro_Total~analfa_Total)
rl1=lm(h1, data = censo2007_wide_Totales)
model1=list('OLS asegurados (I)'=rl1)
modelsummary(model1, title = "Resumen de Regresion Lineal",
stars = TRUE,
output = "kableExtra")| OLS asegurados (I) | |
|---|---|
| (Intercept) | 12397.369*** |
| (748.564) | |
| analfa_Total | −444.967*** |
| (45.970) | |
| Num.Obs. | 1831 |
| R2 | 0.049 |
| R2 Adj. | 0.048 |
| AIC | 40986.3 |
| BIC | 41002.8 |
| Log.Lik. | −20490.143 |
| F | 93.694 |
| RMSE | 17531.34 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Como vemos en la Tabla 1.1, el covariado (conocido como predictor o variable independiente) salió con un valor absoluto alto, negativo, y significativo (es muy poco probable - menos de 0.1% - que no tenga efecto), pero con un R-2 ajustado muy bajo. Como el modelo no nos da buen ajuste, es muy probable que la evaluación del modelo no sea satisfactoria. Así, vemos en la Figura 1.1 que difícilmente este modelo puede ser útil.
par(mfrow = c(2, 2))
plot(rl1, 1,caption = '');title(main="Linealidad")
plot(rl1, 2, caption = '');title(main="Normalidad")
plot(rl1, 3, caption = '');title(main="Homocedasticidad")
plot(rl1, 5, caption = '');title(main="Influyentes")
Figure 1.1: Diagnósticos para el modelo OLS asegurados (I)
Podríamos mejorar este modelo si controlásemos el tamaño de la población en el modelo. Veamos la Tabla 1.2.
library(knitr)
library(modelsummary)
h1control=formula(conSeguro_Total~analfa_Total + poblacion_Total )
rl2=lm(h1control, data = censo2007_wide_Totales)
modelslm=list('OLS asegurados (I)'=rl1,'OLS asegurados (II)'=rl2)
modelsummary(modelslm, title = "Regresiones Lineales",
stars = TRUE,
output = "kableExtra")| OLS asegurados (I) | OLS asegurados (II) | |
|---|---|---|
| (Intercept) | 12397.369*** | 710.956*** |
| (748.564) | (189.771) | |
| analfa_Total | −444.967*** | −12.293 |
| (45.970) | (11.186) | |
| poblacion_Total | 0.387*** | |
| (0.002) | ||
| Num.Obs. | 1831 | 1831 |
| R2 | 0.049 | 0.946 |
| R2 Adj. | 0.048 | 0.946 |
| AIC | 40986.3 | 35719.9 |
| BIC | 41002.8 | 35742.0 |
| Log.Lik. | −20490.143 | −17855.965 |
| F | 93.694 | 16156.363 |
| RMSE | 17531.34 | 4159.25 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Se ve una gran mejora en el R2 ajustado con el modelo II de la Tabla 1.2. Veamos sus gráficas de diagnóstico en la Figura 1.2.
par(mfrow = c(2, 2))
plot(rl2, 1,caption = '');title(main="Linealidad")
plot(rl2, 2, caption = '');title(main="Normalidad")
plot(rl2, 3, caption = '');title(main="Homocedasticidad")
plot(rl2, 5, caption = '');title(main="Influyentes")
Figure 1.2: Diagnósticos para el modelo OLS asegurados (II)
Las gráficas de la Figura 1.2 nos muestran un mejor escenario, pero nota que nuestro predictor dejó de ser significativo ante la presencia de la variable de control. Quizá debimos dar un paso más sencillo: analizar la naturaleza de la variable dependiente. Veamos la Figura 1.3 para entenderla mejor.
library(ggplot2)
VarDep=censo2007_wide_Totales$conSeguro_Total
descris=list(min=min(VarDep),
max=max(VarDep),
media=round(mean(VarDep),2),
var=round(var(VarDep),2),
asim=round(e1071::skewness(VarDep),2),
kurt=round(e1071::kurtosis(VarDep),2))
base=ggplot(data=censo2007_wide_Totales, aes(x=conSeguro_Total)) + theme_classic()
hist=base + geom_histogram(bins=50)
histInfo=hist + annotate("text", x = 100000, y = 1050,
color='grey50',
label = paste0("Minimo: ",descris$min))
histInfo = histInfo + annotate("text", x = 100000, y = 900,
color='grey50',
label = paste0("Máximo: ",descris$max))
histInfo = histInfo + annotate("text", x = 100000, y = 650,
color='grey50',
label = paste0("Media: ",descris$media))
histInfo = histInfo + annotate("text", x = 100000, y = 500,
color='grey50',
label = paste0("Varianza: ",descris$var))
histInfo = histInfo + annotate("text", x = 100000, y = 350,
color='grey50',
label = paste0("Asimetría: ",descris$asim))
histInfo = histInfo + annotate("text", x = 100000, y = 200,
color='grey50',
label = paste0("Curtosis: ",descris$kurt))
histInfo
Figure 1.3: Descripción de la Variable Dependiente
El histograma de la Figura 1.3 nos muestra una distribución con asimetría positiva (cola a la derecha) y super leptocurtica. Ello nos hace reflexionar que nuestra variable dependiente representa conteos, valores enteros positivos. La regresión lineal tendrá problemas pues asume que la variable dependiente tiene valores reales y no acotados. Así mismo, el sesgo presente lo aleja de una ‘campana’ de Gauss, por lo que debilitaría los cálculos de significancia si los datos no siguen una tendencia lineal; claro está que uno puede transformar la variable dependiente, pero ello a la vez complicará la interpretación.
2 Regresión Poisson
La regresión Poisson tiene sus supuestos (Glen 2016):
- Variable Respuesta Es un conteo (Y) por unidad de tiempo o espacio, que puede ser descrita por la distribución Poisson (por ejemplo asaltos por dia). Puede además ser un ratio (\(\lambda\)) cuando la unidad de tiempo o espacio varía para cada conteo (por ejemplo votos a favor del total de votos) .
- Independencia Las observaciones (filas) no deben tener relación entre sí.
- Media=Varianza Por definición, la media de una variable que se distribuye como Poisson debe ser igual a su varianza (equidispersión). Si la varianza supera significativamente a la media hablamos de sobredispersión; pero si la media fuera mucho mayor que la varianza tendríamos subdispersión.
- Linealidad El logaritmo de la variable dependiente (los conteos) debe ser una función lineal de los covariados.
Reimplementemos nuestra hipótesis original usando la regresión Poisson. Los resultados los vemos en la Tabla 2.1, comparándolos con el resultado de la regresión lineal controlada por la población.
rp1=glm(h1, data = censo2007_wide_Totales,
offset=log(poblacion_Total), #exposure
family = poisson(link = "log"))
# displaying results
modelslmpoi=list('OLS asegurados (II)'=rl2,
'POISSON asegurados'=rp1)
modelsummary(modelslmpoi, title = "Regresiones OLS y Poisson",
stars = TRUE,
output = "kableExtra")| OLS asegurados (II) | POISSON asegurados | |
|---|---|---|
| (Intercept) | 710.956*** | −0.900*** |
| (189.771) | (0.000) | |
| analfa_Total | −12.293 | 0.005*** |
| (11.186) | (0.000) | |
| poblacion_Total | 0.387*** | |
| (0.002) | ||
| Num.Obs. | 1831 | 1831 |
| R2 | 0.946 | |
| R2 Adj. | 0.946 | |
| AIC | 35719.9 | 895923.3 |
| BIC | 35742.0 | 895934.4 |
| Log.Lik. | −17855.965 | −447959.669 |
| F | 16156.363 | 21585.036 |
| RMSE | 4159.25 | 4353.10 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Ahora que tenemos ambas regresiones en la Tabla 2.1, vemos que el modelo Poisson le devuelve efecto a la independiente. Nótese que la Poisson está modelando los conteos, teniendo en cuenta la exposición (exposure), añadida usando offset. Esto no siempre es necesario, pero en este caso sí lo era, pues necesitamos controlar la exposure de manera explícita (la población) 1. Sin embargo, el AIC de la Poisson es aun mayor que el de la Gaussiana 2 (veremos luego cómo ello mejora en la Tabla 3.5).
2.1 Interpretación
Ya que sabemos cuándo usar una regresión Poisson, alteremos la hipótesis anterior:
A nivel distrital, la cantidad de personas con algun seguro de salud está afectada por el nivel analfabetismo y por la presencia de trabajadores independientes.
h2=formula(conSeguro_Total~analfa_Total + indep_Total)
rp2=glm(h2, data = censo2007_wide_Totales,
offset=log(poblacion_Total),
family = poisson(link = "log"))
modelsPois=list('POISSON asegurados (I)'=rp1,
'POISSON asegurados (II)'=rp2)
modelsummary(modelsPois,
title = "Regresiones Poisson anidadas",
stars = TRUE,
output = "kableExtra")| POISSON asegurados (I) | POISSON asegurados (II) | |
|---|---|---|
| (Intercept) | −0.900*** | −0.615*** |
| (0.000) | (0.001) | |
| analfa_Total | 0.005*** | 0.016*** |
| (0.000) | (0.000) | |
| indep_Total | −0.011*** | |
| (0.000) | ||
| Num.Obs. | 1831 | 1831 |
| AIC | 895923.3 | 771399.1 |
| BIC | 895934.4 | 771415.6 |
| Log.Lik. | −447959.669 | −385696.535 |
| F | 21585.036 | 75214.880 |
| RMSE | 4353.10 | 3748.51 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
La interpretación de la Tabla 2.2 NO es tan directa como lo era en la regresión lineal. Demos una interpretación simple, sin hacer ningún cálculo, de los resultados del modelo II de la Tabla 2.2:
- En el modelo II ambos predictores son significativos;
- En el modelo II, a mayor analfabetismo, mayor cantidad de asegurados;
- En el modelo II, a mayor cantidad de trabajadores independientes, menor cantidad de asegurados.
Sin embargo, la interpretación precisa de los coeficientes de la Tabla 2.2 requiere más cálculos; tengamos primero en cuenta lo que la regresión Poisson ha calculado usando la Ecuación (2.1):
\[\begin{equation} \log(Y) = \log(\lambda) =\alpha + \beta \cdot X + \epsilon \tag{2.1} \end{equation}\]Cuando el exposure es constante modelamos conteos (Y); cuando no lo es modelamos ratios (\(\lambda\))3. Pero, como vemos en la Ecuación (2.1), los coeficiente necesitan ser exponenciados para saber el efecto sobre Y. Veamos la Tabla 2.3.
# formula para limitar a 4 digitos decimales,
# sin que se muestre notación científica:
formatoNum <- function(x) format(x, digits = 4, scientific = FALSE)
modelsummary(modelsPois,
fmt=formatoNum, # uso mi formula
exponentiate = T, # exponenciar!!!!!
statistic = 'conf.int',
title = "Regresión Poisson - coeficientes exponenciados",
stars = TRUE,
output = "kableExtra")| POISSON asegurados (I) | POISSON asegurados (II) | |
|---|---|---|
| (Intercept) | 0.4066*** | 0.5405*** |
| [0.4063, 0.4069] | [0.5396, 0.5415] | |
| analfa_Total | 1.0051*** | 1.0166*** |
| [1.0050, 1.0051] | [1.0165, 1.0167] | |
| indep_Total | 0.9893*** | |
| [0.9893, 0.9894] | ||
| Num.Obs. | 1831 | 1831 |
| AIC | 895923.3 | 771399.1 |
| BIC | 895934.4 | 771415.6 |
| Log.Lik. | −447959.669 | −385696.535 |
| F | 21585.036 | 75214.880 |
| RMSE | 4353.10 | 3748.51 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
# mas simple:
#cbind(exp(coef(rp2)),exp(confint(rp2)))La Tabla 2.3 tiene los coeficientes exponenciados, así mismo, muestra los intervalos de confianza (exponenciados) en vez de los errores típicos. Nota que mientras en la regresión lineal no deseábamos que nuestro coeficiente esté cerca al cero, es decir, que su intervalo de confianza no incluya al cero, aquí no deseamos que el intervalo de confianza incluya al uno.
Una vez exponenciado, podemos interpretar los coeficientes de manera más sencilla. Así, para el modelo II se ve que:
por cada unidad que aumente el analfabetismo la cantidad esperada de asegurados se multiplica por 1.016, es decir, aumentaría en un 1.6% (100 x |1-1.016|) en promedio.
por cada unidad que aumente los trabajadores independientes, la cantidad esperada de asegurados se multiplica por 0.99, es decir, disminuiría en 1% (100 x |1-0.99|) (Choueiry 2022) en promedio.
Nótese que la regresión propone un efecto multiplicativo sobre el promedio de la respuesta (la regresión OLS o Gaussiana propone un efecto aditivo).
3 Equidispersión
Uno de los supuestos en la Regresión Poisson es que la media y la varianza sean iguales. De los estadigrafos de la Figura 1.3 se mostró que estos valores no están cercanos, de hecho la razón varianza - media es:
round(descris$var/descris$media,2)## [1] 51032.13Aquí es clara la sobredispersión, pero en caso haya dudas podemos poner a prueba la hipótesis de equidispersion (Hilbe 2017). Veamos la Tabla 3.1.
library(magrittr)
overdispersion=AER::dispersiontest(rp2,alternative='greater')$ p.value<0.05
underdispersion=AER::dispersiontest(rp2,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 |
La Tabla 3.1 muestra que es altamente improbable que la varianza sea igual a la media, por lo que se opta por aceptar que lo más probable es que tengamos sobredispersión.
3.1 La Quasi Poisson
La presencia de sobredispersión puede tratarse con la quasipoisson; veamos la Tabla 3.3.
rqp = glm(h2, data = censo2007_wide_Totales,
offset=log(poblacion_Total),
family = quasipoisson(link = "log"))
modelsPQP=list('POISSON asegurados (II)'=rp2,'QUASIPOISSON asegurados (II)'=rqp)
modelsummary(modelsPQP,
title = "Regresiones Poisson y QuasiPoisson",
stars = TRUE,
output = "kableExtra")| POISSON asegurados (II) | QUASIPOISSON asegurados (II) | |
|---|---|---|
| (Intercept) | −0.615*** | −0.615*** |
| (0.001) | (0.018) | |
| analfa_Total | 0.016*** | 0.016*** |
| (0.000) | (0.001) | |
| indep_Total | −0.011*** | −0.011*** |
| (0.000) | (0.001) | |
| Num.Obs. | 1831 | 1831 |
| AIC | 771399.1 | |
| BIC | 771415.6 | |
| Log.Lik. | −385696.535 | |
| F | 75214.880 | 179.984 |
| RMSE | 3748.51 | 3748.51 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
La Tabla 3.3 nos muestra cosas interesantes:
- Los coeficientes son los mismos para ambos modelos:
cbind(coef_Poi=coef(rp2),coef_QuasiPoi=coef(rqp))## coef_Poi coef_QuasiPoi
## (Intercept) -0.61523114 -0.61523114
## analfa_Total 0.01649913 0.01649913
## indep_Total -0.01073999 -0.01073999- Pero no los errores típicos:
library(arm)
cbind(se_Poi=se.coef(rp2),se_QuasiPoi=se.coef(rqp))## se_Poi se_QuasiPoi
## (Intercept) 8.910449e-04 0.0182152135
## analfa_Total 4.609011e-05 0.0009421985
## indep_Total 3.037210e-05 0.0006208826La regresión quasipoisson lidia con la sobredispersión al recalcular los errores típicos, lo que afectaría la significancia de los predictores; de ahí que calcula nuevos intervalos de confianza:
modelsQPexp=list('QuasiPoisson asegurados (II) exponenciado'=rqp)
modelsummary(modelsQPexp,fmt=formatoNum,
exponentiate = T,
statistic = 'conf.int',
title = "EXP() de la Regresión Quasi Poisson (II) para Interpretación",
stars = TRUE,
output = "kableExtra")| QuasiPoisson asegurados (II) exponenciado | |
|---|---|
| (Intercept) | 0.5405*** |
| [0.5215, 0.5601] | |
| analfa_Total | 1.0166*** |
| [1.0148, 1.0185] | |
| indep_Total | 0.9893*** |
| [0.9881, 0.9905] | |
| Num.Obs. | 1831 |
| F | 179.984 |
| RMSE | 3748.51 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
3.2 La Binomial Negativa
Una buena alternativa ante la sobredispersión es usar la regresión binomial negativa. Comparemos todas la regresiones exponenciadas, como se ve en la Tabla 3.5.
h2off=formula(conSeguro_Total ~ analfa_Total + indep_Total + offset(log(poblacion_Total)))
rbn=glm.nb(h2off,data=censo2007_wide_Totales)
modelsQP_BN=list('Poisson asegurados (II)'=rp2,
'QuasiPoisson asegurados (II)'=rqp,
'Binomial Negativa asegurados (II)'=rbn)
modelsummary(modelsQP_BN,fmt=formatoNum,
exponentiate = T,
statistic = 'conf.int',
title = "EXP() de la Regresiones Poisson, Quasi Poisson y Binomial Negativa",
stars = TRUE,
output = "kableExtra")| Poisson asegurados (II) | QuasiPoisson asegurados (II) | Binomial Negativa asegurados (II) | |
|---|---|---|---|
| (Intercept) | 0.5405*** | 0.5405*** | 0.4485*** |
| [0.5396, 0.5415] | [0.5215, 0.5601] | [0.4271, 0.4711] | |
| analfa_Total | 1.0166*** | 1.0166*** | 1.0142*** |
| [1.0165, 1.0167] | [1.0148, 1.0185] | [1.0122, 1.0162] | |
| indep_Total | 0.9893*** | 0.9893*** | 0.9951*** |
| [0.9893, 0.9894] | [0.9881, 0.9905] | [0.9940, 0.9963] | |
| Num.Obs. | 1831 | 1831 | 1831 |
| AIC | 771399.1 | 29080.4 | |
| BIC | 771415.6 | 29102.5 | |
| Log.Lik. | −385696.535 | −14536.206 | |
| F | 75214.880 | 179.984 | 106.748 |
| RMSE | 3748.51 | 3748.51 | 4011.51 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
Nótese en la Tabla 3.5 que los coeficientes obtenidos en la regresión binomial negativa son diferentes a los demas. Además, los AIC son mucho mejores para caso de la Binomial Negativa (recuerda los AIC de la Tabla 2.1); los AIC/BIC podemos verlos así de manera alternativa:
data.frame(Model = c("poisson", "quasipoisson", "negative-binomial"),
AIC = c(AIC(rp2), AIC(rqp), AIC(rbn)),
BIC = c(BIC(rp2), BIC(rqp), BIC(rbn)),stringsAsFactors = FALSE
)%>%
kable(caption = "AIC / BIC para los modelos de conteos")%>%kableExtra::kable_styling(full_width = FALSE)| Model | AIC | BIC |
|---|---|---|
| poisson | 771399.07 | 771415.61 |
| quasipoisson | NA | NA |
| negative-binomial | 29080.41 | 29102.46 |
4 Comparación de modelos
La Tabla 3.6 es la mejor alternativa antes modelos no anidados. Una alternativa adicional sería verificar la sobredispersión producida en los modelos:
# poisson case
performance::check_overdispersion(rp2)## # Overdispersion test
##
## dispersion ratio = 417.897
## Pearson's Chi-Squared = 763915.368
## p-value = < 0.001## Overdispersion detected.# quasipoisson case
performance::check_overdispersion(rqp)## # Overdispersion test
##
## dispersion ratio = 417.897
## Pearson's Chi-Squared = 763915.368
## p-value = < 0.001## Overdispersion detected.# negative binomial case
performance::check_overdispersion(rbn)## # Overdispersion test
##
## dispersion ratio = 0.387
## p-value = 0.008## Underdispersion detected.Como estos modelos NO son anidados, usaremos con cuidado la tabla anova, esta vez pidiendo un test chi-cuadrado; veamos el resultado en la Tabla 4.1.
anova(rp2,rqp,rbn,test = "Chisq") %>%
kable(caption = "Tabla ANOVA para comparar modelos")%>%kableExtra::kable_styling(full_width = FALSE)| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 1828 | 754118.559 | NA | NA | NA |
| 1828 | 754118.559 | 0 | 0.0 | NA |
| 1828 | 1868.019 | 0 | 752250.5 | NA |
La caída del Deviance es tanta para el último caso (de 754119 a 1868) que la mejor opción es la binomial negativa. Sin embargo, la quasipoisson, como se vió en la Tabla 3.5 no tiene ni AIC/BIC ni enlace Log, por lo que se produce p-valores como perdidos. Veamos la alternativa del loglikelihood ratio test como se recomienda en Bobbitt (2021) para casos no anidados, pero sin usar la quasipoisson:
lmtest::lrtest(rp2,rbn)%>%
kable(caption = "loglikelihood ratio test")%>%kableExtra::kable_styling(full_width = FALSE)| #Df | LogLik | Df | Chisq | Pr(>Chisq) |
|---|---|---|---|---|
| 3 | -385696.53 | NA | NA | NA |
| 4 | -14536.21 | 1 | 742320.7 | 0 |
La Tabla 4.3 sugiere la binomial negativa. Por lo general, la binomial negativa es más utilizada que la quasipoisson, pero la binomial negativa no es apropiada para la subdispersión, mientras que la quasipoisson sí se usa para ese caso. Una manera adicional de comparar es la gráfica. Así, la Tabla 3.5 se puede ver de manera gráfica en la Figura 4.1.
library(ggplot2)
dotwhisker::dwplot(list(Poisson=rp2,CuasiPoisso=rqp,BinomialNegativa=rbn),exp=T) + scale_y_discrete(labels=c("trabajo\nindependiente","analfabetismo")) + scale_color_discrete(name="Modelos para:\nCantidad de Asegurados") + geom_vline(
xintercept = 1,
colour = "grey60",
linetype = 2
)
Figure 4.1: Comparación visual de modelos
Finalmente, podemos calcular los coeficientes estandarizados (IBM SPSS 2020) para saber cuál de los predictores tiene mayor efecto (ver Tabla 4.5).
sdVD=sd(censo2007_wide_Totales$conSeguro_Total)
sdVIs=apply(censo2007_wide_Totales[,c("analfa_Total","indep_Total")],2,sd)
DF=list(Poisson=sdVIs*coef(rp2)[c(2,3)]/sdVD,
CuasiPoisson=sdVIs*coef(rqp)[c(2,3)]/sdVD,
BinomNegativa=sdVIs*coef(rbn)[c(2,3)]/sdVD)%>%
data.frame()
DF%>% kable(caption = "Coeficientes Standarizados (ordenar vía valores absolutos)")%>%
kableExtra::kable_styling(full_width = F)| Poisson | CuasiPoisson | BinomNegativa | |
|---|---|---|---|
| analfa_Total | 8.2e-06 | 8.2e-06 | 7.0e-06 |
| indep_Total | -8.4e-06 | -8.4e-06 | -3.8e-06 |