Introducción

Los métodos de reducción de dimensionalidad son técnicas estadísticas que mapean el conjunto de los datos a subespacios derivados del espacio original, de menor dimensión, que permiten hacer una descripción de los datos a un menor costo.

Estas técnicas cobran importancia ya que muchos algoritmos de diversos campos tales como análisis numérico, aprendizaje automático o minería de datos suelen degradar su rendimiento cuando se usan con datos de alta dimensionalidad. En los casos extemos, el algoritmo deja de ser útil para el fin que fue diseñado. La maldición de la dimensión, se refiere a los diversos fenómenos que surgen al analizar y organizar datos de espacios de múltiples dimensiones.

Descripción del dataset

Este conjunto de datos es el resultado de combinar dos conjuntos: uno compuesto exclusivamente por variables cuantitativas y otro por variables cualitativas. La unificación se lleva a cabo a través de la columna del nombre del país, la cual está estandarizada según las normas ISO.

Indices mundiales cuantitativos

Se recopilan 65 indices de 188 países del mundo. Las fuentes principales son:

Las variables de interés seleccionadas para el análisis son:

Es importante aclarar que todas estas variables son cuantitativas.

Categorías mundiales cualitativas

Este es un conjunto de datos alojado por el Banco Mundial. Esta base de datos es mantenida por el Gloabl Poverty Working Group (GPWG), un equipo de expertos en medición de la pobreza del Poverty Reduction and Equity Network, el Development Research Group y el Development Data Group. La organización tiene una plataforma de datos abierta en https://databank.worldbank.org/.

Además, se añadieron manualmente otras 27 variables, las cuales se construyen con diversas fuentes de información en internet.

Si bien este dataset no son todas las categorías disponibles en la base de datos, es un resumen útil, constituido por 169 observaciones en 27 variables, las cuales son:

Preguntas de interés

Definiciones de interés

Los países desarrollados son naciones que han alcanzado niveles elevados de desarrollo humano, caracterizados por estándares de vida elevados y un crecimiento sostenido en las áreas económica, humanística e industrial. En contraste, el término “países en vías de desarrollo” se refiere al progreso económico de un país, afectando diversos aspectos como los políticos y sociales. La economía de estas naciones se encuentra en un estado de transición entre el subdesarrollo y las economías plenamente desarrolladas. El subdesarrollo, por otro lado, es una clasificación para naciones que carecen de los bienes, servicios y mecanismos productivos necesarios para generar su riqueza de manera sostenible. Aunque comúnmente se asocia con la pobreza, implica deficiencias más amplias, como calidad de vida, igualdad social e independencia financiera. No hay consenso sobre el criterio utilizado para distinguir entre países desarrollados, “en vías de desarrollo” o “subdesarrollados”, pero el término denota un rezago en la construcción social, política y económica, con implicaciones negativas para el país.

Carga de datos

Para comenzar, se requiere la carga de los datos alojados en una carpeta designada. Cada conjunto de datos será cargado mediante una sentencia particular y asignado a una variable específica. Posteriormente, se llevará a cabo su combinación.

# Cargamos la librería
library(readr)
## Warning: package 'readr' was built under R version 4.2.3
# Escribimos el directorio en donde está el archivo de las variables numéricas
dir1 <- paste0(dirname(rstudioapi::getActiveDocumentContext()$path),"/Dataset/World_Indexes.csv")
Data_1 <- read.csv(dir1)

# Escribimos el directorio en donde está el archivo de las variables categóricas
dir2 <- paste0(dirname(rstudioapi::getActiveDocumentContext()$path),"/Dataset/World_Categories.csv")
Data_2 <- read.csv2(dir2, na="")

Es correcto convertir todas las variables del dataset de categorías en factores, ya que esto permite un análisis omitiendo ambiguedad en, por ejemplo, las variables binarias.

# Convertir Data_2 a caracteres
Data_2 <- as.data.frame(as.matrix(Data_2))

# Combinar ambos dataframes
Data <- merge(Data_1, Data_2, by.x="Id", by.y="Name")

El resultado final de combinar ambos conjuntos de datos da como resultado un conjunto combinado de 157 observaciones con 91 variables. La disparidad en tamaños se debe a la falta de disponibilidad de todas las variables de interés para todos los países.

Separación de sub-datasets

Luego, procederemos a dividir este marco de datos en tres subconjuntos, cada uno conformado por variables específicas destinadas a técnicas distintas. Esta división se realiza tomando en cuenta las variables que contienen cadenas de texto (categóricas), aquellas variables numéricas (cuantitativas), y las variables lógicas o binarias (categóricas dicotómicas).

# Librería
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Filtrar variables categóricas
Data_categoric <- Data %>%
  select(where(is.character))

# Filtrar variables booleanas
Data_logic <- Data_categoric[,18:29] %>%
  mutate_all(as.factor)

# Eliminar las columnas binarias de el subconjunto categórico
Data_categoric <- Data_categoric[,-c(18:29)]

# Filtrar variables numéricas
Data_numeric <- Data %>%
  select(where(is.numeric))

Análisis Exploratorio

Un primer acercamiento a los datos podría ser el intento de describir el comportamiento de los mismos, mostrando sus valores máximos, mínimos, promedios, etc.

# Descripción de las primeras 10 variables.
summary(Data_numeric)[,1:10]
##  Human.Development.Index.HDI.2014 Gini.coefficient.2005.2013
##  Min.   :0.3483                   Min.   :24.82             
##  1st Qu.:0.5384                   1st Qu.:32.35             
##  Median :0.7150                   Median :37.01             
##  Mean   :0.6810                   Mean   :38.75             
##  3rd Qu.:0.7979                   3rd Qu.:43.30             
##  Max.   :0.9439                   Max.   :65.77             
##  Adolescent.birth.rate.15.19.per.100k.20102015
##  Min.   :  0.617                              
##  1st Qu.: 15.297                              
##  Median : 41.601                              
##  Mean   : 51.247                              
##  3rd Qu.: 75.979                              
##  Max.   :204.789                              
##  Birth.registration.funder.age.5.2005.2013
##  Min.   :  2.00                           
##  1st Qu.: 75.00                           
##  Median : 97.00                           
##  Mean   : 82.92                           
##  3rd Qu.:100.00                           
##  Max.   :100.00                           
##  Carbon.dioxide.emissionsAverage.annual.growth
##  Min.   :-4.8584                              
##  1st Qu.: 0.4147                              
##  Median : 1.5399                              
##  Mean   : 1.8266                              
##  3rd Qu.: 3.2532                              
##  Max.   :14.8995                              
##  Carbon.dioxide.emissions.per.capita.2011.Tones
##  Min.   : 0.02191                              
##  1st Qu.: 0.58949                              
##  Median : 2.25339                              
##  Mean   : 4.00040                              
##  3rd Qu.: 6.36755                              
##  Max.   :37.18764                              
##  Change.forest.percentable.1900.to.2012 Change.mobile.usage.2009.2014
##  Min.   :-81.6667                       Min.   :-21.97               
##  1st Qu.:-13.2982                       1st Qu.: 11.04               
##  Median :  0.0000                       Median : 35.20               
##  Mean   : -0.7578                       Mean   : 61.28               
##  3rd Qu.:  5.2136                       3rd Qu.: 78.81               
##  Max.   :264.3678                       Max.   :577.24               
##  Consumer.price.index.2013 Domestic.credit.provided.by.financial.sector.2013
##  Min.   : 98.24            Min.   :-27.16                                   
##  1st Qu.:107.95            1st Qu.: 31.49                                   
##  Median :113.65            Median : 51.05                                   
##  Mean   :119.48            Mean   : 70.75                                   
##  3rd Qu.:121.96            3rd Qu.: 86.96                                   
##  Max.   :288.65            Max.   :366.53

En este caso vemos que visualizar un resúmen estadístico de las 65 variables numéricas sería una tarea muy engorrosa.

Además, podemos presentar si es que las variables del dataset presentan datos faltantes. Aunque esta tarea también se vuelve complicada de visualizar, presentando solo los valores que sí presentan valores faltantes:

# Mostrar columnas con valores faltantes
colSums(is.na(Data))[colSums(is.na(Data))>0]
##                       Regime.type           First.official.language 
##                                12                                 2 
##                Principal.Religion                  Lending.category 
##                                 2                                33 
##       System.of.National.Accounts Balance.of.Payments.Manual.in.use 
##                                 2                                 1 
##                   System.of.trade     Government.Accounting.concept 
##                                 3                                28

Lo más importante a concluir de este pequeño análisis exploratorio es que, dado que el contenido del dataset es grande, y las variables son demasiadas como para ser analizadas una por una, lo más conveniente sería optar por una técnica descriptiva que intente reducir la dimensionalidad del conjunto de datos, y así poder caracterizar los individuos por estas nuevas variables.

Región

La variable Región presenta, utilizando 7 categorías, la región a la cual pertenece el país registrado. Estas categorías son las siguientes:

summary(as.factor(Data_categoric$Region))
##        East Asia & Pacific      Europe & Central Asia 
##                         22                         46 
##  Latin America & Caribbean Middle East & North Africa 
##                         24                         13 
##              North America                 South Asia 
##                          2                          8 
##         Sub-Saharan Africa 
##                         42

Grupo de Ingreso

La variable Income.Group utiliza cuatro categorías distintas para clasificar los ingresos del país. Estos pueden ser ingresos altos, medio-altos, medio-bajos, o bajos.

summary(as.factor(Data_categoric$Income.Group))
##         High income          Low income Lower middle income Upper middle income 
##                  39                  28                  45                  45

Índice de Desarrollo Humano

Quizá una de las variables más importantes del dataset es el Índice de Desarrollo Humano (IDH). El IDH constituye un paradigma para la observación del progreso o regresión del desarrollo humano de un área geográfica determinada. Implica una forma de examinar el avance de una comunidad y, por consiguiente, evaluar los resultados de las políticas públicas implementadas.

El Índice está compuesto por tres dimensiones: salud, educación y economía. Tradicionalmente, la dimensión salud se calcula a partir de la esperanza de vida; la dimensión de la educación a partir de tasa de alfabetización de los adultos y la tasa bruta de matriculación combinada de educación primaria, secundaria y terciaria; finalmente, la dimensión economía utiliza como insumo el Producto Bruto Interno (PBI).

Esta variable en el dataset se distribuye de la siguiente manera:

summary(Data_numeric$Human.Development.Index.HDI.2014)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3483  0.5384  0.7150  0.6810  0.7979  0.9439

En este resultado observamos una diferencia entre la media y la mediana, lo cual indica dispersión en los datos.

# Librerías
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
# Gráfico
ggplot(Data_numeric, aes(Human.Development.Index.HDI.2014)) +
  geom_histogram(bins=9, color="black", fill="pink") +
  labs(x="Índice de Desarrollo Humano (2014)", y="Frecuencia", title="Histograma del IDH (2014)")

Podemos observar que una gran porción de países tienen un IDH superior al 0.7. Un boxplot de estos datos es el siguiente, utilizando como variable auxiliar la región y la categoría de ingresos de cada país.

# Librerías
library(ggplot2)

# Gráfico
ggplot(Data, aes(x=as.factor(Region), y=Human.Development.Index.HDI.2014)) +
  geom_boxplot() +
  labs(x="Región", y="Índice de Desarrollo Humano", title="IDH separado por región") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

En este gráfico observamos que países de la región del áfrica del sub-Sahara tienen una media de IDH menor a 0.5. En contraste, las regiones de Asia del Este, Latinoamérica, y el Medio Oeste presentan una media de IDH alrededor de 0.7. Tanto Europa como NorteAmérica tienen una mediana de IDH cercana a 0.9 puntos.

# Librerías
library(ggplot2)

# Gráfico
ggplot(Data, aes(x=as.factor(Income.Group), y=Human.Development.Index.HDI.2014)) +
  geom_boxplot() +
  labs(x="Grupo de Ingresos", y="Índice de Desarrollo Humano", title="IDH separado por ingresos") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Como se describió anteriormente, este gráfico apoya la observación de que un país con alto índice de desarrollo humano implica que es un país con una economía desarrollada.

Variables numéricas

Todas las variables numéricas del dataset original van a ser utilizadas en un Análisis de Componentes Principales, el cual es un método estadístico que permite simplificar la complejidad de espacios muestrales con muchas dimensiones a la vez que conserva su información. Por ejemplo, al suponer que existe una muestra con \(n\) individuos cada uno con \(p\) variables \((X_1,X_2,...,X_p)\), es decir, el espacio muestral tiene \(p\) dimensiones. El Análisis de Componentes Principales (PCA) permite encontrar un número de factores subyacentes (\(z<p\)) que explican aproximadamente lo mismo que las \(p\) variables originales. Donde antes se necesitaban \(p\) valores para caracterizar a cada individuo, ahora bastan \(z\) variables. Cada una de estas \(z\) nuevas variables recibe el nombre de componente principal.

El PCA basa su estrategia en la existencia de las denominadas variables latentes, que son variables complicadas de medir pero que capturan un fenómeno subyacente del sistema que se está investigando, y se pueden utilizar estas variables en lugar de las originales debido a su correlación. Un estado determinado de un individuo en el sistema podría considerarse como una variable latente pero, sin embargo, no existe una única medida que determina este estado, sino que existe una combinación de medidas que podrían llegar a describir con cierto nivel de precisión este estado.

Análisis de Correlaciones

Técnicamente, el PCA se fundamenta en la búsqueda de una proyección que optimice la representación de datos mediante mínimos cuadrados. Este proceso transforma un conjunto de observaciones de variables, potencialmente correlacionadas, en un nuevo conjunto de valores de variables no correlacionadas linealmente, conocidas como componentes principales.

El primer componente principal en un conjunto de \(p\) variables se define como una variable derivada que constituye una combinación lineal de las variables originales, explicando la mayor parte de la varianza total de los datos. El segundo componente principal se define como la variable que explica la máxima varianza una vez que se ha eliminado el efecto del primer componente. Este procedimiento se repite a lo largo de \(p\) iteraciones hasta que se logre explicar la totalidad de la varianza. El PCA encuentra su aplicación más común en situaciones donde existe una alta correlación entre las variables y se busca reducir su número a un conjunto independiente.

Estas correlaciones mencionadas se pueden presentar claramente en un mapa de calor de correlaciones (heatmap) entre las variables

# Librería
library(ggplot2)
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.2.3
library(RColorBrewer)

# Modificar los datos numéricos
Data_numeric_alt <- scale(Data_numeric)
colnames(Data_numeric_alt) <- paste0("V", 1:65)

# Correlación
Data_cor <- round(cor(Data_numeric_alt),2) # redondeo a 2 decimales

# Fundir los datos
Data_cor <- melt(Data_cor)

