Lo primero que hacemos es cargar las librerias necesarias para manipular los datos como es: dplyr. Otra llamada ‘haven’ para leer archivos SAS que es el formato donde tenemos los datos de PISA. y luego para el PCA utilizamos Factorminer, factorExtra y missMDA.

# Cargar las librerías necesarias
library(FactoMineR)      # Para análisis de componentes principales (PCA)
library(factoextra)      # Para visualizar resultados del PCA
## Cargando paquete requerido: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(dplyr)           # Para manipulación de datos
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(haven)           # Para leer archivos SAS
library(knitr)           # para tablas y visualizacion

Luego cargamos los datos de escuelas, estudiantes y los del banco mundial

# Leer datos de escuelas desde archivo SAS
escuelas <- read_sas("escuelas.SAS7BDAT", NULL)

# Cargar datos de estudiantes de archivo RData y asignar al objeto 'estudiantes'
estudiantes = load('estudiantesPISAfiltrados.R')
estudiantes <- datosPisaEstudiantes_filtrado

# Cargar datos del Banco Mundial de archivo RData y asignar al objeto 'bancomundial'
bancomundial = load('bancomundial.R')
bancomundial <- datos_filtrado

De la base de datos de estudiantes que es donde estan las notas, obtener las columnas de estas y las guardamos en una variable. Por otro lado calculamos la media de cada conocimiento, es decir, matematicas, ciencia y lectura. Agrupamos por pais y calculamos la media del pais

# Filtrar columnas que contienen las notas plausibles (PV), el país (CNT) y el ID de escuela (CNTSCHID)
notas = estudiantes %>%
  select(matches("^PV|^CNT$|^CNTSCHID$"))

# Calcular la media de las notas  por estudiante en matemáticas, lectura y ciencias
NotasMedias = notas %>% transmute(
  CNT = notas$CNT,
  CNTSCHID = notas$CNTSCHID,
  MathMean = rowMeans(notas[, 3:12]),          # Promedio de PVs de matemáticas
  LectureMean = rowMeans(notas[, 13:22]),      # Promedio de PVs de lectura
  ScienceMean = rowMeans(notas[, 23:32])       # Promedio de PVs de ciencias
)

# Agrupar por país y calcular la media nacional en cada área
NotasMediasPais = NotasMedias %>%
  group_by(CNT) %>%
  summarise(
    MathMean_Country = mean(MathMean, na.rm = TRUE),
    LectureMean_Country = mean(LectureMean, na.rm = TRUE),
    ScienceMean_Country = mean(ScienceMean, na.rm = TRUE)
  )

Ahora unimos las dos bases de datos, lo que tenemos de las notas y la base de datos del banco mundial por el codigo del pais.

# Unir las medias de desempeño por país con los datos del Banco Mundial
# Usar la columna de códigos de país (CNT y Country.Code)


bm_media_res <- merge(
  NotasMediasPais, 
  bancomundial, 
  by.x = "CNT",               # Columna de país en los datos PISA
  by.y = "Country.Code",      # Columna de país en los datos del Banco Mundial
  all.x = TRUE                # Conservar todos los países de PISA aunque falten datos del BM
)

# Exportar el resultado combinado a un archivo CSV
#write.csv(bm_media_res, file = "bm_media_res.csv", row.names = FALSE)

# Instalar y cargar paquete para leer archivos Excel
library(readxl)
## Warning: package 'readxl' was built under R version 4.4.3
# Leer archivo Excel que contiene solo variables numéricas
datos_numericos <- read_excel("datos_bm_res.xlsx")

# Ver los datos numéricos cargados
View(datos_numericos)

# Quitar la columna CNT (nombre o código del país)
datos = subset(datos_numericos, select = -c(CNT))
View(datos)
#summary(datos)  # Ver resumen estadístico

Seleccionamos las columnas que no queremos meter en el PCA, como país y calculamos la media conjunta de los resultados.

# Calcular el promedio de desempeño en las tres áreas y guardarlo en nueva columna
datos$PromedioAreas = rowMeans(datos[, c("MathMean_Country", "LectureMean_Country", "ScienceMean_Country")], na.rm = TRUE)

# Eliminar las columnas individuales de desempeño
datos_media = subset(datos, select = -c(MathMean_Country, LectureMean_Country, ScienceMean_Country))
View(datos_media)

# Realizar Análisis de Componentes Principales (PCA) con escalado de variables
res.pca = PCA(datos_media, scale.unit = TRUE, graph = FALSE, ncp = 10)
## Warning in PCA(datos_media, scale.unit = TRUE, graph = FALSE, ncp = 10):
## Missing values are imputed by the mean of the variable: you should use the
## imputePCA function of the missMDA package
# Obtener los valores propios (eigenvalues)
eig.val <- get_eigenvalue(res.pca)

