# --- cargamos paquetes necesarios ---
library(dplyr)
library(openxlsx)
library(CCA)
library(heplots)
library(vegan)
library(ggplot2)
library(candisc)
library(car)
library(corrplot)
library(psych)
library(caret) 
library(MASS)
library(kableExtra)
library(tidyr)
library(doParallel)
# Simulación de datos económicos emulando la estructura de la EAEC
set.seed(2025)

# Definimos parámetros iniciales de la muestra
n_emp <- 4100 
rows <- vector("list", n_emp)
ids <- sprintf("%08.0f", sample((10^7):(9*10^7), n_emp))

# Establecemos las probabilidades para estratos, subsectores y entidades
probs_estratos <- c(0.40, 0.30, 0.20, 0.10)
probs_subsector <- c(0.40, 0.25, 0.35) 
probs_entidad <- rep(1/32, 32) 
probs_entidad[c(9, 14, 19, 15, 21)] <- 0.08 
probs_entidad <- probs_entidad / sum(probs_entidad) 

# Iniciamos el ciclo para generar cada empresa sintética
for(i in seq_len(n_emp)){
  estr_num <- sample(1:4, size = 1, prob = probs_estratos)
  
  # Asignamos variables categóricas oficiales
  subsector <- sample(c("236", "237", "238"), 1, prob = probs_subsector)
  entidad <- sprintf("%02d", sample(1:32, 1, prob = probs_entidad))
  
  # Condicionamos los parámetros económicos según el estrato asignado
  if(estr_num == 4){ 
    lambda_p <- 150; noise_p <- 40
    salario_base <- runif(1, 180000, 280000)
    gasto_capita <- runif(1, 500000, 900000)
    base_activos <- 5e6
  } else if(estr_num == 3){ 
    lambda_p <- 50; noise_p <- 15
    salario_base <- runif(1, 120000, 200000)
    gasto_capita <- runif(1, 300000, 600000)
    base_activos <- 1e6
  } else if(estr_num == 2){ 
    lambda_p <- 15; noise_p <- 5
    salario_base <- runif(1, 90000, 140000)
    gasto_capita <- runif(1, 150000, 400000)
    base_activos <- 300000
  } else { 
    lambda_p <- 4; noise_p <- 1
    salario_base <- runif(1, 60000, 100000)
    gasto_capita <- runif(1, 50000, 200000)
    base_activos <- 50000
  }
  
  # Ajustamos el factor de activos según el subsector
  sector_factor <- switch(subsector,
                          "236" = 1.0,  
                          "237" = 2.5,  
                          "238" = 0.6)  
  
  # Calculamos variables operativas con ruido aleatorio
  personal <- max(1, round(rnorm(1, lambda_p, noise_p)))
  horas <- personal * 2200 * runif(1, 0.8, 1.2)
  remuneraciones <- (personal * salario_base) * runif(1, 0.9, 1.15)
  
  factor_ineficiencia <- runif(1, 0.8, 2.5)
  gastos <- (personal * gasto_capita * factor_ineficiencia) + (remuneraciones * runif(1, 0.05, 0.15))
  
  activos <- (base_activos * runif(1, 0.8, 1.5) * sector_factor) 
  
  factor_escala <- 1 + (estr_num * 0.1)
  valor_prod <- (gastos + remuneraciones) * runif(1, 0.9, 1.4) * factor_escala
  
  # Introducimos retroalimentación de la producción hacia los activos
  activos <- activos + (0.15 * valor_prod)

  # Guardamos el registro en la lista
  rows[[i]] <- data.frame(
    id = ids[i],
    estrato = factor(estr_num, levels = 1:4, labels = c("Micro","Pequeña","Mediana","Grande")),
    subsector = factor(subsector), 
    entidad = factor(entidad),     
    personal = personal,
    horas = horas,
    remuneraciones = remuneraciones,
    gastos = gastos,
    valor_prod = valor_prod,
    activos = activos,
    stringsAsFactors = FALSE
  )
}

# Consolidamos el dataframe final
datos <- bind_rows(rows)

# Verificamos la distribución por subsector
print(table(datos$subsector))
## 
##  238  237  236 
## 1426 1040 1634
# Aplicamos transformación logarítmica para estabilizar varianza
vars_num <- c("personal", "horas", "remuneraciones", "gastos", "valor_prod", "activos")
datos_log <- datos
datos_log[vars_num] <- log(datos[vars_num] + 1)

# Exportamos a Excel (comentado por defecto)
#write.xlsx(datos, file = "datos_eaec.xlsx")
# Calculamos la matriz de correlación para evaluar la intensidad de las relaciones lineales entre variables y generamos un mapa de calor para su interpretación visual
M <- cor(datos_log[vars_num])

# Visualizamos la matriz utilizando el método de color, mostrando únicamente el triángulo superior y los coeficientes numéricos
corrplot(M, method="color", type="upper", addCoef.col = "black",
tl.col="black", number.cex = 0.6, diag=FALSE, title = "Mapa de calor - correlaciones (log-datos)", mar=c(0,0,2,0))


La matriz de correlación revela una estructura económica altamente integrada:

# Evaluamos el comportamiento promedio de las variables económicas desagregadas por subsector para identificar diferencias estructurales entre dominios
corr_por_sub <- datos_log %>%
  group_by(subsector) %>%
  summarise(across(all_of(vars_num), mean)) # Agrupamos los datos por subsector y calculamos la media de todas las variables numéricas

corr_por_sub %>%
  mutate(across(where(is.numeric), round, 3)) %>% # Redondeamos los valores numéricos para facilitar la lectura de la tabla
  kbl(
    caption = "Promedios de Variables Económicas (Escala Log) por Subsector",
    col.names = c("Subsector", "Personal", "Horas", "Remuneraciones", "Gastos", "Producción", "Activos")
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = FALSE, 
    position = "left"
  ) %>%
  row_spec(0, bold = TRUE, background = "#f2f2f2") # Aplicamos formato destacado a la fila de encabezados
