6. Transversal_analysis.Rmd

```{r include=FALSE}
# Importation des données
mca_data <- load_parquet("mca_data.parquet")
detailed_mca_data <- load_parquet("detailed_mca_data.parquet")
map_data <- load_shapefile("FRANCE.shx")

codes_regions <- fread("data_other/codes_regions_Insee.csv", header = TRUE)
codes_regions$ABRV <- sapply(codes_regions$LIBELLE, get_first_letters)
```

```{r}
# Définition des variables
## Outcomes
mca_relative_outcomes <- c("part_T1", "GD_ratio_T1", "EC_ratio_T1")
mca_relative_outcomes_2T <- c("part_T1", "part_T2", "GD_ratio_T1", "EC_ratio_T1", "EC_ratio_T2")
mca_absolute_outcomes <- c("GIns_ratio_T1", "DIns_ratio_T1", "EIns_ratio_T1", "CIns_ratio_T1")
mca_absolute_outcomes_2T <- c("GIns_ratio_T1", "DIns_ratio_T1", "EIns_ratio_T1", "EIns_ratio_T2", "CIns_ratio_T1", "CIns_ratio_T2")
mca_detailed_outcomes <- grep("^pvoix", names(detailed_mca_data), value = TRUE)
mca_detailed_outcomes <- c(mca_detailed_outcomes, "part_T1", "part_T2")

## Régresseurs
### Abandonnés pour complaire aux économètres du groupe : prefract1791, pclerge1791, pclerge1856, pmessalisants1950
## Modalités supprimées pour éviter la colinéarité : prop014, propnodip, pouvr
mca_regressors <- c("ppropri", "popcomm", "popagglo", "age", "prop1539", "prop4059", "prop60p", "propf", "propbac", "propsup", "pagri", "pindp", "pempl", "ppint", "pcadr", "pchom", "petranger", "revmoy", "pcrimesdelits", "mmoyfortune", "pisf", "prsa", "type_comm")
full_mca_regressors <- c("ppropri", "popcomm", "popagglo", "age", "prop014", "prop1539", "prop4059", "prop60p", "propf", "propnodip", "propbac", "propsup", "pagri", "pindp", "pouvr", "pempl", "ppint", "pcadr", "pchom", "prsa", "petranger", "pcrimesdelits", "revmoy", "pisf", "mmoyfortune", "type_comm")
```

# Analyses descriptives

## Tableau récapitulatif

```{r}
nb_communes <- uniqueN(mca_data, by = "codecommune")

outcomes <- c("part_T1", "part_T2", "GD_ratio_T1", "EC_ratio_T1", "EC_ratio_T2","DIns_ratio_T1", "GIns_ratio_T1", "EIns_ratio_T1", "CIns_ratio_T1", "CIns_ratio_T2", "EIns_ratio_T2")

stat_desc <- data.table()
for (col in c(outcomes)) {
  summary <- mca_data[, {
    fivenum_values <- fivenum(get(col))
    nan_count <- sum(is.na(get(col)))
    nan_count_prop <- round(100 * nan_count / nb_communes, 2)
    nan_pop <- sum(popcomm[is.na(get(col))], na.rm = TRUE)
    data.table(
      Variable = col,
      "Min" = round(fivenum_values[1], 0),
      "1er quart." = fivenum_values[2],
      "Médiane" = fivenum_values[3],
      "3e quart." = fivenum_values[4],
      "Max" = fivenum_values[5],
      "Σ" = nan_count,
      "%" = nan_count_prop
    )
  }]
  stat_desc <- rbind(stat_desc, summary)
}

stat_desc_tbl <- stat_desc %>%
  gt() %>%
  fmt_number(decimals = 2, drop_trailing_zeros = TRUE, suffixing = TRUE) %>%
  tab_header(title = "Statistiques descriptives - 2017", subtitle = NULL, preheader = NULL) %>%
  tab_spanner(label = "Valeurs", columns = c("Min", "1er quart.", "Médiane", "3e quart.", "Max"), id = "values") %>%
  tab_spanner(label = "NaN", columns = c("Σ", "%"), id = "nan") %>%
  cols_label(Variable = "") %>%
  tab_row_group(label = "Approche séquentielle de la participation", rows = 1:5, id = "rel") %>%
  tab_row_group(label = "Approche simultanée de la participation", rows = 6:11, id = "abs") %>%
  row_group_order(groups = c("rel", "abs")) %>%
  tab_style(
    style = list(
      cell_fill(color = "lightgrey"),
      cell_text(style = "italic", align = "center")),
    locations = cells_row_groups(groups = c("rel", "abs"))) %>%
  tab_style(
    style = list(
      cell_borders(sides = "right", color = "lightgrey", style = "solid", weight = px(1))),
    locations = list(
      cells_body(column = c("Variable", "Max")))) %>%
  tab_style(
    style = cell_text(align = "center"),
    locations = cells_body(columns = !matches("Variable")))

stat_desc_tbl
gtsave(stat_desc_tbl, filename = "stat_desc_outcomes.tex", path = "output_transversal")

rm(col, nb_communes, stat_desc, stat_desc_tbl, summary, outcomes)
```

## Distribution du taux de propriétaires

