En este ejemplo vamos a ver la relación que existe entre el desarrollo económico y nivel de democracia en Latinoamérica. Esta idea viene de un trabajo seminal de Przewroski et al (2000), pero nos vamos a limitar en analizar L.A. Los datos los vamos a generar usando APIs del Banco Mundial y Polity Project.

Obviamente con este script puedes analizar dicha relación usando una base de datos mucho más extensiva (i.e. agregando otras regiones o países).

Cargamos las librerías que necesitaremos

library(WDI) # para bajar datos de World Development Indicators
library(countrycode) # para hacer merge entre datos de WDI y Polity
library(dplyr) # para manipular datos
library(tidyverse) # para manipular datos
library (lmtest) # para hacer diagnósticos de nuestra regresión
library(ggplot2) # para graficar
library(stargazer) # para hacer tablas bonitas
library(foreign) # para leer un archivo en SPSS con los datos de Polity
library(scales) # para re-escalar variables
library(sandwich) # para esimar errores estándar robustos

Creamos nuestra base de datos

Para crear la base que necesitaremos en este ejercicio usaremos APIs del Banco Mundial y Polity Project. Del Banco Mundial bajaremos datos de World Development Indicators Database sobre PIB per capita y crecimiento del PIB. Del Polity Project vamos a llamar al score de democracia para cada país-año en nuestra muestra.

wdi <- WDI(country=c('AR', 'BO', 'BR', 'CO', 'CL', 'CR', 'EC', 'GT', 'HN', 'MX', 'NI', 'PA', 'PE', 'PY', 'SV', 'UY','VE'), 
           indicator=c('NY.GDP.MKTP.KD.ZG', 'NY.GDP.PCAP.KD'), 
           start = 1970, end = 2015, 
           extra = TRUE)
wdi$cowcode <- countrycode(wdi$iso2c, "iso2c", "cown") # para hacer compatible con datos de POLITY (códigos de país)
wdi <- rename(wdi, gdppc = NY.GDP.PCAP.KD, gdpgrowth = NY.GDP.MKTP.KD.ZG) # renombramos variable de GDPPC

polity <- read.spss("http://www.systemicpeace.org/inscr/p4v2016.sav", to.data.frame = TRUE) # para bajar los datos de POLITY
polity <- subset(polity, year >=1970 & year<=2015) # limitamos periodo de 1970 a 2015
variables <- c("ccode", "scode", "country", "year", "polity2") # seleccionamos las variables que necesitamos
polity <-polity[variables]  # creamos data.frame únicamente con variables que necestiamos
polity <- rename(polity, cowcode = ccode) # renombramos variable para poder hacer merge con wdi
polity$polity2.r <- rescale(polity$polity2, to = c(0, 10)) # reescalamos el índice de polity de rango (-10 a 10) a rabgo (0 a 1) 

## Ahora hacemos un merge de ambas bases

wdi_polity <- merge(wdi, polity, by = c("cowcode", "year")) # merge en cowcode y year
wdi_polity <- select(wdi_polity, gdpgrowth, gdppc, polity2, polity2.r, year, country.x, cowcode) # seleccionamos variables que usaremos
wdi_polity$loggdppc <- log(wdi_polity$gdppc) # logaritmo de GDPPC

El data.frame que usaremos se llama wdi_polity. Estamos listos para iniciar nuestro análisis.

Análisis Preliminar