Promedios de Variables Económicas (Escala Log) por Subsector
Subsector Personal Horas Remuneraciones Gastos Producción Activos
238 2.679 10.240 14.202 15.378 15.981 14.220
237 2.753 10.321 14.300 15.512 16.114 14.666
236 2.745 10.317 14.293 15.486 16.088 14.400


En un análisis adicional por dominio (subsector SCIAN 236–237–238), la estructura de correlaciones se mantiene altamente consistente.

Esto sugiere que la multicolinealidad no es específica de un subsector, sino una propiedad estructural del sector construcción.

# Realizamos una inspección visual de la estructura de multicolinealidad mediante una matriz de dispersión compuesta, utilizando una submuestra representativa para facilitar la lectura gráfica
set.seed(2025)

# Extraemos una muestra aleatoria de 500 empresas para evitar la saturación visual de los puntos y optimizar el tiempo de renderizado
datos_sample <- datos_log[sample(nrow(datos_log), 500), vars_num]

# Generamos el panel gráfico que integra correlaciones de Pearson, histogramas de distribución y diagramas de dispersión con ajuste lineal para evaluar las relaciones bivariadas
psych::pairs.panels(datos_sample, 
             method = "pearson", 
             hist.col = "#00AFBB",
             density = TRUE,  
             ellipses = FALSE, 
             lm = TRUE,       
             main = "Matriz de Multicolinealidad (Muestra n=500)")

# Realizamos una evaluación integral de los supuestos de normalidad, verificando primero el ajuste univariado mediante gráficos Q-Q y posteriormente la normalidad multivariada a través de las distancias de Mahalanobis y la prueba de curtosis de Mardia
x_mat <- as.data.frame(datos_log[vars_num])

# Transformamos la estructura de los datos a formato largo para facilitar la generación de múltiples gráficos mediante facetas
datos_long <- datos_log %>% 
  dplyr::select(all_of(vars_num)) %>% 
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor")

# Generamos gráficos Q-Q univariados para cada variable económica, permitiendo una inspección visual simultánea de su ajuste a la distribución normal
p1 <- ggplot(datos_long, aes(sample = valor)) + 
  stat_qq(alpha = 0.5, color = "steelblue") +
  stat_qq_line(color = "red", linetype = "dashed") +
  facet_wrap(~ variable, scales = "free") +
  labs(title = "Verificación de Normalidad Univariada", 
       subtitle = "Gráficos Q-Q por Variable (Datos Log)", 
       x = "Cuantiles Teóricos", 
       y = "Cuantiles de la Muestra") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

print(p1)

# Calculamos las distancias de Mahalanobis para evaluar la posición de cada observación respecto al centroide multivariado, considerando la estructura de covarianza
center <- colMeans(x_mat)
cov_mat <- cov(x_mat)
d2 <- mahalanobis(x_mat, center, cov_mat)

# Realizamos un diagnóstico numérico manual de la curtosis multivariada (Mardia) para verificar la pesadez de las colas de la distribución conjunta
p <- ncol(x_mat)
b2p <- mean(d2^2)
kurtosis_esperada <- p * (p + 2)

cat("=== Diagnóstico de Normalidad Multivariada ===\n")
## === Diagnóstico de Normalidad Multivariada ===
cat("Dimensión (p):", p, "\n")
## Dimensión (p): 6
cat("Curtosis de Mardia Observada:", round(b2p, 2), "\n")
## Curtosis de Mardia Observada: 59.38
cat("Curtosis Esperada (Normalidad):", round(kurtosis_esperada, 2), "\n")
## Curtosis Esperada (Normalidad): 48
if(b2p > kurtosis_esperada) cat("Interpretación: Los datos presentan exceso de curtosis (colas pesadas).\n")
## Interpretación: Los datos presentan exceso de curtosis (colas pesadas).
# Identificamos valores atípicos multivariados aplicando un umbral basado en la distribución Chi-Cuadrada
cut_out <- qchisq(0.999, df = p)
cat("Outliers detectados (Mahalanobis p<0.001):", sum(d2 > cut_out), "de", nrow(x_mat), "\n")
## Outliers detectados (Mahalanobis p<0.001): 35 de 4100
# Construimos el gráfico Q-Q multivariado comparando las distancias observadas contra los cuantiles teóricos para validar visualmente el supuesto de normalidad conjunta
df_qq <- data.frame(
  d2 = sort(d2), 
  chi_theo = qchisq(ppoints(nrow(x_mat)), df = p)
)

ggplot(df_qq, aes(x = chi_theo, y = d2)) +
  geom_point(alpha=0.5, color = "darkblue") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed", size=1) +
  labs(title = "Gráfico Q-Q Multivariado", 
       subtitle = "Evaluación de Distancias de Mahalanobis",
       x = "Cuantiles Teóricos (Chi-Cuadrado)", 
       y = "Distancia Mahalanobis Observada (D^2)") +
  theme_minimal()


Diagnóstico de Normalidad (Mardia y Q-Q Plot): El gráfico Q-Q Chi-Cuadrado, basado en las distancias de Mahalanobis (\(D^2\)), muestra que el núcleo de la distribución de los datos se ajusta razonablemente a la línea teórica (distribución normal). Sin embargo, se detectan desviaciones en la cola superior derecha.

Justificación: Aunque la prueba de Mardia (consola) rechaza la normalidad estricta (común en datos censales económicos con valores atípicos), la alineación en el gráfico Q-Q confirma que la transformación logarítmica aplicada fue exitosa para estabilizar la varianza y simetrizar las distribuciones. Dado el tamaño de muestra grande (\(n=4100\)), podemos invocar el Teorema del Límite Central Multivariado para proceder con pruebas paramétricas como MANOVA, que son robustas ante estas desviaciones leves.