# Calcular el porcentaje de varianza explicada por componente si todas fueran iguales
VPmedio = 100 * (1 / nrow(eig.val))

# Visualizar la varianza explicada por componente con línea roja de referencia
fviz_eig(res.pca, addlabels = TRUE) +
  geom_hline(yintercept = VPmedio, linetype = 2, color = "red")

kable(eig.val[1:9,])
eigenvalue variance.percent cumulative.variance.percent
Dim.1 10.634092 25.936809 25.93681
Dim.2 6.429081 15.680686 41.61749
Dim.3 4.192929 10.226655 51.84415
Dim.4 2.478459 6.045021 57.88917
Dim.5 2.374284 5.790938 63.68011
Dim.6 2.123457 5.179163 68.85927
Dim.7 1.688118 4.117361 72.97663
Dim.8 1.444185 3.522404 76.49904
Dim.9 1.039050 2.534268 79.03330

Tras observar la grafica de varianza explicada, dibujar la linea roja que representa la varianza media esperada si todos los componentes fueran igual de importantes. Nos dice el número óptimo de clúster según el criterio de la media,la tabla de la la varianza acumulativa segun aumentan las dimensiones y el creiterio del codo son 4 componentes principales, que representan el 58% de la variabilidad.

Una vez seleccionadas las 4 componentes principales (K = 4) mediante PCA, se calcula para cada país su posición en el nuevo espacio reducido (misScores). A partir de estas coordenadas, se calcula el estadístico Hotelling T² (miT2), que mide cuán lejos está cada país del centro del conjunto de datos en el espacio de componentes principales. Este valor se compara con dos umbrales teóricos de control (F95 y F99), que corresponden a los niveles de confianza del 95% y 99%. Así dectetamos los paises que son anómalos

K = 4
res.pca = PCA(datos_media, scale.unit = TRUE, graph = FALSE, ncp = K, quali.sup = 1)
## Warning in PCA(datos_media, scale.unit = TRUE, graph = FALSE, ncp = K,
## quali.sup = 1): Missing values are imputed by the mean of the variable: you
## should use the imputePCA function of the missMDA package
misScores = res.pca$ind$coord[,1:K]
miT2 = colSums(t(misScores**2)/eig.val[1:K,1])
I = nrow(datos_numericos)
F95 = K*(I**2 - 1)/(I*(I - K)) * qf(0.95, K, I-K)
F99 = K*(I**2 - 1)/(I*(I - K)) * qf(0.99, K, I-K)

plot(1:length(miT2), miT2, type = "p", xlab = 'datos', ylab = "T2")
abline(h = F95, col = "orange", lty = 2, lwd = 2)
abline(h = F99, col = "red3", lty = 2, lwd = 2)

anomalas = which(miT2 > F99)
anomalas
## 29 56 75 
## 29 56 75
datos_numericos[anomalas, "CNT"]

Primero hacemos el PCA con los anólamos, para ver qué influye en que estos lo sean.

library(ggrepel)
library(ggplot2)

etiquetas_paises <- datos_numericos$CNT
media_variable <- rowMeans(datos[, 1:3], na.rm = TRUE)
# Paso 1: Extraer coordenadas de individuos del objeto PCA
pca_df <- as.data.frame(res.pca$ind$coord)

# Paso 2: Agregar nombres de países y grupo según T2 de Hotelling
pca_df$Pais <- etiquetas_paises  # Asegúrate que tenga longitud igual a nrow(pca_df)
pca_df$Grupo <- factor(miT2 > F99)

# Paso 3: Gráfico Dim1 vs Dim2
ggplot(pca_df, aes(x = Dim.1, y = Dim.2, label = Pais, color = Grupo)) +
  geom_text_repel(size = 3, show.legend = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +  # Eje horizontal
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray40") + 
  labs(title = "Individuals - PCA (Dim1 vs Dim2)", x = "Dim1", y = "Dim2") +
  theme_minimal()
## Warning: ggrepel: 39 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

# Paso 4: Gráfico Dim1 vs Dim3
ggplot(pca_df, aes(x = Dim.3, y = Dim.4, label = Pais, color = Grupo)) +
  geom_text_repel(size = 3, show.legend = TRUE) +
  labs(title = "Individuals - PCA (Dim3 vs Dim4)", x = "Dim3", y = "Dim4") +
  theme_minimal()
## Warning: ggrepel: 22 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

Para ver porqué estos son anómalos veremos las contibuciones de la componente 1 para Filipinas e India, y la 2 para Estados Unidos.

Estados Unidos destaca claramente en la segunda componente, principalmente por su alto nivel de desarrollo tecnológico y económico: tiene muchos servidores seguros, alta conectividad, un PIB per cápita elevado y una gran población, lo que lo posiciona como un caso atípico en el análisis.

Por otro lado, India y Filipinas se ubican en el extremo opuesto de la primera componente, influenciada por factores como ingreso, educación y acceso a Internet. Ambos países presentan niveles bajos en estos indicadores y una gran parte de su población trabaja en la agricultura, lo que refleja una economía menos desarrollada. Estas diferencias explican por qué se alejan tanto del resto de países en el gráfico.

library(gridExtra)
## 
## Adjuntando el paquete: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
# Para la Componente 1
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)  # para hablar de phl e ind