#Mapa de Calor
ggplot(Data_cor, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile() +
  scale_fill_gradientn(name= "Correlación", colors=brewer.pal(n=10, name="RdYlBu"), limits=c(-1,1)) +
  geom_text(aes(Var2, Var1, label = value), 
          color = "white", size = 2) +
  theme(axis.text.x = element_text(angle=90, vjust=0.5, hjust=1))

A pesar que es cierto que es complicado observar un gráfico que representa \(65 \times 65\) correlaciones, podemos observar claramente que no todas las variables están correlacionadas, hay claras franjas amarillas o blancas que representan una correlación prácticamente nula (cercana a cero). Podríamos suponer que estas variables, al no estar correlacionadas con ninguna otra, aportan nada más que ruido blanco al conjunto de datos, e incluso podrían entorpecer el análisis.

De todos modos, procedemos a realizar el PCA.

Análisis de Componentes Principales

En primer lugar, el cálculo de los componentes principales depende de las unidades de medida empleadas en las variables. Es por lo tanto importante, antes de aplicar el PCA, estandarizar las variables para que tengan media 0 y desviación estándar 1, ya que, de lo contrario, las variables con mayor varianza dominarían al resto. La estandarización se lleva a cabo restando a cada observación la media y dividiendo entre la desviación estándar de la variable a la que pertenece: \[ \frac{x_i - media(x)}{sd(x)} \]

De todos modos, esta operación se puede realizar fácilmente con el argumento scale.unit=TRUE, de la función PCA (para realizar un análisis de componentes principales), del paquete FactoMineR.

# Librería
library(FactoMineR)
## Warning: package 'FactoMineR' was built under R version 4.2.3
Data_PCA <- PCA(Data_numeric, scale.unit=TRUE, graph=FALSE)

Como se mencionó previamente, los componentes principales representan combinaciones lineales normalizadas de las variables originales dentro de un conjunto de datos. Al realizar el cálculo utilizando variables estandarizadas, los componentes principales se convierten en autovectores obtenidos de la matriz de correlaciones, donde los elementos diagonales tienen un valor de 1. En este contexto, no se utilizan las covarianzas, ya que al estandarizar las variables, ambas matrices se vuelven equivalentes.

Teóricamente, es factible obtener un número de componentes principales igual al número de variables disponibles en el dataset. Sin embargo, la elección de estos componentes se basa en que el primero represente la mayor varianza presente en los datos. El segundo componente captura la variabilidad restante que no está explicada por el primero, y así sucesivamente. La determinación de la cantidad de componentes a retener se realiza considerando un porcentaje significativo de la varianza total, permitiendo una explicación visual a través de gráficos que respalde dicha elección.

El objetivo es identificar las combinaciones lineales que mejor representan las variables \(X_1, ..., X_p\). Sean \((Z_1,Z_2,...,Z_M)\) \(M<p\) combinaciones lineales de las \(p\) variables originales, es decir \[ Z_m = \sum_{j=1}^{p} \phi_{jm} X_j \]

donde \(phi_{1m}, \phi_{2m}, ...\) son las cargas o “loadings” de los componentes principales. Los loadings dan idea sobre qué peso tiene cada variable en cada componente. Cada vector de loadings, de longitud igual a \(p\), define además la dirección en el espacio sobre el cual la varianza de los datos es mayor.

La primera componente principal (\(Z_1\)) es aquella cuya dirección refleja o contiene la mayor variabilidad en los datos (por lo que esta componente será la que más información contenga). Este vector define la línea lo más próxima posible a los datos y que minimiza la suma de las distancias perpendiculares entre cada dato y la línea representada por la componente (usando como medida el promedio de la distancia euclídea al cuadrado). \[ z_{i1} = \phi_{11} x_{i1} + \phi_{21} x_{i2} + ... + \phi_{p1} x_{ip} \]

donde \(\phi_{11}\) corresponde al primer loading de la primera componente principal.

El vector de loadings de la primera componente principal resuelve el problema de optimización \[ \text{maximizar}_{\phi_{11}, ..., \phi_{p1}} \Biggl\{ \frac{1}{n} \sum_{i=1}^{n} \biggl( \sum_{j=1}^{p} \phi_{j1} x_{ij} \biggr)^2 \Biggr\} \qquad \text{sujeto a} \qquad \sum_{j=1}^{p} \phi_{j1}^2 = 1 \]

La segunda componente principal (\(Z_2\)) será una combinación lineal de las variables, que recoja la segunda dirección con mayor varianza de los datos, pero que no esté correlacionada con \(Z_1\). Esta condición es equivalente a decir que la dirección de \(Z_2\) (vector de \(\phi_2\)) ha de ser perpendicular u ortogonal respecto a \(Z_1\) (vector \(\phi_1\)).

Posteriores componentes principales se calculan iterativamente de la misma forma.

Una vez generados los componentes principales, se pueden representar uno frente a otro para obtener visualizaciones de los datos.

# Eigenvalues
round(Data_PCA$eig[1:10,], 3)
##         eigenvalue percentage of variance cumulative percentage of variance
## comp 1      21.720                 33.415                            33.415
## comp 2       3.616                  5.563                            38.978
## comp 3       2.863                  4.405                            43.383
## comp 4       2.550                  3.922                            47.305
## comp 5       2.392                  3.681                            50.986
## comp 6       2.222                  3.418                            54.404
## comp 7       1.990                  3.061                            57.465
## comp 8       1.824                  2.806                            60.271
## comp 9       1.528                  2.351                            62.622
## comp 10      1.469                  2.260                            64.882

En nuestro caso observamos que la primera componente principal explica el 33.4% de la varianza, la segunda componente principal explica el 5.6% de la varianza total, y en conjunto estas dos representan el 38.98% de la varianza total. Sumando una tercer componente principal (4.4% de la varianza), se alcanza a explicar tan solo un 43.4% de la varianza total.

Este resultado indica que nuestro dataset contiene mucho ruido que, si bien no aporta a la construcción de los componentes principales, sí aporta a la varianza total, y por ende los autovalores disminuirán de manera lenta y amortiguada.

# Librería
library(factoextra)
## Warning: package 'factoextra' was built under R version 4.2.3
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
# Gráfico de scree plot
fviz_eig(Data_PCA, addlabels = TRUE) + labs(title="Scree plot de los autovalores del PCA",
                                            x="Dimensiones",
                                            y="Porcentaje de varianza explicada")

La conclusión de este primer intento es que la presencia de numerosas variables contribuye principalmente con ruido a los datos, sugiriendo la posibilidad de reducir el análisis a variables con correlación significativa.

Reducción de las variables correlacionadas

Como se sugirió anteriormente, es posible reducir el conjunto de variables, que en total son 65. Para realizar esta tarea nos podemos basar en la suposición descrita en un principio, que un PCA es idóneo cuando las variables presentan alta correlación entre sí. La reducción de las variables puede ser realizada teniendo en cuenta la suma de las correlaciones por filas o columnas (es decir, la correlación total de cada variable con todas las otras).

# Armar matriz de correlación
Correlation <- round(cor(Data_numeric),2)
colnames(Correlation) <- rownames(Correlation) <- paste0("V", 1:65)

# Presentar las sumas por columnas o por filas
hist(rowSums(abs(Correlation)), xlab="Suma de correlaciones", ylab="Frecuencia", main="Histograma de la correlación total")

# Presentar la media y mediana de las correlaciones totales absolutas
mean(rowSums(abs(Correlation)))
## [1] 17.10585
median(rowSums(abs(Correlation)))
## [1] 14.95

Sabiendo que en un caso de correlación ideal, la suma de las correlaciones por filas o columnas sería de 65, la correlación absoluta total media que presentan las variables es de 17.11. Podemos intentar reducir las variables a solo aquellas que presenten una correlación absoluta total mayor a 15.

# Filtrar los datos según una correlación baja
Data_numeric_fixed <- Data_numeric[, which(colSums(abs(Correlation)) > 15)]

Presentando nuevamente la matriz de correlación con un heatmap:

# Librería
library(ggplot2)
library(reshape2)
library(RColorBrewer)

# Modificar los datos numéricos
Data_numeric_alt <- scale(Data_numeric_fixed)
colnames(Data_numeric_alt) <- paste0("V", 1:dim(Data_numeric_alt)[2])

# Correlación
Data_cor <- round(cor(Data_numeric_alt),2) # redondeo a 2 decimales

# Fundir los datos
Data_cor <- melt(Data_cor)

#Mapa de Calor
ggplot(Data_cor, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile() +
  scale_fill_gradientn(colors=brewer.pal(n=10, name="RdYlBu"), limits=c(-1,1)) +
  geom_text(aes(Var2, Var1, label = value), 
          color = "white", size = 2) +
  theme(axis.text.x = element_text(angle=90, vjust=0.5, hjust=1)) +
  labs(title="Heatmap de las correlaciones", x="Variable 1", y="Variable 2")

Ahora vemos que todas las variables presentan alta correlación, y esto es muy útil para el análisis de componentes principales.

dim(Data_numeric_fixed)
## [1] 157  32

Además nuestro espacio de trabajo se redujo a un total de 32 variables. Podemos realizar nuevamente un PCA con este conjunto de variables.

Nuevo Análisis de Componentes Principales

Se realiza un nuevo PCA con un dataset arreglado, obviando todas las variables que no aportan gran correlación al sistema.

# Librería
library(FactoMineR)

# PCA
Data_PCA <- PCA(Data_numeric_fixed, scale.unit=TRUE, graph=FALSE)

# Modificar nombres de las variables
rownames(Data_PCA$var$coord) <- paste0("V", 1:dim(Data_PCA$var$contrib)[1])
round(Data_PCA$eig[1:10,], 3)
##         eigenvalue percentage of variance cumulative percentage of variance
## comp 1      19.699                 61.559                            61.559
## comp 2       2.189                  6.840                            68.398
## comp 3       1.254                  3.917                            72.316
## comp 4       1.041                  3.253                            75.568
## comp 5       0.849                  2.653                            78.221
## comp 6       0.754                  2.356                            80.577
## comp 7       0.681                  2.129                            82.706
## comp 8       0.610                  1.907                            84.613
## comp 9       0.544                  1.700                            86.313
## comp 10      0.470                  1.468                            87.781

La primera componente principal representa el 61.6% de la varianza total, mientras que la segunda representa un 6.8%. Vemos que dos componentes principales representan el 68.4% de la variabilidad total de los datos, lo cual indica que es correcto continuar con este PCA, decidiendo en un principio trabajar con dos componentes.

# Librería
library(factoextra)

# Gráfico de scree plot
fviz_eig(Data_PCA, addlabels = TRUE) + labs(title="Scree plot de los autovalores del PCA",
                                            x="Dimensiones",
                                            y="Porcentaje de varianza explicada")

Observar el scree plot del nuevo PCA apoya la idea de utilizar solamente dos componentes, ya que existe una disminución muy lenta a partir de la tercer componente.

Análisis de las variables

En esta etapa podemos intentar agrupar las variables e intentar explicar cómo estas aportan individualmente a cada componente principal, para luego intentar explicar la posición de los individuos.

fviz_pca_var(Data_PCA,
             col.var = "cos2", # Color por la calidad de la representación
             gradient.cols = c("darkorchid4", "gold", "darkorange"),
             repel = TRUE
             )

# Librerías
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.2.3
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
# Contribución de las variables a la CP1
grafico_1 <- fviz_contrib(Data_PCA, choice = "var", axes = 1, top = 20)

# Contribución de las variables a la CP2
grafico_2 <- fviz_contrib(Data_PCA, choice = "var", axes = 2, top = 20)

# Gráfico
grid.arrange(grafico_1, grafico_2, ncol=2, top='Contribución de las variables a las primeras tres componentes principales')

En un principio vemos que alrededor de 20 variables aportan significativamente a la construcción de la componente principal. Los nombres de esas variables son:

# Crear dataframe
Var_contrib <- as.data.frame(Data_PCA$var$contrib)

# Mostrar principales contribuciones a la primera dimensión
rownames(Var_contrib %>% arrange(desc(Dim.1)))[1:20]
##  [1] "Human.Development.Index.HDI.2014"                                   
##  [2] "Internet.users.percentage.of.population.2014"                       
##  [3] "Life.expectancy.at.birth..years"                                    
##  [4] "Mean.years.of.schooling...Years"                                    
##  [5] "Secondary.2008.2014"                                                
##  [6] "Infant.Mortality.2013.per.thousands"                                
##  [7] "Expected.years.of.schooling...Years"                                
##  [8] "Under.five.Mortality.2013.thousands"                                
##  [9] "Pupil.teacher.ratio.primary.school.pupils.per.teacher.2008.2014"    
## [10] "Gender.Inequality.Index.2014"                                       
## [11] "Electrification.rate.or.population"                                 
## [12] "Gross.national.income.GNI.per.capita...2011..Dollars"               
## [13] "Population.with.at.least.some.secondary.education.percent.2005.2013"
## [14] "Physicians.per.10k.people"                                          
## [15] "Tertiary..2008.2014"                                                
## [16] "Maternal.mortality.ratio.deaths.per.100.live.births.2013"           
## [17] "Primary.school.dropout.rate.2008.2014"                              
## [18] "Gross.domestic.product.GDP.percapta"                                
## [19] "Domestic.food.price.level.2009.2014.index"                          
## [20] "Adolescent.birth.rate.15.19.per.100k.20102015"

Es interesante observar que el índice de desarrollo humano, el porcentaje de usuarios en internet, la esperanza de vida al nacer, y el promedio de años de escolaridad, son las cuatro variabels que más aportan a la construcción de la primera componente. Siendo estas variables importantes a nivel mundial.

Por otra parte, para la segunda componente observamos que alrededor de 10 variables aportan a su construcción, estas son:

# Mostrar principales contribuciones a la segunda dimensión
rownames(Var_contrib %>% arrange(desc(Dim.2)))[1:10]
##  [1] "Research.and.development.expenditure..2005.2012"         
##  [2] "Gross.domestic.product.GDP.percapta"                     
##  [3] "Gross.national.income.GNI.per.capita...2011..Dollars"    
##  [4] "Renewable.sources.percentage.of.total.2012"              
##  [5] "Fossil.fuels.percentage.of.total.2012"                   
##  [6] "Maternal.mortality.ratio.deaths.per.100.live.births.2013"
##  [7] "Domestic.credit.provided.by.financial.sector.2013"       
##  [8] "Electrification.rate.or.population"                      
##  [9] "Birth.registration.funder.age.5.2005.2013"               
## [10] "Stock.of.immigrants.percentage.of.population.2013"

En este caso, la tres más importantes resultan ser el gasto relacionado a la investigación y desarrollo, el producto bruto doméstico (GDP) per capita, y renta nacional bruta (GNI) per capita. Estas tres variables son exclusivamente económicas.

Análisis de gradiente

Para definir en profundidad la explicación que aporta cada variable individualmente a las componentes, es correcto realizar un análisis de gradiente.

Partiendo de la lista de coordenadas de las variables, es necesario seleccionar el máximo de cada uno de los componentes, y multiplicarlo por \(\frac{2}{3}\), para luego mostrar únicamente los valores superiores a este resultado, indicando cuáles variables aportan a la construcción de cada uno de los dos componentes principales.

# PCA
Data_PCA <- PCA(Data_numeric_fixed, scale.unit=TRUE, graph=FALSE)

# Crear matriz de coordenadas
Coord_var_matrix <- round(Data_PCA$var$coord[,1:2],3)

# Identificar valores menores a 2/3 el máximo de la coordenada de cada componente
for (j in 1:2){
  for (i in 1:dim(Coord_var_matrix)[1]){
    if (abs(Coord_var_matrix[i,j]) < max(abs(Coord_var_matrix[,j]))*(2/3)){
      Coord_var_matrix[i,j] = 0
    }
  }
}

# Presentar tabla
print.table(Coord_var_matrix, zero.print="")
##                                                                      Dim.1
## Human.Development.Index.HDI.2014                                     0.982
## Adolescent.birth.rate.15.19.per.100k.20102015                       -0.778
## Birth.registration.funder.age.5.2005.2013                            0.742
## Carbon.dioxide.emissions.per.capita.2011.Tones                       0.675
## Change.mobile.usage.2009.2014                                             
## Domestic.credit.provided.by.financial.sector.2013                         
## Domestic.food.price.level.2009.2014.index                           -0.804
## Electrification.rate.or.population                                   0.863
## Expected.years.of.schooling...Years                                  0.892
## Fossil.fuels.percentage.of.total.2012                                     
## Gender.Inequality.Index.2014                                        -0.864
## Gross.domestic.product.GDP.percapta                                  0.823
## Gross.national.income.GNI.per.capita...2011..Dollars                 0.849
## Infant.Mortality.2013.per.thousands                                 -0.892
## International.inbound.tourists.thausands.2013                             
## Internet.users.percentage.of.population.2014                         0.927
## Life.expectancy.at.birth..years                                      0.903
## Maternal.mortality.ratio.deaths.per.100.live.births.2013            -0.843
## Mean.years.of.schooling...Years                                      0.901
## Mobile.phone.subscriptions.per.100.people.2014                            
## Physicians.per.10k.people                                            0.844
## Population.with.at.least.some.secondary.education.percent.2005.2013  0.846
## Pre.primary.2008.2014                                                0.724
## Primary.school.dropout.rate.2008.2014                               -0.826
## Pupil.teacher.ratio.primary.school.pupils.per.teacher.2008.2014     -0.871
## Renewable.sources.percentage.of.total.2012                          -0.679
## Research.and.development.expenditure..2005.2012                           
## Secondary.2008.2014                                                  0.898
## Stock.of.immigrants.percentage.of.population.2013                         
## Tertiary..2008.2014                                                  0.844
## Tuberculosis.rate.per.thousands.2012                                      
## Under.five.Mortality.2013.thousands                                 -0.875
##                                                                      Dim.2
## Human.Development.Index.HDI.2014                                          
## Adolescent.birth.rate.15.19.per.100k.20102015                             
## Birth.registration.funder.age.5.2005.2013                                 
## Carbon.dioxide.emissions.per.capita.2011.Tones                            
## Change.mobile.usage.2009.2014                                             
## Domestic.credit.provided.by.financial.sector.2013                         
## Domestic.food.price.level.2009.2014.index                                 
## Electrification.rate.or.population                                        
## Expected.years.of.schooling...Years                                       
## Fossil.fuels.percentage.of.total.2012                                     
## Gender.Inequality.Index.2014                                              
## Gross.domestic.product.GDP.percapta                                  0.441
## Gross.national.income.GNI.per.capita...2011..Dollars                 0.440
## Infant.Mortality.2013.per.thousands                                       
## International.inbound.tourists.thausands.2013                             
## Internet.users.percentage.of.population.2014                              
## Life.expectancy.at.birth..years                                           
## Maternal.mortality.ratio.deaths.per.100.live.births.2013                  
## Mean.years.of.schooling...Years                                           
## Mobile.phone.subscriptions.per.100.people.2014                            
## Physicians.per.10k.people                                                 
## Population.with.at.least.some.secondary.education.percent.2005.2013       
## Pre.primary.2008.2014                                                     
## Primary.school.dropout.rate.2008.2014                                     
## Pupil.teacher.ratio.primary.school.pupils.per.teacher.2008.2014           
## Renewable.sources.percentage.of.total.2012                                
## Research.and.development.expenditure..2005.2012                      0.601
## Secondary.2008.2014                                                       
## Stock.of.immigrants.percentage.of.population.2013                         
## Tertiary..2008.2014                                                       
## Tuberculosis.rate.per.thousands.2012                                      
## Under.five.Mortality.2013.thousands

A partir de estas variables resaltadas podemos construir gradientes, de forma manual, como sigue:

Podemos observar que por ejemplo, las variables relacionadas principalmente al desarrollo humano y educación, como ser el IDH, los años de escolaridad, la posibilidad de tener internet, la cantidad de médicos, etc., son las variables que más aportan a los valores positivos en la construcción de la primer componente principal (ubicada en el eje horizontal). En menor medida se presentan variables no relacionadas, como el registro de nacimientos, las emisiones de dióxido de carbono, y las fuentes de energía renovables, también en los valores positivos del eje horizontal de la primera componente principal. En oposición se presentan variables consideradas como “malas”, es decir, índices que representan factores negativos con su aumento, como ser la tasa de embarazo adolescente, tasas de mortalidad, desigualdad de género, deserción escolar, etc. Valores que se encuentren más hacia la izquierda implican que tienen valores grandes de estos índices.

En la componente principal 2, representada en el eje vertical, está presente una única variable significativamente, esta es el gasto relacionado con la investigación y el desarrollo. Por otra parte, en dirección diagonal hacia el extremo superior derecho existen dos variables exclusivamente económicas, que son el Producto Bruto Interno (PBI) y el Ingrecio Bruto Nacional (IBN). Estas dos últimas variables son de gran importancia en el ámbito económico, y son índices desarrollados para intentar explicar la posición económica de los Estados. De todos modos, estas dos variables tienen más presencia en el eje horizontal (CP 1) que en el eje vertical (CP 2).

Descripción de los individuos

El siguiente paso del análisis es intentar identificar los individuos, observar su distribución entre las dos variables, e intentar inferir sobre posibles grupos formados. El gráfico de individuos es el siguiente.

# Librerías
library(factoextra)
library(gridExtra)
library(RColorBrewer)
library(ggforce)
## Warning: package 'ggforce' was built under R version 4.2.3
library(ggrepel)

# Paleta
Paleta <- brewer.pal(7, "Set1")

# Filtrar labels
Labels_filt <- cbind(Data, Data_PCA$ind$coord) %>%
  select(Id, Dim.1, Dim.2)

# Gráfico
fviz_pca_ind(Data_PCA, geom="point",
             axes=c(1,2)) +
  labs(title="PCA - Individuos") +
  geom_mark_hull(concavity=3, fill="grey", radius=unit(8,"mm"))

Este gráfico muestra la proyección de los individuos en el nuevo sistema de coordenadas, compuesto por las dos componentes principales obtenidas por el PCA. Una primera observación sería que existe una forma más o menos definida de cómo se distribuyen los datos. Esta forma es alargada sobre el eje de la primera componente principal, y tiene una marcada concavidad sobre el eje de la segunda dimensión, teniendo un mínimo sobre la mitad inferior del plano.

Es posible intentar agrupar los individuos a simple vista, teniendo en cuenta que esta observación no será del todo acertada, y es casi completamente subjetiva, pero es útil para intentar inferir sobre las observaciones de los individuos.

En este gráfico intentamos agrupar los individuos en tres grupos notables. Teniendo en cuenta los gradientes descritos anteriormente podemos inferir que, el grupo A estaría compuesto por países desarrollados, con altos índices económicos y de educación, además de una alta inversión en investigación y desarrollo. El grupo B estaŕia compuesto por países denominados emergentes, es decir países que tienen algunas características de un mercado desarrollado, como índices de desarrollo medios, pero que no cumplen estrictamente con los estándares para ser países desarrollados. Estos países pueden convertirse en mercados desarrollados en el futuro, o que lo fueron en el pasado. Finalmente, el grupo C estaría compuesto casi en su totalidad por países sub-desarrollados, es decir países que tienen no solo índices económicos bajos, sino también índices de desarrollo muy bajos.

A continuación podemos mostrar en el gráfico ciertos países, en este caso los pertenecientes al foro internacional de gobernantes y presidentes de bancos centrales G20, para intentar inferir a partir de la observación de las posiciones de estos países en el nuevo sistema de coordenadas, y tomando en cuenta su posición con los grupos generados de forma manual.

A continuación podemos mostrar en el gráfico ciertos países, en este caso es una muestra que intenta ser representativa de todos los países, y se compone de los dos países que presentan mayor Índice de Desarrollo Humano separado por región, los dos países que presentan menor índice de desarrollo humano separado por región, y los 20 países que componen el bloque intergubernamental económico G20. La elección de este indicador como base para la división se basa en la idea que el IDH es un índice que no solo engloba características de educación por ejemplo, sino que también incluye características económicas, como el ingreso bruto nacional.

# Librerías
library(factoextra)
library(gridExtra)
library(RColorBrewer)
library(ggforce)
library(ggrepel)

# Lista de países representativos por región
Representative_countries <- rbind(cbind(Data, Data_PCA$ind$coord) %>%
  group_by(Region) %>%
  slice_max(order_by=Human.Development.Index.HDI.2014, n=2) %>%
  select(Id, Dim.1, Dim.2, Region),
  cbind(Data, Data_PCA$ind$coord) %>%
  group_by(Region) %>%
  slice_min(order_by=Human.Development.Index.HDI.2014, n=2) %>%
  select(Id, Dim.1, Dim.2, Region),
  cbind(Data, Data_PCA$ind$coord) %>%
  mutate(Id = ifelse(G20==1, Id, NA)) %>%
  select(Id, Dim.1, Dim.2, Region) %>%
  na.omit()) %>%
  distinct()

# Paleta
Paleta <- brewer.pal(7, "Set1")

# Gráfico
fviz_pca_ind(Data_PCA, geom="point",
             axes=c(1,2), alpha=0.2) +
  labs(title="PCA - Individuos") +
  geom_mark_hull(concavity=3, fill="grey", radius=unit(8,"mm")) +
  geom_label_repel(data=Representative_countries, aes(x=Dim.1, y=Dim.2, label=Id), alpha=0.7) +
  geom_point(data=Representative_countries, aes(x=Dim.1, y=Dim.2), color="#FF5B22")

Este resultado nos permite realizar observaciones interesantes. Por ejemplo, Estados Unidos es el país que se encuentra en el extremo derecho y, teniendo en cuenta los gradientes descritos anteriormente, diríamos que este país tiene un alto índice de desarrollo y escolaridad, así como también una gran inversión en investigación, junto con índices económicos altos como el PBI y el INB. Países europeos como Suiza, Francia, Alemania, Italia, y el Reino Unido se encuentran también en una zona densa del extremo derecho, indicando observaciones similares a las de Estados Unidos, pero quizá en menor medida. También en este extremo se encuentran países como Australia, Corea del Sur, y Japón, que si bien son países orientales, pertenecen a grupos como la OTAN, y son denominados países desarrollados.

En una zona inferior pero aún en el área derecha de la distribución encontramos a países denominados “emergentes”, que son países con economías que pueden ser potencialmente importantes en una proyección a futuro, como ser Rusia, China, Brasil, México, Turquía, Argentina, Indonesia, Sudáfrica e India. Estos países se caracterizan por tener una economía regular, índices de desarrollo también medios-altos, pero poca capacidad de investigación y desarrollo.

Finalmente podemos identificar todo el brazo izquierdo, compuesto casi en su totalidad por países de medio oriente, países centro y sud africanos, y otros países asiáticos del sur y este. Estos países presentan muy bajos índices de desarrollo económico y humanos, además de altas tasas de mortalidad infantil, embarazo adolescentes, etc, que son características negativas para los países.

Por otra parte, es posible utilizar ciertas variables categóricas presentes en el dataset para poder añadir descripciones a la proyección de individuos en el nuevo sistema con las componentes principales.

# Librerías
library(factoextra)
library(gridExtra)
library(RColorBrewer)

# Filtrar datos
Data_filt <- Data %>%
  mutate(Major.Economic.Group = ifelse(Major.Economic.Group == "None", NA, Major.Economic.Group)) %>%
  select(Region, Income.Group, Major.Economic.Group, Type)
  

# Gráfico
Ind_region <- fviz_pca_ind(Data_PCA,
                           geom="point", axes=c(1,2),
                           col.ind=as.factor(Data_filt$Region),
                           addEllipses = TRUE,
                           ellipse.type="convex") +
  labs(title="PCA - Individuos agrupados por región")

Ind_grupo_ingreso <- fviz_pca_ind(Data_PCA,
                                  geom="point", axes=c(1,2),
                                  col.ind=as.factor(Data_filt$Income.Group),
                                  addEllipses = TRUE,
                                  ellipse.type="convex") +
  labs(title="PCA - Individuos agrupados por grupo de ingreso")

Ind_economic <- fviz_pca_ind(Data_PCA,
                           geom="point", axes=c(1,2),
                           col.ind=as.factor(Data_filt$Major.Economic.Group),
                           addEllipses = TRUE,
                           ellipse.type="convex") +
  labs(title="PCA - Individuos agrupados por grupo económico")

Ind_type <- fviz_pca_ind(Data_PCA,
                         geom="point", axes=c(1,2),
                         col.ind=as.factor(Data_filt$Type),
                         addEllipses = TRUE,
                         ellipse.type="convex") +
  labs(title="PCA - Individuos agrupados por tipo de estado")

grid.arrange(Ind_region, Ind_grupo_ingreso, Ind_economic, Ind_type, ncol=2, nrow=2,
             top="PCA - individuos agrupados por distintas variables categóricas")
## Warning: Removed 71 rows containing missing values (`geom_point()`).
## Warning: Removed 1 rows containing missing values (`geom_point()`).

El primero de estos gráficos, los individuos agrupados por región, muestra que los países africanos están desplazados hacia la izquierda, mientras que casi todos los países europeos y norte-americanos están hacia la derecha. Esto se explica claramente teniendo en cuenta el gradiente del desarrollo humano. En el centro del gráfico existen diversos grupos de países.

El gráfico agrupado por grupo de ingreso es mucho más conciso, y muestra claramente el porqué de las dimensiones. De izquierda a derecha observamos países de bajos ingresos (países sub-desarrollados), seguidos de países de ingresos medio-bajos, medio-altos (países emergentes), y finalmente hacia la derecha los países de ingresos altos (países desarrollados).

La agrupación por grupos económicos indica que hacia la derecha se encuentran grupos poderosos como la unión europea, el NAFTA (North American Free Trade Agreement) o el EFTA (European Free Trade Associatio). Mientras que hacia abajo y el centro existen grupos emergentes como el MERCOSUR, el ASEAN, o el CIS (Commonwealth of Independent States).

Finalmente, el cuarto gráfico permite realizar una comparativa entre las categorías de “País”, “Dependencia” y “País soberano”. Claramente este grupo no aporta demasiada información.

Variables dicotómicas

El siguiente paso del trabajo será explorar y analizar las variables lógicas o dicotómicas, es decir, aquellas que toman únicamente valores binarios como 1/0, VERDADERO/FALSO, SÍ/NO.

El objetivo general de esta sección es similar a la anterior, intentar observar agrupamientos de los individuos, pero ahora tomando como referencia las variables dicotómicas del dataset. Para cumplir este objetivo será útil emplear técnicas de escalado multidimensional como la generación de Coordenadas Principales, las cuales tienen una analogía con las Componentes Principales, pero utilizando otro tipo de distancias en vez de la euclídea.

Preparación de los datos

En primer lugar es necesario convertir los datos lógicos del dataset que figuran como caracteres, a valores binarios de 0 y 1, siendo 0 la negación y 1 la afirmación de la variable a tomar en cuenta.

# Transformar los datos lógicos de caracteres a números
Data_logic <- as.data.frame(apply(Data_logic, 2, function(x) as.numeric(x)))

Matriz de distancias

El siguiente paso es generar una matriz de distancias, la cual se construye a partir de coeficientes de similaridad, que miden el grado de similaridad entre cada par de individuos. Supongamos que se tienen \(p\) variables binarias \(X_1, X_2, ..., X_p\), donde cada \(X_i\) toma valores 0 o 1. Para cada par de individuos se genera una tabla de frecuencias de coincidencias entre pares de valores, donde \(a, b, c, d\) son las frecuencias de \((1,1)\), \((1,0)\), \((0,1)\), y \((0,0)\), respectivamente, con \(p = a + b + c + d\). Un coeficiente de similaridad debería ser función de \(a, b, c, d\). Son conocidos los coeficientes de similaridad: \[ \begin{gather} s_{ij} = \frac{a+d}{p} \qquad \text{Simple Matching}\\ s_{ij} = \frac{a}{a+b+c} \qquad \text{Jaccard} \end{gather} \]

Se puede transformar una disimilaridad en distancia aplicando la fórmula \(d = \sqrt{1 - s}\), con \(s\) como un coeficiente de similaridad, aunque esta transformación la realiza internamenta la función dist.binary del paquete ade4, y al elegir el método “2” se elige utilizar el Simple Matching como coeficiente de similaridad.

# Librerías
library(ade4)
## Warning: package 'ade4' was built under R version 4.2.3
## 
## Attaching package: 'ade4'
## The following object is masked from 'package:FactoMineR':
## 
##     reconst
# Matriz de distancias
Data_logic_dist <- round(dist.binary(Data_logic, method=2, diag=TRUE),2)

Análisis de Coordenadas Principales

Ahora es posible proceder a realizar un análisis de Coordenadas Principales, un método para representar gráficamente, de manera análoga a un PCA, objetos descritos por una matriz que contiene las distancias, construidas a partir de la similitud entre estos individuos. Este método proporciona coordenadas cartesianas compatibles con la matriz de distancias, a diferencia del PCA que utiliza la matriz de covarianzas/correlaciones.

#Librerías
library(factoextra)

# Escalado muldimensional
Data_logic_MDS <- cmdscale(Data_logic_dist, eig=TRUE, k=2)
Data_logic_eig <- data.frame(Eigenvalues=Data_logic_MDS$eig[1:20])

# Gráfico de scree plot para los eigenvalues
ggplot(Data_logic_eig, aes(x=1:20, y=Eigenvalues)) +
  geom_bar(stat="identity") +
  labs(title="PCoA - Scree Plot", x="", y="Autovalores")

A pesar de este resultado, es correcto seleccionar las dos primeras coordenadas principales. A continuación es posible graficar los individuos para realizar algunas observaciones e inferencias sobre la existencia de grupos de individuos. Es importante destacar que, a diferencia del PCA, en el PCoA, las coordenadas dejan de tener sentido cuando se refieren a las varianzas representadas por cada una.

Descripción de los individuos

# Librería
library(ggplot2)
library(ggrepel)

# Generar dataframe
Data_logic_MDS_coord <- data.frame(Name=Data$Id,
                                   Coord1=Data_logic_MDS$points[,1],
                                   Coord2=Data_logic_MDS$points[,2],
                                   Region=Data$Region)

Data_logic_MDS_coord_alt <- cbind(Data_logic_MDS_coord, Data$G20) %>%
  mutate(G20 = ifelse(Data$G20 == 1, 1, NA)) %>%
  na.omit()

# Gráfico
ggplot(Data_logic_MDS_coord, aes(x=Coord1, y=Coord2)) +
  geom_count() +
  scale_size(range=c(2,10)) +
  geom_text_repel(data=Data_logic_MDS_coord_alt, aes(x=Coord1, y=Coord2, label=Name), max.overlaps=20) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="PCoA - Individuos", x="Coordenada 1", y="Coordenada 2")

Una primera observación que podemos realizar es que los datos están muy dispersos alrededor del plano con las nuevas coordenadas generadas. A propósito se utilizó la función geom_count que realiza un conteo de los puntos superpuestos y los representa como puntos más grandes mientras más puntos se acumulen. Teniendo esto en cuenta vemos que al menos cinco puntos sobresalen por su tamaño, hacia la derecha, pero además existen varios puntos individuales dispersos en todo el plano del nuevo sistema de coordenadas.

Además a modo ilustrativo, se colocan los nombres de los países pertenecientes al grupo intergubernamental G20. Con estas indicaciones vemos que hacia la izquierda se encuentran muchos países europeos como el Reino Unido, Francia, Alemania, Italia, y además los Estados Unidos, Australia, Canadá, etc. Sabiendo que los países mencionados son denominados países desarrollados podemos suponer que progresivamente hacia la derecha irán apareciendo países subdesarrollados, e incluso países que no cuentan con los derechos mencionados en las variables tales como el matrimonio igualitario y una cobertura de salud universal. Además podríamos inferir que ciertos conglomerados surgen además por la pertenencia de los países a ciertos bloques económicos y/o militares, como el BRICS, la OTAN, etc, y además si estos países fueron colonias o si fueron imperios en algún momento de su historia. Por esa razón no es sorprendente observar la cercanía de países europeos como el Reino Unido y Francia, que prácticamente están superpuestos.

Utilizando variables auxiliares podemos generar grupos en el gráfico de individuos para observar su comportamiento. Por ejemplo utilizando la variable Región:

# Librería
library(ggplot2)
library(ggrepel)

# Generar dataframe
Data_logic_MDS_coord <- data.frame(Name=Data$Id,
                                   Coord1=Data_logic_MDS$points[,1],
                                   Coord2=Data_logic_MDS$points[,2],
                                   Region=Data$Region)

Data_logic_MDS_coord_alt <- cbind(Data_logic_MDS_coord, Data$G20) %>%
  mutate(G20 = ifelse(Data$G20 == 1, 1, NA)) %>%
  na.omit()

# Gráfico
ggplot(Data_logic_MDS_coord, aes(x=Coord1, y=Coord2, fill=Region)) +
  geom_count(aes(color=Region)) +
  geom_mark_hull(alpha=0.1, concavity=15) +
  scale_size(range=c(3,10)) +
  geom_label_repel(data=Data_logic_MDS_coord_alt, aes(x=Coord1, y=Coord2, label=Name), max.overlaps=20) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="PCoA - Individuos", x="Coordenada 1", y="Coordenada 2")