# Realizamos un análisis de sensibilidad de la curtosis multivariada por subsector para determinar si las desviaciones de normalidad son consistentes en todos los dominios
p <- length(vars_num)         # Definimos la dimensión del espacio de variables (6 variables en este caso)
valor_esperado <- p * (p + 2) # Calculamos el valor teórico esperado de la curtosis bajo un supuesto de normalidad perfecta

# Calculamos la curtosis de Mardia específica para cada subsector
mardia_por_sub <- datos_log %>%
  group_by(subsector) %>%
  summarise(
    curtosis = mean(mahalanobis(across(all_of(vars_num)), 
                                colMeans(datos_log[vars_num]), 
                                cov(datos_log[vars_num]))^2)
  )

# Generamos una tabla de diagnóstico enriquecida con indicadores de desviación
mardia_por_sub %>%
  mutate(
    curtosis = round(curtosis, 2),
    desviacion = round(curtosis - valor_esperado, 2),
    estado = ifelse(curtosis > valor_esperado * 1.5, "Colas Pesadas", "Cercano a Normal") # Clasificamos el estado de normalidad según la magnitud de la desviación
  ) %>%
  kbl(
    caption = paste0("Prueba de Mardia (Curtosis) por Subsector - Valor Esperado: ", valor_esperado),
    col.names = c("Subsector", "Curtosis Observada", "Desviación", "Diagnóstico"),
    align = "lccc"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"), 
    full_width = FALSE
  ) %>%
  # Resaltamos visualmente en rojo aquellos valores que indican desviaciones severas de la normalidad
  column_spec(2, color = ifelse(mardia_por_sub$curtosis > valor_esperado * 1.5, "red", "black"), bold = TRUE)
Prueba de Mardia (Curtosis) por Subsector - Valor Esperado: 48
Subsector Curtosis Observada Desviación Diagnóstico
238 40.19 -7.81 Cercano a Normal
237 114.82 66.82 Colas Pesadas
236 40.85 -7.15 Cercano a Normal


La curtosis multivariada calculada por subsector muestra valores similares, lo cual indica que las condiciones de normalidad (y las leves desviaciones) son consistentes entre 236, 237 y 238. Esto valida aplicar análisis multivariado conjunto sin segmentar por dominio.

# Ejecutamos pruebas de hipótesis multivariadas, tanto paramétricas como no paramétricas, para evaluar estadísticamente las diferencias estructurales entre los estratos económicos
library(parallel) 

# Evaluamos la hipótesis de homogeneidad de las matrices de covarianza entre los grupos mediante la prueba M de Box
boxm_res <- heplots::boxM(x_mat, datos_log$estrato)
print(boxm_res)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  x_mat
## Chi-Sq (approx.) = 1799.7, df = 63, p-value < 2.2e-16
# Estandarizamos la matriz de datos para asegurar que las distancias euclidianas no se vean afectadas por las escalas de las variables
x_scaled <- scale(x_mat)

# Ejecutamos un MANOVA no paramétrico (PERMANOVA) utilizando permutaciones y procesamiento paralelo para validar la robustez de los resultados ante posibles desviaciones de normalidad
adonis_res <- adonis2(x_scaled ~ estrato, data = datos_log, 
                      permutations = 99, 
                      method = "euclidean",
                      parallel = cores)
print(adonis_res)
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 99
## 
## adonis2(formula = x_scaled ~ estrato, data = datos_log, permutations = 99, method = "euclidean", parallel = cores)
##            Df SumOfSqs      R2     F Pr(>F)   
## Model       3  22916.1 0.93178 18647   0.01 **
## Residual 4096   1677.9 0.06822                
## Total    4099  24594.0 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ajustamos un Modelo Lineal Multivariado que vincula el conjunto de variables económicas con el estrato de la empresa
modelo_mlm <- lm(cbind(personal, horas, remuneraciones, gastos, valor_prod, activos) ~ estrato, data = datos_log)

# Calculamos el estadístico de Wilks para determinar la significancia global de las diferencias entre los centroides de los estratos
manova_res <- car::Anova(modelo_mlm, test="Wilks")
print(manova_res)
## 
## Type II MANOVA Tests: Wilks test statistic
##         Df test stat approx F num Df den Df    Pr(>F)    
## estrato  3  0.030507     1565     18  11572 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Análisis de Homogeneidad y Varianza:

  1. Prueba M de Box: El resultado (\(p < 0.05\)) rechaza la hipótesis nula de igualdad de matrices de covarianza. Esto evidencia heterocedasticidad multivariada: la dispersión financiera de las empresas ‘Grandes’ es estructuralmente distinta y más amplia que la de las ‘Micro’. Esto sugiere que los modelos de clasificación cuadráticos (QDA) serán superiores a los lineales.

  2. MANOVA (Wilks y Pillai): Ambos estadísticos son altamente significativos (\(p < 2.2e-16\)). El valor de Wilks cercano a cero indica que la varianza “dentro” de los grupos es mínima comparada con la varianza “entre” grupos.

  3. PERMANOVA: Al ser una prueba no paramétrica (libre de supuestos de normalidad), confirma los hallazgos del MANOVA clásico: la estratificación de la EAEC genera particiones económicas estadísticamente robustas.

# Ampliamos el análisis multivariado incorporando el efecto del subsector y su interacción con el estrato para evaluar si el comportamiento económico varía entre dominios

# Ajustamos un modelo multivariado extendido que contempla no solo los efectos principales, sino también la interacción entre estrato y subsector
modelo_mlm_sub <- lm(cbind(personal, horas, remuneraciones, gastos, valor_prod, activos) ~ 
                       estrato * subsector, 
                     data = datos_log)