# Extraer contribuciones a la Dimensión 1
contribuciones <- get_pca_var(res.pca)$contrib[,1]  # 1 = Dim-1

# Convertir a data frame y ordenar
contrib_df <- data.frame(
  variable = names(contribuciones),
  contrib = contribuciones
)

# Ordenar de mayor a menor contribución
contrib_df_ordenado <- contrib_df[order(-contrib_df$contrib), ]

# Ver las 10 que más contribuyen
head(contrib_df_ordenado, 10)
# Para la Componente 2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10) # para hablar de eeuu

# Extraer contribuciones a la Dimensión 2
contrib_dim2 <- get_pca_var(res.pca)$contrib[, 2]

# Convertir a data frame
contrib_df_dim2 <- data.frame(
  variable = names(contrib_dim2),
  contrib = contrib_dim2
)

# Ordenar de mayor a menor contribución
top10_dim2 <- contrib_df_dim2[order(-contrib_df_dim2$contrib), ][1:10, ]

# Mostrar
top10_dim2

Ahora, eliminamos estos paises y rehacemos el PCA. Representamos los graficos de variables y de individuos.

etiquetas_paises <- datos_numericos$CNT
media_variable <- rowMeans(datos[, 1:3], na.rm = TRUE)
media_intervalos <- cut(media_variable, 
                        breaks = seq(0, max(media_variable, na.rm = TRUE) + 100, by = 100),
                        include.lowest = TRUE, right = FALSE)

datos_filtrados <- datos_media[miT2 <= F99, ]
etiquetas_filtradas <- etiquetas_paises[miT2 <= F99]
media_intervalos_filtrado <- media_intervalos[miT2 <= F99]

# Recalcular el PCA con los datos filtrados
res.pca_filtrado <- PCA(datos_filtrados, scale.unit = TRUE, graph = FALSE, ncp = 10, 
              quali.sup = 1)