A pesar de que el gráfico es muy disperso, sí observamos que todas las áreas convergen en algunos de estos puntos grandes mencionados, y esto es porque hay muchos países en el dataset que no cumplen con ninguno de los requisitos impuestos por las variables, y es lógico que se acumulen en cierta región. Como se había comentado, hacia la derecha se observan principalmente países europeos occidentales, países de nortaeamérica, y algunos países orientales como Japón y Australia. Hacia el centro del gráfico observamos países mas bien emergentes, como Brasil, Argentina, México, Sudáfrica, y Turquía, con la excepción de Surcorea. Es importante recordar que en este caso las variables no son exclusivamente económicas, sino mas bien relacionadas con ciertas características de los países, y más frecuentemente su membrecía a bloques intergubernamentales.

También podemos utilizar la variable Income.Group, la cual diferencia a los países con su pertenencia a una de las cuatro categorías de ingresos. El gráfico es el siguiente.

# Librería
library(ggplot2)
library(ggrepel)

# Generar dataframe
Data_logic_MDS_coord <- data.frame(Name=Data$Id,
                                   Coord1=Data_logic_MDS$points[,1],
                                   Coord2=Data_logic_MDS$points[,2],
                                   Group=Data$Income.Group)

Data_logic_MDS_coord_alt <- cbind(Data_logic_MDS_coord, Data$G20) %>%
  mutate(G20 = ifelse(Data$G20 == 1, 1, NA)) %>%
  na.omit()