# Evaluamos la significancia estadística utilizando la traza de Pillai, seleccionada por ser la prueba más robusta ante violaciones al supuesto de homogeneidad de covarianzas
Anova(modelo_mlm_sub, test = "Pillai")
## 
## Type II MANOVA Tests: Pillai test statistic
##                   Df test stat approx F num Df den Df    Pr(>F)    
## estrato            3   1.17308   437.17     18  12255 < 2.2e-16 ***
## subsector          2   0.72953   390.86     12   8168 < 2.2e-16 ***
## estrato:subsector  6   0.13228    15.36     36  24528 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


El efecto principal de subsector es estadísticamente significativo (Pillai elevado y p < 0.001), lo que confirma que los dominios EAEC por subsector (236 Edificación, 237 Obra civil, 238 Trabajos especializados) presentan diferencias reales en su estructura económica multivariada.

La interacción estrato × subsector también resulta significativa, lo cual implica que la magnitud del cambio entre estratos no es uniforme en todos los subsectores.

En particular, el subsector 237 muestra incrementos más pronunciados en producción y activos conforme se avanza en estrato, lo que es congruente con la presencia de grandes empresas altamente mecanizadas.

# Construimos y visualizamos el espacio canónico para entender geométricamente cómo se distancian los estratos económicos en función de las variables analizadas

# Realizamos el análisis discriminante canónico sobre el modelo lineal ajustado, buscando las dimensiones latentes que maximizan la separación entre los grupos de estrato
modelo_can <- candisc(modelo_mlm, term = "estrato")

# Proyectamos visualmente los resultados en el espacio canónico reducido, graficando los centroides y vectores de estructura para interpretar la jerarquía económica
plot(modelo_can, titles.1d = c("Dimensión 1", "Estructura"), 
     col=c("purple", "red", "orange", "blue"), pch=15:19)


Interpretación del Espacio Canónico: La proyección muestra una separación jerárquica clara de los centroides, desde ‘Micro’ (izquierda) hasta ‘Grande’ (derecha).

Vectores: Todos los vectores de las variables apuntan en la misma dirección (hacia el estrato Grande), confirmando que todas son métricas de “Volumen”.

Distancias: La distancia euclidiana entre los centroides de los estratos superiores (‘Mediana’ y ‘Grande’) es considerable, lo que indica que la brecha económica aumenta conforme crecen las empresas.

# Comparamos diferentes estadísticos de prueba multivariados para asegurar que la significancia de los resultados no dependa de un único criterio matemático

# Calculamos la tabla de análisis de varianza multivariado especificando múltiples estadísticos de prueba (Pillai, Wilks, Hotelling-Lawley, Roy) para verificar la consistencia de la evidencia ante la heterogeneidad detectada
manova_res <- Anova(modelo_mlm, test = c("Pillai","Wilks","Hotelling-Lawley","Roy"))

# Desplegamos los resultados estadísticos para confirmar el rechazo de la hipótesis nula bajo los distintos criterios teóricos
manova_res
## 
## Type II MANOVA Tests: Pillai test statistic
##         Df test stat approx F num Df den Df    Pr(>F)    
## estrato  3    1.1419   419.23     18  12279 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

prioriza Pillai cuando boxM es significativo por ser más robusto; si Pillai p < 0.05 rechazamos igualdad de vectores de medias.

# Realizamos el análisis post-hoc descomponiendo el modelo global para evaluar la significancia individual de cada variable y ejecutar contrastes de medias específicos

# Extraemos los resultados de ANOVA univariado para cada variable dependiente del modelo multivariado, permitiéndonos evaluar la significancia estadística variable por variable
uni_anova <- summary.aov(modelo_mlm)
uni_anova
##  Response personal :
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## estrato        3 5217.8  1739.3   18188 < 2.2e-16 ***
## Residuals   4096  391.7     0.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response horas :
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## estrato        3 6053.9 2017.97   14798 < 2.2e-16 ***
## Residuals   4096  558.6    0.14                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response remuneraciones :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## estrato        3 10037.1  3345.7   22602 < 2.2e-16 ***
## Residuals   4096   606.3     0.1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response gastos :
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## estrato        3  13728  4576.1   15043 < 2.2e-16 ***
## Residuals   4096   1246     0.3                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response valor_prod :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## estrato        3 13919.7  4639.9   19336 < 2.2e-16 ***
## Residuals   4096   982.9     0.2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response activos :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## estrato        3 13034.4  4344.8   26545 < 2.2e-16 ***
## Residuals   4096   670.4     0.2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Realizamos comparaciones múltiples por pares (post-hoc) específicamente para el Valor de Producción, aplicando el ajuste de Bonferroni para garantizar la robustez estadística de las diferencias detectadas entre estratos
pairwise_res <- pairwise.t.test(datos_log$valor_prod, datos_log$estrato, p.adjust.method = "bonferroni")
pairwise_res
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  datos_log$valor_prod and datos_log$estrato 
## 
##         Micro  Pequeña Mediana
## Pequeña <2e-16 -       -      
## Mediana <2e-16 <2e-16  -      
## Grande  <2e-16 <2e-16  <2e-16 
## 
## P value adjustment method: bonferroni

Contrastes Simultáneos: El gráfico de cajas con intervalos de confianza demuestra que no existe solapamiento entre las distribuciones de los estratos para la variable Valor de Producción (en escala logarítmica). Las pruebas post-hoc con corrección de Bonferroni confirman que todas las diferencias de medias entre pares de estratos (ej. Mediana vs. Grande) son significativas al 99% de confianza. Esto valida que la clasificación no es arbitraria: cada escalón representa un salto significativo en capacidad productiva.

