objetivo 4
Marta
2025-05-25
El dataset tenía muchos problemas, como esta con comas, al pasar a tipo numerico los numeros o cambiaban completamente o se sustituian por NA, el porcentaje tp estaba en un formato incorrecto desde el momento de exportarlo
data = read.csv("espcalob.csv")
data1 = data.frame("variable" = colnames(data),
"tipo" = c(rep("categorical",3), rep("numerical", 10),
rep("categorical",2), rep("numerical",4)), stringsAsFactors = FALSE)El siguiente código es una función para pasar las columnas al tipo de datos que especifico. Cuando es de tipo numerica se sustituyen las comas por los puntos, porque anteriormente solo las variables con valores de numeros redondos eran reconocidas como variables numericas.
## Pais Año Estatus Esperanza.de.vida
## Austria : 48 2000 : 78 Developed :1056 Min. :69.90
## Belgium : 48 2001 : 78 Developing: 192 1st Qu.:75.30
## Bulgaria: 48 2002 : 78 Median :78.50
## Croatia : 48 2003 : 78 Mean :78.45
## Cyprus : 48 2004 : 78 3rd Qu.:81.03
## Czechia : 48 2005 : 78 Max. :89.00
## (Other) :960 (Other):780
## Mortalidad.adulta Alcohol..L.año. Porcentaje.gastos IMC
## Min. : 1.00 Min. : 0.01 Min. : 0.0 Min. : 5.10
## 1st Qu.: 63.00 1st Qu.: 9.56 1st Qu.: 143.7 1st Qu.:54.90
## Median : 81.00 Median :10.98 Median : 703.0 Median :57.60
## Mean : 89.23 Mean :10.55 Mean : 2231.0 Mean :53.55
## 3rd Qu.:124.00 3rd Qu.:12.12 3rd Qu.: 3298.0 3rd Qu.:59.90
## Max. :229.00 Max. :17.87 Max. :18961.3 Max. :69.60
## NA's :78
## Gasto.total PIB Poblacion Delgadez.1.19.años.
## Min. : 1.100 Min. : 12.28 Min. : 123 Min. :0.300
## 1st Qu.: 6.320 1st Qu.: 3349.99 1st Qu.: 246354 1st Qu.:0.800
## Median : 7.535 Median : 12102.12 Median : 1468719 Median :1.200
## Mean : 7.250 Mean : 18208.58 Mean : 7826012 Mean :1.423
## 3rd Qu.: 8.920 3rd Qu.: 27019.43 3rd Qu.: 5759450 3rd Qu.:2.000
## Max. :11.970 Max. :119172.74 Max. :82534176 Max. :4.000
## NA's :78 NA's :96 NA's :96
## IDH Region Sexo Media.obes
## Min. :0.703 Europe:1248 Both sexes:416 Min. : 9.43
## 1st Qu.:0.805 Men :416 1st Qu.:15.56
## Median :0.844 Women :416 Median :18.42
## Mean :0.838 Mean :18.50
## 3rd Qu.:0.876 3rd Qu.:20.95
## Max. :0.926 Max. :32.88
## NA's :48
## Animal.kcal Total.kcal Vegetal.kcal
## Min. : 540.0 Min. :2556 Min. :1840
## 1st Qu.: 837.2 1st Qu.:3142 1st Qu.:2177
## Median : 929.5 Median :3342 Median :2350
## Mean : 955.9 Mean :3297 Mean :2341
## 3rd Qu.:1076.8 3rd Qu.:3505 3rd Qu.:2497
## Max. :1427.0 Max. :3804 Max. :2887
## Hay 396 datos faltantes. Proceden de variables como Alchol, Gasto, PIB, Poblacion y IDH. De estos, solo nos interesa estudiar en general el IDH y como se relaciona con la obesidad y las calorias. Asi que los que no tienen IDH se eliminan. Ahora ya no hay NA.
cat("Total de valores faltantes en el conjunto de datos:", sum(is.na(data)), "\n")## Total de valores faltantes en el conjunto de datos: 0cat("Variables con valores faltantes:\n")## Variables con valores faltantes:print(colSums(is.na(data))[colSums(is.na(data)) > 0])## named numeric(0)168 filas se ven afectadas por los datos faltantes, pero solo interesan los respectivos al IDH.
print(data[!complete.cases(data), ])## [1] Pais Año Estatus
## [4] Esperanza.de.vida Mortalidad.adulta Porcentaje.gastos
## [7] IMC Delgadez.1.19.años. IDH
## [10] Region Sexo Media.obes
## [13] Animal.kcal Total.kcal Vegetal.kcal
## <0 rows> (o 0- extensión row.names)Se eliminan las filas sin el IDH, ya que nos centramos en estudiar la obesidad y las calorias en base a eso.
#ahora tiene 1200 observaciónes
dataset = data[!is.na(data$IDH), ]
# Boxplot de variables escaladas
datos_escalados <- dataset %>%
select(where(is.numeric)) %>%
scale() %>%
as.data.frame() %>%
pivot_longer(everything(), names_to = "variable", values_to = "valor")
ggplot(datos_escalados, aes(x = variable, y = valor)) +
geom_boxplot(fill = "lightblue", outlier.color = "red") +
labs(title = "Boxplots", x = "", y = "Valores") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
#Selecciono solo 2015 y de ambos sexos. Como hay datos para hombres, para mujeres y por años, concreto unos valores.
dataset2015 <- dataset %>%
filter(Sexo == "Both sexes", Año == 2015)CON DATASET2015:
Cluster 1: Países con obesidad y dietas normales en calorías (dietas altas en calorías y obesidad, pero IDH no alto). Cluster 2: Corresponde a países desarrollados (alto IDH, buena alimentación, baja obesidad).
datos_numericos <- scale(dataset2015[, c("IDH", "Total.kcal", "Media.obes")])
set.seed(123)
kmeans_result <- kmeans(datos_numericos, centers = 2)
dataset2015$cluster_kmeans <- as.factor(kmeans_result$cluster)
ggplot(dataset2015, aes(x = Media.obes, y = IDH, color = cluster_kmeans)) +
geom_point() + theme_minimal()
kmeans_result$centers## IDH Total.kcal Media.obes
## 1 -0.8053475 -0.4546968 0.6412679
## 2 0.8724598 0.4925882 -0.6947069clusters_2015 <- dataset2015 %>%select(Pais, cluster_kmeans)
dataset <- dataset %>% left_join(clusters_2015, by = "Pais") # Asigna el cluster de 2015 a todos los años
#En dataset solo escojo los datos de ambos sexos
dataset <- dataset %>%
filter(Sexo == "Both sexes")Ahora se realiza un gráfico para ver como han evolucionado los valores de los paises a lo largo de los años. Como hay varios paises, se estudiara según los clusters anteriores de cada pais, usando la media de obesidad por cluster y por año como valor a representar (k means).
ggplot(dataset, aes(
x = Año,
y = Media.obes,
color = as.factor(cluster_kmeans),
group = interaction(Pais, cluster_kmeans))) +
geom_line(alpha = 0.4) +geom_smooth(aes(group = cluster_kmeans), se = FALSE, linewidth = 1.5) +
labs(title = "Evolución del Media.obes por país", x = "Año",
y = "Media.obes", color = "Cluster") + scale_color_manual(values =
c("#1b9e77", "#7570b3"),labels = c("Cluster 1", "Cluster 2")) +theme_minimal()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(dataset, aes(
x = Año,
y = IDH,
color = as.factor(cluster_kmeans),
group = interaction(Pais, cluster_kmeans) )) +
geom_line(alpha = 0.4) +geom_smooth(aes(group = cluster_kmeans), se = FALSE, linewidth = 1.5) +
labs(title = "Evolución del IDH por país", x = "Año",
y = "IDH", color = "Cluster") + scale_color_manual(values =
c("#1b9e77","#7570b3"),labels = c("Cluster 1", "Cluster 2")) +theme_minimal()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(dataset, aes(
x = Año,
y = Total.kcal,
color = as.factor(cluster_kmeans),
group = interaction(Pais, cluster_kmeans))) +
geom_line(alpha = 0.4) +geom_smooth(aes(group = cluster_kmeans), se = FALSE, linewidth = 1.5) +
labs(title = "Evolución del Total.kcal por país", x = "Año",
y = "Total.kcal", color = "Cluster") + scale_color_manual(values =
c("#1b9e77", "#7570b3"),labels = c("Cluster 1", "Cluster 2")) +theme_minimal()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
cluster1_data <- dataset %>%
filter(cluster_kmeans == 1) %>% group_by(Año) %>%
summarise(Total = mean(Total.kcal, na.rm = TRUE), Animal = mean(Animal.kcal, na.rm = TRUE),
Vegetal = mean(Vegetal.kcal, na.rm = TRUE)) %>% tidyr::pivot_longer(cols = c(Total, Animal, Vegetal),
names_to = "Tipo", values_to = "kcal")
ggplot(cluster1_data, aes(x = Año, y = kcal, color = Tipo, linetype = Tipo)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
scale_color_manual( values = c("Total" = "black", "Animal" = "#d90f09", "Vegetal" = "#1b9e77")) + scale_linetype_manual( values = c("Total" = "solid", "Animal" = "dashed", "Vegetal" = "dotted") ) +
labs( title = "Consumo Calórico en Cluster 1 (2000-2015)",
subtitle = "Total.kcal = Animal.kcal + Vegetal.kcal", y = "kcal per cápita/día", x = "Año" ) + theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
cluster1_data <- dataset %>%
filter(cluster_kmeans == 2) %>% group_by(Año) %>%
summarise(Total = mean(Total.kcal, na.rm = TRUE), Animal = mean(Animal.kcal, na.rm = TRUE),
Vegetal = mean(Vegetal.kcal, na.rm = TRUE)) %>% tidyr::pivot_longer(cols = c(Total, Animal, Vegetal),
names_to = "Tipo", values_to = "kcal")
ggplot(cluster1_data, aes(x = Año, y = kcal, color = Tipo, linetype = Tipo)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
scale_color_manual( values = c("Total" = "black", "Animal" = "#d90f09", "Vegetal" = "#1b9e77")) + scale_linetype_manual( values = c("Total" = "solid", "Animal" = "dashed", "Vegetal" = "dotted") ) +
labs( title = "Consumo Calórico en Cluster 2 (2000-2015)",
subtitle = "Total.kcal = Animal.kcal + Vegetal.kcal", y = "kcal per cápita/día", x = "Año" ) + theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
#TESTS
Cual es el mejor metodo para estudiar la relación? Spearman
ggplot(dataset2015, aes(x = Media.obes, y = IDH, color = as.factor(cluster_kmeans))) +
geom_point(size = 3, alpha = 0.7) +
labs(title = "Asignación de Clusters: IDH vs Obesidad")
dist_mat <- dist(scale(dataset2015[, c("Media.obes", "IDH")]))
silhouette_score <- silhouette(as.numeric(dataset2015$cluster_kmeans), dist_mat)
summary(silhouette_score)## Silhouette of 25 units in 2 clusters from silhouette.default(x = as.numeric(dataset2015$cluster_kmeans), dist = dist_mat) :
## Cluster sizes and average silhouette widths:
## 13 12
## 0.5029029 0.5224075
## Individual silhouette widths:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1045 0.4567 0.5645 0.5123 0.6083 0.6421plot(silhouette_score)
#El coeficiente de silhouette no es muy alto, pero es mejor con 2 clusters que con 3 o 4Correlación (r): -0.5138462 -> Relación negativa moderada. A mayor IDH, menor obesidad.
cor.test(dataset2015$IDH, dataset2015$Media.obes, method = "spearman")##
## Spearman's rank correlation rho
##
## data: dataset2015$IDH and dataset2015$Media.obes
## S = 3936, p-value = 0.009439
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.5138462variables <- dataset[, c("Total.kcal", "Media.obes", "IDH", "Vegetal.kcal", "Animal.kcal")]
cor_pearson <- round(cor(variables, method = "pearson"), 2)
cor_spearman <- round(cor(variables, method = "spearman"), 2)
diff_matrix <- abs(cor_pearson - cor_spearman)
print(diff_matrix) #spearman es mejor porque hay relaciones no liniales, por eso es mas alto su valor## Total.kcal Media.obes IDH Vegetal.kcal Animal.kcal
## Total.kcal 0.00 0.05 0.03 0.01 0.05
## Media.obes 0.05 0.00 0.06 0.02 0.03
## IDH 0.03 0.06 0.00 0.01 0.00
## Vegetal.kcal 0.01 0.02 0.01 0.00 0.01
## Animal.kcal 0.05 0.03 0.00 0.01 0.00Estudio las relaciones con el dataset que contiene información sobre diversos años:
El consumo calórico total (Total.kcal) no está asociado significativamente con la obesidad A mayor Kcal, mayor IDH, pero su poca correlación es lineal.
cor_matrix <- cor(variables, method = "spearman")
corrplot(cor_matrix, method = "number")
cor.test(dataset2015$IDH, dataset2015$Media.obes, method = "spearman")##
## Spearman's rank correlation rho
##
## data: dataset2015$IDH and dataset2015$Media.obes
## S = 3936, p-value = 0.009439
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.5138462cor.test(dataset2015$Total.kcal, dataset2015$Media.obes, method = "spearman")##
## Spearman's rank correlation rho
##
## data: dataset2015$Total.kcal and dataset2015$Media.obes
## S = 2866, p-value = 0.6253
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.1023077cor.test(dataset2015$IDH, dataset2015$Total.kcal, method = "pearson")##
## Pearson's product-moment correlation
##
## data: dataset2015$IDH and dataset2015$Total.kcal
## t = 3.0224, df = 23, p-value = 0.006062
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1748815 0.7667649
## sample estimates:
## cor
## 0.5331698cor.test(dataset2015$IDH, dataset2015$Animal.kcal, method = "pearson")##
## Pearson's product-moment correlation
##
## data: dataset2015$IDH and dataset2015$Animal.kcal
## t = 3.5867, df = 23, p-value = 0.00156
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.2669538 0.8038192
## sample estimates:
## cor
## 0.5989106Significancia estadística (p=0.418): El consumo de calorías vegetales no está asociado significativamente con la obesidad en tu muestra.
Coeficiente de correlación (r=-0.312): A mayor consumo de calorías animales, menor obesidad, pero la asociación no es fuerte.
# Correlación entre obesidad y calorías vegetales
cor.test(dataset2015$Vegetal.kcal, dataset2015$Media.obes, method = "pearson")##
## Pearson's product-moment correlation
##
## data: dataset2015$Vegetal.kcal and dataset2015$Media.obes
## t = 0.82428, df = 23, p-value = 0.4182
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2419330 0.5291078
## sample estimates:
## cor
## 0.1693909# Correlación entre obesidad y calorías animales
cor.test(dataset2015$Animal.kcal, dataset2015$Media.obes, method = "pearson")##
## Pearson's product-moment correlation
##
## data: dataset2015$Animal.kcal and dataset2015$Media.obes
## t = -1.5734, df = 23, p-value = 0.1293
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.62934448 0.09511309
## sample estimates:
## cor
## -0.3117334#REPRESENTACIONES
# Países de Europa a dibujar (no se tiene info de todos)
europe_countries <- c("Albania", "Austria", "Belgium", "Bulgaria", "Croatia", "Cyprus", "Czech Republic",
"Denmark", "Estonia", "Finland", "France", "Germany", "Greece", "Hungary",
"Iceland", "Ireland", "Italy", "Latvia", "Lithuania", "Luxembourg", "Malta",
"Netherlands", "Norway", "Poland", "Portugal", "Romania", "Slovakia",
"Slovenia", "Spain", "Sweden", "Switzerland", "UK", "Turkey")
# Se cambian nombres para que coincidan con map_data("world")
dataset2015_europa <- dataset2015 %>%
mutate(
Pais = case_when(
Pais == "Czechia" ~ "Czech Republic",
Pais == "Russia" ~ "Russian Federation",
Pais == "United Kingdom" ~ "UK",
TRUE ~ Pais)) %>%
filter(Pais %in% europe_countries)
europe_map <- map_data("world") %>% filter(region %in% europe_countries)
europe_data <- europe_map %>%
left_join(dataset2015_europa, by = c("region" = "Pais"))
ggplot(europe_data, aes(x = long, y = lat, group = group, fill = Media.obes)) +
geom_polygon(color = "white", size = 0.2) +
scale_fill_gradient(
name = "Obesidad (%)",
low = "#f7fcb9",
high = "#d95f0e",
na.value = "gray80") +
labs(
title = "Prevalencia de Obesidad en Europa (2015)",
subtitle = "Datos para ambos sexos",
caption = "Fuente: FAO/OMS") +
theme_void() +
theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
legend.position = "right", legend.title = element_text(size = 10)) + coord_fixed(1.3)## 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.
paises_europa <- c(
"Albania", "Austria", "Belgium", "Bulgaria", "Croatia", "Cyprus", "Czech Republic",
"Denmark", "Estonia", "Finland", "France", "Germany", "Greece", "Hungary",
"Iceland", "Ireland", "Italy", "Latvia", "Lithuania", "Luxembourg", "Malta",
"Netherlands", "Norway", "Poland", "Portugal", "Romania", "Slovakia",
"Slovenia", "Spain", "Sweden", "Switzerland", "UK", "Turkey")
dataset_europe <- dataset %>%
mutate(
Pais = case_when(
Pais == "Czechia" ~ "Czech Republic",
Pais == "Russia" ~ "Russian Federation",
TRUE ~ Pais ) ) %>%
filter(Pais %in% paises_europa)
europe_map <- map_data("world") %>%
filter(region %in% paises_europa)
europe_clusters <- europe_map %>%
left_join(dataset_europe %>% distinct(Pais, cluster_kmeans),
by = c("region" = "Pais"))
ggplot(europe_clusters, aes(x = long, y = lat, group = group, fill = factor(cluster_kmeans))) +
geom_polygon(color = "white", size = 0.2) + scale_fill_manual( name = "Cluster",
values = c("1" = "#4e79a7", "2" = "#e15759"),
labels = c("Cluster 1", "Cluster 2"),
na.value = "gray90") +labs( title = "Países de Europa por Cluster",
caption = "Clusters basados en IDH, Calorías y Obesidad") +
theme_void() + theme( plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
legend.position = "right") + coord_fixed(1.3) 
Finalmente se va crear un modelo para predecir el IDH en base a las
variables que más relación tienen con él. Usando Esperanza.de.vida,
Mortalidad.adulta, Delgadez.1.19.años., Media.obes y Animal.kcal.
datos_modelo <- dataset[, c("IDH", "Esperanza.de.vida", "Mortalidad.adulta", "Delgadez.1.19.años.", "Animal.kcal")]
cor_matrix <- cor(datos_modelo, method = "spearman")
round(cor_matrix, 2)## IDH Esperanza.de.vida Mortalidad.adulta
## IDH 1.00 0.75 -0.42
## Esperanza.de.vida 0.75 1.00 -0.54
## Mortalidad.adulta -0.42 -0.54 1.00
## Delgadez.1.19.años. -0.50 -0.70 0.50
## Animal.kcal 0.66 0.41 -0.19
## Delgadez.1.19.años. Animal.kcal
## IDH -0.50 0.66
## Esperanza.de.vida -0.70 0.41
## Mortalidad.adulta 0.50 -0.19
## Delgadez.1.19.años. 1.00 -0.35
## Animal.kcal -0.35 1.00corrplot(cor_matrix, method = "number", type = "upper", mar = c(0, 0, 1, 0)) 
set.seed(123)
sample <- sample(c(TRUE, FALSE), nrow(datos_modelo), replace = TRUE, prob = c(0.8, 0.2))
train <- datos_modelo[sample, ]
test <- datos_modelo[!sample, ]
modelo <- lm(IDH ~ Esperanza.de.vida + Delgadez.1.19.años. +I(Esperanza.de.vida^2)+ Animal.kcal, data = train) #quito mortalidad
summary(modelo)##
## Call:
## lm(formula = IDH ~ Esperanza.de.vida + Delgadez.1.19.años. +
## I(Esperanza.de.vida^2) + Animal.kcal, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.098758 -0.015570 0.003588 0.019793 0.051760
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.763e+00 4.938e-01 -3.571 0.00041 ***
## Esperanza.de.vida 5.740e-02 1.225e-02 4.686 4.13e-06 ***
## Delgadez.1.19.años. -2.800e-03 3.019e-03 -0.927 0.35444
## I(Esperanza.de.vida^2) -3.241e-04 7.592e-05 -4.270 2.58e-05 ***
## Animal.kcal 1.058e-04 9.479e-06 11.167 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02664 on 320 degrees of freedom
## Multiple R-squared: 0.7005, Adjusted R-squared: 0.6967
## F-statistic: 187.1 on 4 and 320 DF, p-value: < 2.2e-16predecir <- predict(modelo, newdata = test)
PRESS <- sum((test$IDH - predecir)^2)
TSS <- sum((test$IDH - mean(train$IDH))^2)
Q2 <- 1 - (PRESS / TSS)
print(paste("Q2(bondad de predicción):", round(Q2, 3))) #Predice bien los datos, 64%## [1] "Q2(bondad de predicción): 0.645"RMSE <- sqrt(mean((test$IDH - predecir)^2))
print(paste("RMSE:", round(RMSE, 3)))## [1] "RMSE: 0.025"dataset$IDH_predicho <- predict(modelo, newdata = dataset)El modelo predice con poco error (0.025 respecto al IDH real) los valores de IDH R-squared: 0.6834, significa que el 70% de la varianza esta explicada por el resto de variables. Quitar obesidad