variables <- df_ccaa[, -1]  
mahal_dist <- mahalanobis(variables, 
                         colMeans(variables), 
                         cov(variables))

p_values <- pchisq(mahal_dist, df = ncol(variables), lower.tail = FALSE)

outliers_table <- data.frame(
  CCAA = df_ccaa$CCAA,
  D2 = round(mahal_dist, 2),
  p_value = round(p_values, 2)
)

kable(outliers_table, caption = "Detección de Outliers mediante Distancia de Mahalanobis")

outliers <- which(p_values < 0.05)
if(length(outliers) > 0) {
  outliers_df <- data.frame(
    "Comunidades Outliers" = df_ccaa[outliers, "CCAA"],
    "Valor p" = round(p_values[outliers], 4)
  )
  kable(outliers_df, caption = "Comunidades identificadas como outliers (p < 0.05)")
} else {
  no_outliers_df <- data.frame(
    Resultado = "No se detectaron outliers significativos (p < 0.05)"
  )
  kable(no_outliers_df, caption = "Resultado de detección de outliers")
}

# Calcular matriz de distancias euclídeas
dist_matrix <- dist(variables, method = "euclidean")
dist_subset <- as.matrix(dist_matrix)[1:5, 1:5]
kable(round(dist_subset, 2), caption = "Matriz de distancias euclídeas (primeras 5 filas y columnas)")

par(mfrow = c(2, 3))

hc_ward <- hclust(dist_matrix, method = "ward.D2")
plot(hc_ward, main = "Método Ward", 
     xlab = "Comunidades Autónomas", ylab = "Distancia")
rect.hclust(hc_ward, k = 2, border = "red")

hc_single <- hclust(dist_matrix, method = "single")
plot(hc_single, main = "Vecino más cercano", 
     xlab = "Comunidades Autónomas", ylab = "Distancia")
rect.hclust(hc_single, k = 2, border = "red")

hc_complete <- hclust(dist_matrix, method = "complete")
plot(hc_complete, main = "Vecino más lejano", 
     xlab = "Comunidades Autónomas", ylab = "Distancia")
rect.hclust(hc_complete, k = 2, border = "red")

hc_average <- hclust(dist_matrix, method = "average")
plot(hc_average, main = "Vinculación promedio", 
     xlab = "Comunidades Autónomas", ylab = "Distancia")
rect.hclust(hc_average, k = 2, border = "red")

hc_centroid <- hclust(dist_matrix, method = "centroid")
plot(hc_centroid, main = "Centroide", 
     xlab = "Comunidades Autónomas", ylab = "Distancia")
rect.hclust(hc_centroid, k = 2, border = "red")

par(mfrow = c(1, 1))

library(NbClust)
set.seed(123)
nbclust_result <- NbClust(variables, 
                         distance = "euclidean",
                         min.nc = 2, 
                         max.nc = 8, 
                         method = "ward.D2",
                         index = "all")

best_nc_table <- data.frame(
  Criterio = names(nbclust_result$Best.nc[1,]),
  Número_Clusters = as.numeric(nbclust_result$Best.nc[1,])
)

kable(best_nc_table, caption = "Resultados de NbClust - Número óptimo de clusters por criterio")

freq_table <- table(nbclust_result$Best.nc[1,])
freq_df <- data.frame(
  Número_Clusters = names(freq_table),
  Frecuencia = as.numeric(freq_table)
)

kable(freq_df, caption = "Frecuencia de clusters recomendados (Regla de la mayoría)")

set.seed(123)
kmeans_result <- kmeans(variables, centers = 2, nstart = 25)

kmeans_summary <- data.frame(
  Métrica = c("Número de clusters", "Total within-cluster sum of squares", 
              "Between-cluster sum of squares", "Total sum of squares"),
  Valor = c(2, 
            round(kmeans_result$tot.withinss, 2),
            round(kmeans_result$betweenss, 2),
            round(kmeans_result$totss, 2))
)

kable(kmeans_summary, caption = "Resumen de resultados de K-means")

df_ccaa$cluster_kmeans <- kmeans_result$cluster

cluster1_comunidades <- df_ccaa$CCAA[df_ccaa$cluster_kmeans == 1]
cluster2_comunidades <- df_ccaa$CCAA[df_ccaa$cluster_kmeans == 2]

max_length <- max(length(cluster1_comunidades), length(cluster2_comunidades))
cluster1_comunidades <- c(cluster1_comunidades, rep("", max_length - length(cluster1_comunidades)))
cluster2_comunidades <- c(cluster2_comunidades, rep("", max_length - length(cluster2_comunidades)))

cluster_assignments <- data.frame(
  Cluster_1 = cluster1_comunidades,
  Cluster_2 = cluster2_comunidades
)

kable(cluster_assignments, caption = "Asignación de comunidades por cluster (K-means)")