```{r}
mean_ppropri <- mean(mca_data$ppropri, na.rm = TRUE)

plot <- ggplot(mca_data[is.finite(mca_data$ppropri), ], aes(x = ppropri)) +
  geom_histogram(fill = "lightblue", color = "skyblue", bins = 125) +
  geom_vline(xintercept = mean_ppropri, color = "darkgrey", linetype = "dashed", linewidth = 0.7) +
  labs(
    x = "Taux de propriétaires",
    y = "Nombre de communes") +
  theme_grey() +
  theme(plot.title = element_text(hjust = 0.5))

ggsave("output_transversal/distribution_ppropri_2017.png", plot, width = 8, height = 6, dpi = 300)

plot <- plot + labs(title = "Distribution du taux de propriétaires en 2017")
print(plot)

rm(mean_ppropri, plot)
```

## Cartes du taux de propriétaires

```{r}
draw_map <- function(region, mapped_var, mapped_var_name, base_width) {
  # Chargement des données
  carto <- merge(map_data, mca_data, by.x = "INSEE_COM", by.y = "codecommune", all.x = TRUE)
  if (region != 0) {
    carto <- subset(carto, INSEE_REG == region)
    reg_name <- codes_regions$LIBELLE[codes_regions$REG == region]
    reg_abrv <- codes_regions$ABRV[codes_regions$REG == region]
  } else {
    reg_name <- "France"
    reg_abrv <- "FR"
  }
  
  # Borne inférieure pour éviter que la coloration ne soit déformée par les outliers
  lbound <- 25
  carto[[mapped_var]] <- pmax(carto[[mapped_var]], lbound)

  # Création des chemins
  title <- paste0(mapped_var_name, "\n", reg_name, ", 2017")
  output_file <- paste0("output_transversal/map_ppropri/", mapped_var, "_2017_", reg_abrv, ".png")
  
  # Exportation d'une carte en haute résolution
  ## Calcul des dimensions
  bbox <- st_bbox(carto)
  aspect_ratio <- diff(c(bbox$xmin, bbox$xmax)) / diff(c(bbox$ymin, bbox$ymax))
  base_height <- base_width / aspect_ratio
  font_size <- base_width * 2.2
  ## Création et enregistrement
  plot_png <- ggplot() +
    geom_sf(data = carto, aes(fill = .data[[mapped_var]])) +
    scale_fill_viridis_c(
      name = mapped_var_name, 
      labels = scales::percent_format(scale = 1), 
      limits = c(lbound, 100)) +
    # ggtitle(paste(title)) +
    theme_minimal() +
    theme(
      # plot.title = element_text(hjust = 0.5, size = font_size),
      axis.text = element_blank(),
      panel.grid = element_blank(),
      plot.title.position = "plot",
      legend.title = element_blank(),
      # legend.title = element_text(hjust = 0.5, size = font_size * 0.8),
      legend.text = element_text(size = font_size * 0.7))
  ggsave(output_file, plot = plot_png, width = base_width, height = base_height, units = "in")
  
  # Affichage d'une carte en résolution normale
  display_plot <- ggplot() +
    geom_sf(data = carto, aes(fill = .data[[mapped_var]])) +
    scale_fill_viridis_c(
      name = mapped_var_name,
      labels = scales::percent_format(scale = 1),
      limits = c(lbound,100)) +
    ggtitle(paste(title)) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text = element_blank(),
      panel.grid = element_blank(),
      plot.title.position = "plot",
      legend.title = element_blank(),
      legend.text = element_text())   
    print(display_plot)
}

for (i in c(11, 24, 27, 28, 32, 44, 52, 53, 75, 76, 84, 93)) {draw_map(i, "ppropri", "Taux de propriétaires", 12)}

rm(i, draw_map)
```

## Corrélogramme

```{r}
corr_variables <- c(mca_regressors, mca_relative_outcomes)
mca_corr_data <- mca_data[, ..corr_variables, with = FALSE]
mca_corr_data <- mca_corr_data[, .SD, .SDcols = sapply(mca_corr_data, is.numeric)]
mca_corr_data <- mca_corr_data[complete.cases(mca_corr_data)]
correlation_matrix <- cor(mca_corr_data)
p_values = cor.mtest(mca_corr_data, conf.level = 0.95)
corrplot(correlation_matrix, p.mat = p_values$p, pch.cex = 0.7, type = "upper", method = "color", COL2("RdYlBu"), order = "FPC", tl.srt = 60, tl.col = "darkgrey", diag = FALSE)

# Enregistrement
png(filename = "output_transversal/corrplot.png", width = 2200, height = 1800)
corrplot(correlation_matrix, type = "upper", method = "color", order = "FPC", tl.cex = 3.5, tl.srt = 60, tl.col = "darkgrey", cl.cex = 3.5, diag = FALSE)
dev.off()

rm(corr_variables, correlation_matrix, mca_corr_data, p_values)
```

## Corrélations de ppropri