# Gráfico
ggplot(Data_logic_MDS_coord, aes(x=Coord1, y=Coord2, fill=Group)) +
  geom_count(aes(color=Group)) +
  geom_mark_hull(alpha=0.1, concavity=15) +
  scale_size(range=c(3,10)) +
  geom_label_repel(data=Data_logic_MDS_coord_alt, aes(x=Coord1, y=Coord2, label=Name), max.overlaps=20) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="PCoA - Individuos", x="Coordenada 1", y="Coordenada 2")

Este gráfico clarifica significativamente la observación del gradiente mencionado hacia la izquierda sobre los países desarrollados. En este caso es esperable que los individuos ubicados más hacia la derecha sean economías avanzadas, con poder nuclear y la obtención de ciertos derechos como el matrimonio igualitario y una cobertura de salud universal. Estos países además pertenecen a grupos económicos como el G20, G7, y a grupos militares como la OTAN, pero además es probable que estos países fueron imperios en algún momento de su historia. En una posición media se encuentran frecuentemente a economías emergentes. Y hacia la izquierda se encuentran países subdesarrollados que probablmente no participan de ningun bloque intergubernamental económico, y que no cuentan con los derechos mencionados anteriormente, pero es probable que muchos de estos países pertenezcan al grupo económico emergente denominado G33.

Variables cualitativas

A continuación, llevaremos a cabo un análisis detallado de las variables categóricas que presentan más de dos estados, a diferencia de las variables dicotómicas previamente examinadas. Para abordar este tipo de datos, emplearemos dos técnicas distintas. En primer lugar, aplicaremos un Análisis de Coordenadas Principales, similar al utilizado para las variables dicotómicas, con el objetivo de proyectar los individuos en un nuevo sistema de coordenadas generado a partir de la información categórica del conjunto de datos. Además, exploraremos la opción de realizar un Análisis Factorial de Correspondencia, tanto simple como múltiple, para comprender cómo interactúan los diversos factores presentes en las variables categóricas.

Modificación de los datos

En primer lugar es correcto modificar los datos de las variables cuando estas presentan muchos factores, pero es probable que estos tengan una cierta frecuencia. Lo correcto sería seleccionar únicamente los factores con frecuencias mayores, y englobar todo el resto en una categoría de “Otras”. Podemos identificar por ejemplo la variable Currency.Unit, la cual hace referencia al tipo de moneda utilizada en el país registrado, la cual presenta demasiados factores.

Data_categoric %>% count(Currency.Unit) %>%
  arrange(desc(n)) %>%
  top_n(5)
## Selecting by n
##               Currency.Unit  n
## 1                      Euro 20
## 2    West African CFA franc  8
## 3               U.S. dollar  7
## 4 Central African CFA franc  5
## 5         Australian dollar  2

En este caso deberíamos seleccionar simplemente los 6 primeros factores y a partir de este englobar todo el resto de factores con frecuencias menores en uno “Other”. El caso es el mismo con la variable First.official.language, la cual hace referencia al primer idioma oficial del país registrado. Los factores que presenta esta variable son:

Data_categoric %>% count(First.official.language) %>%
  arrange(desc(n)) %>%
  top_n(8)
## Selecting by n
##    First.official.language  n
## 1                  English 26
## 2                  Spanish 19
## 3                   French 16
## 4                   Arabic 12
## 5               Portuguese  4
## 6                    Dutch  3
## 7                   German  3
## 8                   Creole  2
## 9                    Greek  2
## 10                 Persian  2
## 11                 Russian  2
## 12                    <NA>  2

De la misma forma, deberíamos elegir simplemente los primeros 7 factores, y al resto englobarlos en un factor de “Other”. Estas tareas se realizan el el siguiente bloque de código.

# Modificar la variable del tipo de moneda
Main_currencies <- Data_categoric %>% count(Currency.Unit) %>%
  arrange(desc(n)) %>%
  top_n(5) %>%
  select(Currency.Unit)
## Selecting by n
Data_categoric <- Data_categoric %>%
  mutate(Currency.Unit = ifelse((Currency.Unit %in% Main_currencies[,1]) == TRUE, as.character(Currency.Unit), "Other"))

# Modificar la variable del primer idioma oficial
Main_languages <- Data_categoric %>% count(First.official.language) %>%
  arrange(desc(n)) %>%
  top_n(8) %>%
  select(First.official.language)
## Selecting by n
Data_categoric <- Data_categoric %>%
  mutate(First.official.language = ifelse((First.official.language %in% Main_languages[,1]) == TRUE, as.character(First.official.language), "Other"))

Data_categoric <- Data_categoric[-c(1:3)] %>%
  mutate_all(as.factor)
Data_categoric <- cbind(Data_categoric, Data_logic %>% mutate_all(function(x) ifelse(x==1, "Yes", "No"))) %>%
  mutate_all(as.factor)

Análisis de Coordenadas Principales

Tal como en el análisis de coordenadas principales para variables dicotómicas, se genera una matriz de distancias a partir de los datos. En este caso se utiliza la función daisy de la librería cluster, la cual genera una matriz de similaridad, utilizando el coeficiente de Gower, el cual calcula la disimilaridad entre dos filas según la media ponderada de las contribuciones de cada variable, específicamente, \[ d_{ij} = d(i,j) = \frac{\sum_{k=1}^p w_k \delta_{ij}^{(k)} d_{ij}^{(k)}}{\sum_{k=1}^p w_k \delta_{ij}^{(k)}} \]

De todos modos, programáticamente se calcula la matriz de distancias utilizando \(d = \sqrt{1 - s}\), donde \(s\) sería la similaridad calculada mediante la función daisy. A continuación se realiza el escalado multidimensional utilizado en el análisis de coordenadas principales para las variables dicotómicas.

# Librerías
library(cluster)

# Matriz de similaridad
Data_categoric_dist <- round(sqrt(1 - daisy(Data_categoric, type=list(factor=26), metric="gower")),3)

# Escalado multidimensional
Data_categoric_MDS <- cmdscale(Data_categoric_dist, eig=TRUE, k=2)

Descripción de los individuos

Un primer gráfico de la proyección de los individuos en el nuevo sistema de coordenadas es el siguiente.

# Generar dataframe
Data_categoric_MDS_coord <- data.frame(Name=Data$Id,
                                   Coord1=Data_categoric_MDS$points[,1],
                                   Coord2=Data_categoric_MDS$points[,2])

Data_categoric_MDS_coord_alt <- cbind(Data_categoric_MDS_coord, Data$G20) %>%
  mutate(G20 = ifelse(Data$G20 == 1, 1, NA)) %>%
  na.omit()

# Gráfico
ggplot(Data_categoric_MDS_coord, aes(x=Coord1, y=Coord2)) +
  geom_point() +
  geom_text_repel(data=Data_categoric_MDS_coord_alt, aes(x=Coord1, y=Coord2, label=Name), max.overlaps=20) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="PCoA - Individuos", x="Coordenada 1", y="Coordenada 2")

Una primera apreciación de este gráfico es que los individuos se colocan alrededor del origen de coordenadas, y se distribuyen de una manera bastante radial. Podemos observar que existen puntos muy dispersos y lejanos al gran conglomerado de puntos. A modo ilustrativo se colocan los nombres de los páises pertenecientes al bloque intergubernamental G20, y como ninguno de estos países se encuentra alejado del conglomerado, podemos intentar reducir el gráfico a tan solo los puntos agrupados en el origen.

# Generar dataframe
Data_categoric_MDS_coord <- data.frame(Name=Data$Id,
                                   Coord1=Data_categoric_MDS$points[,1],
                                   Coord2=Data_categoric_MDS$points[,2],
                                   Region=Data$Region)

Data_categoric_MDS_coord_alt <- cbind(Data_categoric_MDS_coord, Data$G20) %>%
  mutate(G20 = ifelse(Data$G20 == 1, 1, NA)) %>%
  na.omit()

# Gráfico
ggplot(Data_categoric_MDS_coord, aes(x=Coord1, y=Coord2, fill=Region)) +
  geom_point(aes(color=Region), size=3) +
  geom_label_repel(data=Data_categoric_MDS_coord_alt, aes(x=Coord1, y=Coord2, label=Name), max.overlaps=20) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="PCoA - Individuos", x="Coordenada 1", y="Coordenada 2") +
  xlim(-0.15, 0.2) + ylim(-0.15, 0.15)
## Warning: Removed 9 rows containing missing values (`geom_point()`).

Es importante destacar nuevamente la alta dispersión de los datos, y la falta de claridad en la generación de grupos a partir de una simple observación. Las regiones no se definen claramente. Lo único que se puede destacar es la cercanía de los países del G20 al centro de coordenadas, con lo cual se podría inferir que una cercanía hacia el centro de coordenadas podría implicar un mayor desarrollo del país.

# Generar dataframe
Data_categoric_MDS_coord <- data.frame(Name=Data$Id,
                                   Coord1=Data_categoric_MDS$points[,1],
                                   Coord2=Data_categoric_MDS$points[,2],
                                   Group=Data$Income.Group)

Data_categoric_MDS_coord_alt <- cbind(Data_categoric_MDS_coord, Data$G20) %>%
  mutate(G20 = ifelse(Data$G20 == 1, 1, NA)) %>%
  na.omit()

# Gráfico
ggplot(Data_categoric_MDS_coord, aes(x=Coord1, y=Coord2, fill=Group)) +
  geom_point(aes(color=Group), size=3) +
  geom_label_repel(data=Data_categoric_MDS_coord_alt, aes(x=Coord1, y=Coord2, label=Name), max.overlaps=20) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="PCoA - Individuos", x="Coordenada 1", y="Coordenada 2") +
  xlim(-0.15, 0.2) + ylim(-0.15, 0.15)
## Warning: Removed 9 rows containing missing values (`geom_point()`).

Un gráfico coloreado por grupo de ingreso tampoco aporta mucha información. La dispersión sigue siendo exagerada, pero podemos distinguir un poco que los países de ingresos altos y medio-altos están más cerca del eje de coordenadas, sin claros grupos visibles.

Análisis de Independencia

Un paso siguiente en el trabajo sería la inclusión de un análisis de independencia entre pares de variables. Este análisis intenta probar la independencia (o la falta de independencia) entre dos variables. En caso de demostrarse la independencia un análisis de correspondencia no sería necesario, y simplemente debería hacerse un análisis bivariado simple para intentar observar alguna potencial correlación entre las variables.

El test de independencia de chi cuadrado examina si es que las filas y columnas de una tabla de contingencia tienen una asociación estadística significativa. Para ello se plantean las siguientes hipótesis:

  • Hipótesis nula (\(H_0\)): las filas y las columnas de una tabla de contingencia son independientes. Es decir, no existe una relación entre ambas variables analizadas.

  • Hipótesis alternativa (\(H_a\)): las filas y las columnas de una tabla de contingencia son dependientes. Es decir, existe una clara relación entre ambas variables analizadas.

Para cada celda, el valor esperado se calcula como sigue: \[ e = \frac{\text{Suma por filas} \cdot \text{Suma por columnas}}{\text{Suma total}} \]

El estadístico \(\chi^2\) observado se calcula como sigue: \[ \chi^2 = \sum \frac{(o - e)^2}{e} \] En donde \(o\) es el valor observado, y \(e\) es el valor esperado.

Al calcular el estadístico \(\chi^2\), este se compara con un valor crítico tabulado con \((r-1)(c-1)\) grados de libertad, donde \(r\) es el número de filas y \(c\) es el número de columntas en la tabla de coningencia, y un nivel de significancia de 0.05.

Por ejemplo, para las primeras dos variables del subconjunto de variables categoricas:

# Generar tabla de contingencia
Tabla <- table(Data_categoric[,c(1,2)])
Tabla
##                    Goverment.type
## Type                Constitutional monarchy Provisional Republic
##   Country                                 1           0        1
##   Dependency                              3           2        5
##   Sovereign country                      19           9      117
# Generar test chi cuadrad
chisq.test(Tabla)
## Warning in chisq.test(Tabla): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  Tabla
## X-squared = 7.5656, df = 4, p-value = 0.1088

El resultado de este test Chi-cuadrado de Pearson indica lo siguiente:

  • X-squared: es el cálculo del estadístico chi cuadrado, este valor debe ser comparado con un valor crítico tabulado al mismo nivel de significancia y grados de libertad.

  • df: indica los grados de libertad, donde \(r = 3\), y \(c=3\), \(gl = (3 - 1) (3 - 1) = 4\)

  • p-value: es el valor p observado asociado a esta prueba en particular. Este valor puede ser más útil que comparar el estadístico chi-cuadrado observado, ya que nos permite comparar con el nivel de significancia establecido de \(\alpha = 0.05\), teniendo en cuenta que estos valores representan el mínimo valor de significancia al cual se rechaza la hipótesis nula con la prueba en particular. Un valor mayor a 0.05 indicaría falta de evidencia para rechazar la hipótesis nula.

Dado que \(0.1 > 0.05\), podemos concluir que, al nivel de significancia establecido, no existe evidencia suficiente como para rechazar la hipótesis nula, es decir, las variables seleccionadas presentan independencia, no están relacionadas, y sería inútil realizar un análisis de correspondencia para intentar observar algún patrón de comportamiento.

Si bien este trabajo puede ser realizado manualmente tomando pares de todas las variables cualitativas categóricas del conjunto de datos, sería útil poder hacerlo programáticamente en conjunto.

Primero podríamos intentar averiguar si es que las variables presentan valores faltantes. La presencia de valores faltantes dificulta la realización de un análisis factorial de correspondencia, si bien “NA” corresponde a un factor aparte, su existencia no aporta información alguna.