Przeworski et al. (2000) analizan la relación entre desarrollo económico (pensemos, por ejemplo, en niveles de PIB per cápita democracia. Específicamente, los autores sugieren que mayores niveles de desarrollo económico resultan en niveles de democracia más altos. Vamos a ver si dicha afirmación encuentra sustento cuando usamos datos exclusivamente de Latinoaméria en el periodo 1970-2015. Para ver la relación entre las variables gdppc y polity2.r usemos primero un scater plot.

ggplot(wdi_polity, aes(x=gdppc, y=polity2.r)) + 
  geom_point( color="#69b3a2")                         # scatter plot básico
## Warning: Removed 46 rows containing missing values (`geom_point()`).

ggplot(wdi_polity, aes(x=gdppc, y=polity2.r)) +
  geom_point() +
  geom_smooth(method=lm , color="red", se=TRUE)  # scatetter plot con tendencia lineal e I.C.
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 46 rows containing non-finite values (`stat_smooth()`).
## Removed 46 rows containing missing values (`geom_point()`).

Hay demasiadas observaciones de cada variable para cada país y el scatter plot se ve bastante desordenado. Pero la relación es positiva.

Una forma alternativa de ver esto es agrupando por país de la siguiente forma.

promedio <- wdi_polity %>% # simplemente estoy agrupando mis datos por país usando las medias d variables
  drop_na() %>%
  group_by(country.x) %>%
  summarise(m_gdppc = mean(gdppc), m_polity = mean(polity2.r))

ggplot(promedio, aes(x = m_gdppc, y = m_polity)) +
  geom_point() +
  geom_smooth(method=lm , color="red", se=TRUE) + # scatetter plot con tendencia lineal e I.C.
  geom_text(aes(label = country.x, hjust = - 0.1)) 
## `geom_smooth()` using formula = 'y ~ x'

Este gráfico esta mucho más ordenado pero la relación positiva no es tan clara (la pendiente es mucho menor). Vamos a ver que nos dice la regresión sobre esta relación.

Estimación de Regresión Lineal Simple

Usamos la función lm para estimar nuestra regresión. Recordemos que nuestra variable dependiente es el nivel de democracia (i.e. polity2.r y la variable independiente es el nivel de desarrollo económico (gdppc).

modelo1 <-lm(polity2.r ~ gdppc, data = wdi_polity) # estimamos modelo de OLS 

stargazer(modelo1, type = "html") # crea tabla con resultados
Dependent variable:
polity2.r
gdppc 0.0002***
(0.00003)
Constant 5.939***
(0.215)
Observations 736
R2 0.050
Adjusted R2 0.049
Residual Std. Error 3.013 (df = 734)
F Statistic 38.754*** (df = 1; 734)
Note: p<0.1; p<0.05; p<0.01

El coeficiente en la variable gdppc es positivo y estadísticamente significativo con un 95% de confianza. Pero, es un efecto muy pequeñito (hay que notar que las variables están en escalas muy diferentes). Para interpretar nuestro coeficiente de regresión: \(\delta\)y = \(\beta\)*\(\delta\)x

A simple vista, sin embargo, el modelo tiene una \(R^2\) muy baja (0.050) lo que indica que nuestro modelo no es muy bueno explicando la variación en nuestra variable dependiente.

Vamos a hacer algunas pruebas de diagnóstico para tener más información sobre la validez de nuestro modelo.

Pruebas de Diagnóstico

Aquí vamos a evaluar la calidad de nuestro modelo. Esta es una regresión lineal simple, por lo que no debemos preocuparnos mucho sobre multicolinealidad. Sin embrago, vamos a evaluar que tan bien (o mal) nuestro modelo cumple con los supuestos de OLS (LINI).

1. Evaluando linealidad

plot(modelo1, which = 1)

Cuando la línea roja es horizontal (o casi horizontal) entonces no tenemos problemas de linealidad.

2. Evaluando Independencia de Errores

bgtest(modelo1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelo1
## LM test = 612.29, df = 1, p-value < 2.2e-16

p es menor a .05. Por lo tanto hay un problema de correlación.

3. Evaluando Normalidad de Errores

Vamos a ver al Q-Q plot

plot(modelo1, which = 2)

Este es un claro ejemplo de la violación del supuesto de normalidad en los errores. Pero no es sorprendente… tenemos una variable dependiente que NO es continua.

4. Homoesquedasticidad de Errores

En este caso lo que estamos buscando es que no haya un patrón claro en la distribución de los residuales. Veamos.

plot(modelo1, which = 3)

Hay un patrón claro en los residuales. Por lo tanto los errores son heteroesqudasticos (no hay homoesquedasticidad).