### ¿VIVIR MÁS O VIVIR MEJOR?
# Instalar y cargar librerías necesarias
library(pacman)
## Warning: package 'pacman' was built under R version 4.4.3
p_load(readxl,tidyverse,tidyr,naniar,ggplot2,dplyr,patchwork,skimr,knitr,visdat,reshape2,corrplot,ggcorrplot,caret,randomForest,quantreg,Matrix,qgam,car)
# Cargar y mostrar base de datos
proyecto <- read_excel("C:/Users/MIRZ/Downloads/Datos_proyecto_1.xlsx")
View(proyecto)
# Conocer el nombre de las columnas
names(proyecto)
## [1] "Pais" "Poblacion" "RPC" "EV" "TM" "TN"
## [7] "IF" "TD" "RP" "IDH" "IG" "GPS"
## [13] "GPS_PIB" "PM" "PH"
# Observar qué tipo de datos son
str(proyecto)
## tibble [168 × 15] (S3: tbl_df/tbl/data.frame)
## $ Pais : chr [1:168] "United States" "China" "Japan" "Germany" ...
## $ Poblacion: num [1:168] 3.25e+08 1.42e+09 1.28e+08 8.27e+07 1.34e+09 ...
## $ RPC : num [1:168] 59939 8612 38214 44680 1980 ...
## $ EV : num [1:168] 78.9 76.7 84.5 81.2 69.4 81.2 82.5 75.7 83.6 82.3 ...
## $ TM : num [1:168] 0.098 0.0737 0.129 0.123 0.0907 0.091 0.092 0.0814 0.112 0.086 ...
## $ TN : chr [1:168] "11" "6,77" "6,3" "8,3" ...
## $ IF : chr [1:168] "1,67" "1,18" "1,26" "1,39" ...
## $ TD : chr [1:168] "4,1" "5,1" "2,5" "3,5" ...
## $ RP : chr [1:168] "11,1" "0,0" "16,0" "15,5" ...
## $ IDH : num [1:168] 0.927 0.788 0.92 0.95 0.644 0.94 0.91 0.76 0.906 0.935 ...
## $ IG : chr [1:168] "41,3" "35,7" "32,9" "29,4" ...
## $ GPS : num [1:168] 0.833 0.549 0.859 0.859 0.391 ...
## $ GPS_PIB : num [1:168] 0.139 0.0952 0.0952 0.1014 0.0455 ...
## $ PM : num [1:168] 0.498 0.49 0.512 0.506 0.484 ...
## $ PH : num [1:168] 0.502 0.51 0.488 0.494 0.516 ...
# Cambiamos los datos chr a numerico
# Lista de columnas que contienen comas
cols <- c("TN", "IF", "TD", "RP", "IG")
# Transformación: reemplazar comas por puntos y convertir a numérico
proyecto[cols] <- lapply(proyecto[cols], function(x) as.numeric(gsub(",", ".", x)))
# Verificar la estructura después de la conversión
str(proyecto)
## tibble [168 × 15] (S3: tbl_df/tbl/data.frame)
## $ Pais : chr [1:168] "United States" "China" "Japan" "Germany" ...
## $ Poblacion: num [1:168] 3.25e+08 1.42e+09 1.28e+08 8.27e+07 1.34e+09 ...
## $ RPC : num [1:168] 59939 8612 38214 44680 1980 ...
## $ EV : num [1:168] 78.9 76.7 84.5 81.2 69.4 81.2 82.5 75.7 83.6 82.3 ...
## $ TM : num [1:168] 0.098 0.0737 0.129 0.123 0.0907 0.091 0.092 0.0814 0.112 0.086 ...
## $ TN : num [1:168] 11 6.77 6.3 8.3 16.27 ...
## $ IF : num [1:168] 1.67 1.18 1.26 1.39 2.01 1.57 1.66 1.63 1.21 1.33 ...
## $ TD : num [1:168] 4.1 5.1 2.5 3.5 4.6 4.4 7.3 13.2 6.4 6.6 ...
## $ RP : num [1:168] 11.1 0 16 15.5 21.9 18.6 15.4 26.5 18.9 9.9 ...
## $ IDH : num [1:168] 0.927 0.788 0.92 0.95 0.644 0.94 0.91 0.76 0.906 0.935 ...
## $ IG : num [1:168] 41.3 35.7 32.9 29.4 32.8 33.5 29.7 52 31.5 31.7 ...
## $ GPS : num [1:168] 0.833 0.549 0.859 0.859 0.391 ...
## $ GPS_PIB : num [1:168] 0.139 0.0952 0.0952 0.1014 0.0455 ...
## $ PM : num [1:168] 0.498 0.49 0.512 0.506 0.484 ...
## $ PH : num [1:168] 0.502 0.51 0.488 0.494 0.516 ...
# Tener una vista de las medidas de tendencia central de cada una de las variables
summary(proyecto)
## Pais Poblacion RPC EV
## Length:168 Min. :1.781e+04 Min. : 293 Min. :52.80
## Class :character 1st Qu.:2.735e+06 1st Qu.: 1996 1st Qu.:67.25
## Mode :character Median :9.469e+06 Median : 5572 Median :74.10
## Mean :4.326e+07 Mean : 13474 Mean :72.54
## 3rd Qu.:3.119e+07 3rd Qu.: 16033 3rd Qu.:77.65
## Max. :1.421e+09 Max. :105280 Max. :84.50
## TM TN IF TD
## Min. :0.01870 Min. : 4.400 Min. :0.700 Min. : 0.100
## 1st Qu.:0.06383 1st Qu.: 9.723 1st Qu.:1.500 1st Qu.: 3.500
## Median :0.07815 Median :16.425 Median :2.070 Median : 5.100
## Mean :0.08223 Mean :18.440 Mean :2.482 Mean : 7.329
## 3rd Qu.:0.09632 3rd Qu.:26.740 3rd Qu.:3.277 3rd Qu.: 8.625
## Max. :0.21400 Max. :45.030 Max. :6.750 Max. :50.000
## RP IDH IG GPS
## Min. : 0.00 Min. :0.3870 Min. :14.50 Min. : 0.0648
## 1st Qu.:14.53 1st Qu.:0.6100 1st Qu.:30.95 1st Qu.: 0.3742
## Median :21.65 Median :0.7445 Median :34.60 Median : 0.5649
## Mean :26.02 Mean :0.7301 Mean :35.86 Mean : 0.7876
## 3rd Qu.:36.83 3rd Qu.:0.8495 3rd Qu.:40.70 3rd Qu.: 0.7174
## Max. :89.00 Max. :0.9660 Max. :59.10 Max. :39.2500
## GPS_PIB PM PH
## Min. :0.00150 Min. :0.3050 Min. :0.4510
## 1st Qu.:0.02042 1st Qu.:0.4950 1st Qu.:0.4880
## Median :0.03955 Median :0.5030 Median :0.4970
## Mean :0.07557 Mean :0.4981 Mean :0.5019
## 3rd Qu.:0.06432 3rd Qu.:0.5120 3rd Qu.:0.5050
## Max. :0.94000 Max. :0.5490 Max. :0.6950
# Calcular estadísticas descriptivas de las variables numéricas
proyecto_n <- proyecto[, -1]
summary_mestadisticas <- data.frame(
Media = sapply(proyecto_n, mean, na.rm = TRUE),
Mediana = sapply(proyecto_n, median, na.rm = TRUE),
Moda = sapply(proyecto_n, function(x) names(sort(table(x), decreasing = TRUE))[1]),
Desviacion_Estandar = sapply(proyecto_n, sd, na.rm = TRUE),
Varianza = sapply(proyecto_n, var, na.rm = TRUE)
)
summary_mestadisticas
## Media Mediana Moda Desviacion_Estandar Varianza
## Poblacion 4.325658e+07 9.468717e+06 17808 1.548580e+08 2.398100e+16
## RPC 1.347382e+04 5.572000e+03 293 1.810622e+04 3.278352e+08
## EV 7.254286e+01 7.410000e+01 76.5 7.575430e+00 5.738713e+01
## TM 8.223333e-02 7.815000e-02 0.066 2.873607e-02 8.257619e-04
## TN 1.843982e+01 1.642500e+01 8 9.841835e+00 9.686172e+01
## IF 2.481548e+00 2.070000e+00 1.26 1.244510e+00 1.548805e+00
## TD 7.329048e+00 5.100000e+00 4.4 6.522115e+00 4.253798e+01
## RP 2.602137e+01 2.165000e+01 21.9 1.680814e+01 2.825137e+02
## IDH 7.301012e-01 7.445000e-01 0.394 1.530660e-01 2.342919e-02
## IG 3.585714e+01 3.460000e+01 35.7 8.152466e+00 6.646270e+01
## GPS 7.875732e-01 5.649500e-01 0.345 3.012643e+00 9.076017e+00
## GPS_PIB 7.556964e-02 3.955000e-02 0.0098 1.457538e-01 2.124416e-02
## PM 4.980649e-01 5.030000e-01 0.503 3.212515e-02 1.032025e-03
## PH 5.019351e-01 4.970000e-01 0.497 3.212515e-02 1.032025e-03
# Histograma de la Distribución de la Esperanza de Vida
ggplot(proyecto, aes(x = EV)) +
geom_histogram(binwidth = 2, fill = "skyblue", color = "black") +
labs(title = "Distribución de la Esperanza de Vida",
x = "Esperanza de vida",
y = "Frecuencia") +
theme_minimal()
# Gráfico de barras de Top 10 Países por Esperanza de Vida
top_10_ev <- proyecto %>% arrange(desc(EV)) %>% head(10)
top_10_ev
## # A tibble: 10 × 15
## Pais Poblacion RPC EV TM TN IF TD RP IDH IG GPS
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Japan 127502725 38214 84.5 0.129 6.3 1.26 2.5 16 0.92 32.9 0.859
## 2 Hong … 7306322 46733 84.1 0.087 4.4 0.7 3.1 23.6 0.956 32.4 0.345
## 3 Singa… 5708041 56746 83.8 0.063 7.9 1.04 1.9 26.6 0.949 32.5 0.566
## 4 Italy 60673701 32038 83.6 0.112 6.4 1.21 6.4 18.9 0.906 31.5 0.740
## 5 Spain 46647428 28175 83.4 0.0896 6.61 1.12 10.6 19.7 0.911 31.5 0.743
## 6 Austr… 24584620 53831 83.3 0.073 11.6 1.63 4 13.4 0.946 34.5 0.719
## 7 Icela… 334393 73233 82.9 0.066 11 1.55 3.5 8.8 0.959 23.7 0.840
## 8 South… 51096415 29958 82.8 0.07 4.7 0.75 2.8 18.3 0.929 32.9 0.632
## 9 Israel 8243848 42852 82.8 0.054 19 2.89 2.6 20.7 0.915 37.9 0.685
## 10 Sweden 9904896 54075 82.7 0.09 9.5 1.45 8.6 14.8 0.952 29.5 0.862
## # ℹ 3 more variables: GPS_PIB <dbl>, PM <dbl>, PH <dbl>
Top10_EV <- ggplot(top_10_ev, aes(x = reorder(Pais, EV), y = EV, fill = Pais)) +
geom_bar(stat = "identity") +
coord_flip() +
labs(title = "Top 10 Países por Esperanza de Vida", y = "Esperanza de Vida (años)", x = "País") +
theme_minimal()
Top10_EV
# Gráfico de dispersión del Gasto en salud vs Esperanza de vida
GPSvsEV <- ggplot(proyecto, aes(x = GPS, y = EV)) +
geom_point(alpha = 0.6, color = "blue") +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "Gasto en Salud vs. Esperanza de Vida", x = "% Gasto en Salud", y = "Esperanza de Vida (años)") +
theme_minimal()
GPSvsEV
## `geom_smooth()` using formula = 'y ~ x'
# Gráfico de dispersión del IDH vs Esperanza de vida
IDHvsEV <- ggplot(proyecto, aes(x = IDH, y = EV)) +
geom_point(alpha = 0.6, color = "darkgreen") +
geom_smooth(method = "lm", color = "red", se = TRUE) +
labs(title = "IDH vs. Esperanza de Vida", x = "Índice de Desarrollo Humano", y = "Esperanza de Vida (años)") +
theme_minimal()
IDHvsEV
## `geom_smooth()` using formula = 'y ~ x'
# Realizar matriz de correlacion para identificar las variables que nos sirven
# Calcular la matriz de correlación excluyendo la columna de países
cor_matrix <- cor(proyecto[ , sapply(proyecto, is.numeric)], use = "pairwise.complete.obs")
# Visualizar la matriz de correlación como un mapa de calor
corrplot(cor_matrix, method = "color", type = "upper",
tl.cex = 0.8, tl.col = "black",
title = "Matriz de correlaciones",
mar = c(0,0,1,0))
# Visualización alternativa
cor_melted <- melt(cor_matrix)
ggplot(data = cor_melted, aes(x = Var1, y = Var2, fill = value)) +
geom_tile(color = "white") +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name = "Correlación") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
ggtitle("Matriz de correlación entre variables")
# Detección de datos faltantes
# Conteo por variable
cat("Cantidad de valores faltantes por variable:\n")
## Cantidad de valores faltantes por variable:
colSums(is.na(proyecto))
## Pais Poblacion RPC EV TM TN IF TD
## 0 0 0 0 0 0 0 0
## RP IDH IG GPS GPS_PIB PM PH
## 0 0 0 0 0 0 0
# Visualización con visdat y naniar
vis_miss(proyecto) + ggtitle("Visualización de valores faltantes")
# Detección de valores atípicos
# Boxplots por variable numérica para detectar outliers
# Convertimos a formato largo para graficar todas las variables numéricas
proyecto_long <- proyecto %>%
select(where(is.numeric)) %>%
pivot_longer(cols = everything(), names_to = "Variable", values_to = "Valor")
ggplot(proyecto_long, aes(x = Variable, y = Valor)) +
geom_boxplot(fill = "orange", outlier.color = "red", outlier.shape = 8) +
theme_minimal() +
labs(title = "Detección de valores atípicos por variable",
x = "Variable", y = "Valor") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Preparar los datos para aplicar un modelo
# Normalizar los datos
# Seleccionar solo las columnas numéricas del dataframe
numeric_cols <- proyecto[, sapply(proyecto, is.numeric)]
# Verificar la estructura de las columnas numéricas
str(numeric_cols)
## tibble [168 × 14] (S3: tbl_df/tbl/data.frame)
## $ Poblacion: num [1:168] 3.25e+08 1.42e+09 1.28e+08 8.27e+07 1.34e+09 ...
## $ RPC : num [1:168] 59939 8612 38214 44680 1980 ...
## $ EV : num [1:168] 78.9 76.7 84.5 81.2 69.4 81.2 82.5 75.7 83.6 82.3 ...
## $ TM : num [1:168] 0.098 0.0737 0.129 0.123 0.0907 0.091 0.092 0.0814 0.112 0.086 ...
## $ TN : num [1:168] 11 6.77 6.3 8.3 16.27 ...
## $ IF : num [1:168] 1.67 1.18 1.26 1.39 2.01 1.57 1.66 1.63 1.21 1.33 ...
## $ TD : num [1:168] 4.1 5.1 2.5 3.5 4.6 4.4 7.3 13.2 6.4 6.6 ...
## $ RP : num [1:168] 11.1 0 16 15.5 21.9 18.6 15.4 26.5 18.9 9.9 ...
## $ IDH : num [1:168] 0.927 0.788 0.92 0.95 0.644 0.94 0.91 0.76 0.906 0.935 ...
## $ IG : num [1:168] 41.3 35.7 32.9 29.4 32.8 33.5 29.7 52 31.5 31.7 ...
## $ GPS : num [1:168] 0.833 0.549 0.859 0.859 0.391 ...
## $ GPS_PIB : num [1:168] 0.139 0.0952 0.0952 0.1014 0.0455 ...
## $ PM : num [1:168] 0.498 0.49 0.512 0.506 0.484 ...
## $ PH : num [1:168] 0.502 0.51 0.488 0.494 0.516 ...
# Normalizar los datos numéricos
datos_normalizados <- as.data.frame(scale(numeric_cols))
# Verificar los datos normalizados
head(datos_normalizados)
## Poblacion RPC EV TM TN IF TD
## 1 1.8199137 2.5662552 0.8391792 0.5486716 -0.7559384 -0.6521022 -0.4950921
## 2 8.8969593 -0.2685163 0.5487666 -0.2969554 -1.1857363 -1.0458315 -0.3417676
## 3 0.5440219 1.3663914 1.5784112 1.6274550 -1.2334916 -0.9815492 -0.7404113
## 4 0.2544384 1.7235062 1.1427923 1.4186583 -1.0302775 -0.8770904 -0.5870868
## 5 8.3652140 -0.6347992 -0.4148751 0.2946355 -0.2204692 -0.3789023 -0.4184299
## 6 0.1515639 1.4391841 1.1427923 0.3050753 -0.8575455 -0.7324551 -0.4490948
## RP IDH IG GPS GPS_PIB PM
## 1 -0.8877465 1.2863656 0.6676332 0.01511191 0.43518845 -0.002019631
## 2 -1.5481407 0.3782605 -0.0192755 -0.07929025 0.13468164 -0.251045688
## 3 -0.5962211 1.2406337 -0.3627299 0.02364262 0.13468164 0.433775969
## 4 -0.6259686 1.4366276 -0.7920478 0.02380859 0.17721914 0.247006426
## 5 -0.2452007 -0.5625103 -0.3749961 -0.13160312 -0.20630439 -0.437815231
## 6 -0.4415342 1.3712963 -0.2891325 0.01033205 0.09283024 0.309262941
## PH
## 1 0.002019631
## 2 0.251045688
## 3 -0.433775969
## 4 -0.247006426
## 5 0.437815231
## 6 -0.309262941
# Calcular la matriz de correlación de los datos normalizados
cor_matrix <- cor(datos_normalizados)
cor_matrix
## Poblacion RPC EV TM TN
## Poblacion 1.000000000 -0.041927905 0.009665429 0.015283903 -0.07706500
## RPC -0.041927905 1.000000000 0.651488267 0.003301292 -0.54768481
## EV 0.009665429 0.651488267 1.000000000 -0.034609205 -0.87207692
## TM 0.015283903 0.003301292 -0.034609205 1.000000000 -0.23402804
## TN -0.077064998 -0.547684814 -0.872076922 -0.234028039 1.00000000
## IF -0.078583134 -0.506007781 -0.844075665 -0.180217363 0.98042467
## TD -0.052572603 -0.170023840 -0.102695462 0.046072117 0.11765910
## RP -0.124548212 -0.369637846 -0.621648920 0.014973882 0.59236093
## IDH -0.006480292 0.722956767 0.905845947 0.100788143 -0.90291069
## IG -0.012638746 -0.335500673 -0.292700352 -0.107245173 0.31790251
## GPS -0.023399654 -0.023840319 -0.039832570 -0.032008708 0.08783549
## GPS_PIB -0.020005525 -0.079222375 -0.188159815 0.058173527 0.25338785
## PM -0.063472145 -0.082677147 0.021521707 0.381253145 -0.04747079
## PH 0.063472145 0.082677147 -0.021521707 -0.381253145 0.04747079
## IF TD RP IDH IG
## Poblacion -0.07858313 -0.05257260 -0.124548212 -0.006480292 -0.01263875
## RPC -0.50600778 -0.17002384 -0.369637846 0.722956767 -0.33550067
## EV -0.84407566 -0.10269546 -0.621648920 0.905845947 -0.29270035
## TM -0.18021736 0.04607212 0.014973882 0.100788143 -0.10724517
## TN 0.98042467 0.11765910 0.592360927 -0.902910687 0.31790251
## IF 1.00000000 0.08964073 0.572007006 -0.866004571 0.24326349
## TD 0.08964073 1.00000000 0.121309155 -0.136827742 0.11899736
## RP 0.57200701 0.12130916 1.000000000 -0.630492231 0.31568179
## IDH -0.86600457 -0.13682774 -0.630492231 1.000000000 -0.33782520
## IG 0.24326349 0.11899736 0.315681791 -0.337825202 1.00000000
## GPS 0.09232714 0.04057730 -0.001932711 -0.059180357 0.01367857
## GPS_PIB 0.30644476 -0.05748814 0.109954623 -0.217761819 -0.08820181
## PM -0.07165410 0.19940539 0.118181026 -0.047220496 0.15091665
## PH 0.07165410 -0.19940539 -0.118181026 0.047220496 -0.15091665
## GPS GPS_PIB PM PH
## Poblacion -0.023399654 -0.020005525 -0.06347215 0.06347215
## RPC -0.023840319 -0.079222375 -0.08267715 0.08267715
## EV -0.039832570 -0.188159815 0.02152171 -0.02152171
## TM -0.032008708 0.058173527 0.38125314 -0.38125314
## TN 0.087835489 0.253387848 -0.04747079 0.04747079
## IF 0.092327141 0.306444756 -0.07165410 0.07165410
## TD 0.040577304 -0.057488142 0.19940539 -0.19940539
## RP -0.001932711 0.109954623 0.11818103 -0.11818103
## IDH -0.059180357 -0.217761819 -0.04722050 0.04722050
## IG 0.013678567 -0.088201813 0.15091665 -0.15091665
## GPS 1.000000000 -0.009452073 0.03460990 -0.03460990
## GPS_PIB -0.009452073 1.000000000 0.02282926 -0.02282926
## PM 0.034609897 0.022829259 1.00000000 -1.00000000
## PH -0.034609897 -0.022829259 -1.00000000 1.00000000
# Extraer la correlación con la variable "EV"
cor_ev <- cor_matrix[, "EV"]
cor_ev
## Poblacion RPC EV TM TN IF
## 0.009665429 0.651488267 1.000000000 -0.034609205 -0.872076922 -0.844075665
## TD RP IDH IG GPS GPS_PIB
## -0.102695462 -0.621648920 0.905845947 -0.292700352 -0.039832570 -0.188159815
## PM PH
## 0.021521707 -0.021521707
# Seleccionar variables con correlación significativa (umbral de 0.3)
high_cor_vars <- names(cor_ev[abs(cor_ev) > 0.3])
high_cor_vars #Variables seleccionadas
## [1] "RPC" "EV" "TN" "IF" "RP" "IDH"
# Filtrar el conjunto de datos solo con variables significativas
datos_reducidos <- datos_normalizados[, high_cor_vars]
# Imprimimos la imagen de la matriz
ggcorrplot(cor_matrix,
method = "circle", # Opción "square" para cuadrados
type = "lower", # Mostrar solo la mitad inferior
lab = TRUE, # Etiquetas con valores numéricos
lab_size = 3, # Tamaño de las etiquetas
title = "Matriz de Correlación de Variables Normalizadas",
colors = c("#6D9EC1", "white", "#E46726")) # Colores de gradiente
# Verificamos multicolinealidad
# Ajustar un modelo lineal con todas las variables significativas
modelo_vif <- lm(EV ~ ., data = datos_reducidos)
vif_values <- vif(modelo_vif)
print(vif_values)
## RPC TN IF RP IDH
## 2.455057 37.121462 27.290880 1.704114 10.179256
# Eliminar TN (por ejemplo) y recalcular modelo y VIF
modelo_vif2 <- lm(EV ~ RPC + IF + RP + IDH, data = datos_reducidos)
vif_values2 <- vif(modelo_vif2)
print(vif_values2)
## RPC IF RP IDH
## 2.435379 4.549623 1.704056 8.024101
# Eliminar IF (por ejemplo) y recalcular modelo y VIF
modelo_vif3 <- lm(EV ~ RPC + TN + RP + IDH, data = datos_reducidos)
vif_values3 <- vif(modelo_vif3)
print(vif_values3)
## RPC TN RP IDH
## 2.446980 6.188464 1.704114 9.992293
# Calcular el VIF para cada variable
vif_values <- vif(modelo_vif)
vif_values
## RPC TN IF RP IDH
## 2.455057 37.121462 27.290880 1.704114 10.179256
# MODELO CUANTIL
# Modelo de regresión cuantil para tau = 0.5 (mediana)
if (!require(quantreg)) install.packages("quantreg")
library(quantreg)
modelo_qr <- rq(EV ~ RPC + IF + RP + IDH, tau = 0.5, data = datos_reducidos)
summary(modelo_qr)
##
## Call: rq(formula = EV ~ RPC + IF + RP + IDH, tau = 0.5, data = datos_reducidos)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 0.02469 -0.02908 0.07225
## RPC 0.09140 0.03045 0.15997
## IF -0.33065 -0.45047 -0.05643
## RP -0.02835 -0.21850 0.07218
## IDH 0.51110 0.41414 0.72364
#Para varios cuantiles
taus <- seq(0.1, 0.9, by = 0.1)
modelos_cuantiles <- lapply(taus, function(tau) {
rq(EV ~ RPC + IF + RP + IDH, tau = tau, data = datos_reducidos)
})
coeficientes <- sapply(modelos_cuantiles, coef)
print(coeficientes)
## [,1] [,2] [,3] [,4] [,5]
## (Intercept) -0.52293441 -0.31998520 -0.20414224 -0.09781509 0.02468586
## RPC 0.06502898 0.07854523 0.06749328 0.04411055 0.09139820
## IF -0.19392417 -0.24906155 -0.18800831 -0.23557384 -0.33064831
## RP -0.23786521 -0.14083430 -0.18800389 -0.09257021 -0.02835079
## IDH 0.66669247 0.65012240 0.65989589 0.66777796 0.51109984
## [,6] [,7] [,8] [,9]
## (Intercept) 0.12192191 0.20444411 0.340686725 0.46510968
## RPC 0.05777073 0.01645561 -0.000635823 0.04786969
## IF -0.28020612 -0.22988308 -0.177711491 -0.20605513
## RP 0.02597823 0.04353982 0.044285475 0.06410422
## IDH 0.58276683 0.65598936 0.656908584 0.58206361
# Ajustar el modelo de regresión cuantil para la mediana (tau = 0.5)
modelo_qr <- rq(EV ~ RPC + IF + RP + IDH, tau = 0.5, data = datos_reducidos)
# Mostrar resumen del modelo con test bootstrap para inferencia
summary(modelo_qr, se = "boot")
##
## Call: rq(formula = EV ~ RPC + IF + RP + IDH, tau = 0.5, data = datos_reducidos)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 0.02469 0.03625 0.68096 0.49686
## RPC 0.09140 0.04699 1.94506 0.05349
## IF -0.33065 0.13575 -2.43579 0.01594
## RP -0.02835 0.08488 -0.33401 0.73880
## IDH 0.51110 0.12764 4.00431 0.00009
# Prueba de bondad entre cuantiles (Wald)
# Ajustar un modelo cuantil para varios cuantiles (ej. 0.25, 0.5 y 0.75)
modelo_wald <- rq(EV ~ RPC + IF + RP + IDH, tau = c(0.25, 0.5, 0.75), data = datos_reducidos)
# Aplicar prueba de Wald (prueba de igualdad de coeficientes entre cuantiles)
anova(modelo_wald)
## Quantile Regression Analysis of Deviance Table
##
## Model: EV ~ RPC + IF + RP + IDH
## Joint Test of Equality of Slopes: tau in { 0.25 0.5 0.75 }
##
## Df Resid Df F value Pr(>F)
## 1 8 496 4.0935 9.695e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Medida de bondad de ajuste
fitnull <- rq(EV ~ 1, tau = c(0.25, 0.5, 0.75), data = datos_reducidos, method = "fn")
fit4 <- rq(EV ~ RPC + IF + RP + IDH, tau = c(0.25, 0.5, 0.75), data = datos_reducidos, method = "fn")
pseudoR2_fit4 <- 1 - fit4$rho / fitnull$rho
pseudoR2_fit4
## [1] 0.6555897 0.6208897 0.5998999
# Gráfico de la prueba de bondad
if (!require(ggplot2)) install.packages("ggplot2")
library(ggplot2)
# Crear un data frame con los resultados
cuantiles <- c(0.25, 0.5, 0.75)
pseudoR2 <- pseudoR2_fit4
df_pseudoR2 <- data.frame(
Cuantil = cuantiles,
Pseudo_R2 = pseudoR2
)
ggplot(df_pseudoR2, aes(x = Cuantil, y = Pseudo_R2)) +
geom_point(color = "blue", size = 4) +
geom_line(color = "blue", size = 1) +
ylim(0,1) +
labs(
title = "Pseudo R-cuadrado para regresión cuantílica en diferentes cuantiles",
x = "Cuantil (τ)",
y = "Pseudo R-cuadrado"
) +
theme_minimal() +
theme(
text = element_text(size = 14)
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Visualización simple para una variable (RPC) con el percentil 90
modelo_rpc <- rq(EV ~ RPC, tau = 0.9, data = datos_reducidos)
summary(modelo_rpc, se = "boot")
##
## Call: rq(formula = EV ~ RPC, tau = 0.9, data = datos_reducidos)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 0.82685 0.04004 20.65009 0.00000
## RPC 0.57902 0.06937 8.34630 0.00000
# Visualización
ggplot(datos_reducidos, aes(x = RPC, y = EV)) +
geom_point() +
geom_abline(intercept = coef(modelo_rpc)[1],
slope = coef(modelo_rpc)[2],
color = "red", size = 1.2) +
labs(title = "Regresión Cuantil (tau = 0.9)",
x = "Renta per cápita (RPC)", y = "Esperanza de vida (EV)")
# Calcular residuos para cada modelo
residuos <- lapply(modelos_cuantiles, function(modelo) modelo$residuals)
# Crear un data frame con los residuos y los cuantiles
df_residuos <- data.frame(
Residuos = unlist(residuos),
Cuantil = factor(rep(cuantiles, each = nrow(datos_reducidos)))
)
# Gráfico de dispersión de los residuos
ggplot(df_residuos, aes(x = 1:nrow(df_residuos), y = Residuos, color = Cuantil)) +
geom_point(alpha = 0.6) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(
title = "Gráfico de Residuos por Cuantil",
x = "Índice",
y = "Residuos"
) +
theme_bw()
# Extraer coeficientes de cada modelo
coeficientes <- lapply(modelos_cuantiles, coef)
coeficientes
## [[1]]
## (Intercept) RPC IF RP IDH
## -0.52293441 0.06502898 -0.19392417 -0.23786521 0.66669247
##
## [[2]]
## (Intercept) RPC IF RP IDH
## -0.31998520 0.07854523 -0.24906155 -0.14083430 0.65012240
##
## [[3]]
## (Intercept) RPC IF RP IDH
## -0.20414224 0.06749328 -0.18800831 -0.18800389 0.65989589
##
## [[4]]
## (Intercept) RPC IF RP IDH
## -0.09781509 0.04411055 -0.23557384 -0.09257021 0.66777796
##
## [[5]]
## (Intercept) RPC IF RP IDH
## 0.02468586 0.09139820 -0.33064831 -0.02835079 0.51109984
##
## [[6]]
## (Intercept) RPC IF RP IDH
## 0.12192191 0.05777073 -0.28020612 0.02597823 0.58276683
##
## [[7]]
## (Intercept) RPC IF RP IDH
## 0.20444411 0.01645561 -0.22988308 0.04353982 0.65598936
##
## [[8]]
## (Intercept) RPC IF RP IDH
## 0.340686725 -0.000635823 -0.177711491 0.044285475 0.656908584
##
## [[9]]
## (Intercept) RPC IF RP IDH
## 0.46510968 0.04786969 -0.20605513 0.06410422 0.58206361
# Crear un data frame con intercepto y pendiente de RPC (primer predictor) y los cuantiles
df_coeficientes <- data.frame(
Intercepto = sapply(coeficientes, function(x) x[1]),
Pendiente_RPC = sapply(coeficientes, function(x) x["RPC"]), # Ajusta el nombre si cambia
Cuantil = cuantiles
)
# Gráfico de coeficientes del intercepto por cuantil
ggplot(df_coeficientes, aes(x = Cuantil, y = Intercepto)) +
geom_line(color = "blue") +
geom_point(color = "blue", size = 3) +
labs(
title = "Coeficientes del Intercepto por Cuantil",
x = "Cuantil",
y = "Intercepto"
) +
theme_bw()
# Gráfico de coeficientes de la pendiente de RPC por cuantil
ggplot(df_coeficientes, aes(x = Cuantil, y = Pendiente_RPC)) +
geom_line(color = "darkgreen") +
geom_point(color = "darkgreen", size = 3) +
labs(
title = "Coeficientes de RPC por Cuantil",
x = "Cuantil",
y = "Pendiente"
) +
theme_bw()
# Distribución de residuos: densidad por cuantil
ggplot(df_residuos, aes(x = Residuos, fill = Cuantil)) +
geom_density(alpha = 0.5) +
labs(
title = "Distribución de los Residuos por Cuantil",
x = "Residuos",
y = "Densidad"
) +
theme_bw()
# Distribución de residuos: histograma por cuantil
ggplot(df_residuos, aes(x = Residuos, fill = Cuantil)) +
geom_histogram(alpha = 0.5, position = "identity", bins = 30) +
labs(
title = "Histograma de los Residuos por Cuantil",
x = "Residuos",
y = "Frecuencia"
) +
theme_bw()
# Obtener predicciones para cada modelo
predicciones <- lapply(modelos_cuantiles, function(modelo) {
predict(modelo, newdata = datos_reducidos)
})
# Crear data frame con predicciones, valores reales y cuantiles
df_pred <- data.frame(
Observado = rep(datos_reducidos$EV, times = length(taus)),
Predicho = unlist(predicciones),
Cuantil = factor(rep(taus, each = nrow(datos_reducidos)))
)
# Gráfico de dispersión predicho vs observado, coloreado por cuantil
ggplot(df_pred, aes(x = Observado, y = Predicho, color = Cuantil)) +
geom_point(alpha = 0.6) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black") +
labs(
title = "Predicciones vs Valores Observados por Cuantil",
x = "Esperanza de Vida Observada",
y = "Esperanza de Vida Predicha"
) +
theme_minimal() +
theme(text = element_text(size = 14))
# Intervalos de confianza
# Definir cuantiles para el intervalo de confianza
taus_ic <- c(0.025, 0.05, 0.5, 0.95, 0.975)
taus_ic
## [1] 0.025 0.050 0.500 0.950 0.975
# Ajustar los modelos para esos cuantiles
modelos_ic <- lapply(taus_ic, function(tau) {
rq(EV ~ RPC + IF + RP + IDH, tau = tau, data = datos_reducidos)
})
# Obtener predicciones para cada modelo
preds_ic <- lapply(modelos_ic, function(mod) predict(mod, newdata = datos_reducidos))
# Crear un data frame con predicciones y observados
df_ic <- data.frame(
Observado = datos_reducidos$EV,
Q2.5 = preds_ic[[1]],
Q5 = preds_ic[[2]],
Mediana = preds_ic[[3]],
Q95 = preds_ic[[4]],
Q97.5 = preds_ic[[5]]
)
# Ordenar por valor observado para que la línea sea clara
df_ic <- df_ic[order(df_ic$Observado), ]
# Graficar. mediana con intervalos de confianza
ggplot(df_ic, aes(x = Observado)) +
geom_line(aes(y = Mediana), color = "blue", size = 1, linetype = "solid") +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "red", alpha = 0.2) + # 95% IC
geom_ribbon(aes(ymin = Q5, ymax = Q95), fill = "orange", alpha = 0.3) + # 90% IC
labs(
title = "Intervalos de Confianza por Cuantiles en Regresión Cuantílica",
x = "Esperanza de Vida Observada",
y = "Esperanza de Vida Predicha"
) +
theme_minimal() +
theme(text = element_text(size = 14))