# Mostrar columnas con valores faltantes
colSums(is.na(Data_categoric))[colSums(is.na(Data_categoric))>0]
##                       Regime.type           First.official.language 
##                                12                                 2 
##                Principal.Religion                  Lending.category 
##                                 2                                33 
##       System.of.National.Accounts Balance.of.Payments.Manual.in.use 
##                                 2                                 1 
##                   System.of.trade     Government.Accounting.concept 
##                                 3                                28

En este caso vemos que son varias las variables que presentan valores faltantes. En este caso, Lending.category y Government.Accounting.concep son las dos que más valores faltantes presentan. No son variables útiles para el siguiente análisis. Además, las variables System.of.Nationa.Accounts, Balance.of.Payments.Manual.in.use, y System.of.trade presentan muy pocos factores, por lo cual se toma la decisión de también abandonar estas variables, y no incluirlas en el análisis. La eliminación de todas estas variables se realiza en la sentencia siguiente.

Data_categoric_fixed <- Data_categoric %>%
  select(-c(Lending.category, System.of.National.Accounts, System.of.trade, Balance.of.Payments.Manual.in.use, Government.Accounting.concept))

A continuación se realizan los tests estadísticos de chi-cuadrado para cada par de variables del conjunto de 9 variables. Finalmente se imprime la matriz de valores p para cada test de hipótesis, ya que se va a intentar realizar el mismo análisis descrito anteriormente comparando el valor-p obtenido.

# Librerías
library(combinat)
## 
## Attaching package: 'combinat'
## The following object is masked from 'package:utils':
## 
##     combn
# Calcular combinaciones posibles de variables
Combinations <- combn(colnames(Data_categoric_fixed), 2)

# Matriz de p-values
p_values_matrix <- matrix(0, nrow=21, ncol=21)
colnames(p_values_matrix) <- colnames(Data_categoric_fixed)
rownames(p_values_matrix) <- colnames(Data_categoric_fixed)

# Ciclo para calcular las tablas de contingencias
for (i in 1:210){
  # Tabla de contingencia
  
  Contingency_table = table(Data_categoric_fixed[,Combinations[,i]])
  
  # Test chi cuadrado
  Chisqt_test = chisq.test(Contingency_table)
  
  # Guardar valor p
  p_values_matrix[Combinations[1,i], Combinations[2,i]] = p_values_matrix[Combinations[2,i], Combinations[1,i]] = Chisqt_test$p.value
}

# Reemplazar diagonal por valores faltantes
for (i in 1:21){
  for (j in 1:21)
    if (i == j){
      p_values_matrix[i,j] = NA
    }
}

p_values_matrix[p_values_matrix <= 0.05] <- NA

# Imprimir matriz de valores p, sin la diagonal
print(p_values_matrix, digits=3, na.print="")
##                          Type Goverment.type Regime.type Currency.Unit Region
## Type                                  0.1088       0.525                0.938
## Goverment.type          0.109                                                
## Regime.type             0.525                                                
## Currency.Unit                                                                
## Region                  0.938                                                
## First.official.language 0.896                        NaN                     
## Principal.Religion      0.993                        NaN         0.274       
## Income.Group            0.459                                                
## Major.Economic.Group    0.722         0.3012                                 
## Advanced.Economies      0.598                                                
## G20                     0.226         0.3294                     0.427       
## G7                      0.654         0.0857                     0.134       
## G33                     0.578         0.0506                     0.135       
## NU.Membership                                                                
## NATO                    0.792                                                
## BRICS                   0.706         0.3119       0.149         0.688  0.441
## Colony                  0.606         0.7252       0.247         0.773  0.833
## Empire                  0.578                                                
## Nuclear.Weapons         0.731         0.7021       0.750         0.845       
## Same.sex.Marriage       0.556                                                
## Universal.healthcare    0.179                                                
##                         First.official.language Principal.Religion Income.Group
## Type                                      0.896             0.9927        0.459
## Goverment.type                                                                 
## Regime.type                                 NaN                NaN             
## Currency.Unit                                               0.2741             
## Region                                                                         
## First.official.language                                        NaN             
## Principal.Religion                          NaN                                
## Income.Group                                                                   
## Major.Economic.Group                                                           
## Advanced.Economies                                                             
## G20                                       0.660                                
## G7                                        0.416                                
## G33                                       0.300             0.2294             
## NU.Membership                                                             0.178
## NATO                                                        0.1835             
## BRICS                                                       0.6091             
## Colony                                    0.644             0.1141        0.690
## Empire                                    0.163             0.1885             
## Nuclear.Weapons                           0.293                           0.292
## Same.sex.Marriage                                                              
## Universal.healthcare                                        0.0625             
##                         Major.Economic.Group Advanced.Economies   G20     G7
## Type                                   0.722              0.598 0.226 0.6540
## Goverment.type                         0.301                    0.329 0.0857
## Regime.type                                                                 
## Currency.Unit                                                   0.427 0.1342
## Region                                                                      
## First.official.language                                         0.660 0.4164
## Principal.Religion                                                          
## Income.Group                                                                
## Major.Economic.Group                                                        
## Advanced.Economies                                                          
## G20                                                                         
## G7                                                                          
## G33                                                             1.000 0.3737
## NU.Membership                                                               
## NATO                                                                        
## BRICS                                  0.101              0.872       1.0000
## Colony                                 0.589              0.291 0.259       
## Empire                                 0.164                                
## Nuclear.Weapons                        0.228                                
## Same.sex.Marriage                                                           
## Universal.healthcare                                                        
##                            G33 NU.Membership  NATO BRICS Colony Empire
## Type                    0.5784               0.792 0.706 0.6056 0.5784
## Goverment.type          0.0506                     0.312 0.7252       
## Regime.type                                        0.149 0.2469       
## Currency.Unit           0.1345                     0.688 0.7727       
## Region                                             0.441 0.8325       
## First.official.language 0.3005                           0.6445 0.1626
## Principal.Religion      0.2294               0.184 0.609 0.1141 0.1885
## Income.Group                           0.178             0.6902       
## Major.Economic.Group                               0.101 0.5893 0.1642
## Advanced.Economies                                 0.872 0.2913       
## G20                     1.0000                           0.2588       
## G7                      0.3737                     1.000              
## G33                                                1.000 1.0000 0.1372
## NU.Membership                                                         
## NATO                                               0.380 0.8496       
## BRICS                   1.0000               0.380       1.0000       
## Colony                  1.0000               0.850 1.000        0.0783
## Empire                  0.1372                           0.0783       
## Nuclear.Weapons         0.4335               0.309       0.0714       
## Same.sex.Marriage                                  0.401 1.0000       
## Universal.healthcare                                     0.9381       
##                         Nuclear.Weapons Same.sex.Marriage Universal.healthcare
## Type                             0.7311             0.556               0.1793
## Goverment.type                   0.7021                                       
## Regime.type                      0.7496                                       
## Currency.Unit                    0.8447                                       
## Region                                                                        
## First.official.language          0.2933                                       
## Principal.Religion                                                      0.0625
## Income.Group                     0.2921                                       
## Major.Economic.Group             0.2276                                       
## Advanced.Economies                                                            
## G20                                                                           
## G7                                                                            
## G33                              0.4335                                       
## NU.Membership                                                                 
## NATO                             0.3089                                       
## BRICS                                               0.401                     
## Colony                           0.0714             1.000               0.9381
## Empire                                                                        
## Nuclear.Weapons                                     0.401               0.2455
## Same.sex.Marriage                0.4014                                       
## Universal.healthcare             0.2455

Ahora, lo importante es comparar cada valor-p obtenido con el valor crítico de 0.05, o 5.00e-02, valores superiores a este demuestran la falta de evidencia estadística para rechazar la hipótesis nula, e implican la independencia entre ese par de variables.

La observación más importante de esta matriz, es que la variable Type, que indica el tipo de Estado, es la variable que más valores significativamente altos presenta. Esto demuestra que falla la intención de rechazar la hipótesis nula para la gran mayoría de las variables del conjunto de datos, es decir, que Type es una variable independiente con casi todas las otras variables, no existe una clara relación entre esta y las otras variables, por lo tanto su inclusión en un análisis factorial de correspondencia sería inútil. Su eliminación se realiza en la siguiente sentencia.

Data_categoric_fixed <- Data_categoric_fixed %>%
  select(-c(Type, BRICS, G33, Colony, Nuclear.Weapons))

Análisis Factorial de Correspondencia Múltiple

El objetivo del AFCM es hallar la mejor representación de la tabla de datos original, es decir encontrar un nuevo sistema de coordenadas preferentemente de dos dimensiones que mejor resuma la información contenida en las tablas de contingencia del conjunto de variables categóricas.

El modelo se desarrolla a partir de la elaboración de la tabla de contingencia en base a lo que se denomina la tabla de códigos que correspondan a las diferentes modalidades de cada una de las características observadas que presentan los individuos. La tabla de códigos es de tamaño \(n \times p\), donde \(n\) representa el número de individuos y \(p\) las variables. Luego, a prtir de esta tabla, es posible generar una tabla que contemple la misma información, denominada tabla disyuntiva completa. Si esto se repite para cada característica observada se contará con una matriz compuesta por 0 y 1 de las \(S\) modalidades de \(p\) características observadas. En donde 1 signifique la presencia de la modalidad en el individuo y 0 la ausencia. La tabla disyuntiva completa asigna a cada individuo una modalidad de la característica observada de modo que la población se divide en tantos subconjuntos como modalidades tenga.

A continuación se generan nuevas dimensiones y se realizan proyecciones de los individuos o de los factores sobre este nuevo sistema de coordenadas. La distancia chi-cuadrado es la medida de la distancia entre dos puntos proyectados en el sistema de coordenadas. Dos individuos que presentaron un gran número de modalidades en común deben ser incluidos en una misma clase de la tipología de individuos.

En un análisis equivalente al de las filas se aplica al de las columnas y luego se realiza un análisis conjunto determinando así un plano factorial de dos dimensiones para las características observadas.

Para la interpretación de los planos factoriales, debemos tener en cuenta que dos modalidades son semejantes cuanto más cerca estén. Es decir, que dos modalidades de respuesta estén cercanas significa que son casi los mismos en que tienen ese tipo de respuesta. Por otra parte debe considerarse que las modalidades cercanas al centro son las respuestas más comunes, mientras que en los extremos de los ejes se ubicarán aquellas más características por alguna razón, expresada en el factor de variabilidad correspondiente.

# Librerías
library(FactoMineR)

# Análisis Factorial de Correspondencia
Data_categoric_MCA <- MCA(Data_categoric_fixed, ncp=2, graph=FALSE)

Descripción de los factores

A continuación mostramos la proyección de las distancias entre los factores en el nuevo sistema de coordenadas.

# Librerías
library(dplyr)
library(stringr)

# Generar grupos de categorías
Variables <- c()
for (i in 1:dim(Data_categoric_fixed)[2]){
  n = as.numeric(count(unique(Data_categoric_fixed[i])))
  name = colnames(Data_categoric_fixed[i])
  Variables = append(Variables, rep(name, n))
}
  
# Filtrar datos
MCA_Categories <- data.frame(Category=rownames(Data_categoric_MCA$var$coord),
                             Data_categoric_MCA$var$coord,
                             Variables,
                             row.names=NULL) %>%
  mutate(Category = str_replace(Category, "First.official.language_", ""),
                          Category = str_replace(Category, "Currency.Unit_", ""))

# Gráfico
ggplot(MCA_Categories, aes(x=Dim.1, y=Dim.2)) +
  geom_point(aes(color=Variables)) +
  geom_text_repel(size=3, aes(label=Category, color=Variables)) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="MCA - Factores", x="Dimensión 1", y="Dimensión 2")

A modo ilustrativo se colorean tanto los puntos como las etiquetas teniendo en cuenta a qué varialbe pertenece cada factor. El resultado es sin duda interesante, presentando en el plano una forma de estrella con tres puntas. En la punta superior podemos observar la relación entre los países latinoamericanos con el idioma español, el principal idioma en este subcontinente, su pertenencia al MERCOSUR, el Central American Parliament, la Andean Community, y su cercanía hacia la religión católica romana, la cual es la más importante en esta región. Llegando hacia el orígen de las coordenadas observamos la aparición de la región de Norteamérica, y el grupo económico NAFTA, compuesto por países norteamericanos, allí se hace más obvia la aparición del dólar como moneda principal de esta rama. Podríamos definir a esta rama superior como compuesta por países principalmente de ingresos medios-altos, aunquen o existe un tipo de régimen claro en esta rama.

La rama derecha observamos que está relacionada con idiomas como el Alemán, el Griego, y el neerlandés, el grupo económico EFTA, la Unión Europea, y la moneda Euro. Es decir, esta rama está completamente relacionada con los países europeos y de Asia central, con la inclusión de Norte América. Tienen asociado un tipo de gobierno monárquico constitucional, un tipo de régimen democrático, y poca religión religiosa. Son claramente países principalmente de altos ingresos y economías avanzadas. Además, con mucha presencia en los grupos G20, G7. Por otra parte, gran cantidad de estos países fueron imperios en su historia. En la actualidad cuentan con derechos universales como el matrimonio igualitario o una cobertura de salud pública.

En la rama inferior izquierda observamos que los individuos tendrían tipos de gobiernos provisionales, y regímenes autoritarios, el islam y el hinduísmo como las religiones principales, y francos africanos como principales monedas, además del Francés, Árabe y Persa como idiomas predominantes. Claramente estos países son de bajos ingresos, y están asociados a países de toda África y el medio oriente.

Al centro del sistema de coordenadas encontramos, como se mencionó, las modalidades más frecuentes, o quizá los inviduos que no se caracterizan por algo en particular. En esta zona se encuentran principalmente países asiáticos, grupos como el ASEAN, el SAARC y el CIS, regímenes poco definidos, pero una alta presencia de gobiernos de tipo república, además de diversos idiomas y monedas poco distinguidos.

Entonces, a partir de estas observaciones podemos generar el siguiente gráfico.

Y siguiendo con los análisis anteriores, podemos presentar los individuos más representativos de cada región, para así poder apoyar nuestras observaciones.

# Generar grupos de categorías
Variables <- c()
for (i in 1:dim(Data_categoric_fixed)[2]){
  n = as.numeric(count(unique(Data_categoric_fixed[i])))
  name = colnames(Data_categoric_fixed[i])
  Variables = append(Variables, rep(name, n))
}
  
# Filtrar datos
MCA_Categories <- data.frame(Category=rownames(Data_categoric_MCA$var$coord),
                             Data_categoric_MCA$var$coord,
                             Variables,
                             row.names=NULL) %>%
  mutate(Category = str_replace(Category, "First.official.language_", ""),
                          Category = str_replace(Category, "Currency.Unit_", ""))

MCA_Ind <- data.frame(Data_categoric_MCA$ind$coord,
                      Region=Data$Region)

# Lista de países representativos por región
Representative_countries <- rbind(cbind(Id=Data$Id,
                                        IDH=Data$Human.Development.Index.HDI.2014,
                                        Region=Data$Region,
                                        data.frame(Data_categoric_MCA$ind$coord)) %>%
  group_by(Region) %>%
  slice_max(order_by=IDH, n=2) %>%
  select(Id, Dim.1, Dim.2, Region),
  cbind(Id=Data$Id,
        IDH=Data$Human.Development.Index.HDI.2014,
        Region=Data$Region,
        data.frame(Data_categoric_MCA$ind$coord)) %>%
  group_by(Region) %>%
  slice_min(order_by=IDH, n=2) %>%
  select(Id, Dim.1, Dim.2, Region),
  cbind(Id=Data$Id,
        G20=Data$G20,
        data.frame(Data_categoric_MCA$ind$coord),
        Region=Data$Region) %>%
  mutate(Id = ifelse(G20==1, Id, NA)) %>%
  select(Id, Dim.1, Dim.2, Region) %>%
  na.omit()) %>%
  distinct()

# Gráfico
ggplot(Representative_countries, aes(x=Dim.1, y=Dim.2)) +
  geom_point(aes(color=Region)) +
  geom_text_repel(size=3, aes(label=Id)) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="MCA - Factores", x="Dimensión 1", y="Dimensión 2")

Este gráfico apoya nuestras observaciones. En la rama superior encontramos países americanos, en la rama inferior izquierda encontramos países africanos y árabes, mientras que hacia la derecha se encuentran los países europeos y Norteamérica.

Métodos de Clasificación

Los métodos de clasificación hacen referencia a un amplio abanico de técnicas no supervisadas cuya finalidad es encontrar patrones o grupos (clusters) dentro de conjunto de observaciones. Las particiones se establecen de forma que, las observaciones que están dentro de un mismo grupo, son similares entre ellas y distintas a las observaciones de otros grupos. Se trata de métodos no supervisados ya que el proceso ignora la variable respuesta que indica a qué grupo pertenece realmente cada observación (si es que existe tal variable).

