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))

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")

datos_ccaa$cluster_kmeans <- kmeans_result$cluster

cluster1_comunidades <- datos_ccaa$CCAA[datos_ccaa$cluster_kmeans == 1]
cluster2_comunidades <- datos_ccaa$CCAA[datos_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)")

library(dplyr)

cluster_means <- datos_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(datos_ccaa)[2:7]

for(var in variables_test) {
  test_result <- t.test(datos_ccaa[[var]] ~ datos_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 <- datos_ccaa$CCAA
pca_df$Cluster <- as.factor(datos_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")