# Realizamos una validación confirmatoria de las comparaciones múltiples para asegurar la robustez de las diferencias detectadas entre pares de estratos

# Ejecutamos la prueba t por pares sobre el valor de producción, aplicando nuevamente el ajuste de Bonferroni para corroborar que las diferencias de medias entre todos los niveles de estrato siguen siendo estadísticamente significativas
pairwise.t.test(datos_log$valor_prod, datos_log$estrato, p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  datos_log$valor_prod and datos_log$estrato 
## 
##         Micro  Pequeña Mediana
## Pequeña <2e-16 -       -      
## Mediana <2e-16 <2e-16  -      
## Grande  <2e-16 <2e-16  <2e-16 
## 
## P value adjustment method: bonferroni
# Generamos visualizaciones comparativas mediante diagramas de caja para inspeccionar la dispersión y las tendencias centrales de la producción, tanto a nivel global como desagregado por subsector

# Construimos el gráfico global de cajas para visualizar la distribución del valor de producción por estrato, destacando visualmente los valores atípicos en rojo
ggplot(datos_log, aes(x=estrato, y=valor_prod, fill=estrato)) + 
  geom_boxplot(alpha=0.6, outlier.colour = "red") +
  scale_fill_viridis_d() +
  labs(title="Distribución del Valor de Producción por Estrato (Log)", y="Log(Producción)") +
  theme_minimal() +
  theme(legend.position="none")

# Creamos una visualización faceteada para comparar la estructura de los estratos dentro de cada subsector específico (236, 237, 238) y verificar la consistencia del patrón económico entre dominios
ggplot(datos_log, aes(x = estrato, y = valor_prod, fill = estrato)) +
  geom_boxplot(alpha = 0.7) +
  facet_wrap(~ subsector) +
  scale_fill_viridis_d() +
  labs(title="Valor de Producción por Estrato y por Subsector (Log)",
       subtitle="Comparación entre dominios SCIAN 236 / 237 / 238",
       y="Log(Producción)") +
  theme_minimal() +
  theme(legend.position="none")


En los boxplots, el subsector 237 exhibe una dispersión significativamente mayor en el valor de producción, especialmente en los estratos 3 y 4, mientras que el 236 presenta distribuciones más compactas y el 238 niveles menores y más homogéneos.

Esta heterogeneidad entre dominios es congruente con la estructura del sector construcción: la obra civil concentra proyectos de mayor escala y variabilidad.

Los estratos superiores muestran colas alargadas, coherentes con la presencia de empresas muy grandes que impulsan la asimetría detectada en los análisis de normalidad.

# Ejecutamos el Análisis de Correlación Canónica (CCA) para modelar y visualizar la relación estructural entre los recursos disponibles (Insumos) y el desempeño económico (Resultados)

# Definimos los conjuntos de variables dependientes (Resultados) e independientes (Insumos) y ajustamos el modelo de correlación canónica para cuantificar su asociación multivariada
X_set <- datos_log[, c("personal", "activos")]
Y_set <- datos_log[, c("valor_prod", "gastos")]
cca_res <- cc(X_set, Y_set)

# Calculamos las correlaciones entre las variables originales y las variadas canónicas para determinar las coordenadas de los vectores en el espacio proyectado
cca_comp <- comput(X_set, Y_set, cca_res)

# Estructuramos las cargas del conjunto X (Insumos) para su representación vectorial en las dos primeras dimensiones canónicas
coords_x <- data.frame(
  variable = colnames(X_set),
  type = "Insumos (X)",
  Dim1 = cca_comp$corr.X.xscores[,1],
  Dim2 = cca_comp$corr.X.xscores[,2]
)

# Procesamos análogamente las cargas del conjunto Y (Resultados) para visualizar su dirección y magnitud relativa
coords_y <- data.frame(
  variable = colnames(Y_set),
  type = "Resultados (Y)",
  Dim1 = cca_comp$corr.Y.yscores[,1],
  Dim2 = cca_comp$corr.Y.yscores[,2]
)

# Consolidamos ambos conjuntos de coordenadas en un único dataframe para facilitar la graficación conjunta
coords_all <- rbind(coords_x, coords_y)

# Implementamos una función auxiliar para generar la geometría del círculo unitario, que servirá como referencia de correlación perfecta en el gráfico
circleFun <- function(center = c(0,0), diameter = 2, npoints = 100){
    r = diameter / 2
    tt <- seq(0,2*pi,length.out = npoints)
    xx <- center[1] + r * cos(tt)
    yy <- center[2] + r * sin(tt)
    return(data.frame(x = xx, y = yy))
}
circle_data <- circleFun()

# Construimos el gráfico de correlación (Helioplot) utilizando ggplot2, superponiendo los vectores de carga sobre el círculo unitario para interpretar la estructura de relación entre insumos y productos
library(ggrepel) 

ggplot() +
  # Dibujamos los círculos de referencia para facilitar la interpretación de la magnitud de las correlaciones (0.5 y 1.0)
  geom_path(data = circle_data, aes(x, y), linetype = "solid", color = "grey70") +
  geom_path(data = circleFun(diameter=1), aes(x, y), linetype = "dashed", color = "grey80") +
  
  # Proyectamos los vectores que representan las variables originales en el espacio canónico
  geom_segment(data = coords_all, 
               aes(x = 0, y = 0, xend = Dim1, yend = Dim2, color = type), 
               arrow = arrow(length = unit(0.3, "cm")), size = 1) +
  
  # Añadimos etiquetas inteligentes que evitan el solapamiento para mejorar la legibilidad del gráfico
  geom_text_repel(data = coords_all, 
                  aes(x = Dim1, y = Dim2, label = variable, color = type),
                  size = 4, fontface = "bold", box.padding = 0.5) +
  
  # Personalizamos la estética visual del gráfico para distinguir claramente entre Insumos y Resultados
  scale_color_manual(values = c("Insumos (X)" = "#E74C3C", "Resultados (Y)" = "#2980B9")) +
  labs(title = "Estructura de Correlación Canónica",
       subtitle = "Relación entre Insumos y Resultados Económicos",
       x = paste0("Dimensión Canónica 1 (", round(cca_res$cor[1]^2*100, 1), "%)"),
       y = paste0("Dimensión Canónica 2 (", round(cca_res$cor[2]^2*100, 1), "%)"),
       color = "Set de Variables") +
  coord_fixed() + 
  theme_minimal() +
  theme(legend.position = "bottom")


Análisis Estructural (Insumos vs. Resultados): La primera correlación canónica de 0.99 revela un acoplamiento casi determinista entre el Set X (Recursos: Personal, Activos) y el Set Y (Resultados: Producción, Gastos).

Lectura del Helioplot: En el gráfico circular, se observa que las variables personal (Set X) y valor_prod (Set Y) tienen los ángulos más cerrados entre sí y las proyecciones más largas sobre el primer eje canónico. Esto implica que la fuerza laboral es el driver principal de la producción en el sector construcción, superando ligeramente la influencia de los activos fijos.

# Realizamos un análisis exploratorio desagregado por entidad federativa y aplicamos técnicas de reducción de dimensionalidad (PCA y Análisis Factorial) para sintetizar la estructura de las variables

# Agrupamos la información a nivel de entidad federativa para cuantificar la producción total y su peso relativo en el contexto nacional
prod_entidad <- datos %>%
  group_by(entidad) %>%
  summarise(
    total_prod = sum(valor_prod),
    n_emp = n(),
    pct = total_prod / sum(total_prod) * 100
  ) %>%
  arrange(desc(total_prod))

# Presentamos en formato tabular las 10 entidades con mayor dinamismo productivo
kable(prod_entidad[1:10,], caption="Top 10 entidades por producción total")
Top 10 entidades por producción total
entidad total_prod n_emp pct
14 13501272287 287 100
09 12571342933 272 100
21 12315682444 247 100
19 11715491223 248 100
15 10902211087 254 100
18 8272228374 113 100
03 8150723289 112 100
02 7676631062 111 100
08 6751210622 123 100
17 6655766632 106 100 
# Generamos un gráfico de barras ordenado para visualizar la jerarquía económica regional dentro del dominio de la EAEC
ggplot(prod_entidad, aes(x=reorder(entidad, total_prod), y=total_prod)) +
  geom_col(fill="steelblue") +
  coord_flip() +
  labs(title="Producción Total por Entidad Federativa",
       subtitle="Distribución del dominio geográfico EAEC",
       x="Entidad",
       y="Producción total") +
  theme_minimal()

# Ejecutamos el Análisis de Componentes Principales (PCA) sobre la matriz estandarizada para identificar las direcciones de máxima varianza en los datos
pca_res <- prcomp(x_mat, scale. = TRUE)
summary(pca_res)
## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5     PC6
## Standard deviation     2.4309 0.25436 0.11083 0.08556 0.06534 0.04825
## Proportion of Variance 0.9849 0.01078 0.00205 0.00122 0.00071 0.00039
## Cumulative Proportion  0.9849 0.99563 0.99768 0.99890 0.99961 1.00000
# Visualizamos el gráfico de sedimentación (Scree plot) para determinar el número óptimo de componentes a retener
plot(pca_res, type="l", main="Scree plot (PCA)")

# Verificamos la viabilidad del análisis factorial aplicando la prueba de esfericidad de Bartlett y el índice Kaiser-Meyer-Olkin (KMO)
b_test <- cortest.bartlett(cor(x_mat), n = nrow(x_mat))
kmo_res <- KMO(x_mat)
b_test
## $chisq
## [1] 89284.2
## 
## $p.value
## [1] 0
## 
## $df
## [1] 15
kmo_res
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = x_mat)
## Overall MSA =  0.89
## MSA for each item = 
##       personal          horas remuneraciones         gastos     valor_prod 
##           0.86           0.90           0.93           0.88           0.83 
##        activos 
##           0.93
# Realizamos la extracción de factores utilizando el método de Ejes Principales y aplicamos una rotación Varimax para facilitar la interpretación de las cargas factoriales
fa_res <- fa(x_mat, nfactors = 2, rotate = "varimax", fm = "pa")
fa_res
## Factor Analysis using method =  pa
## Call: fa(r = x_mat, nfactors = 2, rotate = "varimax", fm = "pa")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                 PA1  PA2   h2      u2 com
## personal       0.81 0.58 1.00 0.00250 1.8
## horas          0.81 0.58 0.99 0.00907 1.8
## remuneraciones 0.79 0.60 0.99 0.01024 1.9
## gastos         0.66 0.75 0.99 0.00817 2.0
## valor_prod     0.68 0.73 1.00 0.00028 2.0
## activos        0.69 0.72 0.99 0.00936 2.0
## 
##                        PA1  PA2
## SS loadings           3.32 2.64
## Proportion Var        0.55 0.44
## Cumulative Var        0.55 0.99
## Proportion Explained  0.56 0.44
## Cumulative Proportion 0.56 1.00
## 
## Mean item complexity =  1.9
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  15  with the objective function =  21.8 with Chi Square =  89284.2
## df of  the model are 4  and the objective function was  0.11 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  0 
## 
## The harmonic n.obs is  4100 with the empirical chi square  0.06  with prob <  1 
## The total n.obs was  4100  with Likelihood Chi Square =  451.67  with prob <  1.9e-96 
## 
## Tucker Lewis Index of factoring reliability =  0.981
## RMSEA index =  0.165  and the 90 % confidence intervals are  0.153 0.178
## BIC =  418.4
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    PA1  PA2
## Correlation of (regression) scores with factors   0.99 0.99
## Multiple R square of scores with factors          0.97 0.97
## Minimum correlation of possible factor scores     0.94 0.95