## Warning in PCA(datos_filtrados, scale.unit = TRUE, graph = FALSE, ncp = 10, :
## Missing values are imputed by the mean of the variable: you should use the
## imputePCA function of the missMDA package
fviz_pca_var(res.pca_filtrado, axes = c(1,2), repel = TRUE, col.var = "contrib", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

Volvemos a hacer el de variables y para que se vea más claro filtramos por las que más contribuyen

contrib <- res.pca_filtrado$var$contrib[, 1:2]  # contribuciones a Dim 1 y Dim 2
contrib_total <- rowSums(contrib)     # suma de la contribución en ambos ejes
variables_filtradas <- which(contrib_total > 5)
fviz_pca_var(res.pca_filtrado,
             axes = c(1, 2),
             select.var = list(name = rownames(res.pca_filtrado$var$coord)[variables_filtradas]),
             repel = TRUE,
             col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

# 1. Extraemos los nombres reales de las variables seleccionadas
variable_names <- rownames(res.pca_filtrado$var$coord)[variables_filtradas]

# 2. Creamos etiquetas numéricas para reemplazar los nombres
etiquetas_numericas <- as.character(seq_along(variable_names))

# 3. Reemplazamos los nombres en una copia del objeto PCA
res.pca.mod <- res.pca_filtrado
rownames(res.pca.mod$var$coord)[variables_filtradas] <- etiquetas_numericas

# 4. Graficamos con etiquetas numéricas
fviz_pca_var(res.pca.mod,
             axes = c(1, 2),
             select.var = list(name = etiquetas_numericas),
             repel = TRUE,
             col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

# 5. Mostramos la leyenda como tabla en consola
leyenda <- data.frame(Número = etiquetas_numericas,
                      Variable = variable_names)
print(leyenda)
##    Número
## 1       1
## 2       2
## 3       3
## 4       4
## 5       5
## 6       6
## 7       7
## 8       8
## 9       9
## 10     10
## 11     11
## 12     12
## 13     13
## 14     14
##                                                                                              Variable
## 1  Educational attainment, at least completed lower secondary, population 25+, total (%) (cumulative)
## 2    Educational attainment, at least completed primary, population 25+ years, total (%) (cumulative)
## 3  Educational attainment, at least completed upper secondary, population 25+, total (%) (cumulative)
## 4                                                          GNI per capita, Atlas method (current US$)
## 5                                                 GNI per capita, PPP (constant 2021 international $)
## 6                                                                       Fixed broadband subscriptions
## 7                                                                       Mobile cellular subscriptions
## 8                                                                             Secure Internet servers
## 9                                                                  GDP per capita (constant 2015 US$)
## 10                                                                        Population ages 0-14, total
## 11                                                                       Population ages 15-64, total
## 12                                                                Population ages 65 and above, total
## 13                                                                                  Population, total
## 14                                                                                      PromedioAreas
pca_data <- as.data.frame(res.pca_filtrado$ind$coord[, 1:2])
colnames(pca_data) <- c("Dim1", "Dim2")
pca_data$label <- etiquetas_filtradas
pca_data$intervalo <- as.factor(media_intervalos_filtrado)

ggplot(pca_data, aes(x = Dim1, y = Dim2, label = label, color = intervalo)) +
  geom_text(size = 3) +
  scale_color_manual(values = c("#1E90FF", "#32CD32", "#DC143C")) +
  labs(title = "PCA - Dimensiones 1 y 2", color = "Intervalo de Media") +
  theme_minimal()

También vemos las variables que más importan en cada componente:

# Para la Componente 1 (res.pca_filtrado)
fviz_contrib(res.pca_filtrado, choice = "var", axes = 1, top = 10)  # top=10 es más informativo

# Extraer contribuciones a la Dimensión 1
contribuciones_filtrado_dim1 <- get_pca_var(res.pca_filtrado)$contrib[, 1]

# Crear y ordenar el data frame
contrib_df_dim1_filtrado <- data.frame(
  variable = names(contribuciones_filtrado_dim1),
  contrib = contribuciones_filtrado_dim1
)
contrib_df_dim1_ordenado <- contrib_df_dim1_filtrado[order(-contrib_df_dim1_filtrado$contrib), ]

# Mostrar las 10 más importantes
head(contrib_df_dim1_ordenado, 10)
# Para la Componente 2 (res.pca_filtrado)
fviz_contrib(res.pca_filtrado, choice = "var", axes = 2, top = 10)

# Extraer contribuciones a la Dimensión 2
contribuciones_filtrado_dim2 <- get_pca_var(res.pca_filtrado)$contrib[, 2]

# Crear y ordenar el data frame
contrib_df_dim2_filtrado <- data.frame(
  variable = names(contribuciones_filtrado_dim2),
  contrib = contribuciones_filtrado_dim2
)
contrib_df_dim2_ordenado <- contrib_df_dim2_filtrado[order(-contrib_df_dim2_filtrado$contrib), ]

# Mostrar las 10 más importantes
head(contrib_df_dim2_ordenado, 10)

Al analizar las variables que más influyen en cada componente del PCA, vemos que la primera dimensión (Dim1) está muy relacionada con el desarrollo económico y educativo. Ahí destacan variables como el PIB per cápita, el GNI per cápita, los niveles de educación alcanzados por la población adulta y el uso de Internet. Básicamente, esta dimensión nos habla del nivel de desarrollo de los países.

La segunda dimensión (Dim2), en cambio, se asocia más con aspectos demográficos y de infraestructura digital. Aquí influyen cosas como la cantidad total de población, la proporción de jóvenes y adultos mayores, las suscripciones móviles y fijas, servidores seguros, y el gasto en educación. Esta dimensión captura cómo están conectados los países y cómo está compuesta su población.

En el gráfico de individuos, donde ya se eliminaron los casos atípicos, los países están coloreados según su media general: azul para medias entre 300 y 400, verde entre 400 y 500, y rojo entre 500 y 600. Ojo, estos colores no representan clústeres, solo rangos de valores.

Los países en azul —como El Salvador, República Dominicana, Panamá o Jamaica— están hacia la izquierda del gráfico. Eso indica que tienen niveles bajos en desarrollo económico y educativo. Algunos también están en la parte baja del gráfico, lo que sugiere menos infraestructura tecnológica o poblaciones más jóvenes.

Los países en verde —como México, Turquía, Colombia, Malasia o Rumanía— están más cerca del centro. Representan niveles intermedios en desarrollo y conectividad.

Y por último, los países en rojo —como Alemania, Suiza, Noruega, Irlanda, Finlandia o Singapur— están a la derecha. Tienen economías sólidas, buen acceso a tecnología y altos niveles educativos. Están cerca del eje horizontal, lo que indica una estructura demográfica más equilibrada.

Este análisis nos da un panorama bastante claro de cómo se posicionan los países según sus características clave. El siguiente paso será aplicar un algoritmo de clustering para ver si estas diferencias también se traducen en grupos definidos de países.