Relación entre variables
Paula Videla
Abril 2021
library(tidyverse)
options(scipen = 99)Cargamos la base de datos. Los datos se obtuvieron de Banco Mundial y se agruparon diversas bases.
El dataset puede descargarse desde aqui.
También puede cargarse desde una query a data.world:
bancomundial <- read.csv(“https://query.data.world/s/owzasudqx4g5a3lqt32ykoudalq3ul”, header=TRUE, stringsAsFactors=TRUE)
bancomundial <- read.csv("data-wb/dataset_bancomundial.csv",
stringsAsFactors = TRUE) En primer lugar vamos a ver que contienen nuestros datos.
head(bancomundial)## Country.Name Country.Code Region IncomeGroup year
## 1 American Samoa ASM East Asia & Pacific Upper middle income 1960
## 2 American Samoa ASM East Asia & Pacific Upper middle income 1961
## 3 American Samoa ASM East Asia & Pacific Upper middle income 1962
## 4 American Samoa ASM East Asia & Pacific Upper middle income 1963
## 5 American Samoa ASM East Asia & Pacific Upper middle income 1964
## 6 American Samoa ASM East Asia & Pacific Upper middle income 1965
## pop_total perc_urban_pop pop_in_largest_city GDP_in_USD
## 1 20123 66.211 NA NA
## 2 20602 66.641 NA NA
## 3 21253 67.068 NA NA
## 4 22034 67.493 NA NA
## 5 22854 67.916 NA NA
## 6 23672 68.334 NA NAsummary(bancomundial)## Country.Name Country.Code Region
## Afghanistan : 61 ABW : 61 East Asia & Pacific :2257
## Albania : 61 AFG : 61 Europe & Central Asia :3538
## Algeria : 61 AGO : 61 Latin America & Caribbean :2562
## American Samoa: 61 ALB : 61 Middle East & North Africa:1281
## Andorra : 61 AND : 61 North America : 183
## Angola : 61 ARE : 61 South Asia : 488
## (Other) :12871 (Other):12871 Sub-Saharan Africa :2928
## IncomeGroup year pop_total
## High income :5002 Min. :1960 Min. : 3893
## Low income :1769 1st Qu.:1975 1st Qu.: 468584
## Lower middle income:3050 Median :1990 Median : 4133322
## Upper middle income:3416 Mean :1990 Mean : 24335531
## 3rd Qu.:2005 3rd Qu.: 13226839
## Max. :2020 Max. :1397715000
## NA's :326
## perc_urban_pop pop_in_largest_city GDP_in_USD
## Min. : 2.077 Min. : 5254 Min. : 8824448
## 1st Qu.: 30.185 1st Qu.: 519169 1st Qu.: 1409873950
## Median : 50.258 Median : 1136880 Median : 7480968858
## Mean : 51.131 Mean : 2618042 Mean : 173953217878
## 3rd Qu.: 71.755 3rd Qu.: 2658965 3rd Qu.: 49745088112
## Max. :100.000 Max. :37468302 Max. :21433226000000
## NA's :345 NA's :4125 NA's :3352¿Qué vemos?
Hay 61 observaciones por país, cada una de ellas corresponde a un año. Los años van desde 1960 a 2020. El país con menos habitantes de la muestra tuvo 3.893hab, y el país más habitado contó con 1.397.715.000hab. El porcentaje de población urbana promedio es 50.258%, mientras existe algún país con 100% de población urbana.
Veamos cuales son las regiones:
table(bancomundial$Region)##
## East Asia & Pacific Europe & Central Asia
## 2257 3538
## Latin America & Caribbean Middle East & North Africa
## 2562 1281
## North America South Asia
## 183 488
## Sub-Saharan Africa
## 2928Ahora, vamos a crear una variable que sea PBI (GDP) per capita (anual)
datos <- bancomundial %>%
mutate(GDP_per_capita = GDP_in_USD / pop_total)
summary(datos$GDP_per_capita)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 34.79 559.00 1954.24 8583.15 8050.54 189422.22 3355El PBI per cápita mínimo es apróximadamente 245 veces más chico que el promedio.
Relación de variables
¿Cómo se relaciona la población con el porcentaje de población urbana?
Para esto, vamos a usar un gráfico de scatterplot (o gráfico de puntos).
datos %>%
ggplot() +
geom_point(aes(y = perc_urban_pop, x = pop_total, color = Region), alpha = 0.6) +
labs(title = "Correlación entre el porcentaje de población urbana y la población total",
x= "Población total (en millones)",
y = "Población urbana (%)") +
scale_x_continuous(labels = c(250, 500,750,1000, 1250),breaks = 10^6 *c(250, 500,750,1000, 1250)) +
theme_light()## Warning: Removed 446 rows containing missing values (geom_point).
Vemos que hay dos países (probablemente sean India y China) que, por su población, distorcionan el gráfico. La gran mayoría de los países tienen menos de 250 millones de habitantes.
Vamos a analizar las regiones como conjunto. Para ello, vamos a generar un nuevo dataset agrupando las regiones.
Previamente, necesitamos agregar nuevos campos para poder obtener los resultados por región.
datos <- datos %>%
mutate(urban_pop = round(perc_urban_pop*pop_total/100, digits = 0))Ahora si, vamos a obtener los datos por región:
datos.region <- datos%>%
group_by(year, Region) %>%
summarise(pop_total = sum(pop_total, na.rm = TRUE), #sumamos población total de la región
urban_pop = sum(urban_pop, na.rm = TRUE), #sumamos población urbana de la región
perc_urban_pop = round(urban_pop/pop_total*100, digits = 2), #obtenemos el % de población urbana
GDP_in_USD = sum(GDP_in_USD, na.rm = TRUE), #sumamos el PBI total de la región
GDP_per_capita = round(GDP_in_USD/pop_total)) #obtenemos el PBI per capita de la región## `summarise()` regrouping output by 'year' (override with `.groups` argument)head(datos.region)## # A tibble: 6 x 7
## # Groups: year [1]
## year Region pop_total urban_pop perc_urban_pop GDP_in_USD GDP_per_capita
## <int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1960 East Asia … 1.03e9 231615684 22.5 1.46e11 142
## 2 1960 Europe & C… 6.60e8 368034505 55.8 2.79e11 422
## 3 1960 Latin Amer… 2.20e8 108735965 49.5 5.74e10 261
## 4 1960 Middle Eas… 1.05e8 36736033 34.9 1.14e10 108
## 5 1960 North Amer… 1.99e8 138875014 69.9 5.84e11 2939
## 6 1960 South Asia 5.73e8 95913995 16.7 4.75e10 83Veamos el scatterplot.
datos.region %>%
ggplot() +
geom_point(aes(y = perc_urban_pop, x = pop_total, color = Region), alpha = 0.6) +
labs(title = "Correlación entre el porcentaje de población urbana y la población total",
subtitle = "Por regiones",
color = "Región",
x= "Población total (en millones)",
y = "Población urbana (%)") +
scale_x_continuous(labels = seq(0,2500,200),breaks = 10^6 *seq(0,2500,200)) +
theme_light()## Warning: Removed 7 rows containing missing values (geom_point).
Vemos que cada región tiene un comportamiento distinto. A continuación queremos evaluar la relación entre la población total y el porcentaje de población urbana en América Latina.
¿Cómo se relaciona la población total de América Latina (variable explicativa) con el porcentaje de población urbana en el continente según cada año?
Vamos a filtrar los datos de datos.region en América Latina.
datos.latam <- datos.region %>%
filter(Region == "Latin America & Caribbean")
head(datos.latam)## # A tibble: 6 x 7
## # Groups: year [6]
## year Region pop_total urban_pop perc_urban_pop GDP_in_USD GDP_per_capita
## <int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1960 Latin Amer… 219825927 108735965 49.5 5.74e10 261
## 2 1961 Latin Amer… 225909659 113558251 50.3 6.08e10 269
## 3 1962 Latin Amer… 232179886 118597846 51.1 9.40e10 405
## 4 1963 Latin Amer… 238608947 123826704 51.9 9.48e10 397
## 5 1964 Latin Amer… 245143683 129227573 52.7 1.06e11 431
## 6 1965 Latin Amer… 251756122 134636970 53.5 1.14e11 453Vamos a realizar nuevamente un gráfico de puntos o scatterplot. En el eje x, vamos a colocar nuestra variable explicativa (población total) y en el eje y colocaremos la variable que queremos explicar (porcentaje de población urbana)
datos.latam %>%
ggplot() +
geom_point(aes(y = perc_urban_pop, x = pop_total)) +
labs(title = "Correlación entre el porcentaje de población urbana y la población total",
subtitle = "En América Latina",
x= "Población total (en millones)",
y = "Población urbana (%)") +
scale_x_continuous(labels = c(200, 400, 600),breaks = 10^6 *c(200, 400, 600) )+
theme_light()## Warning: Removed 1 rows containing missing values (geom_point).
Vemos que existe una relación positiva entre la población total en América Latina y la tasa de población urbana en la región.
Modelo estádistico
Regresión lineal simple
Con el objetivo de analizar la relación entre nuestras variables, vamos a armar un modelo estádistico a partir de una regresión lineal. En R, vamos a usar la función lm() (por “linear model”).
lm1 <- lm(data = datos.latam, perc_urban_pop ~ pop_total )Con el símbolo ~ estamos pidiendo a la función que evalúe perc_urban_pop, nuestra variable dependiente o a explicar, a partir (o versus) pop_total, nuestra variable explicativa o independiente.
Veamos el efecto que tiene la población total sobre el porcentaje de población urbana en América Latina.
lm1##
## Call:
## lm(formula = perc_urban_pop ~ pop_total, data = datos.latam)
##
## Coefficients:
## (Intercept) pop_total
## 37.28021838891 0.00000007151Recordemos como se expresaba la formula de una regresión lineal simple.
\(y_{i} = \beta_{0} + \beta_{1} X_{1} + \epsilon\)
Entonces, de acuerdo a nuestros coeficientes, tendríamos:
\(y_{i} = 37.28021838891 + 0.00000007151 X_{pop.total} + \epsilon\)
¿Qué significa esto?
Que por cada incremento de 1 en la población total, hay un incremento de 0.00000007151 en el porcentaje de población urbana.
Dicho de otra manera, por cada incremento de 1.000.000 de personas en la población total, hay un incremento de 0.07151 en el porcentaje de población urbana.
Vemos como se grafica nuestro modelo:
Recordemos que la regresión lineal puede ser modelada en este caso con una línea
datos.latam %>%
ggplot() +
geom_abline(aes(intercept= 37.28021838891 , slope = 0.00000007151), color = "red", linetype = "dashed") +
geom_point(aes(y = perc_urban_pop, x = pop_total)) +
labs(title = "Correlación entre el porcentaje de población urbana y el PBI per capita",
x= "Población total (en millones)",
y = "Población urbana (%)") +
scale_x_continuous(labels = c(200, 400, 600),breaks = 10^6 *c(200, 400, 600) )+
theme_light()## Warning: Removed 1 rows containing missing values (geom_point).
Podemos observar que si bien nuestro modelo está relacionado con los datos, no coincide perfectamente. ¿Qué tan bien se desempeña nuestro modelo?
Para ello vamos a usar la función summary().
summary(lm1)##
## Call:
## lm(formula = perc_urban_pop ~ pop_total, data = datos.latam)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5398 -1.2617 0.3349 1.4382 1.8408
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.28021838891 0.71371004550 52.23 <0.0000000000000002 ***
## pop_total 0.00000007151 0.00000000157 45.55 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.583 on 58 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9728, Adjusted R-squared: 0.9723
## F-statistic: 2075 on 1 and 58 DF, p-value: < 0.00000000000000022¿Qué podemos deducir?
El r-squared (o R2) tiene un valor de 0.9728. Esta métrica varía de 0 a 1, y permite entender que tan bien nuestro modelo ajusta a las observaciones. En este caso, tenemos un R2 muy alto, con lo cual es mucho lo que nuestro modelo explica. Además el pvalor para la variable pop_total es muy bajo, es decir que la variable es significativamente estadística.
Como conclusión, hay una fuerta correlación entre las variables pop_total y perc_urban_pop (según el p-valor), y además nuestro modelo ajusta bien a los datos (según R2).
Regresión múltiple
Vimos que si bien nuestro modelo es bastante explicativo, también creemos que el PBI per capita de América Latina influye en el porcentaje de población urbana.
Queremos incorporar la variable GDP_per_capita a nuestro modelo de regresión lineal. ¿Cómo se hace? ¡Es muy fácil! Simplemente tenemos que agregar la variable al modelo.
lm2 <- lm(data = datos.latam, perc_urban_pop ~ pop_total + GDP_per_capita)Con el símbolo + agregamos una variable explicativa más, en este caso GDP_per_capita
Vamos a ver los resultados de nuestro modelo.
lm2##
## Call:
## lm(formula = perc_urban_pop ~ pop_total + GDP_per_capita, data = datos.latam)
##
## Coefficients:
## (Intercept) pop_total GDP_per_capita
## 32.21658144028 0.00000009047 -0.00087414338La fórmula para una regresión lineal múltiple es la siguiente:
\(y_{i} = \beta_{0} + \beta_{1} X_{1} + \beta_{2} X_{2} + \epsilon\)
Entonces, de acuerdo a nuestros coeficientes, tendríamos:
\(y_{i} = 32.21658144028 + 0.00000009047 X_{pop.total} -0.00087414338 X_{GDP.per.capita } + \epsilon\)
¿Qué significa esto?
La población sigue teniendo un efecto positivo: por cada incremento de 1.000.000 de personas en América Latina, el porcentaje de población urbana aumenta de 0.09047.
El PBI per capita parecería tener un efecto negativo: por cada incremento de 1000 USD en el PBI per capita, el porcentaje de población urbana disminuye de 0.8741.
Veamos que tan bien se desempeñó nuestro modelo.
summary(lm2)##
## Call:
## lm(formula = perc_urban_pop ~ pop_total + GDP_per_capita, data = datos.latam)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8091 -0.9458 0.5530 0.8425 1.8199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.216581440279 0.932952868194 34.532 < 0.0000000000000002 ***
## pop_total 0.000000090466 0.000000003086 29.316 < 0.0000000000000002 ***
## GDP_per_capita -0.000874143379 0.000131331421 -6.656 0.0000000119 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.198 on 57 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9847, Adjusted R-squared: 0.9842
## F-statistic: 1834 on 2 and 57 DF, p-value: < 0.00000000000000022Al incorporar una segunda variable observamos:
- Las dos variables: pop_total y GDP_per_capita son significativas.
- El R2 es aún mayor, nuestro modelo se ajusta mejora a la realidad.
Bonus track: regresión lineal con variable categórica
En el segundo gráfico pudimos observar que la relación entre la población total y el porcentaje de población urbana estaba fuertemente diferenciado por región.
La regresión lineal permite utilizar como variable explicativa una variable categórica. ¿Cómo funciona? Se crea una variable binaria (dummy) que indica si pertenece o no a esa categoría, y cual es el efecto de pertenecer a esa categoría.
Como vimos antes, las regiones están compuestas por 7 categorías. Entonces nuestro modelo de regresión quedaría así:
- EAP = East Asia & Pacific
- ECA = Europe & Central Asia
- LAC = Latin America & Caribbean
- MENA = Middle East & North Africa
- NA = North America
- SA = South Asia
- SSA = Sub-Saharan Africa
\(y_{i} = \beta_{0} + \beta_{EAP} X_{EAP} + \beta_{ECA} X_{ECA} + \beta_{LAC} X_{LAC} + \beta_{MENA} X_{MENA} + \beta_{NA} X_{NA} + \beta_{SA} X_{SA}+ \beta_{SSA} X_{SSA} +\epsilon\)
Donde los \(\beta_{i}\) toman el valor del efecto y los \(X_{i}\) adoptan valor 1 o 0 si se encuentran en esa región.
Ahora corramos nuestro modelo para evaluar la relación entre porcentaje de población urbana y la región donde se encuentran.
lm3 <- lm(data = datos.region, perc_urban_pop ~ Region)
lm3##
## Call:
## lm(formula = perc_urban_pop ~ Region, data = datos.region)
##
## Coefficients:
## (Intercept) RegionEurope & Central Asia
## 36.41 29.65
## RegionLatin America & Caribbean RegionMiddle East & North Africa
## 32.02 16.54
## RegionNorth America RegionSouth Asia
## 40.05 -11.59
## RegionSub-Saharan Africa
## -9.46¿Qué significa?
Lo primero que observamos es que no están nuestras 7 regiones, falta East Asia & Pacific. Esto ocurre porque el modelo toma la primer región (alfabéticamente) como línea de base. Es decir, si la categoría Region es East Asia & Pacific entonces el porcentaje de población urbana sería 36.41%. En caso de corresponder la categoría South Asia tendrías un porcentaje de población urbana de 26.95 (36.41-9.46).
Volviendo a nuestro modelo, la regresión tendría la siguiente ecuación:
\(y_{i} = 36.41 + 29.65 X_{ECA} + 32.02 X_{LAC} + 16.54X_{MENA} + 40.05X_{NA} -11.592X_{SA} -9.46X_{SSA} +\epsilon\)
Las variables X (dummies) solo se activan 1 vez (en la categoría de corresponde) o ninguna (para la categoría East Asia & Pacific). Para las demás \(X_{i}\), el valor que adopta es 0.
Nota: Con la función relevel(factor, ref) podríamos modificar el nivel de referencia y remplazar East Asia & Pacific por otra categoría.
Veamos ahora los resultados del modelo.
summary(lm3)##
## Call:
## lm(formula = perc_urban_pop ~ Region, data = datos.region)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.9688 -5.4098 0.4153 5.1152 23.2945
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 36.406 1.007 36.165
## RegionEurope & Central Asia 29.649 1.424 20.826
## RegionLatin America & Caribbean 32.023 1.424 22.494
## RegionMiddle East & North Africa 16.543 1.424 11.620
## RegionNorth America 40.051 1.424 28.133
## RegionSouth Asia -11.595 1.424 -8.144
## RegionSub-Saharan Africa -9.460 1.424 -6.645
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## RegionEurope & Central Asia < 0.0000000000000002 ***
## RegionLatin America & Caribbean < 0.0000000000000002 ***
## RegionMiddle East & North Africa < 0.0000000000000002 ***
## RegionNorth America < 0.0000000000000002 ***
## RegionSouth Asia 0.00000000000000456 ***
## RegionSub-Saharan Africa 0.00000000009582287 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.798 on 413 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.8639, Adjusted R-squared: 0.8619
## F-statistic: 436.9 on 6 and 413 DF, p-value: < 0.00000000000000022Con la variable categórica Region observamos que:
- Todas las regiones son estadísticamente significativas.
- El R2 es de 0.8639, nuestro modelo ajusta muy bien, pero los modelos anteriores ajustaban mejor.
Sería interesante realizar un modelo de regresión lineal múltiple incorporando otras variables de nuestro dataset, como el año o la población de la ciudad más grande o el GroupIncome (categoría del Banco Mundial) al que pertenecen. ¡Prueben hacerlo!