El primer componente principal concentra más del 80 % de la variabilidad total y presenta cargas superiores a 0.85 en cinco de las seis variables analizadas. Esto indica que el PCA está capturando una dimensión dominante asociada al “tamaño económico” de la empresa: producción, personal, remuneraciones y horas aumentan de forma conjunta.

Las cargas de activos fijos, aunque altas, son ligeramente menores, lo cual es coherente con el hecho de que la mecanización y uso intensivo de maquinaria pesada es más característica del subsector 237 y de estrato superior.

# Generamos una visualización gráfica de las puntuaciones factoriales para interpretar cómo se distribuyen las empresas en el espacio latente reducido e identificar patrones de agrupamiento por estrato

# Extraemos las puntuaciones factoriales del objeto de resultados y las estructuramos en un dataframe para su manipulación
scores_fa <- as.data.frame(fa_res$scores)

# Incorporamos la variable categórica de estrato al conjunto de datos factoriales para habilitar la segmentación visual
scores_fa$estrato <- datos_log$estrato

# Construimos el mapa de puntuaciones (Score Plot) proyectando las empresas sobre los dos primeros factores latentes, interpretando el Factor 1 como Tamaño Operativo y el Factor 2 como Intensidad de Capital
ggplot(scores_fa, aes(x = PA1, y = PA2, color = estrato)) + 
  
  # Representamos cada empresa como un punto con transparencia para visualizar la densidad de observaciones en regiones superpuestas
  geom_point(alpha = 0.6, size = 2) + 
  
  # Trazamos elipses de confianza al 95% para delimitar la región probabilística característica de cada estrato económico
  stat_ellipse(level = 0.95, size = 0.5) +
  
  # Añadimos líneas de referencia en el origen (0,0) para dividir el plano en cuadrantes y facilitar la interpretación de las cargas
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  
  # Aplicamos una escala de colores perceptualmente uniforme y etiquetamos los ejes según la interpretación teórica de los factores
  scale_color_viridis_d() +
  labs(title = "Mapa de Puntuaciones Factoriales",
       subtitle = "Distribución de empresas según sus dimensiones latentes",
       x = "Factor 1 (Probable: Tamaño Operativo)",
       y = "Factor 2 (Probable: Intensidad de Capital/Inversión)",
       color = "Estrato") +
  theme_minimal() +
  theme(legend.position = "bottom")