cluster_hierarchical <- cutree(hc_ward, k = 2)
df_ccaa$cluster_hierarchical <- cluster_hierarchical

comparison_table <- data.frame(
  CCAA = df_ccaa$CCAA,
  Kmeans = df_ccaa$cluster_kmeans,
  Jerarquico = df_ccaa$cluster_hierarchical,
  Coincide = df_ccaa$cluster_kmeans == df_ccaa$cluster_hierarchical
)

kable(comparison_table, caption = "Comparación entre métodos K-means y Jerárquico")

coincidence_rate <- mean(comparison_table$Coincide) * 100
coincidence_df <- data.frame(
  Métrica = "Porcentaje de coincidencia entre métodos",
  Valor = paste(round(coincidence_rate, 2), "%")
)

kable(coincidence_df, caption = "Coincidencia entre métodos de clustering")

library(dplyr)

cluster_means <- df_ccaa %>%
  group_by(cluster_kmeans) %>%
  summarise(
    n = n(),
    across(automovil:telefono, mean, .names = "mean_{.col}")
  )

kable(round(cluster_means, 2), caption = "Medias de variables por cluster")

t_test_results <- data.frame(
  Variable = character(),
  Estadístico_t = numeric(),
  Valor_p = numeric(),
  Significancia = character()
)

variables_test <- names(df_ccaa)[2:7]

for(var in variables_test) {
  test_result <- t.test(df_ccaa[[var]] ~ df_ccaa$cluster_kmeans)
  significancia <- ifelse(test_result$p.value < 0.01, "***",
                         ifelse(test_result$p.value < 0.05, "**",
                               ifelse(test_result$p.value < 0.1, "*", "")))
  
  t_test_results <- rbind(t_test_results, data.frame(
    Variable = var,
    Estadístico_t = round(test_result$statistic, 3),
    Valor_p = round(test_result$p.value, 4),
    Significancia = significancia
  ))
}

kable(t_test_results, caption = "Pruebas t para diferencias entre clusters")

library(ggplot2)

pca_result <- prcomp(variables, scale. = TRUE)
pca_df <- as.data.frame(pca_result$x[, 1:2])
pca_df$CCAA <- df_ccaa$CCAA
pca_df$Cluster <- as.factor(df_ccaa$cluster_kmeans)

ggplot(pca_df, aes(x = PC1, y = PC2, color = Cluster, label = CCAA)) +
  geom_point(size = 3) +
  geom_text(size = 3, vjust = -0.5, hjust = 0.5) +
  scale_color_manual(values = c("1" = "blue", "2" = "red")) +
  labs(title = "Visualización de Conglomerados en Espacio PCA",
       x = paste("Componente Principal 1 (", 
                 round(summary(pca_result)$importance[2,1]*100, 1), "%)", sep = ""),
       y = paste("Componente Principal 2 (", 
                 round(summary(pca_result)$importance[2,2]*100, 1), "%)", sep = "")) +
  theme_minimal() +
  theme(legend.position = "bottom")

interpretacion_clusters <- data.frame()

for(i in 1:2) {
  cluster_data <- df_ccaa[df_ccaa$cluster_kmeans == i, ]
  comunidades <- paste(cluster_data$CCAA, collapse = ", ")
  
  medias <- colMeans(cluster_data[, 2:7])
  
  for(j in 1:length(medias)) {
    interpretacion_clusters <- rbind(interpretacion_clusters, 
      data.frame(
        Cluster = i,
        Variable = names(medias)[j],
        Media = round(medias[j], 1),
        Comunidades = ifelse(j == 1, comunidades, "")
      ))
  }
}

kable(interpretacion_clusters, caption = "Caracterización detallada de los clusters")

variables_significativas <- t_test_results %>%
  filter(Valor_p < 0.05) %>%
  select(Variable, Valor_p)

if(nrow(variables_significativas) > 0) {
  kable(variables_significativas, caption = "Variables con diferencias estadísticamente significativas (p < 0.05)")
} else {
  no_significativas <- data.frame(
    Mensaje = "No hay variables con diferencias significativas al 5%"
  )
  kable(no_significativas, caption = "Resultado de pruebas de significancia")
}

resumen_ejecutivo <- data.frame(
  Cluster = c(1, 2),
  N_Comunidades = c(sum(df_ccaa$cluster_kmeans == 1), 
                    sum(df_ccaa$cluster_kmeans == 2)),
  Característica_Principal = c("Menor nivel de equipamiento", 
                               "Mayor nivel de equipamiento"),
  Variables_Significativas = c(
    paste(variables_significativas$Variable, collapse = ", "),
    paste(variables_significativas$Variable, collapse = ", ")
  )
)

kable(resumen_ejecutivo, caption = "Resumen ejecutivo de los clusters identificados")