Todos los métodos de clustering tienen una cosa en común, para poder llevar a cabo las agrupaciones necesitan definir y cuantificar la similitud entre las observaciones. El término distancia se emplea dentro del contexto del clustering como cuantificación de la similitud o diferencia entre observaciones. Si se representan las observaciones en un espacio \(p\) dimensional, siendo \(p\) el número de variables asociadas a cada observación, cuando más se asemejen dos observaciones, más próximas estarán, de ahí que se emplee el término distancia. Al igual que en otros métodos estadísticos, la escala en la que se miden las variables y la magnitud de su varianza pueden afectar en gran medida a los resultados obtenidos por clustering. Si una variable tiene una escala mucho mayor que el resto, determinará en gran medida el valor de distancia/similitud obtenido al comparar las observaciones, dirigiendo así la agrupación final. Escalar y centrar las variables de forma que todas ellas tengan media 0 y desviación estándar 1 antes de calcular la matriz de distancias asegura que todas las variables tengan el mismo peso cuando se realice el clustering.

Clustering jerárquico

El clustering jerárquico es una alternativa a los métodos de clustering que no requiere que se pre-especifique el número de clusters. Los métodos que engloba el clustering jerárquico se subdividen en dos tipos dependiendo de la estrategia seguida para crear los grupos:

  • Clustering aglomerativo: el agrupamiento se inicia en la base del árbol, donde cada observación forma un cluster individual. Los clusters se van combinando a medida que la estructura crece hasta converger en una única “rama” central.

  • Clustering divisivo: es la estrategia opuesta al aglomerativo, se inicia con todas la observaciones contenidas en un mismo cluster y se suceden divisiones hasta que cada observación forma un cluster individual.

En cambos casos, los resultados pueden representarse de forma muy intuitiva en una estructura de árbol llamada dendrograma.

En el caso a analizar se opta por utilizar el clustering aglomerativo, ya que es la técnica más explorada. La estructura resultante del clustering aglomerativo se obtiene mediante un algoritmo sencillo.

  1. El proceso se inicia considerando cada una de las observaciones como un cluster individual.

  2. Se inicia un proceso iterativo hasta que todas la sobservaciones pertenecen a un único cluster:

  • Se calcula la distancia entre cada posible par de los \(n\) clusters. El investigador determina el tipo de medida que se emplea para cuantificar la similitud entre observaciones o grupos.

  • Los dos clusters más similares se fusionan, de forma que quedan \(n-1\) clusters.

  1. Determinar dónde cortar la estructura de árbol generada (dendrograma).

Para que el proceso de agrupamiento pueda llevarse a cabo tal como lo indica el algoritmo anterior, es necesario definir cómo se cuantifica la similitud entre dos clusters. Es decir, se tiene que extender el concepto de distancia entre pares de observaciones para que sea aplicable a pares de grupos, cada uno formado por varias observaciones. Los tipos de agrupamiento más empleados son:

  • Simple o mínimo: se calcula la distancia entre todos los posibles pares formados por una observación del cluster A y una del cluster B. La menor de todas ellas se selecciona como la distancia entre los dos clusters. Se trata de la medida menos conservadora.

  • Completo o máximo: se calcula la distancia entre todos los posibles pares formados por una observación del cluster A y una del cluster B. La mayor de todas ellas se selecciona como la distancia entre los dos clusters. Se trata de la medida más conservadora.

  • Average: se calcula la distancia entre todos los posibles pares formados por una observación del cluster A y una del cluster B. El valor promedio de todas ellas se selecciona como la distancia entre los dos clusters.

  • Centroide: se calcula el centroide de cada uno de los clusters y se selecciona la distancia entre ellos como la distancia entre los dos clusters.

  • Ward: se trata de un método general. La selección del par de clusters que se combinan en cada paso del clustering aglomerativo se basa en el valor óptimo de una función objetivo, pudiendo ser esta última cualquier función. El conocido método “Ward’s minimum variance” es un caso particular en el que el objetivo es minimizar la suma total de varianza intra-cluster.

Comparación de técnicas

Una vez creado el modelo de clustering hay que evaluar hasta qué punto su estructura refleja las distancias originales entre observaciones. Una forma de hacerlo es empleando el coeficiente de correlación entre las distancias cofenéticas del dendrograma (las distancias seleccionadas por el método de agrupamiento) y la matriz de distancias original. Cuanto más cercano es el valor a 1, mejor refleja el dendrograma la verdadera similitud entre las observaciones.

# Matriz de distancias
Data_dist <- dist(scale(Data_numeric_fixed), method="euclidean", diag=TRUE)

# Clusters
Data_clust_min <- hclust(Data_dist, method="single")
Data_clust_max <- hclust(Data_dist, method="complete")
Data_clust_upgma <- hclust(Data_dist, method="average")
Data_clust_ward2 <- hclust(Data_dist, method="ward.D2")

# Matrices cofenéticas de los clusters
Data_min_cof <- cophenetic(Data_clust_min)
Data_max_cof <- cophenetic(Data_clust_max)
Data_upgma_cof <- cophenetic(Data_clust_upgma)
Data_ward2_cof <- cophenetic(Data_clust_ward2)

# Correlaciones
paste0("El agrupamiento mínimo resulta en una correlación: ", cor(Data_dist, Data_min_cof))
## [1] "El agrupamiento mínimo resulta en una correlación: 0.439547348372085"
paste0("El agrupamiento máximo resulta en una correlación: ", cor(Data_dist, Data_max_cof))
## [1] "El agrupamiento máximo resulta en una correlación: 0.720744421185822"
paste0("El agrupamiento UPGMA resulta en una correlación: ", cor(Data_dist, Data_upgma_cof))
## [1] "El agrupamiento UPGMA resulta en una correlación: 0.736749798461531"
paste0("El agrupamiento Ward resulta en una correlación: ", cor(Data_dist, Data_ward2_cof))
## [1] "El agrupamiento Ward resulta en una correlación: 0.701317734722266"

Dado que el agrupamiento por el método UPGMA (average) resulta en una correlación igual a 0.74, podemos concluir que este es el método que logra una mejor representación de las distancias originales en el dendrograma, es decir, una menor pérdida de la representación de los datos.

Dendrograma

A continuación podemos representar el dendrograma del modelo de clustering generado utilizando el método de agrupamiento UPGMA.

# Librerías
library(ggdendro)
## Warning: package 'ggdendro' was built under R version 4.2.3
# Matriz de distancias
Data_dist <- dist(scale(Data_numeric_fixed), method="euclidean", diag=TRUE)

# Cluster
Data_clust <- hclust(Data_dist, method="average")
Data_clust$labels <- Data$Id

# Dendrograma
ggdendrogram(Data_clust, rotate=TRUE)

# Librerías
library(factoextra)

# Matriz de distancias
Data_dist <- dist(scale(Data_numeric_fixed), method="euclidean", diag=TRUE)

# Cluster
Data_clust <- hclust(Data_dist, method="average")
Data_clust$labels <- Data$Id

# Dendrograma
fviz_dend(Data_clust,
          cex=0.5, k=5,
          horiz=TRUE, color_labels_by_k=TRUE) +
  labs(title="Dendrograma - clustering jerárquico - UPGMA", y="Distancia")
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <]8;;https://github.com/kassambara/factoextra/issueshttps://github.com/kassambara/factoextra/issues]8;;>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

En este dendrograma podemos observar que existen al menos tres grandes clusters. A pesar que al colocar \(k = 5\) clusters automáticos, se generan dos que no tienen sentido ya que contienen un solo individuo. Esta elección puede resultar arbitraria y se recomienda probar distintos valores y observar sus diferencias. De hecho el aumento del número de clusters esperados no cambia en absoluto el contenido de los tres grupos observados, sino que genera más grupos con uno o solo dos individuos, apoyando nuestra observación de tres grupos mayoritarios.

Podemos utilizar la generación de estos clusters en la proyección de los individuos sobre el plano de coordenadas generado con el análisis de componentes principales.

# Librerías
library(factoextra)
library(gridExtra)
library(RColorBrewer)
library(ggforce)
library(ggrepel)

# Paleta
Paleta <- brewer.pal(7, "Set1")

# Cortar dendrograma
Cutree <- cutree(Data_clust, k=5)

# Filtrar labels
Data_PCA_clust <- cbind(Data, Data_PCA$ind$coord) %>%
  select(Id, Dim.1, Dim.2) %>%
  arrange(Id)
Cutree <- data.frame(Id=names(Cutree), Group=as.factor(Cutree), row.names=NULL)
Data_PCA_clust <- merge(Data_PCA_clust, Cutree, by="Id")

# Lista de países representativos por región
Representative_countries <- rbind(cbind(Id=Data$Id,
                                        IDH=Data$Human.Development.Index.HDI.2014,
                                        Region=Data$Region,
                                        data.frame(Data_PCA$ind$coord)) %>%
  group_by(Region) %>%
  slice_max(order_by=IDH, n=2) %>%
  select(Id, Dim.1, Dim.2, Region),
  cbind(Id=Data$Id,
        IDH=Data$Human.Development.Index.HDI.2014,
        Region=Data$Region,
        data.frame(Data_PCA$ind$coord)) %>%
  group_by(Region) %>%
  slice_min(order_by=IDH, n=2) %>%
  select(Id, Dim.1, Dim.2, Region),
  cbind(Id=Data$Id,
        G20=Data$G20,
        data.frame(Data_PCA$ind$coord),
        Region=Data$Region) %>%
  mutate(Id = ifelse(G20==1, Id, NA)) %>%
  select(Id, Dim.1, Dim.2, Region) %>%
  na.omit()) %>%
  distinct()

# Gráfico
ggplot(Data_PCA_clust, aes(x=Dim.1, y=Dim.2)) +
  geom_point(aes(color=Group)) +
  geom_mark_ellipse(aes(fill=Group)) +
  geom_text_repel(data=Representative_countries, aes(x=Dim.1, y=Dim.2, label=Id)) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="PCA - Clustering", x="Dim. 1 (61.6%)", y="Dim. 2 (6.8%)")

Ahora observamos que de hecho los grupos generados apoyan nuestra observación inicial de la existencia de tres grupos sobre estas componentes principales. Un grupo estaría formado con países africanos (de bajos ingresos), otro grupo con países de economías emergentes, y un tercer gran grupo con países desarrollados. Nuevamente se comenta sobre la existencia de dos grupos con un solo individuo, los cuales no tienen sentido en el análisis.

Número óptimo de clusters

Determinar el número óptimo de clusters es quizá el paso más complicado a la hora de aplicar métodos de clustering. No existe una forma única de averiguar el número adecuado de clusters. Es un proceso bastante subjetivo que depende en gran medida del tipo de clustering empleado y de si se dispone de información previa sobre los datos con los que se está trabajando, por ejemplo, estudios anteriores pueden sugerir o acotar las posibilidades. Sin embargo se desarrollaron varias estrategias que ayudan en el proceso.

Un primer método, llamado “método Elbow”, sigue una estrategia comúnmente empleada para encontrar el valor óptimo de un hiperparámetro. La idea general es probar un rango de valores del hiperparámetro en cuestión, representar gráficamente los resultados obtenidos con cada uno e identificar aquel punto de la curva a partir de la cual la mejora deja de ser sustancial (principio de verosimilitud). En los casos de clustering de fraccionamiento, como por ejemplo K-means, las observaciones se agrupan de una forma tal que se minimiza la varianza intra-cluster. El método Elbow calcula la varianza total intra-cluster en función del número de clusters y elige como óptimo aquel valor a partir del cual añadir más clusters apenas consigue mejoría. Sumado a esto es posible graficar la varianza total entre-clusters, en función del número de clusters. La intersección entre ambos gráficos indica el valor mínimo de clusters del cual se debe partir. Además, observando ambas curvas es posible definir con mayor claridad el número óptimo de clusters.

SS <- data.frame(Cluster=as.factor(1:10), SCD=rep(0,10), SCE=rep(0,10))
for (i in 1:10) {SS$SCD[i] <- sum(kmeans(scale(Data_PCA$ind$coord[,1:2]), centers=i)$withinss)}
for (i in 1:10) {SS$SCE[i] <- sum(kmeans(scale(Data_PCA$ind$coord[,1:2]), centers=i)$betweenss)}

ggplot(SS, aes(x=Cluster)) +
  geom_point(aes(y=SCD), show.legend = F, size = 3, color="red") +
  geom_line(aes(y=SCD), show.legend = F, group = 1, linewidth=1, color="red") +
  geom_point(aes(y=SCE), show.legend = F, size = 3, color="blue") +
  geom_line(aes(y=SCE), show.legend = F, group = 1, linewidth=1, color="blue") +
  labs(x="Número de clusters", y="Suma de cuadrados", title="Comparación SCD-SCE para obtener el número óptimo de clusters")

Este resultado indica, a partir de la comparación de SCD-SCE, que como mínimo se debe optar por 2 clusters. Un tercer cluster también podría resultar beneficioso, porque es claro que a partir de la adición de un cuarto o quinto cluster, la mejoría resulta redundante e innecesaria.

Un segundo método, denominado “average silhouette” es muy similar al de Elbow, con la diferencia de que, en lugar de minimizar el total de la suma de cuadrados intra-clusters, se maximiza la media de los silouette coeficient o índices silueta (\(s_i\)). Este coeficiente cuantifica qué tan buena es la asignación que se ha hecho de una observación comparando su similitud con el resto de observaciones de su cluster frente a las de los otros clusters. Su valor puede estar entre -1 y 1, siendo valores altos un indicativo de que la observación se ha asignado al cluster correcto.

#Librerías
library(factoextra)

# Gráfico
fviz_nbclust(scale(Data_PCA$ind$coord[,1:2]), kmeans, method = "silhouette", k.max = 10) + theme_minimal() + ggtitle("The Silhouette Plot")

El resultado del Silhouette Plot indica que 3 clusters es el número óptimo de clusters, aquel que maximiza la media del silhouette coeficient.

Los métodos mencionados no tiene por qué coincidir exactamente en su estimación del número óptimo de clusters, pero tienden a acotar el rango de posibles valores. Por esta razón es recomendable calcular los tres y en función de los resultados decidir. Además de estos métodos, la función NbClust() del paquete NbClust incorpora 30 índices distintos, dando la posibilidad de calcularlos todos en un único paso. Esto último es muy útil ya que permite identificar el valor en el que coinciden más índices, aportando seguridad de que se está haciendo una buena elección.

# Librerías
library(NbClust)