```{r}
mca_data_complete <- na.omit(mca_data[, c("ppropri", full_mca_regressors), with = FALSE])
mca_data_complete$type_comm <- as.factor(mca_data_complete$type_comm)
dummies <- model.matrix(~ type_comm - 1, data = mca_data_complete)
colnames(dummies) <- paste0("type_comm_", c("village", "bourgs", "banlieues", "centre_metrop"))
mca_data_complete <- cbind(mca_data_complete, dummies)
mca_data_complete$type_comm <- NULL

correlation_results <- sapply(mca_data_complete[, -1, drop = FALSE], function(x) cor(mca_data_complete$ppropri, x))
correlation_results <- correlation_results * 100
correlation_df <- data.frame(Variables = names(correlation_results), Corrélation = correlation_results)
correlation_df <- subset(correlation_df, Variables != "ppropri")
correlation_df <- correlation_df[order(correlation_df$Corrélation, decreasing = TRUE), ]
rownames(correlation_df) <- NULL

correlation_df <- correlation_df %>%
  gt() %>%
  fmt_number(decimals = 2, drop_trailing_zeros = TRUE) %>%
  tab_header("Corrélation des régresseurs avec le taux de propriétaires", subtitle = NULL, preheader = NULL) %>%
  data_color(method = "numeric", columns = !matches("Variables"), palette ="Spectral", alpha = 0.6) %>%
  tab_style(
    style = cell_text(align = "center"),
    locations = list(
      cells_body(column = c("Corrélation")),
      cells_column_labels()))

correlation_df
gtsave(correlation_df, filename = "correl_ppropri_2017.tex", path = "output_transversal")

rm(mca_data_complete, correlation_results, correlation_df, dummies)
```

# Régressions linéaires simples

```{r}
mca_data$type_comm <- factor(mca_data$type_comm)
relative_results <- list()
for (outcome in mca_relative_outcomes_2T) {
  regression_formula <- as.formula(paste(outcome, "~", paste(mca_regressors, collapse = "+")))
  model_name <- paste0(outcome, "_model")
  relative_results[[model_name]] <- lm(regression_formula, data = mca_data)
}
suppressWarnings(stargazer(relative_results, type = "text"))
```

```{r}
mca_data$type_comm <- factor(mca_data$type_comm)
absolute_results <- list()
for (outcome in mca_absolute_outcomes_2T) {
  regression_formula <- as.formula(paste(outcome, "~", paste(mca_regressors, collapse = "+")))
  model_name <- paste0(outcome, "_model")
  absolute_results[[model_name]] <- lm(regression_formula, data = mca_data)
}
suppressWarnings(stargazer(absolute_results, type = "text"))
```

```{r include=FALSE}
suppressWarnings(stargazer(relative_results, type = "latex", out = "output_transversal/2017_relative_results.tex"))
suppressWarnings(stargazer(relative_results, type = "latex", out = "output_transversal/2017_absolute_results.tex"))
rm(model_name, outcome, regression_formula, relative_results,absolute_results)
```

# ACP

```{r}
# Fonction pour tracer les corrplots
draw_corrplot <- function(data, title, subtitle, filename, range) {
  corrplot <- ggcorrplot(
    data,
    method = "square",
    ggtheme = theme_grey,
    show.legend = TRUE,
    outline.color = "darkgrey",
    tl.cex = 8) +
    
    if (identical(range, c(-1,1))) {
      scale_fill_gradient2(limit = c(-1,1), low = "#6D9EC1", mid = "white", high =  "#E46726", midpoint = 0)}
    else if (identical(range, c(0,1))) {
      scale_fill_gradient2(limit = c(0, 1), low = "white", high =  "#E46726")}

  corrplot$theme$plot.title$hjust <- 0.5
  corrplot$theme$plot.title$size <- 12
  corrplot$theme$plot.subtitle$hjust <- 0.5
  corrplot$theme$plot.subtitle$size <- 12
  corrplot$theme$legend.key.width <- unit(0.3, "cm")
  corrplot$theme$legend.key.height <- unit(1, "null")
  corrplot$theme$legend.title <- element_blank()

  ggsave(filename, plot = corrplot)
  
  corrplot$labels$title <- title
  corrplot$labels$subtitle <- subtitle
  
  print(corrplot)
}
```

## Sur les comportements électoraux

On la réalise sur les variables électorales, pour obtenir un espace des comportements politiques sur lequel on projette ensuite les variables socio-démographiques.

