Creación de un nuevo dataset con las nuevas variables retenidas y realizar un análisis estadístico descriptivo de las dos variables principales obtenidas.

# NUEVO DATASET CON LOS DOS COMPONENTES RETENIDOS

# Extraer puntuaciones de los dos primeros componentes
componentes_2011 <- as.data.frame(acp_2011$x[, 1:2])

# Renombrar variables para facilitar la interpretación
colnames(componentes_2011) <- c("PC1_Envejecimiento", "PC2_Extranjeria")

# Añadir identificador de sección censal
componentes_2011$COD_SECC <- datos_2011$COD_SECC

# Reordenar columnas
componentes_2011 <- componentes_2011[, c("COD_SECC", 
                                         "PC1_Envejecimiento", 
                                         "PC2_Extranjeria")]

# Mostrar primeras filas
formattable(head(componentes_2011))

En la anterior tabla puede verse una muestra con las primeras seis secciones censales del nuevo dataset generado con las dos variables basasdas en los dos componentes principales extraídos del dataset original.

Las puntuaciones de las seis primeras secciones censales en los dos componentes principales permiten caracterizar sus perfiles sociodemográficos relativos. El primer componente representa principalmente una dimensión de envejecimiento y hogares unipersonales de personas mayores, mientras que el segundo componente se asocia a la presencia de población extranjera, especialmente extranjera infantil. Las secciones 3907501002, 3907501003 y especialmente 3907501004 presentan puntuaciones positivas en el primer componente y negativas en el segundo, lo que sugiere un perfil más envejecido y con menor presencia extranjera. Por el contrario, la sección 3907501006 presenta baja puntuación en envejecimiento y alta puntuación en extranjería. La sección 3907501005 destaca por obtener puntuaciones elevadas en ambos componentes, lo que apunta a un perfil más complejo, con coexistencia de envejecimiento y presencia extranjera. Por tanto, puede decirse que. en conjunto, las puntuaciones de los componentes permiten sintetizar y comparar los perfiles de las secciones censales de forma más clara que las variables originales.

Seguidamente, se realiza un análisis descriptivo de las dos nuevas variables que componen el nuevo dataset.

Análisis descriptivo de las dos nuevas variables

# Seleccionar solo los componentes
componentes_numericos <- componentes_2011[, c("PC1_Envejecimiento", 
                                              "PC2_Extranjeria")]

# Tabla descriptiva
tabla_desc_componentes <- data.frame(
  Variable = names(componentes_numericos),
  Media = round(sapply(componentes_numericos, mean), 4),
  Desv_tipica = round(sapply(componentes_numericos, sd), 4),
  Minimo = round(sapply(componentes_numericos, min), 4),
  Q1 = round(sapply(componentes_numericos, quantile, probs = 0.25), 4),
  Mediana = round(sapply(componentes_numericos, median), 4),
  Q3 = round(sapply(componentes_numericos, quantile, probs = 0.75), 4),
  Maximo = round(sapply(componentes_numericos, max), 4)
)

# Eliminar nombres de fila para evitar que se repitan las variables
rownames(tabla_desc_componentes) <- NULL

