Regresion usando datos del Barómetro de las Américas (1)
Arturo Maldonado
10/12/2020
En este documento veremos los aspectos más importantes del análisis de regresión lineal simple. Se partirá de un ejemplo similar al documento de correlación que puede ver aquí, analizando el gráfico de dispersión y la línea de regresión, pero se añadirá el análisis sobre los resultados de un modelo de regresión lineal simple. La idea de este documento es trabajar con los 18 países que se analizaron en el documento de correlación. Esta vez usaremos los datos sobre la letalidad del COVID-19 como variable dependiente y buscaremos en la base de datos del Barómetro de las Américas, variables independientes que puedan explicar las diferencias en la tasa de letalidad (número de muertos por cada millón de habitantes) entre países.
SOBRE LOS DATOS
Los datos de afectados por el COVID-19 en América Latina los hemos encontrado en Google. Se realizó una búsqueda “síntomas de Coronavirus en América Latina” y nos llevó a una página web donde se podía encontrar el número de casos y muertes, en total y por cada millón de personas. Esta búsqueda se puede encontrar aquí. En esta web se buscaron los datos de cada uno de los 18 países analizados en el reporte “El pulso de la democracia” actualizados al 9 de diciembre de 2020. Estos datos se registraron en un archivo .xlsx, donde si incluye además el código de país que usa LAPOP para luego integrar los datos del COVID-19 con los datos de cultura política de LAPOP. En primer lugar cargaremos los datos sobre COVID-19.
library(rio)
covid <- import("covid.xlsx")
covid## country pais covid
## 1 Mexico 1 880
## 2 Guatemala 2 292
## 3 El Salvador 3 182
## 4 Honduras 4 321
## 5 Nicaragua 5 28
## 6 Costa Rica 6 366
## 7 Panama 7 879
## 8 Colombia 8 807
## 9 Ecuador 9 833
## 10 Bolivia 10 812
## 11 Peru 11 1157
## 12 Paraguay 12 272
## 13 Chile 13 819
## 14 Uruguay 14 26
## 15 Brasil 15 864
## 16 Argentina 17 880
## 17 Rep. Dom. 21 221
## 18 Jamaica 23 88Estos datos se cargan en un dataframe. Este dataframe incluye tres variables “country”, una variable de tipo “character” con los nombres de los países, la variable “pais” con los códigos LAPOP de cada país, y la variable “covid” con los datos de letalidad (muertes por cada millón de habitantes). Luego cargamos la base de datos del Barómetro de las Américas y seleccionamos la última ronda 2018/19 y los 18 países de la muestra.
library(haven)
lapop <- read_dta("/Users/Arturo/OneDrive - Vanderbilt/C LAPOP/Data/LAPOP_Merge_2004_2018.dta")
lapop <- subset(lapop, wave==2018)
lapop <- subset(lapop, pais<=35)Una vez que hemos leído la base de datos, trabajamos con la variable país, la convertimos en factor y la etiquetamos.
lapop$pais = as.factor(lapop$pais)
levels(lapop$pais) <- c("México", "Guatemala", "El Salvador", "Honduras",
"Nicaragua","Costa Rica", "Panamá", "Colombia",
"Ecuador", "Bolivia", "Perú", "Paraguay",
"Chile", "Uruguay", "Brasil",
"Argentina", "Rep. Dom.", "Jamaica")
levels(lapop$pais)## [1] "México" "Guatemala" "El Salvador" "Honduras" "Nicaragua"
## [6] "Costa Rica" "Panamá" "Colombia" "Ecuador" "Bolivia"
## [11] "Perú" "Paraguay" "Chile" "Uruguay" "Brasil"
## [16] "Argentina" "Rep. Dom." "Jamaica"table(lapop$pais)##
## México Guatemala El Salvador Honduras Nicaragua Costa Rica
## 1580 1596 1511 1560 1547 1501
## Panamá Colombia Ecuador Bolivia Perú Paraguay
## 1559 1663 1533 1682 1521 1515
## Chile Uruguay Brasil Argentina Rep. Dom. Jamaica
## 1638 1581 1498 1528 1516 1513AGREGANDO LOS DATOS DE LA VARIABLE INDEPENDIENTE
El cuestionario del Barómetro de las Américas, que se puede ver aquí, incluye la variable SD6NEW2. ¿Y con la calidad de los servicios médicos y de salud públicos? ¿Está usted… 1. Muy satisfecho(a) 2. Satisfecho(a) 3. Insatisfecho(a) 4. Muy insatisfecho(a)
En los siguientes resultados se puede ver la distribución de respuestas a esta variable.
round(prop.table(table(lapop$sd6new2)), 3)*100##
## 1 2 3 4
## 6.1 37.0 39.8 17.0Vamos a usar esa pregunta, agregada por país, como variable independiente explicativa de la tasa de letalidad del COVID-19 en los 18 países analizados. La hipótesis que manejamos es la evaluación del sistema de sanidad público de los ciudadanos está fuertemente correlacionada con el estado de la sanidad en un país. Es decir, si los ciudadanos tienen una opinión insatisfactoria, es porque el sistema no trabaja bien y, entonces, sería muy probable que este sistema no haya respondido bien a la crisis originada por la pandemia, originando que un país tenga una tasa de letalidad más alta.
Lo primero que tenemos que hacer es recodificar la variable SD6NEW2 de tal manera de recoger la insatisfacción con el sistema de salud, de la siguiente manera:
- Valores 1 y 2 de SD6NEW2 serán 0 en la nueva variable
- Valores 3 y 4 de SD6NEW2 serán 100 en la nueva variable
library(dplyr)
library(car)
lapop$sd6new2 <- as.numeric(lapop$sd6new2)
lapop$vi <- recode(lapop$sd6new2, "1:2=0; 3:4=100")
table(lapop$vi)##
## 0 100
## 5879 7743Con esta nueva variable “vi”, ahora tenemos que agregar los datos de esta variable por país y guardar esta información en un nuevo dataframe “df”.
library(Rmisc) #para poder utilizar el comando summarySE
df <- summarySE(data=lapop, measurevar="vi", groupvar="pais", na.rm=T)
df## pais N vi sd se ci
## 1 México 1539 55.75049 49.68436 1.266486 2.484222
## 2 Guatemala 1540 56.75325 49.55793 1.262853 2.477095
## 3 El Salvador 1481 54.49021 49.81479 1.294437 2.539126
## 4 Honduras 1504 57.71277 49.41798 1.274269 2.499534
## 5 Nicaragua 0 NaN NA NA NaN
## 6 Costa Rica 0 NaN NA NA NaN
## 7 Panamá 0 NaN NA NA NaN
## 8 Colombia 1642 70.40195 45.66213 1.126859 2.210233
## 9 Ecuador 0 NaN NA NA NaN
## 10 Bolivia 0 NaN NA NA NaN
## 11 Perú 1490 68.59060 46.43098 1.202859 2.359478
## 12 Paraguay 1488 54.09946 49.84841 1.292260 2.534846
## 13 Chile 0 NaN NA NA NaN
## 14 Uruguay 0 NaN NA NA NaN
## 15 Brasil 0 NaN NA NA NaN
## 16 Argentina 0 NaN NA NA NaN
## 17 Rep. Dom. 1483 37.82873 48.51234 1.259742 2.471067
## 18 Jamaica 1455 54.43299 49.82022 1.306093 2.562028Como se observa, lamentablemente esta variable solo ha sido recogida en 9 de los 18 países incluidos en el análisis. De todas maneras, como vamos a empezar por hacer un análisis de regresión lineal simple, podemos trabajar solo con 9 casos. El siguiente paso implica incluir los datos de la variable “vi” en un dataframe con los datos de letalidad del COVID-19.
covid <- cbind(covid, df$vi)
colnames(covid)[4] <- "vi" #nombramos la columna de manera más manejable. Por defecto se nombra df$vi haciendo que se tenga que usar ´´ para usar este nombre de variable en los códigos
covid## country pais covid vi
## 1 Mexico 1 880 55.75049
## 2 Guatemala 2 292 56.75325
## 3 El Salvador 3 182 54.49021
## 4 Honduras 4 321 57.71277
## 5 Nicaragua 5 28 NaN
## 6 Costa Rica 6 366 NaN
## 7 Panama 7 879 NaN
## 8 Colombia 8 807 70.40195
## 9 Ecuador 9 833 NaN
## 10 Bolivia 10 812 NaN
## 11 Peru 11 1157 68.59060
## 12 Paraguay 12 272 54.09946
## 13 Chile 13 819 NaN
## 14 Uruguay 14 26 NaN
## 15 Brasil 15 864 NaN
## 16 Argentina 17 880 NaN
## 17 Rep. Dom. 21 221 37.82873
## 18 Jamaica 23 88 54.43299DIAGRAMA DE DISPERSIÓN Y MODELO DE REGRESIÓN
En este dataframe tenemos la variable “covid” que usaremos como variable dependiente y la variable “vi” (insatisfacción con los servicios de salud) que usaremos como variable independiente. En primer lugar graficamos el diagrama de dispersión.
plot(covid$vi, covid$covid,
xlab="% de Insatisfacción con los Servicios de Salud",
ylab="Tasa de Letalidad del COVID-19",
pch=19, xlim=c(30, 80), ylim=c(10, 1200))
text(covid$vi, covid$covid, labels=covid$country, cex=0.5, pos=3)
En este gráfico se observa que podría existir una relación directa, es decir, que países cuyos ciudadanos son más críticos de los servicios de salud tienen, en promedio, una mayor tasa de letalidad. Esto lo analizaremos en más detalle con un modelo de regresión lineal simple.
modelo <- lm(covid~vi, data=covid, na.action = na.exclude)
summary(modelo)##
## Call:
## lm(formula = covid ~ vi, data = covid, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -320.6 -179.0 -127.6 259.5 436.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1056.92 647.18 -1.633 0.1465
## vi 26.92 11.28 2.386 0.0484 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 299.4 on 7 degrees of freedom
## (9 observations deleted due to missingness)
## Multiple R-squared: 0.4485, Adjusted R-squared: 0.3698
## F-statistic: 5.694 on 1 and 7 DF, p-value: 0.04844cor.test(x = covid$vi, y = covid$covid, method = "pearson")##
## Pearson's product-moment correlation
##
## data: covid$vi and covid$covid
## t = 2.3862, df = 7, p-value = 0.04844
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.01011462 0.92322190
## sample estimates:
## cor
## 0.6697375De acuerdo con el test de ANOVA, se obtiene un p-value < 0.05 (0.048), que es el mismo valor que el del coeficiente de la variable independiente. Esto indica que el modelo es válido y que podemos rechazar la hipótesis nula de que la recta de regresión tiene pendiente cero y aceptar la hipótesis alternativa que la pendiente es diferente de cero, indicando una relación entre la variable “insatisfacción con los servicios de salud” y “tasa de letalidad por COVID-19”. De acuerdo al signo positivo del coeficiente de la variable independiente podemos concluir que efectivamente la relación es directa. De acuerdo a nuestra hipótesis, podemos brindar evidencia que países con mayor insatisfacción con el sistema de salud muestran, en promedio, tasas más altas de letalidad por COVID-19. El coeficiente de determinación R nos indica que mediante la variable de insatisfacción con los servicios de salud se explica el 45% de la variación total de la variable de letalidad por COVID-19. El coeficiente r de Pearson, a su vez, indica que las observaciones se ajustan medianamente bien a la recta de predicción. La ecuación de la recta de predicción es:
\[\begin{align*} \hat{y} = \alpha + \beta X \\ \hat{y} = -1059.92 + 26.92*X \end{align*}\]
(OJO: por notación distinguimos \(\hat{y}\) que es el valor predicho por el modelo de y que es el valor observado en la realidad)
Con esta ecuación, si calculáramos el valor predicho para Perú sería 789.5 (=-1059.92+26.92*65.6). El valor observado para Perú es 1157, por lo que el error o residuo del modelo para esta obervación es 367.5.
Finalmente, una vez calculado el modelo, podemos también incluir la recta de predicción en el gráfico de dispersión.
plot(covid$vi, covid$covid,
xlab="% de Insatisfacción con los Servicios de Salud",
ylab="Tasa de Letalidad del COVID-19",
pch=19, xlim=c(30, 80), ylim=c(10, 1200))
text(covid$vi, covid$covid, labels=covid$country, cex=0.5, pos=3)
abline(modelo)
PRESENTACIÓN DE RESULTADOS
Para tener un formato más adecuado de presentación de resultados, podemos usar la librería jtools.
library(jtools)
summ(modelo)## MODEL INFO:
## Observations: 9 (9 missing obs. deleted)
## Dependent Variable: covid
## Type: OLS linear regression
##
## MODEL FIT:
## F(1,7) = 5.69, p = 0.05
## R² = 0.45
## Adj. R² = 0.37
##
## Standard errors: OLS
## -----------------------------------------------------
## Est. S.E. t val. p
## ----------------- ---------- -------- -------- ------
## (Intercept) -1056.92 647.18 -1.63 0.15
## vi 26.92 11.28 2.39 0.05
## -----------------------------------------------------Tenemos que tomar en cuenta que el objeto “modelo” que se crea con el comando lm guarda los vectores de los valores predichos (fitted.values), los residuos (residuals). Para ver los valores predichos, por ejemplo, podemos usar el siguiente código.
modelo$fitted.values## 1 2 3 4 8 11 12 17
## 444.0419 471.0391 410.1116 496.8721 838.5013 789.7347 399.5916 -38.4633
## 18
## 408.5711Y para los residuos:
modelo$residuals## 1 2 3 4 8 11 12
## 435.95809 -179.03907 -228.11161 -175.87208 -31.50127 367.26533 -127.59159
## 17 18
## 259.46330 -320.57110Con estos datos comprobamos que la diferencia entre el valor observado y el valor predicho para Perú (país 11) es 367.3.
TRABAJANDO CON OTRA VARIABLE INDEPENDIENTE
El cuestionario del Barómetro de las Américas también incluye la pregunta M1. Hablando en general acerca del gobierno actual, ¿diría usted que el trabajo que está realizando el Presidnete [Nombre] es…?
- Muy bueno
- Bueno
- Ni bueno, ni malo (regular)
- Malo
- Muy malo (pésimo)
Esta variable se distribuye de la siguiente manera.
table(lapop$m1)##
## 1 2 3 4 5
## 2169 7879 10455 3535 3509Podemos argumentar que aquellos que una peor calificación del trabajo del presidente (incluso si ha sido un presidente anterior) significaría que el Ejecutivo no ha trbajado adecuadamente y, por lo tanto, luego esto podría tener consecuencias en la respuesta del Estado frente a la pandemia. Esta variable, entonces, la tenemos que recodificar para capturar la proporción de ciudadnos críticos con el trabajo del presidente en cada país. Para esto usaremos la siguiente regla de recodificación:
- Puntajes del 1 al 3 en la variable m1 se les asignará el valor de 0 en la nueva variable
- Puntajes de 4 y 5 en m1 se les asignará el valor de 100 en la nueva variable
lapop$m1 <- as.numeric(lapop$m1)
lapop$vi2 <- recode(lapop$m1, "1:3=0; 4:5=100")
table(lapop$vi2)##
## 0 100
## 20503 7044Con esta nueva variable, calculamos el porcentaje de ciudadanos críticos con el trabajo del presidente por país.
df1 <- summarySE(data=lapop, measurevar="vi2", groupvar="pais", na.rm=T)
df1## pais N vi2 sd se ci
## 1 México 1548 3.682171 18.83847 0.4788066 0.9391785
## 2 Guatemala 1555 36.527331 48.16616 1.2214532 2.3958704
## 3 El Salvador 1486 30.215343 45.93460 1.1916000 2.3373983
## 4 Honduras 1521 36.949375 48.28265 1.2380166 2.4284016
## 5 Nicaragua 1498 29.839786 45.77078 1.1825851 2.3196997
## 6 Costa Rica 1489 38.146407 48.59093 1.2592380 2.4700703
## 7 Panamá 1553 48.293625 49.98697 1.2684434 2.4880437
## 8 Colombia 1617 19.913420 39.94727 0.9934182 1.9485233
## 9 Ecuador 1522 27.201051 44.51415 1.1410135 2.2381264
## 10 Bolivia 1674 13.739546 34.43672 0.8416743 1.6508457
## 11 Perú 1504 10.771277 31.01204 0.7996620 1.5685718
## 12 Paraguay 1478 14.073072 34.78613 0.9048339 1.7748962
## 13 Chile 1614 22.862454 42.00770 1.0456279 2.0509319
## 14 Uruguay 1556 28.920308 45.35382 1.1497649 2.2552531
## 15 Brasil 1450 11.448276 31.85065 0.8364392 1.6407612
## 16 Argentina 1514 61.294584 48.72372 1.2522108 2.4562531
## 17 Rep. Dom. 1501 15.789474 36.47638 0.9415023 1.8468008
## 18 Jamaica 1467 10.702113 30.92455 0.8073993 1.5837812Ahora, se tiene que añadir los datos agregados de esta nueva variable en la base de datos del covid.
covid <- cbind(covid, df1$vi2)
colnames(covid)[5] <- "vi2" #nombramos la columna de manera más manejable. Por defecto se nombra df$vi haciendo que se tenga que usar ´´ para usar este nombre de variable en los códigos
covid## country pais covid vi vi2
## 1 Mexico 1 880 55.75049 3.682171
## 2 Guatemala 2 292 56.75325 36.527331
## 3 El Salvador 3 182 54.49021 30.215343
## 4 Honduras 4 321 57.71277 36.949375
## 5 Nicaragua 5 28 NaN 29.839786
## 6 Costa Rica 6 366 NaN 38.146407
## 7 Panama 7 879 NaN 48.293625
## 8 Colombia 8 807 70.40195 19.913420
## 9 Ecuador 9 833 NaN 27.201051
## 10 Bolivia 10 812 NaN 13.739546
## 11 Peru 11 1157 68.59060 10.771277
## 12 Paraguay 12 272 54.09946 14.073072
## 13 Chile 13 819 NaN 22.862454
## 14 Uruguay 14 26 NaN 28.920308
## 15 Brasil 15 864 NaN 11.448276
## 16 Argentina 17 880 NaN 61.294584
## 17 Rep. Dom. 21 221 37.82873 15.789474
## 18 Jamaica 23 88 54.43299 10.702113En la tabla se observa que, contrario a la variable sobre evaluación del sistema de salud, donde solo se tenía información agregada para 9 países, en este caso se cuenta con información para los 18 países evaluados.
DIAGRAMA DE DISPERSIÓN Y MODELO DE REGRESIÓN DE SEGUNDA VARIABLE INDEPENDIENTE
Con este dataframe “covid” tenemos los datos para calcular el modelo de regresión lineal y el gráfico de dispersión.
plot(covid$vi2, covid$covid,
xlab="% de Mala Evaluación del Presidente",
ylab="Tasa de Letalidad del COVID-19",
pch=19, xlim=c(0, 80), ylim=c(10, 1200))
text(covid$vi2, covid$covid, labels=covid$country, cex=0.5, pos=3)
A diferencia de la primera variable independiente, en este caso no es tan claro cuál es la relación entre las variables observando el gráfico. Para asegurarnos acerca de la relación, se calcula el modelo de regresión lineal.
modelo2 <- lm(covid~vi2, data=covid, na.action = na.exclude)
summary(modelo2)##
## Call:
## lm(formula = covid ~ vi2, data = covid, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -508.77 -323.81 49.23 298.74 591.73
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 583.37 181.37 3.216 0.00539 **
## vi2 -1.68 6.17 -0.272 0.78882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 379.3 on 16 degrees of freedom
## Multiple R-squared: 0.004615, Adjusted R-squared: -0.0576
## F-statistic: 0.07418 on 1 and 16 DF, p-value: 0.7888cor.test(x = covid$vi2, y = covid$covid, method = "pearson")##
## Pearson's product-moment correlation
##
## data: covid$vi2 and covid$covid
## t = -0.27237, df = 16, p-value = 0.7888
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5183636 0.4120032
## sample estimates:
## cor
## -0.06793416En este modelo, el p-value > 0.05 (0.79, tanto el de ANOVA como el del coeficiente) nos lleva a concluir que el modelo no es válido, es decir, que la variable independiente, mala evaluación del presidente, no explica la tasa de letalidad por COVID-19. Incluso, mirando el gráfico de dispersión con la recta de predicción, se observa que esta línea tiene una pendiente marginalmente negativa, lo que iría en sentido contrario a la hipótesis planteada. En todo caso, de acuerdo al p-value del coeficiente (es decir, de la pendiente) de la variable independiente, no se puede rechazar la hipótesis nula, que indica que la pendiente es cero.
plot(covid$vi2, covid$covid,
xlab="% de Mala Evaluación del Presidente",
ylab="Tasa de Letalidad del COVID-19",
pch=19, xlim=c(0, 80), ylim=c(10, 1200))
text(covid$vi2, covid$covid, labels=covid$country, cex=0.5, pos=3)
abline(modelo2)
Visto de manera más ordenada en la siguiente tabla, llegamos a la conclusión que la mala evaluación del presidente no expplica la variación de la tasa de letalidad por país.
summ(modelo2)## MODEL INFO:
## Observations: 18
## Dependent Variable: covid
## Type: OLS linear regression
##
## MODEL FIT:
## F(1,16) = 0.07, p = 0.79
## R² = 0.00
## Adj. R² = -0.06
##
## Standard errors: OLS
## ---------------------------------------------------
## Est. S.E. t val. p
## ----------------- -------- -------- -------- ------
## (Intercept) 583.37 181.37 3.22 0.01
## vi2 -1.68 6.17 -0.27 0.79
## ---------------------------------------------------Se puede especular que esta ausencia de relación entre las variable pueda ser debido a que la evaluación del presidente no necesariamente es acerca del presidente que está enfrentando la pandemia. A pesar que los ciudadanos podrían asumir que la buena a mala respuesta del Estado a la pandemia es debida a las gestiones presidenciales recientes, y no solo la última, al parecer no siempre es así. El trabajo de campo de las encuestas del Barómetro de las Américas de 2018/19 recogió información a inicios del 2018 (incluso a finales de 2017), por lo que podría haber pasado hasta 3 años hasta el inicio de la pandemia, pudiendo haber cambiado la administración en ese periodo.