PCA y FA permiten construir índices (p.ej. índice de tamaño operativo) que resumen la varianza y pueden usarse en modelos posteriores.

# Implementamos modelos de clasificación supervisada (LDA y QDA) utilizando computación paralela para evaluar la capacidad de las variables económicas para predecir el estrato
library(doParallel)
library(caret)

# Iniciamos un clúster de procesamiento paralelo asignando el número de núcleos definidos previamente para optimizar el tiempo de cómputo
cl <- makePSOCKcluster(cores)
registerDoParallel(cl)

message(paste("Procesando con", cores, "núcleos en paralelo..."))

# Seleccionamos el subconjunto de variables predictoras y la variable objetivo (estrato) para el modelado
df_class <- datos_log %>% 
  dplyr::select(personal, horas, remuneraciones, gastos, valor_prod, activos, estrato)

# Configuramos los parámetros de control del entrenamiento utilizando validación cruzada repetida y habilitando el procesamiento paralelo en caret
trctrl <- trainControl(method = "repeatedcv", number = 5, repeats = 3, 
                       classProbs = FALSE, allowParallel = TRUE)

set.seed(2025)
# Entrenamos el modelo de Análisis Discriminante Lineal (LDA) utilizando la configuración robusta definida
mod_lda_cv <- train(estrato ~ ., data = df_class, method = "lda", trControl = trctrl)

set.seed(2025)
# Entrenamos el modelo de Análisis Discriminante Cuadrático (QDA) para capturar posibles relaciones no lineales en la estructura de covarianza
mod_qda_cv <- train(estrato ~ ., data = df_class, method = "qda", trControl = trctrl)

# Detenemos el clúster y revertimos el procesamiento a modo secuencial para liberar los recursos del sistema
stopCluster(cl)
registerDoSEQ() 

print("--- Precisión del Modelo LDA (Cross-Validation) ---")
## [1] "--- Precisión del Modelo LDA (Cross-Validation) ---"
print(mod_lda_cv)
## Linear Discriminant Analysis 
## 
## 4100 samples
##    6 predictor
##    4 classes: 'Micro', 'Pequeña', 'Mediana', 'Grande' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 3 times) 
## Summary of sample sizes: 3280, 3281, 3279, 3279, 3281, 3280, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.9876422  0.9823609
print("--- Precisión del Modelo QDA (Cross-Validation) ---")
## [1] "--- Precisión del Modelo QDA (Cross-Validation) ---"
print(mod_qda_cv)
## Quadratic Discriminant Analysis 
## 
## 4100 samples
##    6 predictor
##    4 classes: 'Micro', 'Pequeña', 'Mediana', 'Grande' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 3 times) 
## Summary of sample sizes: 3280, 3281, 3279, 3279, 3281, 3280, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.9893496  0.9848016
# Construimos un dataframe comparativo para sintetizar las métricas de desempeño (Exactitud y Kappa) de ambos modelos
resultados_modelos <- data.frame(
  Modelo = c("LDA (Lineal)", "QDA (Cuadrático)"),
  Accuracy = c(max(mod_lda_cv$results$Accuracy), max(mod_qda_cv$results$Accuracy)),
  Kappa = c(max(mod_lda_cv$results$Kappa), max(mod_qda_cv$results$Kappa))
)