kable(
  tabla_desc_componentes,
  caption = "Análisis descriptivo de los dos componentes principales retenidos",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

# Histogramas
par(mfrow = c(2, 3))
for (v in names(componentes_numericos)) {
  hist(componentes_numericos[[v]],
       main = paste("Histograma de", v),
       xlab = v,
       col = "blue",
       border = "white")
}
par(mfrow = c(1, 1))

# Boxplots para detectar dispersión y posibles valores extremos
boxplot(componentes_numericos,
        main = "Distribución de variables sociodemográficas 2011",
        las = 2,
        col = "blue")

El análisis descriptivo de los dos componentes principales muestra que ambas variables tienen medias próximas a cero, como consecuencia de haberse obtenido a partir de variables estandarizadas. El primer componente, asociado al envejecimiento y a los hogares unipersonales de personas mayores, presenta una mayor dispersión que el segundo, lo que es coherente con su mayor capacidad explicativa dentro del ACP. Su distribución es relativamente equilibrada, aunque con secciones situadas en ambos extremos del eje de envejecimiento. En cambio, el segundo componente, vinculado a la presencia de población extranjera y extranjera infantil, muestra una distribución más asimétrica hacia la derecha, con muchas secciones en valores bajos o moderados y algunos casos con puntuaciones elevadas. Esto sugiere que la dimensión de extranjería está más concentrada territorialmente en un número menor de secciones censales.

En conjunto, los dos componentes permiten sintetizar dos patrones diferenciados: uno más extendido y gradual, vinculado al envejecimiento, y otro más concentrado espacialmente, vinculado a la presencia de población extranjera.

Correlación entre los dos componentes

# Matriz de correlación entre los componentes retenidos
cor_componentes <- cor(componentes_numericos) 

# Redondeo de los componentes para simplificar
tabla_cor_componentes <- round(cor_componentes, 4)

# Generamos una tabla con los resultados de correlación
kable(
  tabla_cor_componentes,
  caption = "Matriz de correlación entre los componentes principales retenidos",
  align = "c"
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

La matriz de correlación entre los dos componentes principales retenidos muestra una correlación nula entre PC1_Envejecimiento y PC2_Extranjeria. Este resultado es coherente con la lógica del análisis de componentes principales, ya que los componentes obtenidos son ortogonales por construcción. Por tanto, las dos nuevas variables sintetizan dimensiones independientes de la variabilidad del dataset: una primera dimensión asociada al envejecimiento y los hogares unipersonales de personas mayores, y una segunda dimensión vinculada a la presencia de población extranjera y extranjera infantil.

Análisis Factorial Exploratorio

Se realiza un análisis factorial exploratorio como una técnica de reducción de la dimensionalidad orientada a identificar posibles factores latentes que expliquen las relaciones observadas entre las variables sociodemográficas seleccionadas para el año 2011. A diferencia del análisis de componentes principales, cuyo objetivo principal es sintetizar la varianza total de las variables originales, el análisis factorial busca explicar la varianza común compartida entre ellas.

En este ejercicio se utiliza la misma base de datos de 2011 y se parte de la matriz de correlaciones entre las variables originales. Esta matriz se obtiene mediante la función cor(), utilizando las correlaciones de Pearson, que son las que R calcula por defecto. Posteriormente, se aplica un análisis factorial exploratorio con la librería psych, empleando el método de ejes principales (fm = "pa"), adecuado cuando el interés se centra en estimar factores comunes a partir de la estructura correlacional de los datos.

Dado que el análisis previo de componentes principales sugería la existencia de dos dimensiones principales —una vinculada al envejecimiento y los hogares unipersonales de personas mayores, y otra asociada a la presencia de población extranjera—, se propone inicialmente una solución de dos factores. No obstante, teniendo en cuenta que el índice KMO obtenido previamente mostraba una adecuación moderada-baja, los resultados del análisis factorial se interpretan con finalidad exploratoria y con la debida cautela.

# Análisis factorial exploratorio con ejes principales
fa_2011 <- fa(
  r = cor_2011,
  nfactors = 2,
  fm = "pa",
  rotate = "varimax",
  n.obs = nrow(X2011)
)


fa_2011

## Factor Analysis using method =  pa
## Call: fa(r = cor_2011, nfactors = 2, n.obs = nrow(X2011), rotate = "varimax", 
##     fm = "pa")
## Standardized loadings (pattern matrix) based upon correlation matrix
##            PA1  PA2    h2   u2 com
## P75       0.86 0.04 0.738 0.26 1.0
## HOG_64    0.73 0.18 0.564 0.44 1.1
## HOG_MONO -0.25 0.05 0.065 0.93 1.1
## EXT       0.11 0.49 0.249 0.75 1.1
## EXT_MEN  -0.12 0.79 0.634 0.37 1.0
## 
##                        PA1  PA2
## SS loadings           1.36 0.90
## Proportion Var        0.27 0.18
## Cumulative Var        0.27 0.45
## Proportion Explained  0.60 0.40
## Cumulative Proportion 0.60 1.00
## 
## Mean item complexity =  1.1
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  10  with the objective function =  0.78 with Chi Square =  109.49
## df of  the model are 1  and the objective function was  0.01 
## 
## The root mean square of the residuals (RMSR) is  0.02 
## The df corrected root mean square of the residuals is  0.06 
## 
## The harmonic n.obs is  144 with the empirical chi square  0.87  with prob <  0.35 
## The total n.obs was  144  with Likelihood Chi Square =  1.86  with prob <  0.17 
## 
## Tucker Lewis Index of factoring reliability =  0.913
## RMSEA index =  0.077  and the 90 % confidence intervals are  0 0.252
## BIC =  -3.11
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    PA1  PA2
## Correlation of (regression) scores with factors   0.90 0.82
## Multiple R square of scores with factors          0.81 0.68
## Minimum correlation of possible factor scores     0.61 0.35

print(fa_2011$loadings, cutoff = 0.30)

## 
## Loadings:
##          PA1    PA2   
## P75       0.858       
## HOG_64    0.728       
## HOG_MONO              
## EXT              0.487
## EXT_MEN          0.788
## 
##                  PA1   PA2
## SS loadings    1.355 0.895
## Proportion Var 0.271 0.179
## Cumulative Var 0.271 0.450

fa_2011$Vaccounted

##                             PA1       PA2
## SS loadings           1.3552668 0.8953015
## Proportion Var        0.2710534 0.1790603
## Cumulative Var        0.2710534 0.4501137
## Proportion Explained  0.6021887 0.3978113
## Cumulative Proportion 0.6021887 1.0000000

fa.diagram(fa_2011)

El análisis factorial exploratorio mediante el método de ejes principales y rotación varimax permite identificar una solución de dos factores con una interpretación sustantiva clara. El primer factor, PA1, presenta cargas elevadas en la población de 75 años o más y en los hogares unipersonales de mayores de 64 años, por lo que puede interpretarse como una dimensión de envejecimiento demográfico y vulnerabilidad residencial de personas mayores. El segundo factor, PA2, muestra cargas más elevadas en la población extranjera infantil y, en menor medida, en la población extranjera total, por lo que puede interpretarse como una dimensión asociada a la presencia de población extranjera.

La solución factorial explica conjuntamente el 45,01 % de la varianza total, correspondiendo el 27,11 % al primer factor y el 17,91 % al segundo. Aunque este porcentaje es inferior al obtenido en el ACP, debe tenerse en cuenta que el análisis factorial trata de explicar la varianza común entre las variables, no la varianza total. Los indicadores de ajuste son razonables: el RMSR es bajo, el TLI se sitúa por encima de 0.90, el BIC es negativo y el contraste de suficiencia de dos factores no resulta significativo, por lo que no se rechaza la adecuación de la solución bifactorial.

No obstante, la interpretación debe realizarse con cautela. La variable HOG_MONO presenta una comunalidad muy baja y no carga de forma relevante en ninguno de los dos factores, lo que indica que los hogares monoparentales no quedan bien representados por la estructura factorial obtenida. En conjunto, el análisis factorial permite resumir la información en dos dimensiones principales —envejecimiento y extranjería—, aunque la solución debe entenderse como exploratoria.

Se indica qué procedimiento se ha utilizado para obtener la matriz de correlación.

La matriz de correlación se obtuvo a partir de las variables sociodemográficas seleccionadas para el año 2011 mediante la función cor() de R. Concretamente, se utilizó el objeto X2011, que contiene únicamente las variables numéricas del análisis, excluyendo el identificador de sección censal. El procedimiento aplicado fue:

cor_2011 <- cor(X2011, use = “pairwise.complete.obs”)

Al no especificar el argumento method, la función cor() calcula por defecto la matriz de correlaciones de Pearson. Por tanto, la matriz utilizada en el análisis factorial recoge las asociaciones lineales entre las variables originales de 2011.

El argumento use = "pairwise.complete.obs" indica que, en caso de existir valores perdidos, cada correlación se calcularía utilizando todos los pares de observaciones disponibles para las dos variables implicadas. No obstante, en este dataset no se detectaron valores perdidos, por lo que la matriz se calculó sobre el conjunto completo de secciones censales.

Esta matriz de correlaciones fue posteriormente utilizada como entrada del análisis factorial exploratorio con la función fa() de la librería psych, aplicando el método de ejes principales.

Obtención de las puntuaciones factoriales

Para obtener puntuaciones factoriales necesitamos aplicar fa() sobre la matriz de datos, no solo sobre la matriz de correlaciones. Como ya tenemos las variables del dataset original estandarizadas en X2011_z, usaremos este objeto como base para esta tarea.

# Comprobar que X2011_z no contiene el código de sección censal
colnames(X2011_z)

## [1] "P75"      "HOG_64"   "HOG_MONO" "EXT"      "EXT_MEN"

# Análisis factorial sobre los datos estandarizados para obtener puntuaciones
fa_2011_scores <- fa(
  r = X2011_z,
  nfactors = 2,
  fm = "pa",
  rotate = "varimax",
  scores = "regression"
)

## maximum iteration exceeded

# Comprobar nombres de las puntuaciones factoriales
colnames(fa_2011_scores$scores)

## [1] "PA1" "PA2"

# Crear directamente el nuevo dataset de puntuaciones
puntuaciones_factores <- data.frame(
  COD_SECC = datos_2011$COD_SECC,
  F1_Envejecimiento = fa_2011_scores$scores[, 1],
  F2_Extranjeria = fa_2011_scores$scores[, 2]
)

head(puntuaciones_factores)

##     COD_SECC F1_Envejecimiento F2_Extranjeria
## 1 3907501001        -0.2127151     0.04248055
## 2 3907501002         0.7120695    -0.32247211
## 3 3907501003         0.4653397    -0.33637468
## 4 3907501004         2.0834454    -0.22251810
## 5 3907501005         0.4270473     1.67110942
## 6 3907501006        -0.4197269     0.99483423

Se realiza un análisis estadístico de las puntuaciones y se indica el método utilizado para obtener los nuevos factores.

El análisis estadístico descriptivo de estas puntuaciones permite conocer la posición relativa de las secciones censales en cada dimensión factorial. Las medias de las puntuaciones suelen situarse próximas a cero, ya que proceden de variables estandarizadas, mientras que la dispersión indica el grado de diferenciación existente entre secciones. Valores altos en el primer factor indican secciones con mayor peso relativo de población de 75 años o más y hogares unipersonales de mayores de 64 años. Valores altos en el segundo factor indican secciones con mayor presencia relativa de población extranjera, especialmente población extranjera infantil.

# Seleccionar únicamente las puntuaciones factoriales
factores_numericos <- puntuaciones_factores[
  , c("F1_Envejecimiento", "F2_Extranjeria")
]

# Tabla descriptiva
tabla_desc_factores <- data.frame(
  Factor = names(factores_numericos),
  Media = round(sapply(factores_numericos, mean), 4),
  Desv_tipica = round(sapply(factores_numericos, sd), 4),
  Minimo = round(sapply(factores_numericos, min), 4),
  Q1 = round(sapply(factores_numericos, quantile, probs = 0.25), 4),
  Mediana = round(sapply(factores_numericos, median), 4),
  Q3 = round(sapply(factores_numericos, quantile, probs = 0.75), 4),
  Maximo = round(sapply(factores_numericos, max), 4)
)

rownames(tabla_desc_factores) <- NULL

kable(
  tabla_desc_factores,
  caption = "Análisis descriptivo de las puntuaciones factoriales por sección censal",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

# Histogramas
par(mfrow = c(2, 3))
for (v in names(factores_numericos)) {
  hist(factores_numericos[[v]],
       main = paste("Histograma de", v),
       xlab = v,
       col = "blue",
       border = "white")
}
par(mfrow = c(1, 1))

# Boxplots para detectar dispersión y posibles valores extremos
boxplot(factores_numericos,
        main = "Distribución de variables sociodemográficas 2011",
        las = 2,
        col = "blue")

Para obtener las puntuaciones factoriales de cada sección censal se aplicó el análisis factorial exploratorio sobre la matriz de datos estandarizados correspondiente al año 2011. La extracción de los factores se realizó mediante el método de ejes principales (fm = "pa"), manteniendo una solución de dos factores con rotación varimax. La rotación se utilizó para facilitar la interpretación de la estructura factorial, permitiendo diferenciar con mayor claridad una dimensión asociada al envejecimiento y otra vinculada a la presencia de población extranjera.

Las puntuaciones factoriales se obtuvieron mediante el método de regresión (scores = "regression"). Este procedimiento estima, para cada sección censal, una puntuación en cada uno de los factores retenidos. De este modo, cada sección queda caracterizada por dos nuevas variables: una puntuación en el factor de envejecimiento y vulnerabilidad residencial de personas mayores, y otra puntuación en el factor de presencia de población extranjera.

En consecuencia, las puntuaciones factoriales permiten sintetizar la información original en dos dimensiones interpretables y comparar las secciones censales de Santander según su perfil sociodemográfico.

Identificación de las 10 secciones más vulnerables demográficamente de Santander.

Para identificar las secciones censales más vulnerables demográficamente se utiliza la puntuación factorial correspondiente al factor de envejecimiento. Este factor presenta cargas elevadas en la población de 75 años o más y en los hogares unipersonales de mayores de 64 años, por lo que valores altos en este factor indican secciones con mayor vulnerabilidad demográfica asociada al envejecimiento y a la posible soledad residencial de personas mayores. A partir de este criterio, se ordenaron las secciones censales de mayor a menor puntuación en F1_Envejecimiento y se seleccionaron las diez primeras.

Estas secciones no deben interpretarse como las más vulnerables en términos generales, sino específicamente como las más vulnerables desde el punto de vista demográfico asociado al envejecimiento y a los hogares unipersonales de personas mayores.

# Ordenar secciones según el factor de envejecimiento
top10_vulnerables <- puntuaciones_factores[
  order(puntuaciones_factores$F1_Envejecimiento, decreasing = TRUE),
][1:10, ]

# Incorporar columna con el puesto en la clasificación
top10_vulnerables <- data.frame(
  Puesto = 1:nrow(top10_vulnerables),
  top10_vulnerables
)

# Eliminar nombres de fila
rownames(top10_vulnerables) <- NULL

# Mostrar tabla
kable(
  top10_vulnerables,
  caption = "Diez secciones censales con mayor vulnerabilidad demográfica considerada como mayor puntuación en F1_Envejecimiento",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

Los resultados muestran que la sección censal 3907506010 ocupa el primer lugar, con una puntuación de 2.59 en F1_Envejecimiento, seguida de las secciones 3907506011´, ´3907502009 y 3907501004, todas ellas con puntuaciones superiores a 2. Estos valores indican que dichas secciones se sitúan claramente por encima de la media del conjunto de Santander en la dimensión de envejecimiento, por lo que pueden considerarse las áreas con mayor vulnerabilidad demográfica desde este criterio.

La segunda puntuación factorial, F2_Extranjeria, permite matizar el perfil de estas secciones. Algunas de ellas, como 3907506010, 3907501004, 3907505001, 3907502015, 3907502012, 3907507003 y 3907502013, presentan valores negativos en el factor de extranjería, lo que sugiere que su vulnerabilidad demográfica se asocia principalmente al envejecimiento y no a una presencia elevada de población extranjera. En contrapunto a las anteriores, las secciones 3907506011, 3907502009 y 3907502026 combinan puntuaciones altas en envejecimiento con valores positivos en extranjería, por lo que presentan un perfil sociodemográfico más mixto.

Se vuelve a recordar que, en conjunto, estas diez secciones no deben interpretarse como las más vulnerables en términos generales, sino específicamente como las más vulnerables desde la perspectiva del envejecimiento demográfico y de la posible soledad residencial de personas mayores. El ranking obtenido permitiría priorizar territorialmente aquellas secciones en las que el peso relativo de la población mayor y de los hogares unipersonales de mayores es más elevado.

Justificación del método de agrupamiento

Se utiliza un análisis de clúster jerárquico aglomerativo a partir de las puntuaciones factoriales obtenidas en el análisis factorial exploratorio. La distancia utilizada es la distancia Manhattan, tal como solicita el enunciado. Para seleccionar el método de agrupamiento más adecuado se comparan distintos métodos de enlace —simple, completo, promedio y McQuitty— mediante la correlación cofenética. Este criterio permite evaluar en qué medida el dendrograma conserva las distancias originales entre observaciones. Se seleccionará el método con mayor correlación cofenética, al ser el que reproduce con mayor fidelidad la estructura de distancias del conjunto de secciones censales.

No se utiliza el método de Ward porque está más directamente vinculado a la distancia euclídea y a la minimización de la varianza interna de los grupos, mientras que, en este ejercicio, el enunciado impone el uso de distancia Manhattan.

# 1. DATOS DE PARTIDA: PUNTUACIONES FACTORIALES

# Seleccionamos las puntuaciones factoriales obtenidas en el ejercicio 2
datos_cluster <- puntuaciones_factores[, c("F1_Envejecimiento", "F2_Extranjeria")]

# Estandarizamos las puntuaciones para que ambos factores tengan el mismo peso en la distancia
datos_cluster_z <- scale(datos_cluster)

# Añadimos los códigos de sección censal como nombres de fila
rownames(datos_cluster_z) <- puntuaciones_factores$COD_SECC


# 2. MATRIZ DE DISTANCIAS: MANHATTAN

dist_manhattan <- dist(datos_cluster_z, method = "manhattan")


# 3. COMPARACIÓN DE MÉTODOS DE AGRUPAMIENTO

metodos <- c("single", "complete", "average", "mcquitty")

modelos_hc <- lapply(
  metodos,
  function(m) hclust(dist_manhattan, method = m)
)

names(modelos_hc) <- metodos

# Correlación cofenética para comparar métodos
cor_cophenetica <- sapply(
  modelos_hc,
  function(modelo) cor(dist_manhattan, cophenetic(modelo))
)

tabla_metodos <- data.frame(
  Metodo = names(cor_cophenetica),
  Correlacion_cofenetica = round(cor_cophenetica, 4)
)

kable(
  tabla_metodos,
  caption = "Comparación de métodos de agrupamiento mediante correlación cofenética",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

# Seleccionamos automáticamente el método con mayor correlación cofenética
metodo_elegido <- names(which.max(cor_cophenetica))
metodo_elegido

## [1] "average"

cluster_hc <- modelos_hc[[metodo_elegido]]

Cómo puede observarse, el resultado del análisis de las correlaciones cofenéticas de los cuatro métodos de enlace testeados da como mejor opción el método del promedio o average, puesto que da el indice más elevado de las cuatro opciones, lo que significa que dicho método garantiza que el dendograma reproduzca con la mayor fidelidad la estructura de distancias del conjunto de secciones censales.

Selección del número de clústeres

Se utiliza el método de la silueta media (o average silhouette method) para seleccionar el número de clústeres que represente mejor la calidad de las diferentes posibilidades de agrupamiento y ayude a identificar el número óptimo de clústers para un conjunto de datos.

# 4. Selección del número de clústeres mediante silueta media

ks <- 2:8

silueta_media <- sapply(
  ks,
  function(k) {
    grupos <- cutree(cluster_hc, k = k)
    mean(silhouette(grupos, dist_manhattan)[, 3])
  }
)

tabla_silueta <- data.frame(
  Numero_clusters = ks,
  Silueta_media = round(silueta_media, 4)
)

kable(
  tabla_silueta,
  caption = "Silueta media según el número de clústeres",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

# Número óptimo de clústeres según la mayor silueta media
k_optimo <- ks[which.max(silueta_media)]
k_optimo

## [1] 2

# Gráfico sencillo de silueta media
ggplot(tabla_silueta, aes(x = Numero_clusters, y = Silueta_media)) +
  geom_line(color = "#333333", linewidth = 0.8) +
  geom_point(color = "#1F4E79", size = 3) +
  scale_x_continuous(breaks = ks) +
  labs(
    title = "Selección del número de clústeres",
    subtitle = "Criterio: mayor silueta media",
    x = "Número de clústeres",
    y = "Silueta media"
  ) +
  theme_minimal(base_size = 12)

La selección inicial de una solución de dos clústeres se justifica a partir del criterio de silueta media, ya que el valor máximo se alcanza para \(k = 2\). Esto indica que, desde el punto de vista de la separación interna de los grupos, la partición en dos conglomerados es la que ofrece mejor equilibrio entre cohesión dentro de los clústeres y diferenciación entre ellos.

Dendrograma con el número de clústeres seleccionado

plot(
  cluster_hc,
  labels = FALSE,
  hang = -1,
  main = paste("Dendrograma de secciones censales - método", metodo_elegido),
  xlab = "Secciones censales",
  sub = ""
)

rect.hclust(
  cluster_hc,
  k = k_optimo,
  border = 2:(k_optimo + 1)
)

El dendrograma obtenido mediante distancia Manhattan y método de enlace average refuerza la primera lectura recibida de la interpretación de la silueta media, ya que muestra una separación temprana entre un grupo muy reducido (sección única) y el conjunto principal de secciones censales. El corte en dos grupos permite, por tanto, identificar una estructura básica de diferenciación: por un lado, la mayoría de secciones con perfiles relativamente próximos entre sí y, por otro, una sección o grupo muy singular que se separa claramente del resto.

No obstante, esta solución debe interpretarse con cautela, ya que la partición en dos clústeres parece estar muy influida por la presencia de un caso extremo o perfil atípico. Por ello, aunque \(k = 2\) constituye la primera opción desde el punto de vista del índice de silueta, más adelante convendría valorar soluciones alternativas con mayor número de clústeres si el objetivo es obtener una segmentación territorial más rica e interpretable.

Inclusión en la base de datos de resultados del cluster al que pertenece cada sección censal.

# 5. ASIGNACIÓN DE CLUSTER A CADA SECCIÓN CENSAL

grupo_cluster <- cutree(cluster_hc, k = k_optimo)

resultados_cluster <- data.frame(
  COD_SECC = puntuaciones_factores$COD_SECC,
  F1_Envejecimiento = puntuaciones_factores$F1_Envejecimiento,
  F2_Extranjeria = puntuaciones_factores$F2_Extranjeria,
  Cluster = factor(grupo_cluster)
)

# Ordenar por clúster y sección censal
resultados_cluster <- resultados_cluster[
  order(resultados_cluster$Cluster, resultados_cluster$COD_SECC),
]

rownames(resultados_cluster) <- NULL


# Generar tabla con los 20 primeros distritos mostrando el clúster al que se le ha asignado
kable(
  head(resultados_cluster, 20),
  caption = "Primeras secciones censales con clúster asignado",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

set.seed(123)  # Opcional: permite reproducir siempre la misma muestra aleatoria

kable(
  resultados_cluster[sample(nrow(resultados_cluster), 20), ],
  caption = "Veinte secciones censales seleccionadas aleatoriamente con clúster asignado",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

Perfil medio de cada clúster

Para ayudar a una mejor interpretación los grupos generados al clusterizar, se analiza el perfil medio de las secciones que componen cada uno de los dos clústers en los que se ha particionado inicialmente el total de secciones censales.

# Número de secciones por clúster
tabla_tamano_cluster <- as.data.frame(table(resultados_cluster$Cluster))
colnames(tabla_tamano_cluster) <- c("Cluster", "Numero_secciones")

# Perfil medio de los clústeres
perfil_cluster <- aggregate(
  resultados_cluster[, c("F1_Envejecimiento", "F2_Extranjeria")],
  by = list(Cluster = resultados_cluster$Cluster),
  FUN = mean
)

perfil_cluster <- merge(
  perfil_cluster,
  tabla_tamano_cluster,
  by = "Cluster"
)

kable(
  perfil_cluster,
  caption = "Perfil medio de los clústeres según las puntuaciones factoriales",
  align = "c",
  row.names = FALSE,
  digits = 3
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

La solución de dos clústeres muestra una partición muy desequilibrada. El clúster 1 agrupa 143 secciones censales y presenta valores medios prácticamente nulos en ambos factores (F1_Envejecimiento = 0.012; F2_Extranjeria = -0.020), por lo que representa el perfil general o mayoritario del municipio. En cambio, el clúster 2 está formado por una única sección censal, con bajo envejecimiento relativo (F1 = -1.699) y una puntuación muy elevada en extranjería (F2 = 2.912).

Por tanto, aunque la solución de dos grupos es la mejor según el criterio de silueta media, su interpretación debe realizarse con cautela, ya que no segmenta el territorio en dos perfiles amplios, sino que identifica una sección claramente atípica frente al resto. Esta solución es útil para detectar un caso singular de alta extranjería y bajo envejecimiento, pero puede resultar limitada si el objetivo del análisis es obtener una tipología territorial más rica y equilibrada.

En resumen, parece que situar el agrupamiento en \(k = 2\) maximiza la silueta porque detecta muy bien un outlier territorial, no porque existan dos grandes grupos homogéneos y comparables.

Gráfico de clústeres sobre los dos factores

ggplot(
  resultados_cluster,
  aes(
    x = F1_Envejecimiento,
    y = F2_Extranjeria,
    color = Cluster
  )
) +
  geom_point(size = 2.5, alpha = 0.85) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray60") +
  labs(
    title = "Clústeres de secciones censales de Santander",
    subtitle = "Agrupación sobre puntuaciones factoriales. Distancia Manhattan",
    x = "Factor 1: Envejecimiento",
    y = "Factor 2: Extranjería",
    color = "Clúster"
  ) +
  theme_minimal(base_size = 12)

Tras los análisis realizados en base a una solución inicial de dos clústeres, se observa que dicha clusterización genera una partición muy desequilibrada, con un primer grupo formado por 143 secciones censales y un segundo grupo compuesto por una única sección. Aunque esta solución puede ser estadísticamente posible, su interpretación debe realizarse con cautela, ya que no refleja propiamente una tipología territorial de dos grupos, sino la separación de una sección censal con un perfil claramente atípico.

El segundo clúster presenta una puntuación negativa en el factor de envejecimiento y una puntuación muy elevada en el factor de extranjería, lo que sugiere que se trata de una sección singular, caracterizada por baja vulnerabilidad asociada al envejecimiento y alta presencia relativa de población extranjera. Por tanto, más que un grupo territorial amplio, este clúster puede interpretarse como un caso extremo.

En consecuencia, para seleccionar el número final de clústeres no debería atenderse únicamente a la silueta media, sino que también debe dársele importancia a la interpretabilidad de la solución y al tamaño de los grupos resultantes. Si una solución alternativa con tres o cuatro clústeres ofrece grupos más equilibrados y perfiles sociodemográficos diferenciados, resultaría preferible para construir una clasificación útil de las secciones censales de Santander.

Clusterización alternativa

Para realizar una exploración de la mejor opción respecto al número de grupos, se analizan soluciones alternativas entre 3 y 8 clústeres, considerando no solo la silueta media, sino también el tamaño de los grupos, la estructura observada en el dendrograma y la interpretabilidad de los perfiles resultantes.

# Comparación de soluciones entre 2 y 8 clústeres

ks <- 2:8

diagnostico_clusters <- data.frame(
  Numero_clusters = integer(),
  Silueta_media = numeric(),
  Minimo_secciones_cluster = integer(),
  Maximo_secciones_cluster = integer()
)

for (k in ks) {
  
  grupos_k <- cutree(cluster_hc, k = k)
  sil_k <- silhouette(grupos_k, dist_manhattan)
  tam_k <- table(grupos_k)
  
  diagnostico_clusters <- rbind(
    diagnostico_clusters,
    data.frame(
      Numero_clusters = k,
      Silueta_media = round(mean(sil_k[, 3]), 4),
      Minimo_secciones_cluster = min(tam_k),
      Maximo_secciones_cluster = max(tam_k)
    )
  )
}

kable(
  diagnostico_clusters,
  caption = "Diagnóstico de soluciones alternativas de clústeres",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

Solución alternativa propuesta: 4 clústeres

La solución de 4 clústeres se considera una alternativa más útil, al permitir diferenciar perfiles sociodemográficos más matizados sin reducir el análisis a la separación de un único caso extremo.

Esta solución mantiene la distancia Manhattan y el método jerárquico aglomerativo con enlace promedio (average). El método average resulta adecuado porque ofrece una agrupación menos sensible a casos extremos que el enlace simple y más equilibrada que otros métodos, conservando una estructura interpretable de distancias entre secciones.

En consecuencia, se adopta una solución alternativa de cuatro clústeres, no porque maximice estrictamente la silueta media, sino porque ofrece una segmentación territorial más informativa y operativamente útil de las secciones censales de Santander.

# SELECCIÓN ALTERNATIVA DE 4 CLÚSTERES

k_alt <- 4

grupo_cluster_alt <- cutree(cluster_hc, k = k_alt)

# Base de resultados con clúster asignado
resultados_cluster_alt <- data.frame(
  COD_SECC = puntuaciones_factores$COD_SECC,
  F1_Envejecimiento = puntuaciones_factores$F1_Envejecimiento,
  F2_Extranjeria = puntuaciones_factores$F2_Extranjeria,
  Cluster = factor(grupo_cluster_alt)
)

# Ordenar por clúster y sección
resultados_cluster_alt <- resultados_cluster_alt[
  order(resultados_cluster_alt$Cluster, resultados_cluster_alt$COD_SECC),
]

rownames(resultados_cluster_alt) <- NULL

kable(
  resultados_cluster_alt[sample(nrow(resultados_cluster_alt), 20), ],
  caption = "Primeras secciones censales con clúster asignado. Solución alternativa de 4 clústeres",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

# PERFIL MEDIO DE CADA UNO DE LOS 4 CLÚSTERES

# Tamaño de cada clúster
tamano_cluster_alt <- as.data.frame(table(resultados_cluster_alt$Cluster))
colnames(tamano_cluster_alt) <- c("Cluster", "Numero_secciones")

# Perfil medio de cada clúster
perfil_cluster_alt <- aggregate(
  resultados_cluster_alt[, c("F1_Envejecimiento", "F2_Extranjeria")],
  by = list(Cluster = resultados_cluster_alt$Cluster),
  FUN = mean
)

perfil_cluster_alt <- merge(
  perfil_cluster_alt,
  tamano_cluster_alt,
  by = "Cluster"
)

kable(
  perfil_cluster_alt,
  caption = "Perfil medio de los clústeres según las puntuaciones factoriales. Solución de 4 clústeres",
  align = "c",
  row.names = FALSE,
  digits = 3
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

Dendrograma con 4 clústeres

plot(
  cluster_hc,
  labels = FALSE,
  hang = -1,
  main = "Dendrograma de secciones censales - método average",
  xlab = "Secciones censales",
  sub = ""
)

rect.hclust(
  cluster_hc,
  k = k_alt,
  border = 2:(k_alt + 1)
)

Gráfico de los clústeres sobre los factores

ggplot(
  resultados_cluster_alt,
  aes(
    x = F1_Envejecimiento,
    y = F2_Extranjeria,
    color = Cluster
  )
) +
  geom_point(size = 2.5, alpha = 0.85) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray60") +
  labs(
    title = "Clústeres de secciones censales de Santander",
    subtitle = "Solución alternativa de 4 clústeres. Distancia Manhattan y método average",
    x = "Factor 1: Envejecimiento",
    y = "Factor 2: Extranjería",
    color = "Clúster"
  ) +
  theme_minimal(base_size = 12)

La solución alternativa de cuatro clústeres permite obtener una segmentación más informativa de las secciones censales de Santander que la solución inicial de dos grupos. Aunque la solución de dos clústeres presentaba una silueta media elevada, generaba una partición muy desequilibrada, separando una única sección censal del resto. Por ello, se optó por explorar una solución de cuatro clústeres, más útil desde el punto de vista interpretativo.

El clúster 2 agrupa la mayoría de las secciones censales y representa un perfil sociodemográfico medio, con puntuaciones bajas en extranjería y sin un envejecimiento especialmente elevado. El clúster 1 reúne secciones con mayor presencia relativa de población extranjera, pero sin destacar por envejecimiento. El clúster 3 identifica las secciones con mayor vulnerabilidad demográfica asociada al envejecimiento, al presentar una media claramente elevada en el factor de envejecimiento. Finalmente, el clúster 4 está formado por una única sección, caracterizada por una puntuación muy alta en extranjería y baja en envejecimiento, por lo que debe interpretarse como un caso singular o atípico.

En conjunto, esta clasificación permite diferenciar cuatro perfiles territoriales: secciones con mayor extranjería, secciones de perfil medio-bajo en ambas dimensiones, secciones envejecidas y una sección atípica con alta extranjería. Esta solución resulta más adecuada para segmentar las secciones censales de Santander, ya que combina criterios estadísticos con interpretabilidad sustantiva.

Validación interna de las solución de 4 grupos

Además de la silueta media, se utilizan otros índices de validación interna para evaluar la calidad de la solución de clústeres. Para ello se calculan índice de Dunn, que mide la calidad de una partición en clusters comparando dos magnitudes: la separación entre clusters (distancia mínima entre puntos de clusters distintos) y la compacidad interna (diámetro máximo, es decir, la mayor distancia entre dos puntos dentro del mismo cluster). Valores altos indican clusters bien separados y compactos, lo que se considera una buena solución. También se mide la versión robusta de este idicador (Dunn2), que sustituye los operadores sensibles a outliers por promedios: en lugar de tomar la distancia mínima entre clusters usa la distancia media entre puntos de clusters distintos, y en lugar del diámetro máximo usa la distancia intra-cluster promedio.Al igual que Dunn, valores más altos indican una mejor solución: clusters compactos y bien separados.

También se calcula el índice Pearson gamma, el cual se usa para medir la calidad de una clusterización jerárquica comparando cuántas relaciones entre pares de elementos quedan mejor explicadas por la partición en clústeres que por una estructura aleatoria. En esencia, evalúa si los grupos formados son más coherentes internamente y más separados entre sí.

# Estadísticos de validación interna
validacion_cluster <- cluster.stats(
  d = dist_manhattan,
  clustering = grupo_cluster_alt
)

# Tabla comparativa de estadísticos

ks <- 2:8

tabla_validacion <- data.frame(
  Numero_clusters = integer(),
  Silueta_media = numeric(),
  Dunn = numeric(),
  Dunn2 = numeric(),
  Pearson_gamma = numeric(),
  Minimo_secciones_cluster = integer(),
  Maximo_secciones_cluster = integer()
)

for (k in ks) {
  
  grupos_k <- cutree(cluster_hc, k = k)
  
  stats_k <- cluster.stats(
    d = dist_manhattan,
    clustering = grupos_k
  )
  
  tam_k <- table(grupos_k)
  
  tabla_validacion <- rbind(
    tabla_validacion,
    data.frame(
      Numero_clusters = k,
      Silueta_media = round(stats_k$avg.silwidth, 4),
      Dunn = round(stats_k$dunn, 4),
      Dunn2 = round(stats_k$dunn2, 4),
      Pearson_gamma = round(stats_k$pearsongamma, 4),
      Minimo_secciones_cluster = min(tam_k),
      Maximo_secciones_cluster = max(tam_k)
    )
  )
}

kable(
  tabla_validacion,
  caption = "Comparación de soluciones de clústeres mediante índices de validación interna",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

La comparación de soluciones de clústeres muestra que la partición en dos grupos obtiene la mayor silueta media y los valores más elevados en los índices de Dunn. Sin embargo, esta solución genera una distribución extremadamente desequilibrada, con un clúster formado por 143 secciones censales y otro compuesto por una única sección. Por tanto, aunque estadísticamente detecta una separación clara, su utilidad para segmentar territorialmente Santander es limitada, ya que refleja principalmente la presencia de un caso atípico.

Por este motivo, la elección final del número de clústeres no se basa únicamente en la maximización de la silueta media, sino en una valoración conjunta de distintos criterios: calidad interna del agrupamiento, tamaño de los grupos, estructura del dendrograma e interpretabilidad sustantiva. La solución de cuatro clústeres presenta una silueta media aceptable y un valor elevado de Pearson_gamma, muy próximo al máximo observado. Además, permite diferenciar perfiles sociodemográficos más útiles: secciones con mayor presencia extranjera, secciones de perfil medio, secciones envejecidas y una sección singular o atípica.

En consecuencia, se selecciona la solución de cuatro clústeres como alternativa más adecuada para la segmentación de las secciones censales de Santander. Esta elección sacrifica parcialmente la maximización estricta de algunos índices, pero mejora la utilidad interpretativa y territorial del análisis.

Preparación de datos:

# Seleccionar solo el año 1986 y las variables solicitadas
Produc_1986 <- subset(
  Produc,
  year == 1986,
  select = c(
    state,
    year,
    pcap,
    hwy,
    water,
    util,
    pc,
    gsp,
    emp,
    unemp
  )
)

# Revisar el resultado
head(Produc_1986)

##           state year      pcap      hwy    water     util        pc    gsp
## 17      ALABAMA 1986  19723.37  8813.24  2308.99  8601.14  61628.88  48409
## 34      ARIZONA 1986  20212.42  6553.40  2923.15 10735.88  44189.19  46058
## 51     ARKANSAS 1986   9787.20  4652.71   960.66  4173.84  32889.30  28168
## 68   CALIFORNIA 1986 139042.68 43350.54 24592.33 71099.80 363779.70 464550
## 85     COLORADO 1986  18889.54  6148.80  4340.24  8400.50  47485.05  51781
## 102 CONNECTICUT 1986  17290.61  6878.14  3314.55  7097.92  40413.22  61750
##         emp unemp
## 17   1463.3   9.8
## 34   1337.8   6.9
## 51    813.8   8.7
## 68  11258.0   6.7
## 85   1408.3   7.4
## 102  1604.2   3.8

str(Produc_1986)

## 'data.frame':    48 obs. of  10 variables:
##  $ state: Factor w/ 48 levels "ALABAMA","ARIZONA",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ year : int  1986 1986 1986 1986 1986 1986 1986 1986 1986 1986 ...
##  $ pcap : num  19723 20212 9787 139043 18890 ...
##  $ hwy  : num  8813 6553 4653 43351 6149 ...
##  $ water: num  2309 2923 961 24592 4340 ...
##  $ util : num  8601 10736 4174 71100 8400 ...
##  $ pc   : num  61629 44189 32889 363780 47485 ...
##  $ gsp  : int  48409 46058 28168 464550 51781 61750 10072 152269 88827 11672 ...
##  $ emp  : num  1463 1338 814 11258 1408 ...
##  $ unemp: num  9.8 6.9 8.7 6.7 7.4 3.8 4.3 5.7 5.9 8.7 ...

summary(Produc_1986)

##          state         year           pcap             hwy       
##  ALABAMA    : 1   Min.   :1986   Min.   :  2936   Min.   : 1830  
##  ARIZONA    : 1   1st Qu.:1986   1st Qu.:  9461   1st Qu.: 4137  
##  ARKANSAS   : 1   Median :1986   Median : 19306   Median : 8179  
##  CALIFORNIA : 1   Mean   :1986   Mean   : 27393   Mean   :10736  
##  COLORADO   : 1   3rd Qu.:1986   3rd Qu.: 29882   3rd Qu.:12292  
##  CONNECTICUT: 1   Max.   :1986   Max.   :139043   Max.   :43351  
##  (Other)    :42                                                  
##      water              util               pc              gsp        
##  Min.   :  335.5   Min.   :  770.8   Min.   :  6939   Min.   :  7585  
##  1st Qu.: 1023.5   1st Qu.: 3042.7   1st Qu.: 28129   1st Qu.: 21183  
##  Median : 2926.7   Median : 8315.6   Median : 50303   Median : 47956  
##  Mean   : 4366.6   Mean   :12290.5   Mean   : 72260   Mean   : 75459  
##  3rd Qu.: 4867.9   3rd Qu.:13173.6   3rd Qu.: 83594   3rd Qu.: 87571  
##  Max.   :24592.3   Max.   :75237.1   Max.   :363780   Max.   :464550  
##                                                                       
##       emp              unemp       
##  Min.   :  196.3   Min.   : 2.800  
##  1st Qu.:  579.6   1st Qu.: 5.300  
##  Median : 1373.2   Median : 6.750  
##  Mean   : 2049.9   Mean   : 6.929  
##  3rd Qu.: 2586.4   3rd Qu.: 8.275  
##  Max.   :11258.0   Max.   :13.000  
## 

dim(Produc_1986)

## [1] 48 10

El dataset filtrado contiene 48 observaciones, correspondientes a distintos estados de Estados Unidos en el año 1986, y 10 variables. La variable state actúa como identificador territorial y la variable year toma un único valor, 1986, por lo que no aporta variabilidad al análisis. Las variables restantes son de tipo económico y recogen información sobre capital público (pcap), infraestructuras viarias (hwy), infraestructuras de agua y saneamiento (water), otras infraestructuras públicas (util), capital privado (pc), producto bruto estatal (gsp), empleo no agrícola (emp) y tasa de desempleo (unemp).

La exploración descriptiva muestra que las variables económicas presentan escalas muy diferentes y una elevada dispersión. Algunas variables, como pcap, pc, gsp y emp, tienen valores máximos muy superiores a sus medianas, lo que sugiere una distribución asimétrica y la presencia de estados con un tamaño económico claramente mayor. En cambio, unemp se expresa como tasa y tiene una escala mucho más reducida. Por este motivo, antes de aplicar el escalado multidimensional resulta necesario transformar y estandarizar las variables numéricas, evitando que las distancias entre estados queden dominadas por las variables de mayor magnitud.

Antes de aplicar el escalado multidimensional se realiza una transformación de los datos para evitar que las distancias entre estados quedaran dominadas por las variables de mayor escala. En primer lugar, se excluirán del análisis las variables state y year, conservando state únicamente como identificador territorial. Dado que las variables económicas de volumen presentan fuertes diferencias de magnitud y asimetría positiva, se aplicará una transformación logarítmica a las variables de capital, infraestructuras, producto bruto estatal y empleo. La tasa de desempleo se mantiene en su escala original al tratarse de una variable porcentual. Finalmente, todas las variables son estandarizadas para que contribuyeran de forma comparable al cálculo de la matriz de distancias. Esta matriz de distancias constituye la entrada del escalado multidimensional.

# Variables económicas seleccionadas
Produc_1986_num <- Produc_1986[, c(
  "pcap",
  "hwy",
  "water",
  "util",
  "pc",
  "gsp",
  "emp",
  "unemp"
)]

# Asignar nombres de los estados como identificadores de fila
rownames(Produc_1986_num) <- Produc_1986$state

# Transformación logarítmica para variables de tamaño o volumen
Produc_1986_transf <- Produc_1986_num

variables_log <- c(
  "pcap",
  "hwy",
  "water",
  "util",
  "pc",
  "gsp",
  "emp"
)

Produc_1986_transf[, variables_log] <- log(
  Produc_1986_transf[, variables_log]
)

# La variable unemp se mantiene en escala original porque ya es una tasa
# Después se estandarizan todas las variables
Produc_1986_scaled <- scale(Produc_1986_transf)

# Comprobación
summary(Produc_1986_scaled)

##       pcap              hwy               water              util        
##  Min.   :-1.9268   Min.   :-1.79322   Min.   :-1.8459   Min.   :-2.0375  
##  1st Qu.:-0.6818   1st Qu.:-0.78623   1st Qu.:-0.8257   1st Qu.:-0.7796  
##  Median : 0.0801   Median : 0.05191   Median : 0.1363   Median : 0.1472  
##  Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.5457   3rd Qu.: 0.55844   3rd Qu.: 0.6018   3rd Qu.: 0.5682  
##  Max.   : 2.1846   Max.   : 2.11510   Max.   : 2.0843   Max.   : 2.1704  
##        pc                gsp                emp               unemp        
##  Min.   :-2.11378   Min.   :-1.76622   Min.   :-1.84821   Min.   :-1.8575  
##  1st Qu.:-0.58863   1st Qu.:-0.75777   1st Qu.:-0.78053   1st Qu.:-0.7329  
##  Median : 0.04485   Median : 0.04445   Median : 0.07235   Median :-0.0806  
##  Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.59865   3rd Qu.: 0.63577   3rd Qu.: 0.69772   3rd Qu.: 0.6054  
##  Max.   : 2.20139   Max.   : 2.27419   Max.   : 2.15033   Max.   : 2.7309

round(apply(Produc_1986_scaled, 2, mean), 3)

##  pcap   hwy water  util    pc   gsp   emp unemp 
##     0     0     0     0     0     0     0     0

round(apply(Produc_1986_scaled, 2, sd), 3)

##  pcap   hwy water  util    pc   gsp   emp unemp 
##     1     1     1     1     1     1     1     1

Una vez aplicadas las transformaciones, se presenta el resumen descriptivo de las variables que se utilizarán en el escalado multidimensional. Las variables económicas de volumen fueron transformadas mediante logaritmos para reducir la asimetría y suavizar el efecto de los estados con mayor tamaño económico. Posteriormente, todas las variables fueron estandarizadas, por lo que sus medias se sitúan próximas a 0 y sus desviaciones típicas son iguales a 1. Esta transformación permite que todas las variables contribuyan de forma comparable al cálculo de las distancias entre estados.

Presentación de una tabla descriptiva de las variables empleadas una vez transformadas.

# Convertir la matriz estandarizada en data frame
Produc_1986_scaled_df <- as.data.frame(Produc_1986_scaled)

# Tabla descriptiva de variables transformadas y estandarizadas
tabla_desc_transformadas <- data.frame(
  Variable = names(Produc_1986_scaled_df),
  Media = round(sapply(Produc_1986_scaled_df, mean), 4),
  Desv_tipica = round(sapply(Produc_1986_scaled_df, sd), 4),
  Minimo = round(sapply(Produc_1986_scaled_df, min), 4),
  Q1 = round(sapply(Produc_1986_scaled_df, quantile, probs = 0.25), 4),
  Mediana = round(sapply(Produc_1986_scaled_df, median), 4),
  Q3 = round(sapply(Produc_1986_scaled_df, quantile, probs = 0.75), 4),
  Maximo = round(sapply(Produc_1986_scaled_df, max), 4)
)

# Eliminar nombres de fila
rownames(tabla_desc_transformadas) <- NULL

# Mostrar tabla formateada
kable(
  tabla_desc_transformadas,
  caption = "Resumen descriptivo de las variables transformadas y estandarizadas",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

Como puede verse en la anterior tabla, tras la transformación y estandarización de las variables, todas presentan media 0 y desviación típica 1, lo que permite compararlas en una escala común. Este paso resulta necesario antes de aplicar el escalado multidimensional, ya que las variables originales tenían magnitudes muy diferentes y podían distorsionar el cálculo de distancias. Los rangos observados muestran que la transformación logarítmica ha suavizado la dispersión de las variables económicas de volumen, situando la mayoría de los valores entre aproximadamente -2 y 2. La tasa de desempleo, aunque no fue transformada logarítmicamente, también se estandarizó y muestra un valor máximo más elevado, lo que sugiere la presencia de algún estado con desempleo relativamente alto. En conjunto, la tabla confirma que las variables están preparadas para construir una matriz de distancias equilibrada entre estados.

Análisis de la estructura de relación entre variables con una representación adecuada.

Para atender a esta cuestión, se considera que lo más adecuado es analizar la matriz de correlaciones entre las variables transformadas y estandarizadas, y representarla mediante un correlograma. Así puedes ver qué variables se mueven juntas y si existen bloques de variables relacionadas.

# Convertir a data frame las variables transformadas y estandarizadas
Produc_1986_scaled_df <- as.data.frame(Produc_1986_scaled)

# Matriz de correlaciones
cor_produc_1986 <- cor(
  Produc_1986_scaled_df,
  use = "pairwise.complete.obs"
)

# Tabla de correlaciones
kable(
  round(cor_produc_1986, 3),
  caption = "Matriz de correlaciones entre variables transformadas y estandarizadas",
  align = "c"
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

# Correlograma
corrplot(
  cor_produc_1986,
  method = "color",
  type = "lower",
  order = "hclust",
  addCoef.col = "black",
  number.cex = 0.70,
  tl.col = "#1F2937",
  tl.srt = 35,
  tl.cex = 0.85,
  diag = FALSE,
  col = colorRampPalette(c("#B2182B", "white", "#2166AC"))(200),
  mar = c(0, 0, 2, 0),
  title = "Estructura de relación entre variables económicas. Año 1986"
)

La matriz de correlaciones muestra una estructura relacional muy marcada. Las variables vinculadas al capital público, infraestructuras, capital privado, producto bruto estatal y empleo presentan correlaciones positivas muy elevadas entre sí. Esto indica que comparten una dimensión común relacionada con el tamaño económico y la dotación infraestructural de los estados. En particular, variables como gsp y emp, o pcap y util, muestran asociaciones prácticamente lineales, lo que evidencia una elevada redundancia informativa entre estos indicadores.

En cambio, la tasa de desempleo presenta correlaciones bajas con el resto de variables, lo que sugiere que mide una dimensión distinta, más vinculada a la situación del mercado laboral que al volumen económico o infraestructural de cada estado. Por tanto, unemp puede aportar una fuente de diferenciación adicional dentro del análisis.

De cara al escalado multidimensional, esta estructura implica que las distancias entre estados estarán muy condicionadas por el bloque de variables económicas altamente correlacionadas. En consecuencia, la representación MDS tenderá a ordenar los estados principalmente según una dimensión de tamaño económico e infraestructural, mientras que la tasa de desempleo podrá introducir diferencias secundarias entre estados que, aun siendo similares en tamaño económico, presenten situaciones laborales distintas.

Matiz metodológico:La alta correlación entre varias variables no impide aplicar el MDS, pero sí debe tenerse en cuenta en la interpretación, ya que las variables redundantes pueden reforzar en exceso una misma dimensión latente dentro del cálculo de distancias.

Se aplican MDS clásico con cmdscale(), MDS no métrico con isoMDS() y Sammon con sammon().

A continuación, se realiza el análisis de escalado multidimensional utilizando tres metodologías distintas: un MDS clásico, un MDS no métrico y, por último un MDS mediante el método de Sammon, una variante del escalado no métrico que busca minimizar una función de error o Stress con una característica particular: otorga un mayor peso a las distancias pequeñas originales entre los datos. Esto permite que el mapa resultante preserve con mayor precisión las relaciones de proximidad local entre los individuos u objetos analizados.

Posteriormente, una vez obtenidos los resultados se mide la bondad de los mismos con distintas medidas de bondad (o maldad):

El GOF y las medidas de stress son indicadores utilizados para evaluar la bondad del ajuste en el escalamiento multidimensional, es decir, para comprobar hasta qué punto el mapa obtenido representa adecuadamente las disimilitudes o distancias originales entre los casos.

El GOF, o Goodness of Fit, mide la bondad del ajuste desde la perspectiva contraria: indica qué proporción de la información original queda bien representada en la solución obtenida. Se utiliza especialmente en el MDS clásico y suele calcularse a partir de los autovalores. En este caso, cuanto más alto es el valor, mejor es el ajuste. Valores cercanos a 1 indican que la representación es muy buena.

Las medidas de stress indican el grado de “maldad”, es decir, de error o desajuste de la solución. Comparan las distancias representadas en el espacio reducido con las disimilitudes originales. Por tanto, cuanto menor es el stress, mejor es el ajuste. Un valor cercano a 0 indica una representación muy buena, mientras que valores elevados indican que el mapa no reproduce bien la estructura original de los datos.

En resumen, el stress mide el error de la solución y debe ser bajo, mientras que el GOF mide la calidad de la representación y debe ser alto. Además, el ajuste también puede evaluarse visualmente mediante el diagrama de Shepard, donde una buena solución se observa cuando los puntos siguen aproximadamente una línea diagonal.

# 1. MATRIZ DE DATOS PARA MDS

# Usamos los estados como nombres de fila
rownames(Produc_1986_scaled_df) <- Produc_1986$state

# Matriz de distancias euclídeas
dist_mds <- dist(Produc_1986_scaled_df, method = "euclidean")

# 2. MDS CLÁSICO CON cmdscale()

mds_clasico <- cmdscale(
  dist_mds,
  k = 2,
  eig = TRUE
)

coord_clasico <- as.data.frame(mds_clasico$points)
colnames(coord_clasico) <- c("Dim1", "Dim2")
coord_clasico$state <- rownames(Produc_1986_scaled_df)
coord_clasico$Metodo <- "MDS clásico"

# Bondad de ajuste del MDS clásico
mds_clasico$GOF

## [1] 0.9756712 0.9756712

# 3. MDS NO MÉTRICO CON isoMDS()

mds_iso <- isoMDS(
  dist_mds,
  y = mds_clasico$points,
  k = 2,
  maxit = 100
)

## initial  value 1.683819 
## iter   5 value 1.457087
## final  value 1.452029 
## converged

coord_iso <- as.data.frame(mds_iso$points)
colnames(coord_iso) <- c("Dim1", "Dim2")
coord_iso$state <- rownames(Produc_1986_scaled_df)
coord_iso$Metodo <- "MDS no métrico"

# Stress del MDS no métrico
mds_iso$stress

## [1] 1.452029

# 4. MDS DE SAMMON CON sammon()

mds_sammon <- sammon(
  dist_mds,
  y = mds_clasico$points,
  k = 2,
  niter = 100
)

## Initial stress        : 0.00371
## stress after  10 iters: 0.00228, magic = 0.500
## stress after  20 iters: 0.00228, magic = 0.500

coord_sammon <- as.data.frame(mds_sammon$points)
colnames(coord_sammon) <- c("Dim1", "Dim2")
coord_sammon$state <- rownames(Produc_1986_scaled_df)
coord_sammon$Metodo <- "Sammon"

# Stress de Sammon
mds_sammon$stress

## [1] 0.002282167

# Tabla comparativa de ajuste

resumen_mds <- data.frame(
  Metodo = c(
    "MDS clásico - cmdscale()",
    "MDS no métrico - isoMDS()",
    "Sammon - sammon()"
  ),
  Indicador = c(
    "GOF",
    "Stress",
    "Stress"
  ),
  Valor = c(
    round(mds_clasico$GOF[1], 4),
    round(mds_iso$stress, 4),
    round(mds_sammon$stress, 4)
  )
)

kable(
  resumen_mds,
  caption = "Indicadores de ajuste de los métodos de escalado multidimensional",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

Los indicadores de ajuste obtenidos muestran que la representación bidimensional mediante escalado multidimensional es adecuada. En el MDS clásico, el valor de bondad de ajuste es muy elevado, \(GOF = 0.9757\), lo que indica que las dos dimensiones retenidas conservan una proporción muy alta de la estructura de distancias original entre los estados. Por su parte, el MDS no métrico presenta un valor de stress bajo, 1.4520, lo que sugiere que la ordenación relativa de las disimilitudes entre estados se reproduce correctamente en el espacio bidimensional.

El método de Sammon obtiene un stress especialmente reducido, 0.0023, lo que indica una representación muy fiel de las distancias, especialmente entre estados próximos. En conjunto, los tres métodos ofrecen resultados consistentes y apoyan la utilización de una solución bidimensional para representar las semejanzas y diferencias económicas entre los estados en 1986.

No obstante, los indicadores no son directamente comparables entre sí, ya que el MDS clásico utiliza una medida de bondad de ajuste, mientras que isoMDS() y sammon() utilizan medidas de stress. En el primer caso, valores más altos indican mejor ajuste; en los otros dos, valores más bajos indican mejor representación.

Se genera un gráfico comparativo de las soluciones obtenidas.

Se realiza un gráfico comparativo que permita observar la estructura de proximidad entre estados bajo los tres métodos de escalado multidimensional previamente ejecutados. La interpretación debe centrarse en la cercanía relativa entre estados y en la formación de agrupaciones, más que en el valor absoluto de los ejes, ya que la orientación y escala de las dimensiones puede variar entre métodos. En conjunto, la comparación permite valorar si las soluciones obtenidas son consistentes en la representación de las similitudes económicas entre estados.

# Gráfico comparativo
coord_todos_mds <- rbind(
  coord_clasico,
  coord_iso,
  coord_sammon
)

ggplot(
  coord_todos_mds,
  aes(
    x = Dim1,
    y = Dim2,
    label = state
  )
) +
  geom_point(
    size = 2.2,
    alpha = 0.85,
    color = "#1F4E79"
  ) +
  geom_text(
    size = 2.3,
    vjust = -0.6,
    check_overlap = TRUE
  ) +
  facet_wrap(
    ~ Metodo,
    scales = "free"
  ) +
  labs(
    title = "Comparación de soluciones de escalado multidimensional",
    subtitle = "MDS clásico, MDS no métrico y método de Sammon",
    x = "Dimensión 1",
    y = "Dimensión 2"
  ) +
  theme_minimal(base_size = 12)

Los tres mapas muestran una estructura muy parecida. Esto es importante porque indica que la solución es estable y no depende excesivamente del método utilizado.

En los tres gráficos se observa que algunos estados aparecen claramente alejados del núcleo central. Por ejemplo, California y New York se sitúan hacia la derecha, lo que probablemente refleja su mayor tamaño económico, mayor empleo, mayor producto bruto estatal y mayor dotación de capital e infraestructuras. Esto es coherente con la matriz de correlaciones, donde las variables económicas de volumen estaban muy fuertemente asociadas entre sí.

También aparecen agrupaciones de estados más pequeños o con menor volumen económico hacia la zona izquierda del gráfico. En la parte inferior aparecen estados como Mississippi, West Virginia o Louisiana, que se diferencian del resto en la segunda dimensión. Esta segunda dimensión podría estar recogiendo diferencias vinculadas a la tasa de desempleo o a una estructura económica menos alineada con el eje principal de tamaño económico.

No obstante, conviene recordar que en MDS los ejes no se interpretan de forma directa como variables originales. Lo relevante es la proximidad relativa entre estados, no tanto la orientación exacta del gráfico. Además, los mapas pueden rotarse o reflejarse sin cambiar su significado.

Se selecciona una solución principal, interpretando económicamente las dos dimensiones del mapa.

Se selecciona como solución principal el MDS clásico mediante cmdscale(), dado que presenta una bondad de ajuste muy elevada (GOF = 0.9757) y ofrece una representación bidimensional clara y consistente con las soluciones obtenidas mediante MDS no métrico y Sammon. Además, desde el criterio de parsimonia, resulta preferible al no requerir procedimientos iterativos de optimización ni transformaciones adicionales de las distancias, proporcionando una solución estable y fácilmente interpretable que reproduce adecuadamente la estructura de proximidades observada en los datos. Aunque el método de Sammon obtiene un valor de stress extremadamente reducido y el MDS no métrico también ofrece un ajuste satisfactorio, ambos generan configuraciones muy similares a la obtenida mediante el MDS clásico, sin aportar una mejora sustancial en la interpretación económica de los resultados. Por ello, se opta por el modelo más sencillo que explica adecuadamente la estructura de los datos.

La primera dimensión del mapa puede interpretarse como un eje de tamaño económico y dotación de capital e infraestructuras, separando a los estados de mayor volumen económico, como California, New York o Texas, de estados de menor dimensión económica. La segunda dimensión introduce una diferenciación adicional vinculada al perfil económico-laboral de los estados, separando algunos estados del sur y de menor dinamismo relativo de otros estados con perfiles económicos más favorables. En conjunto, el mapa permite representar de forma sintética las semejanzas y diferencias económicas entre los estados en 1986.

La interpretación de los ejes en MDS debe hacerse con prudencia, ya que las dimensiones no son variables observadas directamente. Lo más relevante es la posición relativa entre estados y la formación de agrupaciones o separaciones dentro del mapa.

Se identifican aquellos estados que aparecen como más próximos y cuáles son los más alejados.

Para contestar esta cuestión de forma rigurosa, lo mejor es no identificar los estados visualmente observando el gráfico, sino calcular las distancias entre estados en la solución MDS seleccionada, que en nuestro caso ha sido el MDS clásico.

# Top 10 pares de estados más cercanos y más alejados en la solución MDS clásica

# Coordenadas de la solución principal
coord_principal <- coord_clasico[, c("Dim1", "Dim2")]
rownames(coord_principal) <- coord_clasico$state

# Matriz de distancias euclídeas en el mapa MDS
dist_mapa <- as.matrix(dist(coord_principal))

# Eliminar la diagonal
diag(dist_mapa) <- NA

# Convertir matriz de distancias a tabla
tabla_distancias <- as.data.frame(as.table(dist_mapa))

colnames(tabla_distancias) <- c("Estado_1", "Estado_2", "Distancia")

# Eliminar NA
tabla_distancias <- tabla_distancias[!is.na(tabla_distancias$Distancia), ]

# Eliminar pares duplicados simétricos
tabla_distancias <- tabla_distancias[
  as.character(tabla_distancias$Estado_1) < as.character(tabla_distancias$Estado_2),
]

# Top 10 estados más cercanos

pares_mas_cercanos <- tabla_distancias[
  order(tabla_distancias$Distancia),
][1:10, ]

# Añadir columna de posición
pares_mas_cercanos <- data.frame(
  Puesto = 1:nrow(pares_mas_cercanos),
  pares_mas_cercanos
)

rownames(pares_mas_cercanos) <- NULL

kable(
  pares_mas_cercanos,
  caption = "Diez pares de estados más cercanos en la solución MDS clásica",
  align = "c",
  row.names = FALSE,
  digits = 3
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

# Top 10 estados más alejados

pares_mas_alejados <- tabla_distancias[
  order(tabla_distancias$Distancia, decreasing = TRUE),
][1:10, ]

# Añadir columna de posición
pares_mas_alejados <- data.frame(
  Puesto = 1:nrow(pares_mas_alejados),
  pares_mas_alejados
)

rownames(pares_mas_alejados) <- NULL

kable(
  pares_mas_alejados,
  caption = "Diez pares de estados más alejados en la solución MDS clásica",
  align = "c",
  row.names = FALSE,
  digits = 3
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

Descripción inicial de ambos bloques (variables previo escalado)

# Función descriptiva
tabla_descriptiva <- function(datos) {
  data.frame(
    Variable = names(datos),
    Media = round(sapply(datos, mean), 4),
    Desv_tipica = round(sapply(datos, sd), 4),
    Minimo = round(sapply(datos, min), 4),
    Q1 = round(sapply(datos, quantile, probs = 0.25), 4),
    Mediana = round(sapply(datos, median), 4),
    Q3 = round(sapply(datos, quantile, probs = 0.75), 4),
    Maximo = round(sapply(datos, max), 4),
    row.names = NULL
  )
}

desc_produccion <- tabla_descriptiva(bloque_produccion)
desc_labfin <- tabla_descriptiva(bloque_laboral_financiero)

kable(
  desc_produccion,
  caption = "Resumen descriptivo del bloque de producción",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1F4E79")

kable(
  desc_labfin,
  caption = "Resumen descriptivo del bloque laboral-financiero",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(0, bold = TRUE, color = "white", background = "#1F4E79")

# GRÁFICOS BLOQUE PRODUCCIÓN

# Histogramas
par(mfrow = c(2, 3))
for (v in names(bloque_produccion)) {
  hist(bloque_produccion[[v]],
       main = paste("Histograma de", v),
       xlab = v,
       col = "blue",
       border = "white")
}
par(mfrow = c(1, 1))

# Boxplots para detectar dispersión y posibles valores extremos
boxplot(bloque_produccion,
        main = "Distribución de variables bloque producción",
        las = 2,
        col = "blue")

# GRÁFICOS BLOQUE LABORAL-FINANCIERO

# Histogramas
par(mfrow = c(2, 3))
for (v in names(bloque_laboral_financiero)) {
  hist(bloque_laboral_financiero[[v]],
       main = paste("Histograma de", v),
       xlab = v,
       col = "blue",
       border = "white")
}
par(mfrow = c(1, 1))

# Boxplots para detectar dispersión y posibles valores extremos
boxplot(bloque_laboral_financiero,
        main = "Distribución de variables bloque laboral-financiero",
        las = 2,
        col = "blue")

La descripción inicial muestra que las variables del bloque de producción presentan medias y medianas relativamente próximas, lo que sugiere distribuciones moderadamente equilibradas, aunque con diferencias apreciables entre empresas. Dado que y, i y k están expresadas en logaritmos, las diferencias reflejan principalmente variaciones relativas en producción, consumos intermedios y capital real. En el bloque laboral-financiero, el empleo y los salarios presentan también una estructura relativamente estable, especialmente en el caso de w, cuya dispersión es menor.

La variable f, correspondiente al flujo de caja real, presenta un comportamiento claramente diferente. Su media es muy superior a la mediana y la desviación típica es elevada, con presencia tanto de valores negativos como de valores máximos muy altos. Esto indica una fuerte asimetría y la existencia de empresas con flujos de caja extremos. Por este motivo, aunque se mantiene en escala original en la descripción inicial, para la correlación canónica se estandarizan todas las variables. Además, debe interpretarse con cautela el papel de f, ya que sus valores extremos pueden influir en las correlaciones y en las combinaciones canónicas.

Como comprobación de robustez, podría considerarse una transformación arco seno hiperbólico como \(asinh(f)\), que permite tratar valores negativos y reduce la influencia de valores extremos, manteniendo una interpretación similar a una transformación logarítmica para valores elevados.

Se analizan las correlaciones cruzadas entre los dos grupos de variables o bloques.

# CORRELACIONES CRUZADAS ENTRE LOS DOS BLOQUES DE VARIABLES

# Bloque de producción: y, i, k
# Bloque laboral-financiero: n, w, f


# Calcular matriz de correlaciones cruzadas
cor_cruzadas <- cor(
  X_prod,
  Y_labfin,
  use = "pairwise.complete.obs"
)


# TABLA FORMATEADA

kable(
  round(cor_cruzadas, 3),
  caption = "Correlaciones cruzadas entre el bloque de producción y el bloque laboral-financiero",
  align = "c"
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  ) %>%
  column_spec(
    1,
    bold = TRUE
  )

# GRÁFICO DE CORRELACIONES CRUZADAS


corrplot(
  cor_cruzadas,
  method = "color",
  is.corr = FALSE,
  addCoef.col = "black",
  number.cex = 0.85,
  tl.col = "#1F2937",
  tl.srt = 35,
  tl.cex = 0.9,
  col = colorRampPalette(
    c("#B2182B", "white", "#2166AC")
  )(200),
  mar = c(0, 0, 2, 0),
  title = "Correlaciones cruzadas entre bloques"
)

La matriz de correlaciones cruzadas muestra asociaciones positivas entre todas las variables del bloque de producción y las variables del bloque laboral-financiero. Las relaciones más intensas se observan entre las variables productivas y el empleo (n), especialmente entre producción real y empleo (r = 0.900), consumos intermedios y empleo (r = 0.855) y capital real y empleo (r = 0.820). Esto indica que las empresas con mayor dimensión productiva tienden también a presentar mayores niveles de empleo.

Las correlaciones con los salarios (w) y con el flujo de caja real (f) son también positivas, aunque de intensidad moderada. Por tanto, las empresas más productivas tienden a mostrar mayores salarios y mayor flujo de caja, pero estas relaciones son menos directas que la existente con el empleo. En conjunto, los resultados anticipan la existencia de una relación multivariante relevante entre ambos bloques, que será sintetizada posteriormente mediante la correlación canónica.

Correlación Canónica clásica con cancor().

Se aplica una correlación canónica clásica mediante cancor() sobre los dos bloques previamente estandarizados. El objetivo es obtener combinaciones lineales del bloque de producción y del bloque laboral-financiero que maximicen la correlación entre ambos conjuntos de variables. Las correlaciones canónicas indican la intensidad de la relación entre cada par de variables canónicas, mientras que las cargas permiten interpretar qué variables contribuyen más a cada dimensión.

# Aplicar correlación canónica clásica
cca_1990 <- cancor(
  x = X_prod,
  y = Y_labfin
)

# Ver resultados principales
cca_1990

## $cor
## [1] 0.970441254 0.034395533 0.006036322
## 
## $xcoef
##           [,1]        [,2]        [,3]
## y -0.093068210  0.07691265 -0.25899607
## i  0.060269887 -0.01226041  0.25278788
## k -0.002497595 -0.07145545  0.02071595
## 
## $ycoef
##           [,1]        [,2]        [,3]
## n -0.030464156 -0.00300201  0.02904165
## w -0.010004692 -0.02618981 -0.02717680
## f -0.003901059  0.03310731 -0.02409966
## 
## $xcenter
##             y             i             k 
##  2.503173e-16  8.988947e-17 -7.223431e-17 
## 
## $ycenter
##             n             w             f 
## -1.696891e-16 -1.852605e-17  1.129686e-17

# Graficamos las correlaciones canónicas
cor_can <- cca_1990$cor

barplot(
  cor_can,
  names.arg = paste0("Dimensión ", seq_along(cor_can)),
  ylim = c(0, 1),
  col = "#1F4E79",
  main = "Correlaciones canónicas",
  ylab = "Correlación canónica"
)

abline(h = 0.5, lty = 2, col = "grey40")

Los resultados de la correlación canónica clásica muestran que la relación entre el bloque productivo (y,i,k) y el bloque laboral-financiero (n,w,f) se concentra casi exclusivamente en la primera dimensión canónica. La primera correlación canónica alcanza un valor muy elevado, \(\rho_1 = 0.970\), lo que indica una asociación lineal muy intensa entre una combinación de las variables productivas y una combinación de las variables laborales-financieras. En cambio, la segunda y tercera correlación canónica, con valores de 0.034 y 0.006, son prácticamente nulas, por lo que no aportan dimensiones interpretativas relevantes adicionales.

En cuanto a los coeficientes canónicos, en el bloque productivo la primera función canónica está principalmente determinada por y e i, mientras que k tiene un peso directo muy reducido. En el bloque laboral-financiero, la primera función canónica aparece dominada por n, seguida de forma mucho más débil por w y f. Esto sugiere que la asociación principal entre ambos bloques se articula especialmente a través de la relación entre la dinámica productiva y la variable laboral n.

No obstante, los coeficientes canónicos deben interpretarse con cautela, ya que pueden verse afectados por la escala de las variables y por la colinealidad interna dentro de cada bloque. Por ello, la interpretación sustantiva debe complementarse con las cargas canónicas y los mapas de variables. En conjunto, el resultado permite concluir que existe una estructura común muy fuerte entre ambos grupos de variables, pero que dicha estructura es esencialmente unidimensional.

Por lo tanto, existe una relación canónica muy fuerte entre producción y variables laborales-financieras, y las otras dos dimensiones son irrelevantes desde el punto de vista interpretativo. El primer eje parece estar asociado sobre todo a la conexión entre y/i en el bloque productivo y n en el bloque laboral-financiero, es decir entre producción real y consumos intermedios ,desde el bloque productivo, con la variable empleo en el bloque laboral-financiero.

Obtención e interpretación de las correlaciones canónicas.

# Correlaciones canónicas

cor_can <- cca_1990$cor

tabla_cor_can <- data.frame(
  Funcion_canonica = paste0("Par canónico ", seq_along(cor_can)),
  Correlacion_canonica = round(cor_can, 4),
  Correlacion_cuadrada = round(cor_can^2, 4)
)

kable(
  tabla_cor_can,
  caption = "Correlaciones canónicas entre el bloque de producción y el bloque laboral-financiero",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

# Coeficientes canónicos

coef_prod <- as.data.frame(round(cca_1990$xcoef, 4))
coef_labfin <- as.data.frame(round(cca_1990$ycoef, 4))

colnames(coef_prod) <- paste0("U", 1:ncol(coef_prod))
colnames(coef_labfin) <- paste0("V", 1:ncol(coef_labfin))

coef_prod$Variable <- rownames(coef_prod)
coef_labfin$Variable <- rownames(coef_labfin)

coef_prod <- coef_prod[, c("Variable", paste0("U", 1:(ncol(coef_prod) - 1)))]
coef_labfin <- coef_labfin[, c("Variable", paste0("V", 1:(ncol(coef_labfin) - 1)))]

rownames(coef_prod) <- NULL
rownames(coef_labfin) <- NULL

kable(
  coef_prod,
  caption = "Coeficientes canónicos del bloque de producción",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )

kable(
  coef_labfin,
  caption = "Coeficientes canónicos del bloque laboral-financiero",
  align = "c",
  row.names = FALSE
) %>%
  kable_styling(
    full_width = FALSE,
    position = "center",
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  row_spec(
    0,
    bold = TRUE,
    color = "white",
    background = "#1F4E79"
  )