# Número óptimo de clusters
NbClust(scale(Data_numeric_fixed), distance="euclidean", min.nc=2, max.nc=10, method="kmeans", index="alllong")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 9 proposed 2 as the best number of clusters 
## * 11 proposed 3 as the best number of clusters 
## * 1 proposed 5 as the best number of clusters 
## * 2 proposed 8 as the best number of clusters 
## * 4 proposed 10 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  3 
##  
##  
## *******************************************************************
## $All.index
##        KL       CH Hartigan     CCC     Scott      Marriot    TrCovW   TraceW
## 2  3.5602 131.6744  47.9308 -0.2283  297.5315 1.319314e+49 27467.722 2699.090
## 3  4.0485 109.4508  14.2862  1.8357  594.2836 4.483869e+48  9071.860 2061.584
## 4  1.2344  83.9495  11.9422  1.9002  781.8221 2.414127e+48  6733.317 1886.572
## 5  3.9943  70.3990   4.2106  2.4807 1011.3167 8.744815e+47  5636.266 1749.980
## 6  0.6961  58.3351   5.0620  1.9091 1132.2208 5.829982e+47  5339.079 1702.810
## 7  0.5500  50.7475   8.0421  1.7090 1304.6973 2.645179e+47  5005.598 1647.579
## 8  4.3886  46.6657   2.6510  2.3826 1430.0023 1.555320e+47  4389.356 1563.740
## 9  0.2486  41.6092   7.4357  1.8404 1527.7666 1.056060e+47  4293.783 1536.404
## 10 1.6712  39.4023   4.8160  2.6116 1703.3057 4.262122e+46  3876.374 1462.906
##    Friedman  Rubin Cindex     DB Silhouette    Duda Pseudot2    Beale Ratkowsky
## 2  387.5779 1.8495 0.4004 1.0782     0.3699  0.8887  13.4025   2.7970    0.4713
## 3  513.8575 2.4214 0.3836 1.2544     0.2848  1.4800 -20.4333  -7.1554    0.4364
## 4  538.2310 2.6461 0.4088 1.6428     0.2108  1.4705 -15.6770  -7.0199    0.3889
## 5  573.0296 2.8526 0.3909 1.9426     0.1560  0.8719   3.5265   3.1861    0.3562
## 6  592.8753 2.9316 0.3903 1.9324     0.1520  0.5161   6.5629  20.0000    0.3276
## 7  598.6264 3.0299 0.3832 1.8620     0.1525  1.4181 -13.5623  -6.4127    0.3064
## 8  618.0161 3.1923 0.3721 1.9651     0.1237 11.7903 -43.0137 -10.3355    0.2903
## 9  656.7206 3.2491 0.3567 1.9829     0.1126  0.8203   4.8193   4.7229    0.2747
## 10 655.1810 3.4124 0.3442 1.8946     0.1220  2.2961 -22.0143 -12.1952    0.2641
##         Ball Ptbiserial     Gap    Frey McClain  Gamma    Gplus      Tau   Dunn
## 2  1349.5449     0.6109 -0.0896  0.8494  0.5918 0.7150 436.2026 2188.683 0.1925
## 3   687.1948     0.5946 -0.6184  1.1551  1.1064 0.7849 295.2659 2154.378 0.1599
## 4   471.6429     0.5536 -1.1387  1.6205  1.4953 0.7901 254.4840 1915.976 0.1979
## 5   349.9960     0.4919 -1.4347  0.6298  2.1154 0.7752 227.0279 1566.121 0.1948
## 6   283.8017     0.4850 -1.7843  0.0155  2.2442 0.7815 212.0305 1517.105 0.1974
## 7   235.3684     0.4898 -1.8760  1.5899  2.2634 0.7977 193.7241 1527.554 0.1974
## 8   195.4675     0.4412 -2.3805  0.3610  2.9637 0.7850 172.3084 1258.523 0.2133
## 9   170.7116     0.4325 -2.5173  0.1888  3.2468 0.8050 144.3385 1191.899 0.1802
## 10  146.2906     0.4311 -2.6383 -0.0806  3.4149 0.8285 120.2708 1162.209 0.1802
##    Hubert SDindex Dindex   SDbw
## 2   4e-04  0.7539 3.9645 0.6854
## 3   5e-04  0.8387 3.4557 0.5811
## 4   5e-04  1.0065 3.3129 0.5620
## 5   6e-04  1.1484 3.1827 0.5341
## 6   6e-04  1.1320 3.1465 0.5640
## 7   6e-04  1.3457 3.1091 0.6612
## 8   6e-04  1.3456 3.0126 0.6314
## 9   6e-04  1.2055 2.9856 0.5904
## 10  6e-04  1.1146 2.9231 0.5447
## 
## $All.CriticalValues
##    CritValue_Duda CritValue_PseudoT2 Fvalue_Beale CritValue_Gap
## 2          0.8964            12.3681            0        0.5350
## 3          0.8610            10.1675            1        0.5285
## 4          0.8485             8.7509            1        0.3062
## 5          0.8252             5.0834            0        0.3617
## 6          0.7989             1.7619            0        0.1046
## 7          0.8308             9.3652            1        0.5208
## 8          0.5241            42.6754            1        0.1548
## 9          0.8154             4.9799            0        0.1403
## 10         0.8189             8.6256            1        0.1493
## 
## $Best.nc
##                     KL       CH Hartigan     CCC    Scott      Marriot   TrCovW
## Number_clusters 8.0000   2.0000   3.0000 10.0000   3.0000 3.000000e+00     3.00
## Value_Index     4.3886 131.6744  33.6446  2.6116 296.7521 6.639524e+48 18395.86
##                   TraceW Friedman   Rubin  Cindex     DB Silhouette Duda
## Number_clusters   3.0000   3.0000  3.0000 10.0000 2.0000     2.0000 3.00
## Value_Index     462.4928 126.2796 -0.3473  0.3442 1.0782     0.3699 1.48
##                 PseudoT2   Beale Ratkowsky     Ball PtBiserial     Gap Frey
## Number_clusters   3.0000  3.0000    2.0000   3.0000     2.0000  2.0000    1
## Value_Index     -20.4333 -7.1554    0.4713 662.3501     0.6109 -0.0896   NA
##                 McClain   Gamma    Gplus      Tau   Dunn Hubert SDindex Dindex
## Number_clusters  2.0000 10.0000  10.0000    2.000 8.0000      0  2.0000      0
## Value_Index      0.5918  0.8285 120.2708 2188.683 0.2133      0  0.7539      0
##                   SDbw
## Number_clusters 5.0000
## Value_Index     0.5341
## 
## $Best.partition
##   [1] 2 3 3 2 1 3 1 1 3 2 1 1 3 2 3 3 3 3 3 1 2 2 2 2 1 2 2 1 3 3 2 3 2 1 1 1 1
##  [38] 2 3 3 3 3 2 2 1 2 3 1 1 3 2 3 1 2 1 3 2 2 3 2 3 1 1 3 3 3 3 1 1 1 3 1 3 1
##  [75] 2 3 2 1 1 2 2 1 1 2 2 1 3 2 1 2 3 3 3 3 3 3 3 2 2 3 2 1 3 2 2 1 2 3 3 2 3
## [112] 3 3 1 1 3 1 2 3 2 3 3 2 1 1 2 3 1 2 1 3 2 3 1 1 3 3 2 3 2 3 3 3 3 3 2 1 1
## [149] 1 3 3 3 3 3 2 2 2

Este resultado indica, a través de una regla de mayoría, que el número óptimo de clusters es 3.

K-means clustering

El método K-means clusterings agrupa las observaciones en \(K\) clusters distintos, donde el número \(K\) lo determina el analista antes de ejecutar el algoritmo. K-means clusterings encuentra los \(K\) mejores clusters, entendiendo como mejor cluster aquel cuya varianza interna (intra-cluster) sea lo más pequeña posible. Se trata por lo tanto de un problema de optimización, en el que se reparten las observaciones en \(K\) clusters de forma que la suma de las varianzas internas de todos ellos sea lo menor posible. Para poder solucionar este problema es necesario definir un modo de cuantificar la varianza interna.

Considérese \(C_1, ..., C_k\) como los sets formados por los índices de las observaciones de cada uno de los clusters. Por ejemplo, el set \(C_1\) contiene los índices de las observaciones agrupadas en el cluster 1. La nomenclatura empleada para indicar que la observación \(i\) pertenece al cluster \(k\) es \(i \in C_k\). Todos los sets satisfacen dos propiedades:

  1. \(C_1 \cup C_2 \cup ... \cup C_k = \{1, ..., n\}\). Significa que toda observación pertenece al menos a uno de los \(k\) clusters.

  2. \(C_k \cap C_{k'} = \emptyset\) para todo \(k \ne k'\). Implica que los clusters no solapan, ninguna observación pertenece a más de un cluster a la vez.

Dos de las medidas más comúnmente empleadas definen la varianza interna de un cluster \(W(C_k)\) como:

  • La suma de las distancias euclídeas al cuadrado entre cada observación \(x_i\) y el centroide \(\mu\) de su cluster. Esto equivale a la suma de cuadrados internos del cluster. \[ W(C_k) = \sum_{x_i \in C_k} (x_i - \mu_k)^2 \]

  • La suma de las distancias euclídeas al cuadrado entre todos los pares de observaciones que forman el cluster, dividida entre el número de observaciones del cluster. \[ W(C_k) = \frac{1}{|C_k|} \sum_{i, i' \in C_k} \sum_{j = 1}^{p} (x_{ij} - x_{i' j})^2 \]

Minimizar la suma total de varianza interna \(\sum_{k=1}^{k} W(C_k)\) de forma exacta (encontrar el mínimo global) es un proceso muy complejo debido a la cantidad de formas en las que \(n\) observaciones se pueden dividir en \(K\) grupos. Sin embargo, es posible obtener una solución que, aun no siendo la mejor de entre todas las posibles, es muy buena (óptimo local). El algoritmo empleado para ello es:

  1. Asignar aleatoriamente un número entre 1 y \(K\) a cada observación. Esto sirve como asignación inicial aleatoria de las observaciones a los clusters.

  2. Iterar los siguientes pasos hasta que la asignación de las observaciones a los clusters no cambie o se alcance un número máximo de iteraciones establecido por el usuario.

  • Para cada uno de los clusters calcular su centroide. Entendiendo por centroide la posición definida por la media de cada una de las dimensiones (variables) de las observaciones que forman el cluster. Aunque no es siempre equivalente, puede entenderse como el centro de gravedad.

  • Asignar cada observación al cluster cuyo centroide está más próximo.

Para el conjunto de datos estudiado se utiliza la información recolectada al analizar el número óptimo de clusters, entendiendo que hay fuerte evidencia para suponer que existen tres clusters.

La función kmeans() del paquete stats realiza K-mean-clustering. Entre sus argumentos destacan: centers, que determina el número K de clusters que se van a generar y nstart, que determina el número de veces que se va a repetir el proceso, cada vez con una asignación aleatoria inicial distinta. Es recomendable que este último valor sea alto, entre 25-50, para no obtener resultados malos debido a una iniciación poco afortunada del proceso.

A continuación se utiliza la construcción de estos tres clusters para graficar los grupos en el sistema de ejes generados a partir del análisis de componentes principales.

# Librerías
library(factoextra)
library(gridExtra)
library(RColorBrewer)
library(ggforce)
library(ggrepel)

# Paleta
Paleta <- brewer.pal(7, "Set1")

# Modelo kmean
Data_kmeans <- kmeans(scale(Data_PCA$ind$coord[,1:2]), centers=3, nstart=30)

# Filtrar labels
Data_kmeans_clust <- cbind(Data, Data_PCA$ind$coord, Group=as.factor(Data_kmeans$cluster)) %>%
  select(Id, Dim.1, Dim.2, Group)

# Lista de países representativos por región
Representative_countries <- rbind(cbind(Id=Data$Id,
                                        IDH=Data$Human.Development.Index.HDI.2014,
                                        Region=Data$Region,
                                        data.frame(Data_PCA$ind$coord)) %>%
  group_by(Region) %>%
  slice_max(order_by=IDH, n=2) %>%
  select(Id, Dim.1, Dim.2, Region),
  cbind(Id=Data$Id,
        IDH=Data$Human.Development.Index.HDI.2014,
        Region=Data$Region,
        data.frame(Data_PCA$ind$coord)) %>%
  group_by(Region) %>%
  slice_min(order_by=IDH, n=2) %>%
  select(Id, Dim.1, Dim.2, Region),
  cbind(Id=Data$Id,
        G20=Data$G20,
        data.frame(Data_PCA$ind$coord),
        Region=Data$Region) %>%
  mutate(Id = ifelse(G20==1, Id, NA)) %>%
  select(Id, Dim.1, Dim.2, Region) %>%
  na.omit()) %>%
  distinct()

# Gráfico
ggplot(Data_kmeans_clust, aes(x=Dim.1, y=Dim.2)) +
  geom_point(aes(color=Group)) +
  geom_mark_hull(aes(fill=Group), concavity=10, radius=unit(5,"mm")) +
  geom_text_repel(data=Representative_countries, aes(x=Dim.1, y=Dim.2, label=Id)) +
  geom_vline(xintercept=0, linetype=2) + geom_hline(yintercept=0, linetype=2) +
  labs(title="PCA - K-means Clustering", x="Dim. 1 (61.6%)", y="Dim. 2 (6.8%)")

Este resultado es muy interesante porque apoya la observación que se hizo al realizar el Análisis de Componentes Principales, quizá con unas pequeñas modificaciones, pero los grupos eran evidentes y ahora son apoyados por los clusters generados, incluso con el óptimo número de clusters.

Clusters de variables

Las funciones del paquete ClustOfVar permiten realizar clusters pero utilizando las variables como elementos.

# Librería
library(ClustOfVar)
## Warning: package 'ClustOfVar' was built under R version 4.2.3
# Cluster
Treevar <- hclustvar(Data_numeric_fixed)

plot(Treevar)

Podemos observar que hay tres grandes grupos de variables.

part <- cutreevar(Treevar, 5)

part_init <- part$cluster
k.means <- kmeansvar(Data_numeric_fixed, init=part_init, matsim=TRUE)

plot(k.means)

## $coord.quanti
## $coord.quanti$Cluster1
##                                                                          dim 1
## Gender.Inequality.Index.2014                                        -0.8802759
## Change.mobile.usage.2009.2014                                       -0.5977124
## Pre.primary.2008.2014                                                0.7521271
## Physicians.per.10k.people                                            0.8684098
## Tertiary..2008.2014                                                  0.8795240
## Population.with.at.least.some.secondary.education.percent.2005.2013  0.8822518
## Secondary.2008.2014                                                  0.9031298
## Expected.years.of.schooling...Years                                  0.9175719
## Mean.years.of.schooling...Years                                      0.9287164
## Internet.users.percentage.of.population.2014                         0.9355143
## Human.Development.Index.HDI.2014                                     0.9750548
## 
## $coord.quanti$Cluster2
##                                                                      dim 1
## Life.expectancy.at.birth..years                                 -0.9257927
## Electrification.rate.or.population                              -0.9094668
## Birth.registration.funder.age.5.2005.2013                       -0.8033899
## Mobile.phone.subscriptions.per.100.people.2014                  -0.6445446
## Tuberculosis.rate.per.thousands.2012                             0.6553876
## Adolescent.birth.rate.15.19.per.100k.20102015                    0.7989590
## Primary.school.dropout.rate.2008.2014                            0.8538186
## Pupil.teacher.ratio.primary.school.pupils.per.teacher.2008.2014  0.8992574
## Maternal.mortality.ratio.deaths.per.100.live.births.2013         0.9297182
## Infant.Mortality.2013.per.thousands                              0.9447315
## Under.five.Mortality.2013.thousands                              0.9483730
## 
## $coord.quanti$Cluster3
##                                                           dim 1
## Domestic.food.price.level.2009.2014.index            -0.8302609
## Stock.of.immigrants.percentage.of.population.2013     0.6634477
## Domestic.credit.provided.by.financial.sector.2013     0.7504138
## Carbon.dioxide.emissions.per.capita.2011.Tones        0.7581560
## Research.and.development.expenditure..2005.2012       0.8100044
## Gross.domestic.product.GDP.percapta                   0.9546004
## Gross.national.income.GNI.per.capita...2011..Dollars  0.9614203
## 
## $coord.quanti$Cluster4
##                                                dim 1
## Fossil.fuels.percentage.of.total.2012      -0.961471
## Renewable.sources.percentage.of.total.2012  0.961471
## 
## $coord.quanti$Cluster5
##                                               dim 1
## International.inbound.tourists.thausands.2013     1
## 
## 
## $coord.levels
## $coord.levels$Cluster1
## [1] "No categorical variables in this cluster"
## 
## $coord.levels$Cluster2
## [1] "No categorical variables in this cluster"
## 
## $coord.levels$Cluster3
## [1] "No categorical variables in this cluster"
## 
## $coord.levels$Cluster4
## [1] "No categorical variables in this cluster"
## 
## $coord.levels$Cluster5
## [1] "No categorical variables in this cluster"

Y este resultado apoya nuestras observaciones del análisis de gradiente. Existe un cluster relacionado con variables de desarrollo humano como el IDH y la educación, otro cluster con variables relacionadas a los nacimientos/muertes, y un tercer cluster se genera con variables económicas, como el PBI, el IBN, y el gasto en investigación y desarrollo.

Métodos de aprendizaje automático

A continuación se intenta construir un modelo que prediga el resultado de una variable a partir de los datos del registro ingresado. Para esta tarea se utiliza la técnica de árboles de decisión, y su versión generalizada de Random Forest.

Muchos métodos predictivos generan modelos globales en los que una única ecuación se aplica a todo el espacio muestral. Cuando el caso de uso implica múltiples predictores, que interaccionan entre ellos de forma compleja y no lineal, es muy difícil encontrar un único modelo global que sea capaz de reflejar la relación entre las variables. Los métodos estadísticos y de machine learning basados en árboles engloban a un conjunto de técnicas supervisadas no paramétricas que consiguen segmentar el espacio de los predictores en regiones simples, dentro de las cuales es más sencillo manejar las interacciones. Es esta característica la que les proporciona gran parte de su potencial.

Árboles de decisión

Los árboles de decisión son modelos predictivos formados por reglas binarias (si/no) con las que se consigue repartir las observaciones en función de sus atributos y predecir así el valor de la variable respuesta.

Los árboles de regresión son el subtipo de árboles de predicción que se aplica cuando la variable respuesta es categórica.

Los modelos Random Forest están formados por un conjunto de árboles de decisión individuales, cada uno entrenado con una muestra ligeramente distinta de los datos de entrenamiento generada mediante bootstrapping. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.

Para construir un árbol de clasificación se emplea el método “recursive binary splitting”, sin embargo, como la variable respuesta es cualitativa, no es posible emplear el RSS como criterio de selección de las divisiones óptimas, sino que se utilizan otras alternativas, todas ellas con el objetivo de encontrar nodos lo más puros/homogéneos posible.