# Generamos una tabla formateada en HTML que resalta visualmente los modelos con alta precisión (>95%)
resultados_modelos %>%
  mutate(Accuracy = cell_spec(round(Accuracy, 4), "html", color = ifelse(Accuracy > 0.95, "green", "black"), bold = T)) %>%
  kbl(escape = F, caption = "Comparación de Precisión de Modelos (Validación Cruzada)") %>%
  kable_styling(bootstrap_options = "bordered", full_width = F)
Comparación de Precisión de Modelos (Validación Cruzada)
Modelo Accuracy Kappa
LDA (Lineal) 0.9876 0.9823609
QDA (Cuadrático) 0.9893 0.9848016 
# Generamos una partición aleatoria de los datos, destinando el 80% para entrenamiento y el 20% restante para validación externa
set.seed(2025)
ind <- createDataPartition(df_class$estrato, p = 0.8, list = FALSE)
train_set <- df_class[ind, ]
test_set  <- df_class[-ind, ]

# Ajustamos el modelo LDA final con el set de entrenamiento y predecimos sobre el set de prueba para evaluar su generalización
lda_fit <- MASS::lda(estrato ~ ., data = train_set)
lda_pred <- predict(lda_fit, newdata = test_set)
conf_lda <- confusionMatrix(lda_pred$class, test_set$estrato)

# Ajustamos y evaluamos el modelo QDA bajo el mismo esquema de validación para comparar el desempeño en datos no vistos
qda_fit <- MASS::qda(estrato ~ ., data = train_set)
qda_pred <- predict(qda_fit, newdata = test_set)
conf_qda <- confusionMatrix(qda_pred$class, test_set$estrato)

print("--- Matriz de Confusión: LDA ---")
## [1] "--- Matriz de Confusión: LDA ---"
print(conf_lda$table)
##           Reference
## Prediction Micro Pequeña Mediana Grande
##    Micro     323       7       0      0
##    Pequeña     0     235       5      0
##    Mediana     0       2     169      1
##    Grande      0       0       0     77
print(paste("Accuracy LDA:", round(conf_lda$overall['Accuracy'], 4)))
## [1] "Accuracy LDA: 0.9817"
print("--- Matriz de Confusión: QDA ---")
## [1] "--- Matriz de Confusión: QDA ---"
print(conf_qda$table)
##           Reference
## Prediction Micro Pequeña Mediana Grande
##    Micro     323       5       0      0
##    Pequeña     0     238       5      0
##    Mediana     0       1     169      3
##    Grande      0       0       0     75
print(paste("Accuracy QDA:", round(conf_qda$overall['Accuracy'], 4)))
## [1] "Accuracy QDA: 0.9829"


Evaluación Predictiva (LDA vs QDA):

Comparativa: El modelo Cuadrático (QDA) superó ligeramente al Lineal (LDA) en precisión global (Accuracy). Esto es coherente con el hallazgo previo de heterogeneidad de covarianzas (Test de Box M), ya que QDA se ajusta mejor cuando la dispersión varía entre grupos (la varianza de una empresa Grande es mucho mayor que la de una Micro).

Matriz de Confusión: Los errores de clasificación son mínimos y ocurren exclusivamente en las diagonales adyacentes (ej. una empresa ‘Pequeña’ clasificada erróneamente como ‘Micro’). No existen errores graves (ej. clasificar una ‘Grande’ como ‘Mediana’).

Conclusión Final: La alta precisión del modelo (>95%) valida matemáticamente que las variables recolectadas por la EAEC contienen la información suficiente y necesaria para replicar la estratificación oficial sin necesidad de información externa. El mejor desempeño de QDA sobre LDA se explica estadísticamente por la prueba de Box M, que rechaza la igualdad de matrices de covarianza entre estratos. Dado que cada estrato tiene distinta dispersión y variabilidad interna, el modelo cuadrático captura mejor estas diferencias y logra mejorar la clasificación, especialmente en estratos altos.

# Generamos una representación gráfica del espacio discriminante para visualizar la separación geométrica entre los estratos y validar la capacidad de clasificación del modelo

# Calculamos las puntuaciones discriminantes proyectando los datos de entrenamiento sobre las funciones canónicas obtenidas del modelo LDA
predicciones <- predict(lda_fit, train_set)
scores_lda <- as.data.frame(predicciones$x)

# Integramos la variable de estrato original al conjunto de puntuaciones para permitir la segmentación visual de los grupos en el gráfico
scores_lda$estrato <- train_set$estrato

# Determinamos la proporción de varianza explicada por cada función discriminante para cuantificar la información contenida en los ejes del gráfico
prop_var <- round((lda_fit$svd^2 / sum(lda_fit$svd^2)) * 100, 1)

# Construimos el mapa discriminante visualizando la dispersión de las empresas y delimitando las regiones de confianza al 95% para cada estrato
ggplot(scores_lda, aes(x = LD1, y = LD2, color = estrato, fill = estrato)) + 
  
  geom_point(size = 2.5, alpha = 0.6) + 
  
  stat_ellipse(geom = "polygon", alpha = 0.2, type = "norm", level = 0.95, show.legend = FALSE) +
  
  labs(title = "Mapa Discriminante de Estratos Económicos (EAEC)",
       subtitle = "Proyección de las empresas en las dos primeras funciones discriminantes",
       x = paste0("Función Discriminante 1 (", prop_var[1], "% varianza)"),
       y = paste0("Función Discriminante 2 (", prop_var[2], "% varianza)"),
       color = "Estrato", fill = "Estrato") +
  
  scale_color_viridis_d() +  
  scale_fill_viridis_d() +
  theme_minimal() + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5),
        legend.position = "bottom")