Aplicación 2.6: Análisis econométrico mixto (datos espaciales, forma funcional, inferencia y predicción)
Medición de la calidad de vida
La renta per cápita es la forma más usual de cuantificar el nivel de bienestar de un país, permitiendo además la comparación entre distintas economías, independientemente de su tamaño. Sin embargo, utilizar sólo el PIB de los habitantes de un territorio como indicador del bienestar social y del progreso del país resulta bastante inapropiado, como se reconoce en el Informe Stiglitz-Sen-Fitoussi del año 2009.
Distintos organismos, como la OCDE, las Naciones Unidas o la Comisión Europea, han reconocido la debilidad del PIB o la renta per cápita para valorar el nivel de desarrollo humano de una sociedad, proponiendo distintas alternativas para determinar la calidad de vida de las personas que residen en un determinado territorio.
Para este ejercicio utilizaremos los datos del año 2021 elaborados por las Naciones Unidas para construir el Índice de Desarrollo Humano (https://hdr.undp.org/data-center/human-development-index#/indicies/HDI), que es un indicador de bienestar agregado de los países (con escala entre 0 y 100) que engloba tanto los niveles de salud [medida por la variable LifeExp, la esperanza de vida al nacer, medida en años] y educación [medida por la variable ExpEduc, el grado de escolarización esperado para cada persona, medido en años] como el nivel económico de cada país [medido por la variable GNIpc, la renta per cápita, medida en dólares constantes en términos de paridad de poder adquisitivo].
Code
# Lectura de librerías
library(tidyverse)
library(readxl)
library(sf)
library(geojsonsf)
library(viridis)
library(car)
library(alr4)
# Lectura de los datos
datos <- read_excel("data/HDI.xlsx")
summary(datos[,2:5]) HDI LifeExp ExpEduc GNIpc
Min. :38.50 Min. :52.53 Min. : 5.543 Min. : 731.8
1st Qu.:59.95 1st Qu.:65.75 1st Qu.:11.601 1st Qu.: 4592.9
Median :73.90 Median :71.69 Median :13.405 Median : 12306.3
Mean :72.06 Mean :71.31 Mean :13.535 Mean : 20249.1
3rd Qu.:83.50 3rd Qu.:76.70 3rd Qu.:15.624 3rd Qu.: 30079.8
Max. :96.20 Max. :85.47 Max. :21.055 Max. :146829.7 Code
# Lectura del mapa
mapa <- geojson_sf("data/WORLDmap.geojson")
# Fusión de datos y mapa
datos_geo <- left_join(mapa, datos, by = "name")
# Gráfica del HDI
ggplot(datos_geo) +
geom_sf(aes(fill=HDI)) +
theme_bw() +
labs(title = "Índice de Desarrollo Humano 2021 (Escala 0-100)") +
scale_fill_viridis(option="magma")
Code
# Regresión para explicar la variable HDI:
# HDI = f(LifeExp, ExpEduc, GNIpc) + e
# Gráficas parciales
# Variables originales
scatterplotMatrix(~ HDI + LifeExp + ExpEduc + GNIpc, data=datos)
Code
# Variables transformadas
scatterplot(log(HDI) ~ LifeExp, data=datos,
smooth=list(smoother=loessLine, var=FALSE, lwd.smooth=3),
col="blue", regLine=list(lwd=3))
Code
scatterplot(log(HDI) ~ ExpEduc, data=datos,
smooth=list(smoother=loessLine, var=FALSE, lwd.smooth=3),
col="blue", regLine=list(lwd=3))
Code
scatterplot(log(HDI) ~ log(GNIpc), data=datos,
smooth=list(smoother=loessLine, var=FALSE, lwd.smooth=3),
col="blue", regLine=list(lwd=3))
Code
# Regresión por MCO del modelo
# log(HDI) = beta1 + beta2*LifeExp + beta3*ExpEduc + beta4*log(GNIpc) + e
# Interpretación beta2 y beta3: semi-elasticidades
# Interpretación beta4: elasticidad
mod_log <- lm(log(HDI) ~ LifeExp + ExpEduc + log(GNIpc), data = datos)
summary(mod_log)
Call:
lm(formula = log(HDI) ~ LifeExp + ExpEduc + log(GNIpc), data = datos)
Residuals:
Min 1Q Median 3Q Max
-0.133214 -0.023027 0.004089 0.032086 0.105939
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.5298416 0.0352821 71.703 < 2e-16 ***
LifeExp 0.0064339 0.0008519 7.552 1.84e-12 ***
ExpEduc 0.0206765 0.0020530 10.072 < 2e-16 ***
log(GNIpc) 0.1056086 0.0059811 17.657 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.046 on 187 degrees of freedom
Multiple R-squared: 0.9587, Adjusted R-squared: 0.958
F-statistic: 1445 on 3 and 187 DF, p-value: < 2.2e-16Code
# Intervalo de confianza del 95% para los parámetros del modelo
confint(mod_log, level=.95) 2.5 % 97.5 %
(Intercept) 2.460239566 2.59944364
LifeExp 0.004753341 0.00811453
ExpEduc 0.016626586 0.02472649
log(GNIpc) 0.093809495 0.11740775Code
# Gráficos de efectos
plot(allEffects(mod_log), grid=TRUE, rug=TRUE)
Code
# Inferencia estadística: contraste de la hipótesis {beta2<=beta3}
# Test t específico: H_0:{beta2-beta3<=0}
# t_0 = (beta2-beta3-0)/se(beta2-beta3)
b2 <- coef(mod_log)[[2]] # estimación del parámetro beta2
b3 <- coef(mod_log)[[3]] # estimación del parámetro beta3
varb2 <- vcov(mod_log)[2,2] # estim. var(beta2)
varb3 <- vcov(mod_log)[3,3] # estim. var(beta3)
covb2b3 <- vcov(mod_log)[2,3] # estim. cov(beta2,beta3)
se_beta2_menos_beta3 <- sqrt(varb2 + varb3 - 2 * covb2b3)
t_0 <- (b2-b3)/se_beta2_menos_beta3
cat("Estadístico t0 = ", round(t_0,3))Estadístico t0 = -5.814Code
df <- df.residual(mod_log)
pval = 1-pt(t_0,df)
cat("p-valor = ", round(pval,3))p-valor = 1Code
# Test t general de la hipótesis {R*beta (>)(=)(<) r}
# Definición del contraste
t_test<-function(model,R,r){
VCOV<-vcov(model)
Rbeta<-as.numeric(t(R) %*% coef(model))
varRbeta<-as.numeric(t(R) %*% VCOV %*% R)
t_0<-(Rbeta-r)/sqrt(varRbeta)
df<-df.residual(model)
pval_unilat_izqda = pt(t_0,df)
pval_bilat = 2*(1-pt(abs(t_0),df))
pval_unilat_dcha = 1-pt(t_0,df)
list<-c(t_0, pval_unilat_izqda, pval_bilat, pval_unilat_dcha)
return(list)
}
# Ejecución del contraste
R <- c(0,1,-1,0)
r <- 0
p <- t_test(mod_log,R,r)
cat(" Estadístico t0 = ",p[[1]],"\n",
" p-valor unilat_izqda (Ho: Rbeta >= r) = ",p[[2]],"\n",
" p-valor bilateral (Ho: Rbeta = r) = ",p[[3]],"\n",
" p-valor unilat_dcha (Ho: Rbeta <= r) = ",p[[4]]," ") Estadístico t0 = -5.814319
p-valor unilat_izqda (Ho: Rbeta >= r) = 1.294589e-08
p-valor bilateral (Ho: Rbeta = r) = 2.589178e-08
p-valor unilat_dcha (Ho: Rbeta <= r) = 1 Code
# Predicción del HDI para un país con los siguientes indicadores parciales:
# LifeExp = 71.7, ExpEduc = 13.4, GNIpc = 12306.3
nuevo_PAIS <- data.frame(LifeExp=c(71.7), ExpEduc=c(13.4), GNIpc=c(12306.3))
nuevo_PAIS LifeExp ExpEduc GNIpc
1 71.7 13.4 12306.3Code
# HDI esperado con intervalo de confianza
pred_log_HDI_IC <- predict(mod_log, nuevo_PAIS,
interval="confidence", level=0.95)
pred_log_HDI_IC fit lwr upr
1 4.262828 4.256127 4.26953Code
pred_HDI_IC=exp(pred_log_HDI_IC)
pred_HDI_IC fit lwr upr
1 71.01054 70.53624 71.48802Code
# Lectura de librerías
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy as sp
import scipy.stats as scs
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Lectura de los datos
datos = pd.read_excel('data/HDI.xlsx')
datos.describe() HDI LifeExp ExpEduc GNIpc
count 191.000000 191.00000 191.000000 191.000000
mean 72.057592 71.31286 13.534658 20249.088223
std 15.066063 7.64596 2.923911 21825.277076
min 38.500000 52.52540 5.542510 731.786709
25% 59.950000 65.74720 11.601258 4592.919612
50% 73.900000 71.69400 13.404920 12306.341000
75% 83.500000 76.69930 15.623665 30079.789725
max 96.200000 85.47340 21.054590 146829.700600Code
# Lectura del mapa
mapa = gpd.read_file("data/WORLDmap.geojson")
# Fusión de datos y mapa
datos_geo = mapa.merge(datos, on = 'name', how = 'left')
# Gráfica del HDI
fig, ax = plt.subplots(1,1)
datos_geo.plot(column="HDI",
legend=True,
cmap='magma',
ax=ax)
plt.title('Índice de Desarrollo Humano 2021 (Escala 0-100)')
plt.tight_layout()
plt.show()
Code
# Regresión para explicar la variable HDI:
# HDI = f(LifeExp, ExpEduc, GNIpc) + e
# Gráficas parciales
# Variables originales
plt.figure(2)
sns.pairplot(data = datos, diag_kind = 'kde');
plt.show()
Code
# Variables tranformadas
plt.figure(3)
sns.regplot(data = datos, x="LifeExp", y="HDI");
plt.yscale('log')
plt.show()
Code
plt.figure(4)
sns.regplot(data = datos, x="ExpEduc", y="HDI");
plt.yscale('log')
plt.show()
Code
plt.figure(5)
sns.regplot(data = datos, x="GNIpc", y="HDI", logx=True);
plt.xscale('log')
plt.yscale('log')
plt.show()
Code
# Regresión por MCO del modelo
# log(HDI) = beta1 + beta2*LifeExp + beta3*ExpEduc + beta4*log(GNIpc) + e
# Interpretación beta2 y beta3: semi-elasticidades
# Interpretación beta4: elasticidad
mod_log = smf.ols(formula = "np.log(HDI) ~ LifeExp + ExpEduc + np.log(GNIpc)",
data = datos).fit()
print(mod_log.summary()) OLS Regression Results
==============================================================================
Dep. Variable: np.log(HDI) R-squared: 0.959
Model: OLS Adj. R-squared: 0.958
Method: Least Squares F-statistic: 1445.
Date: Sun, 09 Feb 2025 Prob (F-statistic): 4.60e-129
Time: 13:38:13 Log-Likelihood: 319.11
No. Observations: 191 AIC: -630.2
Df Residuals: 187 BIC: -617.2
Df Model: 3
Covariance Type: nonrobust
=================================================================================
coef std err t P>|t| [0.025 0.975]
---------------------------------------------------------------------------------
Intercept 2.5298 0.035 71.703 0.000 2.460 2.599
LifeExp 0.0064 0.001 7.552 0.000 0.005 0.008
ExpEduc 0.0207 0.002 10.072 0.000 0.017 0.025
np.log(GNIpc) 0.1056 0.006 17.657 0.000 0.094 0.117
==============================================================================
Omnibus: 3.315 Durbin-Watson: 2.023
Prob(Omnibus): 0.191 Jarque-Bera (JB): 2.977
Skew: -0.298 Prob(JB): 0.226
Kurtosis: 3.138 Cond. No. 781.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.Code
# Intervalo de confianza del 95% para los parámetros del modelo
mod_log.conf_int() 0 1
Intercept 2.460240 2.599444
LifeExp 0.004753 0.008115
ExpEduc 0.016627 0.024726
np.log(GNIpc) 0.093809 0.117408Code
# Inferencia estadística: contraste de la hipótesis {beta2<=beta3}
# Test t específico: H_0:{beta2-beta3<=0}
# t_0 = (beta2-beta3-0)/se(beta2-beta3)
b2 = mod_log.params['LifeExp'] # estimación del parámetro beta2
b3 = mod_log.params['ExpEduc'] # estimación del parámetro beta3
varb2 = (mod_log.bse['LifeExp'])**2 # estim. var(beta2)
varb3 = (mod_log.bse['ExpEduc'])**2 # estim. var(beta3)
covb2b3 = mod_log.cov_params().loc['LifeExp', 'ExpEduc'] # estim. cov(beta2,beta3)
se_beta2_menos_beta3 = np.sqrt(varb2 + varb3 - 2 * covb2b3)
t_0 = (b2-b3)/se_beta2_menos_beta3
print("Estadístico t0 = ", t_0.round(3))Estadístico t0 = -5.814Code
df = mod_log.df_resid
pval = 1-sp.stats.t.cdf(t_0, df)
print("p-valor = ", pval.round(3))p-valor = 1.0Code
# Test t general de la hipótesis {R*beta (>)(=)(<) r}
# Definición del contraste
def t_test(model, R, r):
VCOV = model.cov_params()
b = model.params
Rbeta = R.T.dot(b)
varRbeta = R.T.dot(VCOV).dot(R)
t_0 = (Rbeta - r)/np.sqrt(varRbeta)
df = model.df_resid
pval_unilat_izqda = scs.t.cdf(t_0, df)
pval_bilat = 2*(1 - scs.t.cdf(np.abs(t_0), df))
pval_unilat_dcha = 1 - scs.t.cdf(t_0, df)
return(t_0, pval_unilat_izqda, pval_bilat, pval_unilat_dcha)
# Ejecución del contraste
R = np.array([0, 1, -1, 0])
r = 0
p = t_test(mod_log, R, r)
print('Resultados:', '\n',
'Estadístico t0 = ', round(p[0], 4), '\n',
' p-valor unilat_izqda (Ho: Rbeta >= r) = ', round(p[1], 4), '\n',
' p-valor bilateral (Ho: Rbeta = r) = ', round(p[2], 4), '\n',
' p-valor unilat_dcha (Ho: Rbeta <= r) = ', round(p[3], 4))Resultados:
Estadístico t0 = -5.8143
p-valor unilat_izqda (Ho: Rbeta >= r) = 0.0
p-valor bilateral (Ho: Rbeta = r) = 0.0
p-valor unilat_dcha (Ho: Rbeta <= r) = 1.0Code
# Predicción del HDI para un país con los siguientes indicadores parciales:
# LifeExp = 71.7, ExpEduc = 13.4, GNIpc = 12306.3
nuevo_PAIS = pd.DataFrame({'LifeExp': [71.7], 'ExpEduc': [13.4], 'GNIpc': [12306.3]})
print(nuevo_PAIS) LifeExp ExpEduc GNIpc
0 71.7 13.4 12306.3Code
# HDI esperado con intervalo de confianza
pred_log_HDI_IC = mod_log.get_prediction(nuevo_PAIS).summary_frame(alpha=0.05)
print(pred_log_HDI_IC[['mean','mean_ci_lower', 'mean_ci_upper']]) mean mean_ci_lower mean_ci_upper
0 4.262828 4.256127 4.26953Code
pred_HDI_IC=np.exp(pred_log_HDI_IC)
print(pred_HDI_IC[['mean','mean_ci_lower', 'mean_ci_upper']]) mean mean_ci_lower mean_ci_upper
0 71.01054 70.536244 71.488024