Podemos generar un primer árbol, para lo cual primero debemos fraccionar los datos en un subconjunto de entrenamiento y un subconjunto de prueba.

# Librerías
library(dplyr)

# Establecer semilla
set.seed(42)

# Filtrar datos
Data_filtered <- cbind(Data_numeric_fixed,
                       Income.Group=Data$Income.Group)

# Fraccionar datos de entrenamiento
Data_training <- Data_filtered %>%
  group_by(Income.Group) %>%
  mutate(num_rows=n()) %>%
  sample_frac(0.7, wight=num_rows) %>%
  select(-num_rows)

# Fraccionar datos de prueba
Data_test <- setdiff(Data_filtered, Data_training)

Y ya que tenemos los dos subconjuntos de datos, armamos el arbol de decisión.

# Librerías
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'tidyr' was built under R version 4.2.3
## Warning: package 'purrr' was built under R version 4.2.3
## Warning: package 'forcats' was built under R version 4.2.3
## Warning: package 'lubridate' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ gridExtra::combine() masks dplyr::combine()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ dplyr::lag()         masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(rpart)
library(rpart.plot)
## Warning: package 'rpart.plot' was built under R version 4.2.3
library(caret)
## Warning: package 'caret' was built under R version 4.2.3
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
# Establecer semilla
set.seed(42)

# Arbol de decisión
Modelo_arbol <- rpart(formula=Income.Group~., data=Data_training, method ='class')

# Presentar modelo
summary(Modelo_arbol)
## Call:
## rpart(formula = Income.Group ~ ., data = Data_training, method = "class")
##   n= 109 
## 
##          CP nsplit rel error    xerror       xstd
## 1 0.3076923      0 1.0000000 1.1923077 0.04736893
## 2 0.2884615      1 0.6923077 0.9743590 0.06149728
## 3 0.0100000      3 0.1153846 0.2692308 0.05278896
## 
## Variable importance
##     Gross.national.income.GNI.per.capita...2011..Dollars 
##                                                       19 
##                      Gross.domestic.product.GDP.percapta 
##                                                       19 
##                         Human.Development.Index.HDI.2014 
##                                                       17 
##             Internet.users.percentage.of.population.2014 
##                                                       11 
##                      Expected.years.of.schooling...Years 
##                                                       10 
##                Domestic.food.price.level.2009.2014.index 
##                                                        6 
## Maternal.mortality.ratio.deaths.per.100.live.births.2013 
##                                                        6 
##                       Electrification.rate.or.population 
##                                                        5 
##           Carbon.dioxide.emissions.per.capita.2011.Tones 
##                                                        4 
##          Research.and.development.expenditure..2005.2012 
##                                                        4 
## 
## Node number 1: 109 observations,    complexity param=0.3076923
##   predicted class=Lower middle income  expected loss=0.7155963  P(node) =1
##     class counts:    27    20    31    31
##    probabilities: 0.248 0.183 0.284 0.284 
##   left son=2 (24 obs) right son=3 (85 obs)
##   Primary splits:
##       Human.Development.Index.HDI.2014                     < 0.8181529 to the right, improve=23.43270, (0 missing)
##       Gross.national.income.GNI.per.capita...2011..Dollars < 21173.05  to the right, improve=22.20791, (0 missing)
##       Gross.domestic.product.GDP.percapta                  < 9265.098  to the right, improve=21.64271, (0 missing)
##       Internet.users.percentage.of.population.2014         < 61.23     to the right, improve=21.54420, (0 missing)
##       Expected.years.of.schooling...Years                  < 15.21803  to the right, improve=20.95542, (0 missing)
##   Surrogate splits:
##       Domestic.food.price.level.2009.2014.index            < 2.885     to the left,  agree=0.972, adj=0.875, (0 split)
##       Gross.national.income.GNI.per.capita...2011..Dollars < 21173.05  to the right, agree=0.972, adj=0.875, (0 split)
##       Internet.users.percentage.of.population.2014         < 71.055    to the right, agree=0.972, adj=0.875, (0 split)
##       Expected.years.of.schooling...Years                  < 15.21803  to the right, agree=0.963, adj=0.833, (0 split)
##       Gross.domestic.product.GDP.percapta                  < 20888.58  to the right, agree=0.963, adj=0.833, (0 split)
## 
## Node number 2: 24 observations
##   predicted class=High income          expected loss=0  P(node) =0.2201835
##     class counts:    24     0     0     0
##    probabilities: 1.000 0.000 0.000 0.000 
## 
## Node number 3: 85 observations,    complexity param=0.2884615
##   predicted class=Lower middle income  expected loss=0.6352941  P(node) =0.7798165
##     class counts:     3    20    31    31
##    probabilities: 0.035 0.235 0.365 0.365 
##   left son=6 (20 obs) right son=7 (65 obs)
##   Primary splits:
##       Gross.national.income.GNI.per.capita...2011..Dollars     < 2422.489  to the left,  improve=22.28416, (0 missing)
##       Gross.domestic.product.GDP.percapta                      < 2302.225  to the left,  improve=20.89046, (0 missing)
##       Human.Development.Index.HDI.2014                         < 0.5511474 to the left,  improve=17.76226, (0 missing)
##       Maternal.mortality.ratio.deaths.per.100.live.births.2013 < 305       to the right, improve=17.25180, (0 missing)
##       Carbon.dioxide.emissions.per.capita.2011.Tones           < 1.401049  to the left,  improve=14.60100, (0 missing)
##   Surrogate splits:
##       Gross.domestic.product.GDP.percapta                      < 2302.225  to the left,  agree=0.988, adj=0.95, (0 split)
##       Human.Development.Index.HDI.2014                         < 0.5511474 to the left,  agree=0.953, adj=0.80, (0 split)
##       Maternal.mortality.ratio.deaths.per.100.live.births.2013 < 305       to the right, agree=0.953, adj=0.80, (0 split)
##       Internet.users.percentage.of.population.2014             < 8         to the left,  agree=0.929, adj=0.70, (0 split)
##       Electrification.rate.or.population                       < 45.15     to the left,  agree=0.918, adj=0.65, (0 split)
## 
## Node number 6: 20 observations
##   predicted class=Low income           expected loss=0  P(node) =0.1834862
##     class counts:     0    20     0     0
##    probabilities: 0.000 1.000 0.000 0.000 
## 
## Node number 7: 65 observations,    complexity param=0.2884615
##   predicted class=Lower middle income  expected loss=0.5230769  P(node) =0.5963303
##     class counts:     3     0    31    31
##    probabilities: 0.046 0.000 0.477 0.477 
##   left son=14 (31 obs) right son=15 (34 obs)
##   Primary splits:
##       Gross.domestic.product.GDP.percapta                  < 9265.098  to the left,  improve=19.46119, (0 missing)
##       Gross.national.income.GNI.per.capita...2011..Dollars < 11702.73  to the left,  improve=18.53792, (0 missing)
##       Domestic.food.price.level.2009.2014.index            < 3.848333  to the right, improve=13.03675, (0 missing)
##       Human.Development.Index.HDI.2014                     < 0.6956107 to the left,  improve=12.55897, (0 missing)
##       Carbon.dioxide.emissions.per.capita.2011.Tones       < 1.246358  to the left,  improve=11.58397, (0 missing)
##   Surrogate splits:
##       Gross.national.income.GNI.per.capita...2011..Dollars < 8797.936  to the left,  agree=0.969, adj=0.935, (0 split)
##       Human.Development.Index.HDI.2014                     < 0.6794895 to the left,  agree=0.815, adj=0.613, (0 split)
##       Carbon.dioxide.emissions.per.capita.2011.Tones       < 2.101573  to the left,  agree=0.815, adj=0.613, (0 split)
##       Expected.years.of.schooling...Years                  < 12.33248  to the left,  agree=0.800, adj=0.581, (0 split)
##       Research.and.development.expenditure..2005.2012      < 0.4058467 to the left,  agree=0.800, adj=0.581, (0 split)
## 
## Node number 14: 31 observations
##   predicted class=Lower middle income  expected loss=0.09677419  P(node) =0.2844037
##     class counts:     0     0    28     3
##    probabilities: 0.000 0.000 0.903 0.097 
## 
## Node number 15: 34 observations
##   predicted class=Upper middle income  expected loss=0.1764706  P(node) =0.3119266
##     class counts:     3     0     3    28
##    probabilities: 0.088 0.000 0.088 0.824
rpart.plot(Modelo_arbol)

Y ahora podemos realizar una predicción con este modelo generado.

# Predicción
Modelo_predict <- predict(Modelo_arbol, Data_test, type="class")

# Matriz de confusión
confusionMatrix(Modelo_predict, as.factor(Data_test[["Income.Group"]]))
## Confusion Matrix and Statistics
## 
##                      Reference
## Prediction            High income Low income Lower middle income
##   High income                  11          0                   0
##   Low income                    0          8                   1
##   Lower middle income           0          0                  12
##   Upper middle income           1          0                   1
##                      Reference
## Prediction            Upper middle income
##   High income                           1
##   Low income                            0
##   Lower middle income                   4
##   Upper middle income                   9
## 
## Overall Statistics
##                                           
##                Accuracy : 0.8333          
##                  95% CI : (0.6978, 0.9252)
##     No Information Rate : 0.2917          
##     P-Value [Acc > NIR] : 1.023e-14       
##                                           
##                   Kappa : 0.7754          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: High income Class: Low income
## Sensitivity                      0.9167            1.0000
## Specificity                      0.9722            0.9750
## Pos Pred Value                   0.9167            0.8889
## Neg Pred Value                   0.9722            1.0000
## Prevalence                       0.2500            0.1667
## Detection Rate                   0.2292            0.1667
## Detection Prevalence             0.2500            0.1875
## Balanced Accuracy                0.9444            0.9875
##                      Class: Lower middle income Class: Upper middle income
## Sensitivity                              0.8571                     0.6429
## Specificity                              0.8824                     0.9412
## Pos Pred Value                           0.7500                     0.8182
## Neg Pred Value                           0.9375                     0.8649
## Prevalence                               0.2917                     0.2917
## Detection Rate                           0.2500                     0.1875
## Detection Prevalence                     0.3333                     0.2292
## Balanced Accuracy                        0.8697                     0.7920

Obtenemos una precisión del 83%.

Random Forest

Un modelo Random Forest está formado por un conjunto (ensemble) de árboles de decisión individuales, cada uno entrenado con una muestra aleatoria extraída de los datos de entrenamiento originales mediante bootstrapping. Esto implica que cada árbol se entrena con unos datos ligeramente distintos. En cada árbol individual, las observaciones se van distribuyendo por bifurcaciones (nodos) generando la estructura del árbol hasta alcanzar un nodo terminal. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.

El algoritmo de Random Forest es una modificación del proceso de bagging que consigue mejorar los resultados gracias a que decorrelaciona aún más los árboles generados en el proceso.

# Librerías
library(randomForest)
## Warning: package 'randomForest' was built under R version 4.2.3
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:gridExtra':
## 
##     combine
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:dplyr':
## 
##     combine
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
# Establecer semilla
set.seed(43)

# Filtrar datos
Data_filtered <- cbind(Data_numeric_fixed,
                       Income.Group=as.factor(Data$Income.Group))

# Fraccionar datos de entrenamiento
Sample_training <- sample(1:nrow(Data_filtered), 0.7*(nrow(Data_filtered)))

# Random forest
Modelo_randomforest <- randomForest(Income.Group ~ ., data=Data_filtered, subset=Sample_training)

# Presentar modelo
plot(Modelo_randomforest)

Y ahora realizar una predicción.

Random_forest_predict <- predict(Modelo_randomforest, Data_filtered[-Sample_training,])

confusionMatrix(Random_forest_predict, as.factor(Data_filtered[-Sample_training,][["Income.Group"]]))
## Confusion Matrix and Statistics
## 
##                      Reference
## Prediction            High income Low income Lower middle income
##   High income                   9          0                   0
##   Low income                    0          9                   1
##   Lower middle income           0          2                  12
##   Upper middle income           0          0                   2
##                      Reference
## Prediction            Upper middle income
##   High income                           4
##   Low income                            0
##   Lower middle income                   4
##   Upper middle income                   5
## 
## Overall Statistics
##                                           
##                Accuracy : 0.7292          
##                  95% CI : (0.5815, 0.8472)
##     No Information Rate : 0.3125          
##     P-Value [Acc > NIR] : 3.681e-09       
##                                           
##                   Kappa : 0.6364          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: High income Class: Low income
## Sensitivity                      1.0000            0.8182
## Specificity                      0.8974            0.9730
## Pos Pred Value                   0.6923            0.9000
## Neg Pred Value                   1.0000            0.9474
## Prevalence                       0.1875            0.2292
## Detection Rate                   0.1875            0.1875
## Detection Prevalence             0.2708            0.2083
## Balanced Accuracy                0.9487            0.8956
##                      Class: Lower middle income Class: Upper middle income
## Sensitivity                              0.8000                     0.3846
## Specificity                              0.8182                     0.9429
## Pos Pred Value                           0.6667                     0.7143
## Neg Pred Value                           0.9000                     0.8049
## Prevalence                               0.3125                     0.2708
## Detection Rate                           0.2500                     0.1042
## Detection Prevalence                     0.3750                     0.1458
## Balanced Accuracy                        0.8091                     0.6637

Obteniendo una precisión de 73%.

varImpPlot(Modelo_randomforest)

Conclusiones

A lo largo de este trabajo, se exploraron diversas técnicas estadísticas para analizar la estructura de los datos. En primer lugar, se aplicaron técnicas de análisis de componentes principales (PCA) utilizando variables numéricas exclusivamente. Este enfoque buscaba reducir la dimensionalidad de los datos y reveló tres grupos distintos al proyectar los individuos en un nuevo sistema de ejes coordenados. Estos grupos se asociaron como países africanos y de Medio Oriente (considerados subdesarrollados), países asiáticos y americanos (denominados emergentes), y países europeos y norteamericanos (identificados como desarrollados). A pesar de la subjetividad inherente a este tipo de agrupamientos, la representación en el nuevo plano proporcionó información valiosa respaldada por medidas de distancia utilizadas en la generación de los componentes principales.

A continuación, se emplearon variables categóricas dicotómicas y, posteriormente, categóricas para realizar un análisis de coordenadas principales (PCoA), un enfoque análogo al PCA que utiliza distintas medidas de distancias. Sin embargo, estos estudios no ofrecieron información tan importante como el PCA. Además, todas estas variables se utilizaron en un Análisis Factorial de Correspondencia (AFC), que se centra en los factores o estados de las variables categóricas en lugar de los individuos. Este análisis reveló una distribución definida en tres grandes ramas, destacando características distintivas de nuevamente tres grupos: países americanos, países africanos y árabes, y países europeos y norteamericanos. Los factores menos definidos o más prominentes se ubicaron en el centro de coordenadas.

Con la información obtenida de estos análisis, se aplicaron técnicas de clasificación, como clustering jerárquico y k-means, para consolidar la observación de tres grupos predominantes en los datos. Estos grupos coincidieron en su mayoría con la clasificación previa de países subdesarrollados, emergentes y desarrollados.

Finalmente, se utilizaron métodos de aprendizaje automático, específicamente la construcción de modelos de árboles de decisión y random forest, para predecir la variable categórica “Grupo de Ingreso” (alto, medio-alto, medio-bajo, bajo) a partir de las variables numéricas del conjunto de datos. Estos modelos exhibieron una notable precisión, cerrando así la investigación con un análisis predictivo efectivo.

La existencia de estos tres grupos de países, es decir, desarrollados, emergentes y subdesarrollados, refleja las disparidades económicas y sociales que existen a nivel global. Estas categorizaciones se basan en indicadores económicos, sociales y de desarrollo humano que permiten clasificar a los países según su nivel de desarrollo relativo. Los países desarrollados generalmente exhiben altos niveles de ingresos, infraestructuras avanzadas y estándares de vida elevados. Los países emergentes están en proceso de industrialización y experimentan un crecimiento económico significativo, aunque aún enfrentan desafíos en términos de desarrollo humano. Por otro lado, los países subdesarrollados confrontan obstáculos considerables en áreas como la pobreza, la salud y la educación.

Esta clasificación, sin embargo, no debe interpretarse como estática ni como una descripción exhaustiva de la complejidad de las realidades nacionales. Existen variaciones significativas dentro de cada categoría, y algunos países han logrado transiciones exitosas de un grupo a otro a lo largo del tiempo. Además, la noción de desarrollo es multidimensional y va más allá de los indicadores económicos, incluyendo aspectos como la equidad, la sostenibilidad ambiental y la calidad de vida. En última instancia, estas categorías son herramientas útiles para comprender las tendencias generales, pero no deben utilizarse de manera rígida al considerar la diversidad y dinámica de las naciones a lo largo del tiempo.