Introducción a las regresiones.
Juan José Echeverry
2023-09-24
Fuente de datos
Primer paso: Cargar e instalar las librerias
- Librerias a usar:
library(car)
library(skimr)
library(ggplot2)
library(ggcorrplot)
library(tidyverse)
library(texreg)
library(prediction)
library(lmtest)
library(sandwich)
library(miceadds)
library(gridExtra)
library(ggpubr)
library(pastecs)
Segundo paso: Cargar y analizar los datos
Estadísticas descriptivas y distribución de las variables en el modelo
Antes de estimar un modelo lineal con mínimos cuadrados ordinarios (MCO) se recomienda identificar primero la distribución de las variables que te interesan: la variable dependiente y (también llamada variable de respuesta) y la variable independiente de interés x (también llamada variable explicativa o regresiva).
El objetivo es prestar atención a las siguientes cuestiones:
Variación en x : que las variables independientes (pero sobre todo la de interés) tengan variación en nuestra muestra. Si no hay variación en x, no podemos estimar cómo esta variación afecta a la variación de y.
Variación en y: si la variable dependiente no varía, no podemos explicar su variación según la variación de x.
Unidad de medida de las variables: aquí es donde evaluamos cómo se miden nuestras variables (además de revisar los libros de códigos que suelen venir con las bases de datos con los que estamos trabajando), para poder entender qué tipos de variables son (nominales, ordinales, continuas), y también para interpretar correctamente los resultados más adelante.
Tipo de variables: en la estimación del MCO, la variable dependiente debe ser, en general, continua, (aunque es posible trabajar con variables dependientes dicotómicas). Por lo tanto, debemos asegurarnos de que la variable dependiente es continua y numérica. Además, es importante conocer el tipo de variable independiente y comprobar que su tipo es coherente con su formato de codificación (es decir, si tenemos una variable independiente de “rangos de edad,” la variable debe ser categórica o factorial, no numérica), para que nuestras interpretaciones de los resultados sean correctas. Los problemas con el formato de codificación son muy comunes, ¡ten cuidado!
Identificar los valores perdidos: si nuestras variables tienen demasiados valores perdidos, debemos comprobar dónde están (en qué variables, para cuales observaciones) y, eventualmente, decidir si una imputación es deseable y posible
#Quitamos notacion cientifica
options(scipen=999)
#Cargamos datos
welfare <- read_excel("Data_set.xlsx")
#Revisamos variable
skimr::skim(welfare)| Name | welfare |
| Number of rows | 1074 |
| Number of columns | 14 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 13 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| country_id | 0 | 1 | 3 | 3 | 0 | 25 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| year | 0 | 1.00 | 1990.98 | 12.41 | 1970.00 | 1980.00 | 1991.00 | 2002.00 | 2012.00 | ▇▇▇▇▇ |
| population | 129 | 0.88 | 36.36 | 7.01 | 18.99 | 30.79 | 36.57 | 42.26 | 48.73 | ▂▆▇▇▇ |
| gini | 649 | 0.40 | 49.54 | 6.73 | 28.90 | 44.00 | 50.60 | 55.00 | 68.30 | ▁▅▇▇▁ |
| sector_dualism | 606 | 0.44 | -7.19 | 11.61 | -36.66 | -14.59 | -9.13 | -0.84 | 18.51 | ▁▅▇▅▃ |
| gdp | 194 | 0.82 | 7345.94 | 4981.75 | 1421.71 | 4046.04 | 6264.84 | 8667.02 | 29343.90 | ▇▅▁▁▁ |
| foreign_inv | 198 | 0.82 | 2.19 | 4.18 | -55.24 | 0.48 | 1.43 | 3.29 | 39.81 | ▁▁▇▆▁ |
| ethnic_diversity | 370 | 0.66 | 0.45 | 0.50 | 0.00 | 0.00 | 0.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| regime_type | 0 | 1.00 | 3.03 | 1.22 | 0.00 | 2.00 | 4.00 | 4.00 | 4.00 | ▁▃▁▃▇ |
| education_budget | 325 | 0.70 | 3.80 | 1.65 | 0.40 | 2.60 | 3.80 | 4.80 | 9.10 | ▃▇▇▂▁ |
| health_budget | 325 | 0.70 | 2.50 | 1.89 | 0.10 | 1.10 | 2.10 | 3.40 | 15.00 | ▇▃▁▁▁ |
| socialsec_budget | 467 | 0.57 | 4.10 | 3.57 | 0.00 | 1.50 | 3.00 | 5.70 | 16.80 | ▇▃▂▁▁ |
| legislative_bal | 2 | 1.00 | 0.00 | 0.27 | -0.82 | -0.15 | 0.00 | 0.11 | 0.92 | ▁▃▇▂▁ |
| repression | 0 | 1.00 | 0.10 | 0.30 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
stat.desc(welfare)## country_id year population gini
## nbr.val NA 1074.000000000 945.0000000 425.0000000
## nbr.null NA 0.000000000 0.0000000 0.0000000
## nbr.na NA 0.000000000 129.0000000 649.0000000
## min NA 1970.000000000 18.9925880 28.9000000
## max NA 2012.000000000 48.7311592 68.3000000
## range NA 42.000000000 29.7385712 39.4000000
## sum NA 2138314.000000000 34356.0089188 21053.9933333
## median NA 1991.000000000 36.5708504 50.6000000
## mean NA 1990.981378026 36.3555650 49.5388078
## SE.mean NA 0.378561808 0.2281924 0.3262479
## CI.mean NA 0.742805393 0.4478231 0.6412646
## var NA 153.913911986 49.2078319 45.2360207
## std.dev NA 12.406204576 7.0148294 6.7257729
## coef.var NA 0.006231201 0.1929506 0.1357678
## sector_dualism gdp foreign_inv ethnic_diversity
## nbr.val 468.0000000 880.000000 876.0000000 704.00000000
## nbr.null 0.0000000 0.000000 0.0000000 384.00000000
## nbr.na 606.0000000 194.000000 198.0000000 370.00000000
## min -36.6597777 1421.710000 -55.2422031 0.00000000
## max 18.5111249 29343.900000 39.8092346 1.00000000
## range 55.1709027 27922.190000 95.0514377 1.00000000
## sum -3365.4469246 6464423.000000 1915.2881060 320.00000000
## median -9.1254131 6264.845000 1.4272784 0.00000000
## mean -7.1911259 7345.935227 2.1864019 0.45454545
## SE.mean 0.5367223 167.934709 0.1410747 0.01877977
## CI.mean 1.0546898 329.599822 0.2768843 0.03687115
## var 134.8171616 24817818.479232 17.4342052 0.24828656
## std.dev 11.6110793 4981.748536 4.1754287 0.49828362
## coef.var -1.6146400 0.678164 1.9097260 1.09622396
## regime_type education_budget health_budget socialsec_budget
## nbr.val 1074.00000000 749.00000000 749.00000000 607.0000000
## nbr.null 22.00000000 0.00000000 0.00000000 13.0000000
## nbr.na 0.00000000 325.00000000 325.00000000 467.0000000
## min 0.00000000 0.40000000 0.10000000 0.0000000
## max 4.00000000 9.10000000 15.00000000 16.8000000
## range 4.00000000 8.70000000 14.90000000 16.8000000
## sum 3259.00000000 2846.30000000 1875.38000000 2486.6000000
## median 4.00000000 3.80000000 2.10000000 3.0000000
## mean 3.03445065 3.80013351 2.50384513 4.0965404
## SE.mean 0.03723480 0.06034266 0.06897698 0.1447749
## CI.mean 0.07306129 0.11846112 0.13541151 0.2843214
## var 1.48902640 2.72728608 3.56361033 12.7225788
## std.dev 1.22025669 1.65144969 1.88775272 3.5668724
## coef.var 0.40213430 0.43457676 0.75394149 0.8707036
## legislative_bal repression
## nbr.val 1072.000000000 1074.000000000
## nbr.null 274.000000000 961.000000000
## nbr.na 2.000000000 0.000000000
## min -0.816666678 0.000000000
## max 0.916149985 1.000000000
## range 1.732816663 1.000000000
## sum -2.198186918 110.000000000
## median 0.000000000 0.000000000
## mean -0.002050547 0.102420857
## SE.mean 0.008248672 0.009185576
## CI.mean 0.016185391 0.018023730
## var 0.072939509 0.090618552
## std.dev 0.270073155 0.301029154
## coef.var -131.707826749 2.939139194- Los resultados están ordenados por el tipo de variable, e indica el número de valores perdidos de cada una de ellas, su respectiva media, su desviación estándar, los valores correspondientes de los percentiles y un pequeño histograma que muestra cómo está distribuida la variable.
Matriz de correlación de variables independientes
Después de identificar todas las variables que incorporaremos al modelo, se recomienda observar cómo se relacionan entre sí. Para ello, crearemos una matriz de correlación de las variables independientes con el comando ggcorrplot() del paquete homónimo, con la que podremos evaluar la correlación de Pearson entre todas las variables.
De todas formas, es importante recordar que la correlación no implica causalidad. Aquí, simplemente queremos entender si las variables del modelo están relacionadas de cierta manera. Este paso no sólo es importante para reconocer nuestros datos y variables, sino también porque queremos evitar tener multicolinealidad perfecta (que haya variables independientes que estén perfectamente correlacionadas) en nuestro modelo, ya que ésta es una suposición central de el MCO.
corr_selected <- welfare %>%
select(gini, education_budget, sector_dualism, foreign_inv, gdp,ethnic_diversity, regime_type, health_budget, socialsec_budget, legislative_bal, population) %>%
# calculate correlation matrix and round to 1 decimal place:
cor(use = "pairwise") %>% round(1)
ggcorrplot(corr_selected, type = "lower", lab = T, show.legend = F)
- Ahora que conocemos todas las variables que se incorporarán al modelo, y cómo se correlacionan entre sí, profundizaremos en las variables clave de interés: las dependientes y las independientes.
Distribución de las variables de interés
- Como hemos mencionado anteriormente, queremos estimar cómo el cambio de una variable independiente (su variación) afecta a la variación de una variable dependiente. Es decir, cómo cambia y cuando cambia x. En este caso, supongamos que nos interesa estimar cómo varían los niveles de desigualdad de un país (medidos como el porcentaje del PIB asignado al presupuesto de educación). Así, nuestra variable de interés independiente es el gasto en educación, mientras que la variable dependiente es la desigualdad.
#Observemos cómo se distribuyen estas variables en nuestra base de datos:
Gini <- ggplot(welfare, aes(x = gini, na.rm = T)) +
geom_histogram(binwidth = 1) +
labs(x = "Gini Index", y = "Frequency",
caption = "Source: Huber et al (2012)")
Education <- ggplot(welfare, aes(x = education_budget, na.rm = T))+
geom_histogram(binwidth = 1) +
labs(x = "Education Expenditure", y = "Frequency", caption = "Source: Huber et al (2012)")
#comparamos graficamente modelo 1, 2 y 3
grid.arrange(Gini, Education, ncol=2) #opcion 1 nombrando las graficas## Warning: Removed 649 rows containing non-finite values (`stat_bin()`).## Warning: Removed 325 rows containing non-finite values (`stat_bin()`).
Relación entre la variable dependiente e independiente
- Después de observar cómo se distribuyen las variables de interés, podemos ver gráficamente cómo se relacionan. Es decir, graficamos la correlación entre estas dos variables: en el eje x (horizontal) ubicamos la variable independiente, mientras que en el eje y (vertical) la variable dependiente. Como resultado, cada “punto” representa una observación de nuestra muestra con un valor particular del gasto en educación (x) y un valor particular en el índice de Gini (y).
ggplot(welfare, aes(education_budget, gini)) +
geom_point() +
labs(x = "Education Expenditure (% of GDP)", y = "Gini",
caption = "Source: Huber and Stephens (2012)")## Warning: Removed 718 rows containing missing values (`geom_point()`).
- Esta es una primera visualización de la relación entre nuestras variables que nos permite observar si hay algún tipo de relación entre ellas. Aquí vemos claramente una relación positiva (cuanto más alto es el gasto en educación, más alto es el Gini). De todos modos, hasta ahora no podemos decir nada concluyente sobre el efecto del gasto en educación en los niveles de desigualdad. Para ello es necesario estimar un modelo. Hasta ahora sólo hemos explorado nuestros datos. ¡Pasemos a las regresiones!
Tercer paso: Estimar Modelo
- El modelo lineal simple asume que una variable de respuesta aleatoria y es una función lineal de una variable independiente x más un término de error u . También decimos que la variable dependiente y es el resultado de un proceso de generación de datos (DGP, por sus siglas en inglés) que puede escribirse como:
\(Y=β0+β1x+u\)
La función lm que forma parte de la base R es la principal herramienta para la estimación de los modelos lineales. A partir de la cual se entiende que un modelo lineal (lm, por sus siglas en inglés) se estima para una variable dependiente Y regresada (~) en una variable independiente X. El “1” no suele incluirse, pero lo añadimos para denotar la intercepción (β0).
Basándonos en la investigación de Huber et al. (2006), nuestro modelo plantea que la desigualdad es una función lineal del gasto en educación, además de un término de error no observado u y una constante β0. Formalment, el supuesto de observaciones independientes y distribuidas de forma idéntica es lo que nos permite escribir el modelo para un individuo de i escogido al azar como:
\(Desigualdadi=β0+β1GastoEucacionali+ui\)
- Dado que los datos se almacenan en un data.frame, necesitamos indicar su nombre en la función para traer los datos de la base de datos correspondiente. Eso es lo que ocurre después de la coma en el siguiente comando
# after the comma we indicate the data.frame that contains the data
model_1 <- lm(gini ~ 1 + education_budget, data = welfare)
class(model_1) # we verify that the object class is "lm" ## [1] "lm"En la primera línea de código creamos un objeto (model_1) que guarda los resultados de la función lm. Esta función crea objetos de clase lm, que son vectores que incluyen los coeficientes estimados, los errores estándar, los residuos, los valores predichos, entre otros resultados de la estimación.
Para ver los componentes del objeto, una forma rápida de hacerlo es usar la función summary.
summary(model_1) ##
## Call:
## lm(formula = gini ~ 1 + education_budget, data = welfare)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.036 -5.616 1.086 4.990 15.571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 44.8052 1.0230 43.798 < 0.0000000000000002 ***
## education_budget 1.2326 0.2502 4.927 0.00000129 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.703 on 354 degrees of freedom
## (718 observations deleted due to missingness)
## Multiple R-squared: 0.06417, Adjusted R-squared: 0.06152
## F-statistic: 24.27 on 1 and 354 DF, p-value: 0.000001286- Se pueden lograr presentaciones más elegantes con la función screenreg del paquete texreg. Esto se utilizará en todos los capítulos de los modelos econométricos. Veamos la presentación de los resultados con la función screenreg
screenreg(model_1) ##
## ============================
## Model 1
## ----------------------------
## (Intercept) 44.81 ***
## (1.02)
## education_budget 1.23 ***
## (0.25)
## ----------------------------
## R^2 0.06
## Adj. R^2 0.06
## Num. obs. 356
## ============================
## *** p < 0.001; ** p < 0.01; * p < 0.05#De esta forma, podemos añadirle nombres a las variables:
screenreg(model_1,
custom.model.names = "Model 1", custom.coef.names = c("Constant", "Education expenditure"))##
## =================================
## Model 1
## ---------------------------------
## Constant 44.81 ***
## (1.02)
## Education expenditure 1.23 ***
## (0.25)
## ---------------------------------
## R^2 0.06
## Adj. R^2 0.06
## Num. obs. 356
## =================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Como puedes notar, los resultados se muestran mejor con texreg que con summary(). Se ve claramente que gasto_educ –el gasto en educación– tiene un efecto positivo de magnitud 1.233. En particular, cuando el gasto en educación como porcentaje del PIB aumenta en una unidad, la desigualdad aumenta en 1,23%. El efecto del gasto en educación es significativo con un nivel de confianza del 99,9%. Lo sabemos porque además del coeficiente aparecen tres estrellas que se refieren a un nivel de significación del 0,01%.
Una aproximación manual de la distancia beta estimada en la distribución del estimador bajo la hipótesis nula β1=0 se obtiene cuando dividimos la estimación por su error estándar:
1.233 / 0.25## [1] 4.932El mismo valor es entregado por la tercera columna de la sección “Coeficientes” en el summary() del modelo_1. El valor t siempre se interpreta como la distancia de la estimación ^β1 de la media de la distribución del estimador bajo H0=β1=0. En este caso, el valor 1,233 está a 4,93 desviaciones estándar de la distribución del estimador cuando H0 es verdadero (la media de la distribución es 0).
En este caso, la estimación de ^β1 es más de 4,93 desviaciones estándar de la media de la distribución bajo H0=β1=0, por lo que es poco probable que se haya observado un efecto de 1,23 si H0fuera cierto.
La probabilidad precisa se puede observar en la cuarta columna del modelo summary, que se puede solicitar a R con el siguiente comando.
coef(summary(model_1))[, "Pr(>|t|)"] ## (Intercept)
## 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000005548453
## education_budget
## 0.0000012863896523289702605305018581760556628523772815242409706115722656250000000000000000000000000000000000000000000000000000000000000000000000000000000- La probabilidad de observar una estimación de 1,23 si el H0 fuera cierto es de 0,00000128. Por lo tanto, podemos rechazar H0 incluso con un nivel de confianza del 99,9%.
Cuarto paso: Representación gráfica
- Como vimos anteriormente, una de las formas más fáciles de presentar la relación entre dos variables es a través de gráficos. Ya se ha dicho que ggplot2 es una herramienta conveniente para generar varios tipos de gráficos. Exploremos su uso para MCO. En el primer código, todas las observaciones se grafican según sus valores de variables independientes y dependientes.
ggplot(data = welfare, # the dataset is selected
aes(x = education_budget, y = gini))+ # indep. and dep. variables
geom_point() + # the observed values are plotted
geom_smooth(method = "lm", # the regression line is overlapped
se = F, # the error area is not plotted at a 95% CI
color = "black") + # color line
labs (x = "Education Expenditure", y = "Inequality")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 718 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 718 rows containing missing values (`geom_point()`).
- Por lo general, también es útil mostrar una representación gráfica del error de predicción de la línea. ggplot2 nos permite editar un área sombreada donde los valores predichos con un cierto nivel de significación podrían haber sido localizados. Aunque el 95% de confianza es el valor por defecto, también podemos editar ese valor. El primer bloque muestra la línea de regresión y su error para un nivel de significancia estadística del 95%. Ten en cuenta que esta línea no representará el coeficiente que obtuviste a través de lm después de incluir controles a tu regresión.
ggplot(data = welfare, aes(x = education_budget, y = gini)) +
geom_point() +
geom_smooth(method = "lm", color = "black",
se = T) + # we add the prediction
labs(x = "Education Expenditure", y = "Inequality") ## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 718 rows containing non-finite values (`stat_smooth()`).## Warning: Removed 718 rows containing missing values (`geom_point()`).