Concernant l'analyse des axes, il faut distinguer : - les coordonnées factorielles des variables = corrélation avec l'axe (tableau à lire verticalement pour comprendre la composition des axes) - les cos2 des variables = coordonnées au carré = qualité de de représentation (tableau à lire horizontalement pour comprendre à quel point chaque variable participe à l'espace finalement retenu) - les contributions à l'inertie des axes = 100 \* cos2 / valeur propre de l'axe

```{r}
# Préparation : sélection des variables, suppression des valeurs manquantes et création d'indicatrices pour les types de communes (car l'argument quali.sup donne des résultats étranges)
variables_of_interest <- c(mca_detailed_outcomes, full_mca_regressors)
mca_data_subset <- detailed_mca_data[, ..variables_of_interest]

mca_data_subset <- na.omit(mca_data_subset)

mca_data_subset$type_comm <- as.factor(mca_data_subset$type_comm)
dummies <- model.matrix(~ type_comm - 1, data = mca_data_subset)
colnames(dummies) <- paste0("type_comm_", 1:ncol(dummies))
mca_data_subset <- cbind(mca_data_subset, dummies)
mca_data_subset$type_comm <- NULL

# Analyse
quant_regressors <- setdiff(full_mca_regressors, "type_comm")
quant_regressors <- c(quant_regressors, "type_comm_1", "type_comm_2", "type_comm_3", "type_comm_4")
result <- PCA(mca_data_subset, scale.unit = TRUE, quanti.sup=quant_regressors, ncp = 10, graph = FALSE)

# Résultats
## Histogramme de la contribution des axes à la variance totale
screeplot <- fviz_screeplot(result, addlabels = TRUE, ylim = c(0, 30), title = "") + theme_gray()
screeplot$layers[[1]]$aes_params$fill <- "lightskyblue"
screeplot$layers[[1]]$aes_params$colour <- "cornflowerblue"
screeplot$layers[[2]]$aes_params$colour <- "darkblue"
screeplot$layers[[3]]$aes_params$colour <- "darkblue"
screeplot$layers[[4]]$aes_params$colour <- "darkblue"
screeplot$labels$x <- "Dimensions"
screeplot$labels$y <- "% de variance expliquée"
ggsave("output_transversal/screeplot_ACP1.png", plot = screeplot, width = 10, height = 6, units = "in", dpi = 300)

screeplot$labels$title = "Histogramme de l'inertie totale"
screeplot$theme$plot.title$hjust <- 0.5
screeplot

## Contribution des variables aux axes (comportements électoraux seulement)
## Correspond à la qualité de représentation des variables (= cos2) normalisée entre 1 et 100
## Présentable à l'aide d'un tableau...
coeffs <- result$var$contrib
labels <- rownames(result$var$contrib)
axis_df <- data.frame(labels, coeffs, row.names = NULL)
colnames(axis_df) <- c("Variables", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10")
axis_tbl <- axis_df %>%
  gt() %>%
  fmt_number(decimals = 2, drop_trailing_zeros = TRUE) %>%
  tab_header("Contribution des variables à l'inertie des axes", subtitle = NULL, preheader = NULL) %>%
  tab_spanner(label = "Axes", columns = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"), id = "axes") %>%
  data_color(method = "numeric", columns = !matches("Variables"), palette ="Oranges", ) %>%
  tab_style(
    style = cell_text(align = "center"),
    locations = list(
      cells_column_labels(columns = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10")),
      cells_body(columns = !matches("Variables")))) %>%
  tab_style(
    style = cell_fill(color = "lightgrey"),
    locations = cells_column_labels(columns = c("1", "3", "5", "7", "9"))) 
axis_tbl
gtsave(axis_tbl, filename = "contrib_axes_acp1.tex", path = "output_transversal")
## ... ou d'un corrplot
draw_corrplot(t(result$var$contrib[, 1:10]/100), "Contribution des variables à l'inertie des axes", "", "output_transversal/contrib_axes_ACP1.png", range = c(0,1))

## Corrélation des variables aux axes
## Avec à nouveau deux représentations possibles...
## Des corrplots
draw_corrplot(t(result$var$coord[, 1:3]), title = "Corrélation des outcomes aux axes", "", "output_transversal/corr_outc_axes_ACP1.png", range = c(-1,1))
draw_corrplot(t(result$quanti.sup$coord[, 1:3]), "Corrélation des régresseurs", "aux axes", "output_transversal/corr_rég_axes_ACP1.png", range = c(-1,1))
## Des projections sur les plans factoriels
plot.PCA(result, axes = c(1,2), choix = "var", select = "cos2 7", title = "Projection sur les axes 1 et 2")
plot.PCA(result, axes=c(1,3), choix="var", select = "cos2 7", title = "Projection sur les axes 1 et 3")
plot.PCA(result, axes=c(2,3), choix="var", select = "cos2 7", title = "Projection sur les axes 2 et 3")

## Qualité de représentation des variables
draw_corrplot(t(result$var$cos2[, 1:3]), "Qualité de représentation", "des variables électorales", "output_transversal/qual_rep_outc_ACP1.png", range = c(0,1))
draw_corrplot(t(result$quanti.sup$cos2[, 1:3]), "Qualité de représentation", "des régresseurs", "output_transversal/qual_rep_régr_ACP1.png", range = c(0,1))

# Nettoyage
rm(axis_df, axis_tbl, coeffs, dummies, labels, mca_data_subset, quant_regressors, result, variables_of_interest, screeplot)
```

## Sur l'ensemble des variables

On refait une ACP sans distinguer les outcomes et les régresseurs. Les corrplots deviennent illisibles, mais comme on ne les utilisera pas...

```{r}
# Préparation : sélection des variables, suppression des valeurs manquantes et création d'indicatrices pour les types de communes (car l'argument quali.sup donne des résultats étranges)
variables_of_interest <- c(mca_detailed_outcomes, full_mca_regressors)
mca_data_subset <- detailed_mca_data[, ..variables_of_interest]

mca_data_subset <- na.omit(mca_data_subset)

mca_data_subset$type_comm <- as.factor(mca_data_subset$type_comm)
dummies <- model.matrix(~ type_comm - 1, data = mca_data_subset)
colnames(dummies) <- paste0("type_comm_", 1:ncol(dummies))
mca_data_subset <- cbind(mca_data_subset, dummies)
mca_data_subset$type_comm <- NULL

# Analyse
result <- PCA(mca_data_subset, scale.unit = TRUE, ncp = 10, graph = FALSE)

# Résultats
## Histogramme de la contribution des axes à la variance totale
screeplot <- fviz_screeplot(result, addlabels = TRUE, ylim = c(0, 30), title = "") + theme_gray()
screeplot$layers[[1]]$aes_params$fill <- "lightskyblue"
screeplot$layers[[1]]$aes_params$colour <- "cornflowerblue"
screeplot$layers[[2]]$aes_params$colour <- "darkblue"
screeplot$layers[[3]]$aes_params$colour <- "darkblue"
screeplot$layers[[4]]$aes_params$colour <- "darkblue"
screeplot$theme$plot.title$hjust <- 0.5
screeplot$labels$x <- "Dimensions"
screeplot$labels$y <- "% de variance expliquée"
ggsave("output_transversal/screeplot_ACP2.png", plot = screeplot, width = 10, height = 6, units = "in", dpi = 300)

screeplot$labels$title = "Histogramme de l'inertie totale"
screeplot$theme$plot.title$hjust <- 0.5
screeplot

## Contribution des variables aux axes (comportements électoraux seulement)
## Correspond à la qualité de représentation des variables (= cos2) normalisée entre 1 et 100
## Présentable à l'aide d'un tableau...
coeffs <- result$var$contrib
labels <- rownames(result$var$contrib)
axis_df <- data.frame(labels, coeffs, row.names = NULL)
colnames(axis_df) <- c("Variables", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10")
axis_tbl <- axis_df %>%
  gt() %>%
  fmt_number(decimals = 2, drop_trailing_zeros = TRUE) %>%
  tab_header("Contribution des variables à l'inertie des axes", subtitle = NULL, preheader = NULL) %>%
  tab_spanner(label = "Axes", columns = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"), id = "axes") %>%
  data_color(method = "numeric", columns = !matches("Variables"), palette ="Oranges", ) %>%
  tab_row_group(label = "Outcomes", rows = 1:15, id = "out") %>%
  tab_row_group(label = "Régresseurs", rows = 16:44, id = "reg") %>%
  row_group_order(groups = c("out", "reg")) %>%
  tab_style(
    style = list(
      cell_fill(color = "lightgrey"),
      cell_text(style = "italic", align = "center")),
    locations = cells_row_groups(groups = c("out", "reg"))) %>%  
  tab_style(
    style = cell_text(align = "center"),
    locations = list(
      cells_column_labels(columns = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10")),
      cells_body(columns = !matches("Variables")))) %>%
  tab_style(
    style = cell_fill(color = "lightgrey"),
    locations = cells_column_labels(columns = c("1", "3", "5", "7", "9"))) 
axis_tbl
gtsave(axis_tbl, filename = "contrib_axes_acp2.tex", path = "output_transversal")
## ... ou d'un corrplot
draw_corrplot(t(result$var$contrib[, 1:10]/100), "Contribution des variables à l'inertie des axes", "", "output_transversal/contrib_axes_ACP2.png", range = c(0,1))

## Corrélation des variables aux axes
## Avec à nouveau deux représentations possibles...
## Des corrplots
draw_corrplot(t(result$var$coord[1:15, 1:3]), title = "Corrélation des outcomes aux axes", "", "output_transversal/corr_outc_axes_ACP2.png", range = c(-1,1))
draw_corrplot(t(result$var$coord[16:44, 1:3]), "Corrélation des régresseurs", "aux axes", "output_transversal/corr_rég_axes_ACP2.png", range = c(-1,1))
## Des projections sur les plans factoriels
plot.PCA(result, axes = c(1,2), choix = "var", select = "cos2 13", title = "Projection sur les axes 1 et 2")
plot.PCA(result, axes=c(1,3), choix="var", select = "cos2 13", title = "Projection sur les axes 1 et 3")
plot.PCA(result, axes=c(2,3), choix="var", select = "cos2 13", title = "Projection sur les axes 2 et 3")

## Qualité de représentation des variables
draw_corrplot(t(result$var$cos2[1:15, 1:3]), "Qualité de représentation", "des variables électorales", "output_transversal/qual_rep_outc.png", range = c(0,1))
draw_corrplot(t(result$var$cos2[16:44, 1:3]), "Qualité de représentation", "des régresseurs", "output_transversal/qual_rep_régr.png", range = c(0,1))

# Nettoyage
rm(axis_df, axis_tbl, coeffs, dummies, labels, mca_data_subset, result, variables_of_interest, screeplot)
```

# K-means

On inverse l'approche : on construit un espace socio-démographique, sur lequel on projette les comportements politiques

## Détermination du nombre de clusters optimal

```{r}
set.seed(789)

## Préparation des données
## On élimine popagglo pour éviter qu'elle ne déforme trop les résultats en région parisienne
regressors <- setdiff(full_mca_regressors, "popagglo")
mca_data_subset <- detailed_mca_data[, ..regressors]
mca_data_subset <- na.omit(mca_data_subset)
mca_data_subset <- scale(mca_data_subset)

## Méthode elbow
wcss <- numeric(length = 20)
for (i in 1:20) {
  kmeans_model <- kmeans(mca_data_subset, centers = i)
  wcss[i] <- kmeans_model$tot.withinss
}
plot(1:20, wcss, type = "b", xlab = "Nombre de clusters", ylab = "Somme des distances intra-clusters au carré")

## Méthode silhouette
silhouette_avg <- numeric(length = 20)
for (i in 2:20) {
  kmeans_model <- kmeans(mca_data_subset, centers = i)
  sil <- silhouette(kmeans_model$cluster, dist(mca_data_subset))
  silhouette_avg[i - 1] <- mean(sil[, 3])
}
plot(2:20, silhouette_avg[2:20], type = "b", xlab = "Nombre de clusters", ylab = "Coefficient silhouette")

## Enregistrement des graphs
png("output_transversal/WCSS.png", width = 1000, height = 500)
plot(1:20, wcss, type = "b", xlab = "Nombre de clusters", ylab = "Somme des distances intra-clusters au carré")
dev.off()
png("output_transversal/Silhouette.png", width = 1000, height = 500)
plot(2:20, silhouette_avg[2:20], type = "b", xlab = "Nombre de clusters", ylab = "Coefficient silhouette")
dev.off()

rm(regressors, mca_data_subset, kmeans_model, wcss, silhouette_avg, sil)
```

## Classification avec k=6

```{r}
set.seed(789)

# Préparation des données
## Sélection des variables
regressors <- c("codecommune", full_mca_regressors)
regressors <- setdiff(regressors, "popagglo")
mca_data_subset <- detailed_mca_data[, ..regressors]
mca_data_subset <- na.omit(mca_data_subset)
## Création d'indicatrices pour les types de communes
mca_data_subset$type_comm <- as.factor(mca_data_subset$type_comm)
dummies <- model.matrix(~ type_comm - 1, data = mca_data_subset)
colnames(dummies) <- paste0("type_comm_", 1:ncol(dummies))
mca_data_subset <- cbind(mca_data_subset, dummies)
mca_data_subset$type_comm <- NULL
## Normalisation
scaling_regressors <- setdiff(full_mca_regressors, c("popagglo", "type_comm"))
regressors <- c(scaling_regressors, paste0("type_comm_", 1:4))
mca_data_subset_scaled <- mca_data_subset[, ..regressors]
mca_data_subset_scaled[, (scaling_regressors) := lapply(.SD, scale), .SDcols = scaling_regressors]

# Détermination des clusters et ajout des variables électorales
kmeans_result <- kmeans(mca_data_subset_scaled, centers = 6)
mca_data_subset$cluster <- kmeans_result$cluster
columns <- c("codecommune", "cluster", paste0("type_comm_", 1:4))
cluster_data <- mca_data_subset[, ..columns, drop = FALSE]
cluster_data <- merge(cluster_data, detailed_mca_data, by = "codecommune")
```

```{r}
# Tableau avec les valeurs moyennes par commune
## Calcul des moyennes par cluster
variables_of_interest <- c(mca_detailed_outcomes, full_mca_regressors, paste0("type_comm_", 1:4))
variables_of_interest <- setdiff(variables_of_interest, "type_comm")
means <- cluster_data[, lapply(.SD, mean, na.rm = TRUE), by = cluster, .SDcols = variables_of_interest]
## Pour calculer plutôt les moyennes pondérées par cluster
# means <- cluster_data[, lapply(.SD, function(x) weighted.mean(x, w = popcomm, na.rm = TRUE)), by = cluster, .SDcols = variables_of_interest]
## Calcul des effectifs des clusters
sum <- cluster_data[, .(Population = sum(popcomm, na.rm = TRUE)), by = cluster]
means <- merge(means, sum, by = "cluster")
means$cluster <- as.factor(means$cluster)
cluster_counts <- as.data.frame(table(cluster_data$cluster))
means <- merge(means, cluster_counts, by.x = "cluster", by.y = "Var1")
## Calcul des moyennes pour la population
global_means <- cluster_data[, lapply(.SD, mean, na.rm = TRUE), .SDcols = variables_of_interest]
global_means$cluster <- "mean"
means <- rbind(means, global_means, fill = TRUE)
## Expression des types de commune en %
means[, paste0("type_comm_", 1:4)] <- lapply(means[, paste0("type_comm_", 1:4)], as.numeric)
means[, paste0("type_comm_", 1:4)] <- means[, paste0("type_comm_", 1:4)] * 100
means$cluster <- as.character(means$cluster)
## Transposition pour afficher les clusters en colonne
t_means <- t(means)
t_means <- data.frame(matrix = rownames(t_means), t_means, row.names = NULL)
colnames(t_means) <- t_means[1, ]
t_means <- t_means[-1, ]
t_means[setdiff(names(t_means), "cluster")] <- lapply(t_means[setdiff(names(t_means), "cluster")], as.numeric)
## Réarrangement des lignes
first_rows <- t_means[-((nrow(t_means) - 1):nrow(t_means)), ]
last_rows <- t_means[(nrow(t_means) - 1):nrow(t_means), ]
t_means <- rbind(last_rows, first_rows)
t_means$cluster[2] <- "Communes"
## Intervalle sur lequel appliquer la coloration (afin qu'elle soit centrée sur la moyenne de la population plutôt que des clusters)
t_means <- as.data.frame(t_means)
t_means$min <- apply(t_means[, 2:7], 1, min)
t_means$max <- apply(t_means[, 2:7], 1, max)
t_means$diff_min <- t_means$mean - t_means$min
t_means$diff_max <- t_means$max - t_means$mean
t_means$spread <- pmax(t_means$diff_min, t_means$diff_max)
t_means$range_inf <- t_means$mean - t_means$spread
t_means$range_sup <- t_means$mean + t_means$spread
t_means <- subset(t_means, select = -c(min, max, diff_min, diff_max, spread))
## Création du tableau avec gt
means_tbl <- t_means %>%
  gt() %>%
  fmt_number(rows = 3:46, decimals = 2, drop_trailing_zeros = TRUE) %>%
  fmt_number(rows = (c(1,19,20)), decimals = 0, scale_by = 1000) %>%
  fmt_number(rows = (c(2,40,42)), decimals = 0) %>%
  cols_hide(c("range_inf", "range_sup")) %>%
  tab_header(title = "Moyennes des variables par cluster", subtitle = "Exprimées selon communes") %>%
  data_color(
    method = "numeric",
    columns = !matches("cluster|mean"),
    rows = 3:46,
    direction = "row",
    palette = "RColorBrewer::RdBu",
    alpha = 0.6) %>%
  tab_spanner(label = "Clusters", columns = c("1", "2", "3", "4", "5", "6"), id = "clust") %>%
  tab_spanner(label = "Total", columns = "mean", id = "moy") %>%
  cols_label(cluster = "", mean = "") %>%
  tab_row_group(label = "Effectifs", rows = 1:2, id = "pop") %>%
  tab_row_group(label = "Outcomes", rows = 3:17, id = "out") %>%
  tab_row_group(label = "Régresseurs", rows = 18:46, id = "reg") %>%
  row_group_order(groups = c("pop", "out", "reg")) %>%
  sub_missing(missing_text = "-") %>%
  tab_style(
    style = list(
      cell_fill(color = "lightgrey"),
      cell_text(style = "italic", align = "center")),
    locations = cells_row_groups(groups = c("pop", "out", "reg"))) %>%
  tab_style(
    style = cell_text(style = "italic", align = "center"),
    locations = list(
      cells_column_labels(columns = "mean"),
      cells_body(columns = "mean"))) %>%
 tab_style(
    style = cell_text(align = "center"),
    locations = list(
      cells_column_labels(columns = c("1", "2", "3", "4", "5", "6", "mean")),
      cells_body(columns = !matches("cluster"))))
means_tbl
gtsave(means_tbl, filename = "k6_means.tex", path = "output_transversal")
```

```{r}
# Tableau avec les valeurs moyennes pondérées par la population
## Calcul des moyennes par cluster
variables_of_interest <- c(mca_detailed_outcomes, full_mca_regressors, paste0("type_comm_", 1:4))
variables_of_interest <- setdiff(variables_of_interest, "type_comm")
means <- cluster_data[, lapply(.SD, function(x) weighted.mean(x, w = popcomm, na.rm = TRUE)), by = cluster, .SDcols = variables_of_interest]
## Calcul des effectifs des clusters
sum <- cluster_data[, .(Population = sum(popcomm, na.rm = TRUE)), by = cluster]
means <- merge(means, sum, by = "cluster")
means$cluster <- as.factor(means$cluster)
cluster_counts <- as.data.frame(table(cluster_data$cluster))
means <- merge(means, cluster_counts, by.x = "cluster", by.y = "Var1")
## Calcul des moyennes pour la population
global_means <- cluster_data[, lapply(.SD, function(x) weighted.mean(x, w = popcomm, na.rm = TRUE)), .SDcols = variables_of_interest]
global_means$cluster <- "mean"
means <- rbind(means, global_means, fill = TRUE)
## Expression des types de commune en %
means[, paste0("type_comm_", 1:4)] <- lapply(means[, paste0("type_comm_", 1:4)], as.numeric)
means[, paste0("type_comm_", 1:4)] <- means[, paste0("type_comm_", 1:4)] * 100
means$cluster <- as.character(means$cluster)
## Transposition pour afficher les clusters en colonne
t_means <- t(means)
t_means <- data.frame(matrix = rownames(t_means), t_means, row.names = NULL)
colnames(t_means) <- t_means[1, ]
t_means <- t_means[-1, ]
t_means[setdiff(names(t_means), "cluster")] <- lapply(t_means[setdiff(names(t_means), "cluster")], as.numeric)
## Réarrangement des lignes
first_rows <- t_means[-((nrow(t_means) - 1):nrow(t_means)), ]
last_rows <- t_means[(nrow(t_means) - 1):nrow(t_means), ]
t_means <- rbind(last_rows, first_rows)
t_means$cluster[2] <- "Communes"
## Intervalle sur lequel appliquer la coloration (afin qu'elle soit centrée sur la moyenne de la population plutôt que des clusters)
t_means <- as.data.frame(t_means)
t_means$min <- apply(t_means[, 2:7], 1, min)
t_means$max <- apply(t_means[, 2:7], 1, max)
t_means$diff_min <- t_means$mean - t_means$min
t_means$diff_max <- t_means$max - t_means$mean
t_means$spread <- pmax(t_means$diff_min, t_means$diff_max)
t_means$range_inf <- t_means$mean - t_means$spread
t_means$range_sup <- t_means$mean + t_means$spread
t_means <- subset(t_means, select = -c(min, max, diff_min, diff_max, spread))
## Création du tableau avec gt
means_tbl <- t_means %>%
  gt() %>%
  fmt_number(rows = 3:46, decimals = 2, drop_trailing_zeros = TRUE) %>%
  fmt_number(rows = (c(1,19,20)), decimals = 0, scale_by = 1000) %>%
  fmt_number(rows = (c(2,40,42)), decimals = 0) %>%
  cols_hide(c("range_inf", "range_sup")) %>%
  tab_header(title = "Moyennes des variables par cluster", subtitle = "Exprimées selon la population") %>%
  data_color(
    method = "numeric",
    columns = !matches("cluster|mean"),
    rows = 3:46,
    direction = "row",
    palette = "RColorBrewer::RdBu",
    alpha = 0.6) %>%
  tab_spanner(label = "Clusters", columns = c("1", "2", "3", "4", "5", "6"), id = "clust") %>%
  tab_spanner(label = "Total", columns = "mean", id = "moy") %>%
  cols_label(cluster = "", mean = "") %>%
  tab_row_group(label = "Effectifs", rows = 1:2, id = "pop") %>%
  tab_row_group(label = "Outcomes", rows = 3:17, id = "out") %>%
  tab_row_group(label = "Régresseurs", rows = 18:46, id = "reg") %>%
  row_group_order(groups = c("pop", "out", "reg")) %>%
  sub_missing(missing_text = "-") %>%
  tab_style(
    style = list(
      cell_fill(color = "lightgrey"),
      cell_text(style = "italic", align = "center")),
    locations = cells_row_groups(groups = c("pop", "out", "reg"))) %>%
  tab_style(
    style = cell_text(style = "italic", align = "center"),
    locations = cells_body(columns = "mean")) %>%
 tab_style(
    style = cell_text(align = "center"),
    locations = list(
      cells_column_labels(columns = c("1", "2", "3", "4", "5", "6", "mean")),
      cells_body(columns = !matches("cluster"))))
means_tbl
gtsave(means_tbl, filename = "k6_weighted_means.tex", path = "output_transversal")
```

```{r}
# Nettoyage
rm(cluster_counts, columns, dummies, first_rows, global_means, kmeans_result, last_rows, mca_data_subset, mca_data_subset_scaled, means, means_tbl, regressors, scaling_regressors, sum, t_means, variables_of_interest)
```

Les résultats sont plus intéressants, avec des clusters qui se différencient nettement.

NB : les données sur le RSA sont manquantes pour la Corse, ainsi que les villes de Lyon et Marseille. Celles-ci sont donc exclues de la classification !

## Cartes

```{r}
draw_map_cluster <- function(region, base_width) {
  carto <- merge(map_data, cluster_data, by.x = "INSEE_COM", by.y = "codecommune", all.x = TRUE)
  carto[["cluster"]] <- factor(carto[["cluster"]])
  if (region != 0) {
    carto <- subset(carto, INSEE_REG == region)
    reg_name <- codes_regions$LIBELLE[codes_regions$REG == region]
    reg_abrv <- codes_regions$ABRV[codes_regions$REG == region]
  } else {
    reg_name <- "France"
    reg_abrv <- "FR"
  }

  output_file <- paste0("output_transversal/map_clusters/clusters_", reg_abrv, ".png")
  title <- paste0("Clusters en ", reg_name)
  legend_labels <- c("1" = "Centre des grandes villes",
                     "2" = "Banlieues et villes pauvres",
                     "5" = "Banlieues et zones périurbaines aisées",
                     "4" = "Périurbain populaire",
                     "3" = "Rural précaire",
                     "6" = "Rural à dominante paysanne")
  legend_colors <- c("1" = "blue3",
                     "2" = "cornflowerblue",
                     "5" = "lightskyblue1",
                     "6" = "darkgreen",
                     "3" = "limegreen",
                     "4" = "lightgreen")

  # Exportation d'une carte en haute résolution
  ## Calcul des dimensions
  bbox <- st_bbox(carto)
  aspect_ratio <- diff(c(bbox$xmin, bbox$xmax)) / diff(c(bbox$ymin, bbox$ymax))
  base_height <- base_width / aspect_ratio
  font_size <- base_width * 2.2
  ## Création et enregistrement
  plot_png <- ggplot() +
    # geom_sf(data = carto, aes(fill = .data[["cluster"]]), show.legend = FALSE) +
    geom_sf(data = carto, aes(fill = .data[["cluster"]])) +
    # ggtitle(paste(title)) +
    theme_minimal() +
    theme(
      # plot.title = element_text(hjust = 0.5, size = font_size),
      axis.text = element_blank(),
      panel.grid = element_blank(),
      plot.title.position = "plot",
      # legend.title = element_text(hjust = 0.5, size = font_size * 0.8),
      legend.text = element_text(size = font_size * 0.5)) +
    scale_fill_manual(
      values = legend_colors,
      breaks = names(legend_labels),
      labels = legend_labels,
      name = NULL)
  ggsave(output_file, plot = plot_png, width = base_width, height = base_height, units = "in")
  
  display_plot <- ggplot() +
    geom_sf(data = carto, aes(fill = .data[["cluster"]])) +
    ggtitle(paste(title)) +
    theme_minimal() +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text = element_blank(),
      panel.grid = element_blank(),
      plot.title.position = "plot",
      legend.title = element_text(hjust = 0.5),
      legend.text = element_text()) +
    scale_fill_manual(
      values = legend_colors,
      breaks = names(legend_labels),
      labels = legend_labels,
      name = NULL)
    print(display_plot)
}

for (i in c(11, 24, 27, 28, 32, 44, 52, 53, 75, 76, 84, 93)) {draw_map_cluster(i, 12)}

rm(cluster_data, i)
```