Annexe G : Code R

setwd("~/Cours_master/Mémoire")

library(tidyverse)
library(ggpubr)
library(readxl)
library(rstatix)
library(coin)
library(readr)
library(dplyr)
library(Hmisc)
library(ggplot2)
library(sf)
library(spdep)
library(spatialreg)
library(lmtest)
library(car)

df <- read_excel("ecoles.xlsx")
df$type_ecole <- as.factor(df$type_ecole)

# variables
indicators <- c("vegetation", "canopee", "purification", 
                "temp_regulation", "pm10", "heatrisk")

df_test <- df %>% filter(type_ecole %in% c("Public", "Libre"))
df_test$type_ecole <- factor(df_test$type_ecole, levels = c("Public", "Libre"))

# Boucle de calcul pour chaque indicateur
results_list <- list()
for (ind in indicators) {
  
  #Mann-Whitney U
  mw_test <- wilcox.test(as.formula(paste(ind, "~ type_ecole")), data = df_test)
  
  #Taille d'effet 
  eff_size <- df_test %>% wilcox_effsize(as.formula(paste(ind, "~ type_ecole")))
  
  #Assemblage de la ligne pour cet indicateur
  row <- data.frame(
    Indicateur = ind,
    Mann_Whitney_p = round(mw_test$p.value, 3),
    Effect_Size_r = round(eff_size$effsize, 3)
  )
  
  results_list[[ind]] <- row
}

# tableau
tab4 <- bind_rows(results_list)
print(tab4)

#exporter

#### statistiques descriptives indicateurs environnementaux ####
#inclut toutes les écoles

niveau_prive <- setdiff(unique(df$type_ecole), c("Public", "Libre"))
df$type_ecole   <- factor(df$type_ecole,   levels = c("Public", "Libre", niveau_prive))
df$LANGUE <- factor(df$LANGUE, levels = c("FR", "NL", "NA"))

# mise en forme
format_stat <- function(x, type = c("median_iqr", "mean_sd")) {
  type <- match.arg(type)
  x <- x[!is.na(x)]
  if (length(x) == 0) return(NA_character_)
  if (type == "median_iqr") sprintf("%.2f (%.2f)", median(x), IQR(x))
  else sprintf("%.2f (%.2f)", mean(x), sd(x))
}

#resumé par groupe
resumer_groupe <- function(data, group_var) {
  data %>%
    group_by(.data[[group_var]]) %>%
    summarise(
      n = n(),
      `Végétation intérieure`      = format_stat(vegetation, "median_iqr"),
      `Canopée de rue`             = format_stat(canopee, "median_iqr"),
      `EVP accessibles*`           = format_stat(EVP, "mean_sd"),
      `Purification de l'air`     = format_stat(purification, "median_iqr"),
      `Régulation thermique`      = format_stat(temp_regulation, "median_iqr"),
      `Pollution de l'air (PM10)` = format_stat(pm10, "median_iqr"),
      `Risque de chaleur`         = format_stat(heatrisk, "median_iqr"),
      .groups = "drop"
    ) %>%
    rename(Groupe = 1) %>%
    mutate(Groupe = paste0(Groupe, " (n=", n, ")")) %>%
    select(-n)
}

#blocs
tab_type   <- resumer_groupe(df, "type_ecole")   %>% mutate(Section = "Type d'enseignement", .before = 1)
tab_reseau <- resumer_groupe(df, "LANGUE") %>% mutate(Section = "Réseau linguistique", .before = 1)

tab_rbc <- df %>%
  summarise(
    n = n(),
    `Végétation intérieure`      = format_stat(vegetation, "median_iqr"),
    `Canopée de rue`             = format_stat(canopee, "median_iqr"),
    `EVP accessibles*`           = format_stat(EVP, "mean_sd"),
    `Purification de l'air`     = format_stat(purification, "median_iqr"),
    `Régulation thermique`      = format_stat(temp_regulation, "median_iqr"),
    `Pollution de l'air (PM10)` = format_stat(pm10, "median_iqr"),
    `Risque de chaleur`         = format_stat(heatrisk, "median_iqr")
  ) %>%
  mutate(Groupe = paste0("Région de Bruxelles-Capitale (n=", n, ")"), Section = "Ensemble de l'échantillon") %>%
  select(Section, Groupe, everything(), -n)

tab_commune <- resumer_groupe(df, "COMMUNE") %>%
  mutate(Section = "Par commune", .before = 1) %>%
  arrange(desc(as.numeric(sub(".*n=(\\d+).*", "\\1", Groupe))))  # triées par n décroissant

#tableau
tableau_final <- bind_rows(tab_type, tab_reseau, tab_rbc, tab_commune)

print(tableau_final, n = Inf, width = Inf)

write.csv(tableau_final, "annexe_F_ecoles_secondaires.csv", row.names = FALSE)

### Classement écoles vulnérables  ###

INDICATEURS_FAVORABLES <- c("vegetation", "canopee",
                            "purification", "temp_regulation")
INDICATEURS_DEFAVORABLES <- c("pm10", "heatrisk")

# calcul du score 
df_classement <- df %>%
  mutate(across(all_of(INDICATEURS_FAVORABLES), ~ 1 - ., .names = "inv_{.col}")) %>%
  mutate(score_vulnerabilite = rowMeans(
    select(., starts_with("inv_"), all_of(INDICATEURS_DEFAVORABLES)),
    na.rm = TRUE
  )) %>%
  filter(type_ecole != setdiff(unique(type_ecole), c("Public", "Libre"))[1]) %>%
  arrange(desc(score_vulnerabilite))

# histogramme distribution score
summary(df_classement$score_vulnerabilite)
sd(df_classement$score_vulnerabilite, na.rm = TRUE)

mediane_score_env <- median(df_classement$score_vulnerabilite, na.rm = TRUE)

hist_score_env <- ggplot(df_classement, aes(x = score_vulnerabilite)) +
  geom_histogram(binwidth = 0.05, fill = "#4472C4", color = "white") +
  geom_vline(aes(xintercept = median(score_vulnerabilite)), linetype = "dashed", color = "red") +
  annotate("text", x = mediane_score_env, y = Inf,
           label = paste0("Médiane = ", round(mediane_score_env, 2)),
           color = "red", vjust = 1.5, hjust = -0.1, size = 3.5) +
  scale_x_continuous(breaks = seq(0, 1, by = 0.2)) +
  scale_y_continuous(labels = scales::label_number(accuracy = 1)) +
  labs(
    x = "Score de vulnérabilité",
    y = "Nombre d'écoles"
  ) +
  theme_minimal()

hist_score_env
ggsave("histogramme_score_vulnerabilite_env.png", hist_score_env, width = 7, height = 5, dpi = 300)

# top 10
top10 <- df_classement %>%
  mutate(rang = row_number()) %>%
  select(rang, NOM, COMMUNE, type_ecole, LANGUE, score_vulnerabilite,
         vegetation, canopee, EVP, purification, temp_regulation, pm10, heatrisk) %>%
  head(10)

print(top10, width = Inf)

write.csv(df_classement %>% select(NOM, COMMUNE, type_ecole, LANGUE, score_vulnerabilite),
          "classement_complet_ecoles.csv", row.names = FALSE)
write.csv(top10, "top10_ecoles_prioritaires.csv", row.names = FALSE)
cat("\nExportés : classement_complet_ecoles.csv (toutes les écoles) et",
    "top10_ecoles_prioritaires.csv (les 10 premières)\n")

# -------------------------------------------------------------------------------------
#
env <- read_excel("ecoles.xlsx")
socioeco <- read_excel("ecoles_socioeco.xlsx")

brut <- env %>%
  inner_join(socioeco %>% select(NOM, COMMUNE, RUE, revenu, education),
             by = c("NOM", "COMMUNE", "RUE"))

# Conversion des coordonnees degres-minutes -> metres (Lambert 72)
dms_to_dd <- function(chaine) {
  chaine <- trimws(chaine)
  signe <- ifelse(grepl("^-", chaine), -1, 1)
  nombres <- as.numeric(regmatches(chaine, gregexpr("[0-9]+\\.?[0-9]*", chaine))[[1]])
  signe * (nombres[1] + nombres[2] / 60)
}
brut <- brut %>% mutate(lon = sapply(x, dms_to_dd), lat = sapply(y, dms_to_dd))
brut_sf <- st_as_sf(brut, coords = c("lon", "lat"), crs = 4326, remove = FALSE) %>% st_transform(31370)
cc <- st_coordinates(brut_sf)
brut$x_m <- cc[, 1]
brut$y_m <- cc[, 2]

set.seed(123)
doublons <- duplicated(brut[, c("x_m", "y_m")]) | duplicated(brut[, c("x_m", "y_m")], fromLast = TRUE)
brut$x_m[doublons] <- brut$x_m[doublons] + runif(sum(doublons), -0.5, 0.5)
brut$y_m[doublons] <- brut$y_m[doublons] + runif(sum(doublons), -0.5, 0.5)

df <- brut %>% filter(!is.na(revenu), !is.na(education))

coords <- as.matrix(df[, c("x_m", "y_m")])


# autocorrélation spatiale du score de vulnérabilité envi

df_classement <- df_classement %>%
  mutate(
    lon = sapply(x, dms_to_dd),
    lat = sapply(y, dms_to_dd)
  )

df_classement_sf <- st_as_sf(df_classement, coords = c("lon", "lat"), crs = 4326, remove = FALSE) %>%
  st_transform(31370)

coords_classement <- st_coordinates(df_classement_sf)
df_classement$x_m <- coords_classement[, 1]
df_classement$y_m <- coords_classement[, 2]

# gestion des doublons de coordonnées (meme logique que pour "brut")
doublons_classement <- duplicated(df_classement[, c("x_m", "y_m")]) |
  duplicated(df_classement[, c("x_m", "y_m")], fromLast = TRUE)
df_classement$x_m[doublons_classement] <- df_classement$x_m[doublons_classement] +
  runif(sum(doublons_classement), -0.5, 0.5)
df_classement$y_m[doublons_classement] <- df_classement$y_m[doublons_classement] +
  runif(sum(doublons_classement), -0.5, 0.5)

# matrice de poids spatiaux k=3, propre a cette population
coords_classement_mat <- as.matrix(df_classement[, c("x_m", "y_m")])
listw_classement <- nb2listw(knn2nb(knearneigh(coords_classement_mat, k = 3)), style = "W")

# indice de Moran sur le score de vulnerabilite
moran_score_vuln <- moran.test(df_classement$score_vulnerabilite, listw_classement)
moran_score_vuln

# version Monte Carlo (plus robuste, ne suppose pas la normalite des donnees)
moran_score_vuln_mc <- moran.mc(df_classement$score_vulnerabilite, listw_classement, nsim = 999)
moran_score_vuln_mc

# Analyses statistiques ##

#variables
IND_ENV <- c("vegetation", "canopee", "purification",
             "temp_regulation", "pm10", "heatrisk")
VAR_SES <- c("revenu", "education")


### Corrélation

indicateurs <- brut %>%
  select(vegetation, canopee, purification, temp_regulation,
         pm10, heatrisk, education, revenu) %>%
  as.matrix()

noms_affiches <- c("Vegetation interieure", "Canopee de rue",
                   "Purification de l'air", "Regulation thermique",
                   "Pollution (PM10)", "Risque de chaleur",
                   "Education", "Revenu")
colnames(indicateurs) <- noms_affiches

res <- rcorr(indicateurs, type = "spearman")
M <- res$r
p_mat <- res$P
vars <- noms_affiches
n <- length(vars)

donnees_plot <- data.frame()
for (i in 2:n) {
  for (j in 1:(i - 1)) {
    r_val <- M[i, j]
    p_val <- p_mat[i, j]
    etoiles <- dplyr::case_when(
      p_val < 0.001 ~ "***",
      p_val < 0.01  ~ "**",
      p_val < 0.05  ~ "*",
      TRUE ~ ""
    )
    donnees_plot <- rbind(donnees_plot, data.frame(
      ligne = vars[i],
      colonne = vars[j],
      r = r_val,
      p = p_val,
      significatif = p_val < 0.05,
      etiquette = paste0(sprintf("%.2f", r_val),
                         ifelse(etoiles != "", paste0("\n", etoiles), ""))
    ))
  }
}
donnees_plot$colonne <- factor(donnees_plot$colonne, levels = vars[1:(n - 1)])
donnees_plot$ligne <- factor(donnees_plot$ligne, levels = rev(vars[2:n]))

p <- ggplot(donnees_plot, aes(x = colonne, y = ligne, fill = r)) +
  geom_tile(color = "white", linewidth = 0.6) +
  geom_text(aes(label = etiquette), size = 3, lineheight = 0.85) +
  scale_fill_gradient2(low = "steelblue4", mid = "white", high = "firebrick3",
                       midpoint = 0, limits = c(-1, 1), name = NULL,
                       guide = guide_colorbar(direction = "horizontal",
                                              barwidth = 15, barheight = 0.8,
                                              title.position = "top")) +
  scale_x_discrete(position = "top") +
  coord_fixed() +
  theme_minimal(base_size = 12) +
  theme(
    axis.title = element_blank(),
    axis.text = element_text(color = "grey20", size = 11),
    axis.text.x.top = element_text(color = "grey20", size = 11, angle = 45, hjust = 0, vjust = 0),
    panel.grid = element_blank(),
    legend.position = "bottom"
  )
p
ggsave("correlation_spearman_ggplot.png", p, width = 9.5, height = 8.5, dpi = 300)

# Corrélation avec correction spatiale 

use_spatialpack <- requireNamespace("SpatialPack", quietly = TRUE)

resultats_bivarie <- data.frame()
for (ind in IND_ENV) {
  for (ses in VAR_SES) {
    paire <- brut %>% filter(!is.na(.data[[ind]]), !is.na(.data[[ses]]))
    
    rho_classique <- cor.test(paire[[ind]], paire[[ses]], method = "spearman")
    
    if (use_spatialpack) {
      test_spatial <- SpatialPack::modified.ttest(
        x = rank(paire[[ind]]), y = rank(paire[[ses]]),
        coords = as.matrix(paire[, c("x_m", "y_m")])
      )
      p_final <- test_spatial$p.value
      n_eff   <- test_spatial$dof + 2
    } else {
      p_final <- rho_classique$p.value
      n_eff   <- NA
    }
    
    resultats_bivarie <- rbind(resultats_bivarie, data.frame(
      indicateur_env = ind, variable_ses = ses,
      rho = round(unname(rho_classique$estimate), 3),
      p_classique = signif(rho_classique$p.value, 3),   
      p_spatial   = signif(p_final, 3),                 
      n_effectif  = round(n_eff, 1),                    
      n = nrow(paire)
    ))
  }
}
print(resultats_bivarie)
write.csv(resultats_bivarie, "resultats_bivaries_detail.csv", row.names = FALSE)

resultats_bivarie <- resultats_bivarie %>%
  mutate(
    etoiles = case_when(
      p_spatial < 0.001 ~ "***",
      p_spatial < 0.01  ~ "**",
      p_spatial < 0.05  ~ "*",
      TRUE ~ ""
    ),
    etiquette = paste0(sprintf("%.2f", rho),
                       ifelse(etoiles != "", paste0("\n", etoiles), ""))
  )

noms_lisibles <- c(
  vegetation = "Végétation intérieure", canopee = "Canopée de rue",
  purification = "Purification de l'air", temp_regulation = "Régulation thermique",
  pm10 = "Pollution (PM10)", heatrisk = "Risque de chaleur"
)
noms_ses <- c(revenu = "Revenu", education = "Education")

resultats_bivarie$indicateur_env <- factor(noms_lisibles[resultats_bivarie$indicateur_env],
                                           levels = noms_lisibles)
resultats_bivarie$variable_ses <- factor(noms_ses[resultats_bivarie$variable_ses],
                                         levels = noms_ses)

graph_bivarie <- ggplot(resultats_bivarie, aes(x = indicateur_env, y = variable_ses, fill = rho)) +
  geom_tile(color = "white", linewidth = 0.6) +
  geom_text(aes(label = etiquette), size = 3.5, lineheight = 0.85) +
  scale_fill_gradient2(low = "steelblue4", mid = "white", high = "firebrick3",
                       midpoint = 0, limits = c(-1, 1), name = NULL,
                       guide = guide_colorbar(direction = "horizontal",
                                              barwidth = 12, barheight = 0.8,
                                              title.position = "top")) +
  scale_x_discrete(position = "top") +
  coord_fixed() +
  theme_minimal(base_size = 12) +
  theme(
    axis.title = element_blank(),
    axis.text = element_text(color = "grey20", size = 11),
    axis.text.x.top = element_text(color = "grey20", size = 11, angle = 45, hjust = 0, vjust = 0),
    panel.grid = element_blank(),
    legend.position = "bottom"
  )

graph_bivarie
ggsave("correlation_bivariee_spatiale.png", graph_bivarie, width = 9, height = 4, dpi = 300)

# -----------------------------------------------------------------------------
### Régressions  ###

listw <- nb2listw(knn2nb(knearneigh(coords, k = 3)), style = "W")

extraire_p <- function(tests, ancien, nouveau) {
  if (!is.null(tests[[ancien]])) return(tests[[ancien]]$p.value)
  if (!is.null(tests[[nouveau]])) return(tests[[nouveau]]$p.value)
  NA
}

### vegetation ###

#Modele OLS
veg_ols <- lm(vegetation ~ revenu + education, data = df)
summary(veg_ols)

#Heteroscedasticite (test de Breusch-Pagan)
bptest(veg_ols)

#Normalite des residus (test de Shapiro-Wilk)
shapiro.test(residuals(veg_ols))

#Multicollinearite
vif(veg_ols)

#Autocorrelation spatiale des residus (I de Moran)
veg_moran <- lm.morantest(veg_ols, listw)
veg_moran

#LM test (pour choisir OLS, SLAG ou SEM)
veg_lm <- lm.LMtests(veg_ols, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
veg_lm

p_lmerr  <- extraire_p(veg_lm, "LMerr", "RSerr")
p_lmlag  <- extraire_p(veg_lm, "LMlag", "RSlag")
p_rlmerr <- extraire_p(veg_lm, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(veg_lm, "RLMlag", "adjRSlag")

#modele spatial
if (veg_moran$p.value >= 0.05) {
  veg_type <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  veg_type <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  veg_type <- "SLAG"
} else if (p_lmerr < 0.05) {
  veg_type <- "SEM"
} else {
  veg_type <- "OLS"
}

if (veg_type == "SLAG") {
  veg_final <- lagsarlm(vegetation ~ revenu + education, data = df, listw = listw)
} else if (veg_type == "SEM") {
  veg_final <- errorsarlm(vegetation ~ revenu + education, data = df, listw = listw)
} else {
  veg_final <- veg_ols
}

cat("\n>>> Modele retenu pour vegetation :", veg_type, "\n")
summary(veg_final)

### canopee ###

# OLS
canopee_ols <- lm(canopee ~ revenu + education, data = df)
summary(canopee_ols)

#Heteroscedasticite
bptest(canopee_ols)

#Normalite des residus
shapiro.test(residuals(canopee_ols))

#Multicollinearite
vif(canopee_ols)

#Autocorrelation spatiale
canopee_moran <- lm.morantest(canopee_ols, listw)
canopee_moran

# LM tests
canopee_lm <- lm.LMtests(canopee_ols, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
canopee_lm

p_lmerr  <- extraire_p(canopee_lm, "LMerr", "RSerr")
p_lmlag  <- extraire_p(canopee_lm, "LMlag", "RSlag")
p_rlmerr <- extraire_p(canopee_lm, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(canopee_lm, "RLMlag", "adjRSlag")

# modele spatial
if (canopee_moran$p.value >= 0.05) {
  canopee_type <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  canopee_type <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  canopee_type <- "SLAG"
} else if (p_lmerr < 0.05) {
  canopee_type <- "SEM"
} else {
  canopee_type <- "OLS"
}

if (canopee_type == "SLAG") {
  canopee_final <- lagsarlm(canopee ~ revenu + education, data = df, listw = listw)
} else if (canopee_type == "SEM") {
  canopee_final <- errorsarlm(canopee ~ revenu + education, data = df, listw = listw)
} else {
  canopee_final <- canopee_ols
}

cat("\n>>> Modele retenu pour canopee :", canopee_type, "\n")
summary(canopee_final)

### EVP ###

#OLS
evp_ols <- lm(EVP ~ revenu + education, data = df)
summary(evp_ols)

bptest(evp_ols)

shapiro.test(residuals(evp_ols))

vif(evp_ols)

# autocorrélation spatiale 
evp_moran <- lm.morantest(evp_ols, listw)
evp_moran

# lm tests
evp_lm <- lm.LMtests(evp_ols, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
evp_lm

p_lmerr  <- extraire_p(evp_lm, "LMerr", "RSerr")
p_lmlag  <- extraire_p(evp_lm, "LMlag", "RSlag")
p_rlmerr <- extraire_p(evp_lm, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(evp_lm, "RLMlag", "adjRSlag")

# modèle spatial
if (evp_moran$p.value >= 0.05) {
  evp_type <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  evp_type <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  evp_type <- "SLAG"
} else if (p_lmerr < 0.05) {
  evp_type <- "SEM"
} else {
  evp_type <- "OLS"
}

if (evp_type == "SLAG") {
  evp_final <- lagsarlm(EVP ~ revenu + education, data = df, listw = listw)
} else if (evp_type == "SEM") {
  evp_final <- errorsarlm(EVP ~ revenu + education, data = df, listw = listw)
} else {
  evp_final <- evp_ols
}

cat("\n>>> Modele retenu pour EVP :", evp_type, "\n")
summary(evp_final)

### purification ###

#OLS
purif_ols <- lm(purification ~ revenu + education, data = df)
summary(purif_ols)

bptest(purif_ols)

shapiro.test(residuals(purif_ols))

vif(purif_ols)

purif_moran <- lm.morantest(purif_ols, listw)
purif_moran

# lm tests
purif_lm <- lm.LMtests(purif_ols, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
purif_lm

p_lmerr  <- extraire_p(purif_lm, "LMerr", "RSerr")
p_lmlag  <- extraire_p(purif_lm, "LMlag", "RSlag")
p_rlmerr <- extraire_p(purif_lm, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(purif_lm, "RLMlag", "adjRSlag")

# modèle spatial
if (purif_moran$p.value >= 0.05) {
  purif_type <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  purif_type <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  purif_type <- "SLAG"
} else if (p_lmerr < 0.05) {
  purif_type <- "SEM"
} else {
  purif_type <- "OLS"
}

if (purif_type == "SLAG") {
  purif_final <- lagsarlm(purification ~ revenu + education, data = df, listw = listw)
} else if (purif_type == "SEM") {
  purif_final <- errorsarlm(purification ~ revenu + education, data = df, listw = listw)
} else {
  purif_final <- purif_ols
}

cat("\n>>> Modele retenu pour purification :", purif_type, "\n")
summary(purif_final)

### temp_regulation ###

# OLS
temp_ols <- lm(temp_regulation ~ revenu + education, data = df)
summary(temp_ols)

bptest(temp_ols)

shapiro.test(residuals(temp_ols))

#autocorrélation spatiale
temp_moran <- lm.morantest(temp_ols, listw)
temp_moran

#lm tests
temp_lm <- lm.LMtests(temp_ols, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
temp_lm

p_lmerr  <- extraire_p(temp_lm, "LMerr", "RSerr")
p_lmlag  <- extraire_p(temp_lm, "LMlag", "RSlag")
p_rlmerr <- extraire_p(temp_lm, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(temp_lm, "RLMlag", "adjRSlag")

# modèle spatial
if (temp_moran$p.value >= 0.05) {
  temp_type <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  temp_type <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  temp_type <- "SLAG"
} else if (p_lmerr < 0.05) {
  temp_type <- "SEM"
} else {
  temp_type <- "OLS"
}

if (temp_type == "SLAG") {
  temp_final <- lagsarlm(temp_regulation ~ revenu + education, data = df, listw = listw)
} else if (temp_type == "SEM") {
  temp_final <- errorsarlm(temp_regulation ~ revenu + education, data = df, listw = listw)
} else {
  temp_final <- temp_ols
}

cat("\n>>> Modele retenu pour temp_regulation :", temp_type, "\n")
summary(temp_final)


### pm10 ###

#OLS 
pm10_ols <- lm(pm10 ~ revenu + education, data = df)
summary(pm10_ols)

bptest(pm10_ols)

shapiro.test(residuals(pm10_ols))

#autocorrélation spatiale
pm10_moran <- lm.morantest(pm10_ols, listw)
pm10_moran

# lm tests
pm10_lm <- lm.LMtests(pm10_ols, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
pm10_lm

p_lmerr  <- extraire_p(pm10_lm, "LMerr", "RSerr")
p_lmlag  <- extraire_p(pm10_lm, "LMlag", "RSlag")
p_rlmerr <- extraire_p(pm10_lm, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(pm10_lm, "RLMlag", "adjRSlag")

# modèle spatial
if (pm10_moran$p.value >= 0.05) {
  pm10_type <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  pm10_type <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  pm10_type <- "SLAG"
} else if (p_lmerr < 0.05) {
  pm10_type <- "SEM"
} else {
  pm10_type <- "OLS"
}

if (pm10_type == "SLAG") {
  pm10_final <- lagsarlm(pm10 ~ revenu + education, data = df, listw = listw)
} else if (pm10_type == "SEM") {
  pm10_final <- errorsarlm(pm10 ~ revenu + education, data = df, listw = listw)
} else {
  pm10_final <- pm10_ols
}

cat("\n>>> Modele retenu pour pm10 :", pm10_type, "\n")
summary(pm10_final)

### risque de chaleur ###

#OLS
heat_ols <- lm(heatrisk ~ revenu + education, data = df)
summary(heat_ols)

bptest(heat_ols)

shapiro.test(residuals(heat_ols))

# autocorrélation spatiale 
heat_moran <- lm.morantest(heat_ols, listw)
heat_moran

# lm tests
heat_lm <- lm.LMtests(heat_ols, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
heat_lm

p_lmerr  <- extraire_p(heat_lm, "LMerr", "RSerr")
p_lmlag  <- extraire_p(heat_lm, "LMlag", "RSlag")
p_rlmerr <- extraire_p(heat_lm, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(heat_lm, "RLMlag", "adjRSlag")

# modèle spatial
if (heat_moran$p.value >= 0.05) {
  heat_type <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  heat_type <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  heat_type <- "SLAG"
} else if (p_lmerr < 0.05) {
  heat_type <- "SEM"
} else {
  heat_type <- "OLS"
}

if (heat_type == "SLAG") {
  heat_final <- lagsarlm(heatrisk ~ revenu + education, data = df, listw = listw)
} else if (heat_type == "SEM") {
  heat_final <- errorsarlm(heatrisk ~ revenu + education, data = df, listw = listw)
} else {
  heat_final <- heat_ols
}

cat("\n>>> Modele retenu pour heatrisk :", heat_type, "\n")
summary(heat_final)

# ---------------------------------------------------------------------
# tableau résultats

extraire_ligne_finale <- function(nom_var, prefix_type, modele) {
  if (class(modele)[1] == "lm") {
    coefs <- summary(modele)$coefficients
    r2 <- summary(modele)$r.squared
    lag_coeff <- NA   # pas de coefficient de lag/erreur pour un OLS
  } else if (prefix_type == "SLAG") {
    coefs <- summary(modele)$Coef
    r2 <- cor(df[[nom_var]], fitted(modele))^2
    lag_coeff <- round(modele$rho, 3)
  } else {
    coefs <- summary(modele)$Coef
    r2 <- cor(df[[nom_var]], fitted(modele))^2
    lag_coeff <- round(modele$lambda, 3)
  }
  
  get_val <- function(terme, colonne) {
    if (terme %in% rownames(coefs)) coefs[terme, colonne] else NA
  }
  
  data.frame(
    variable = nom_var,
    modele = prefix_type,
    intercept = round(get_val("(Intercept)", 1), 3),
    p_intercept = signif(get_val("(Intercept)", 4), 3),
    coef_revenu = round(get_val("revenu", 1), 3),
    p_revenu = signif(get_val("revenu", 4), 3),
    coef_education = round(get_val("education", 1), 3),
    p_education = signif(get_val("education", 4), 3),
    R2 = round(r2, 3),
    lag_coeff = lag_coeff,
    AIC = round(AIC(modele), 1),
    logLik = round(as.numeric(logLik(modele)), 1)
  )
}

tableau_resultats <- rbind(
  extraire_ligne_finale("vegetation", veg_type, veg_final),
  extraire_ligne_finale("canopee", canopee_type, canopee_final),
  extraire_ligne_finale("EVP", evp_type, evp_final),
  extraire_ligne_finale("purification", purif_type, purif_final),
  extraire_ligne_finale("temp_regulation", temp_type, temp_final),
  extraire_ligne_finale("pm10", pm10_type, pm10_final),
  extraire_ligne_finale("heatrisk", heat_type, heat_final)
)

print(tableau_resultats, row.names = FALSE)

library(openxlsx)
write.xlsx(tableau_resultats, "tableau_resultats_normalise.xlsx", overwrite = TRUE)

# --------------------------------------------------------------------------
# tableau annexe 

extraire_ligne_modele <- function(nom_var, nom_modele, modele) {
  if (class(modele)[1] == "lm") {
    coefs <- summary(modele)$coefficients
    r2 <- summary(modele)$r.squared
  } else {
    coefs <- summary(modele)$Coef
    r2 <- cor(df[[nom_var]], fitted(modele))^2
  }
  data.frame(
    variable = nom_var, modele = nom_modele,
    intercept = round(coefs["(Intercept)", 1], 3),
    p_intercept = signif(coefs["(Intercept)", 4], 3),
    coef_revenu = round(coefs["revenu", 1], 3),
    p_revenu = signif(coefs["revenu", 4], 3),
    coef_education = round(coefs["education", 1], 3),
    p_education = signif(coefs["education", 4], 3),
    R2 = round(r2, 3),
    AIC = round(AIC(modele), 1),
    logLik = round(as.numeric(logLik(modele)), 1)
  )
}

tableau_annexe <- rbind(
  extraire_ligne_modele("vegetation", "OLS", veg_ols),
  extraire_ligne_modele("canopee", "OLS", canopee_ols),
  extraire_ligne_modele("canopee", canopee_type, canopee_final),
  extraire_ligne_modele("EVP", "OLS", evp_ols),
  extraire_ligne_modele("EVP", evp_type, evp_final),
  extraire_ligne_modele("purification", "OLS", purif_ols),
  extraire_ligne_modele("purification", purif_type, purif_final),
  extraire_ligne_modele("temp_regulation", "OLS", temp_ols),
  extraire_ligne_modele("temp_regulation", temp_type, temp_final),
  extraire_ligne_modele("pm10", "OLS", pm10_ols),
  extraire_ligne_modele("pm10", pm10_type, pm10_final),
  extraire_ligne_modele("heatrisk", "OLS", heat_ols),
  extraire_ligne_modele("heatrisk", heat_type, heat_final)
)

print(tableau_annexe, row.names = FALSE)

write.xlsx(tableau_annexe, "tableau_annexe_normalise.xlsx", overwrite = TRUE)

 #### Régression - revenu seul  #### 

#OLS
veg_ols_rs <- lm(vegetation ~ revenu, data = df)
summary(veg_ols_rs)

bptest(veg_ols_rs)

shapiro.test(residuals(veg_ols_rs))

#autocorrélation spatiale
veg_moran_rs <- lm.morantest(veg_ols_rs, listw)
veg_moran_rs

#lm tests
veg_lm_rs <- lm.LMtests(veg_ols_rs, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
veg_lm_rs

p_lmerr  <- extraire_p(veg_lm_rs, "LMerr", "RSerr")
p_lmlag  <- extraire_p(veg_lm_rs, "LMlag", "RSlag")
p_rlmerr <- extraire_p(veg_lm_rs, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(veg_lm_rs, "RLMlag", "adjRSlag")

#modèle spatial
if (veg_moran_rs$p.value >= 0.05) {
  veg_type_rs <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  veg_type_rs <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  veg_type_rs <- "SLAG"
} else if (p_lmerr < 0.05) {
  veg_type_rs <- "SEM"
} else {
  veg_type_rs <- "OLS"
}

if (veg_type_rs == "SLAG") {
  veg_final_rs <- lagsarlm(vegetation ~ revenu, data = df, listw = listw)
} else if (veg_type_rs == "SEM") {
  veg_final_rs <- errorsarlm(vegetation ~ revenu, data = df, listw = listw)
} else {
  veg_final_rs <- veg_ols_rs
}

cat("\n>>> Modele retenu pour vegetation (revenu seul) :", veg_type_rs, "\n")
summary(veg_final_rs)

### canopee -- revenu seul

#OLS
canopee_ols_rs <- lm(canopee ~ revenu, data = df)
summary(canopee_ols_rs)

# 2. Heteroscedasticite (test de Breusch-Pagan)
bptest(canopee_ols_rs)

# 3. Normalite des residus (test de Shapiro-Wilk)
shapiro.test(residuals(canopee_ols_rs))

# 4. Autocorrelation spatiale des residus (I de Moran)
canopee_moran_rs <- lm.morantest(canopee_ols_rs, listw)
canopee_moran_rs

# 5. Tests du multiplicateur de Lagrange (pour choisir OLS, SLAG ou SEM)
canopee_lm_rs <- lm.LMtests(canopee_ols_rs, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
canopee_lm_rs

p_lmerr  <- extraire_p(canopee_lm_rs, "LMerr", "RSerr")
p_lmlag  <- extraire_p(canopee_lm_rs, "LMlag", "RSlag")
p_rlmerr <- extraire_p(canopee_lm_rs, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(canopee_lm_rs, "RLMlag", "adjRSlag")

# 6. Choix et ajustement du modele spatial final
if (canopee_moran_rs$p.value >= 0.05) {
  canopee_type_rs <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  canopee_type_rs <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  canopee_type_rs <- "SLAG"
} else if (p_lmerr < 0.05) {
  canopee_type_rs <- "SEM"
} else {
  canopee_type_rs <- "OLS"
}

if (canopee_type_rs == "SLAG") {
  canopee_final_rs <- lagsarlm(canopee ~ revenu, data = df, listw = listw)
} else if (canopee_type_rs == "SEM") {
  canopee_final_rs <- errorsarlm(canopee ~ revenu, data = df, listw = listw)
} else {
  canopee_final_rs <- canopee_ols_rs
}

cat("\n>>> Modele retenu pour canopee (revenu seul) :", canopee_type_rs, "\n")
summary(canopee_final_rs)


# =======================================================================
# EVP -- revenu seul
# =======================================================================

# 1. Modele OLS de base
evp_ols_rs <- lm(EVP ~ revenu, data = df)
summary(evp_ols_rs)

# 2. Heteroscedasticite (test de Breusch-Pagan)
bptest(evp_ols_rs)

# 3. Normalite des residus (test de Shapiro-Wilk)
shapiro.test(residuals(evp_ols_rs))

# 4. Autocorrelation spatiale des residus (I de Moran)
evp_moran_rs <- lm.morantest(evp_ols_rs, listw)
evp_moran_rs

# 5. Tests du multiplicateur de Lagrange (pour choisir OLS, SLAG ou SEM)
evp_lm_rs <- lm.LMtests(evp_ols_rs, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
evp_lm_rs

p_lmerr  <- extraire_p(evp_lm_rs, "LMerr", "RSerr")
p_lmlag  <- extraire_p(evp_lm_rs, "LMlag", "RSlag")
p_rlmerr <- extraire_p(evp_lm_rs, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(evp_lm_rs, "RLMlag", "adjRSlag")

# 6. Choix et ajustement du modele spatial final
if (evp_moran_rs$p.value >= 0.05) {
  evp_type_rs <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  evp_type_rs <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  evp_type_rs <- "SLAG"
} else if (p_lmerr < 0.05) {
  evp_type_rs <- "SEM"
} else {
  evp_type_rs <- "OLS"
}

if (evp_type_rs == "SLAG") {
  evp_final_rs <- lagsarlm(EVP ~ revenu, data = df, listw = listw)
} else if (evp_type_rs == "SEM") {
  evp_final_rs <- errorsarlm(EVP ~ revenu, data = df, listw = listw)
} else {
  evp_final_rs <- evp_ols_rs
}

cat("\n>>> Modele retenu pour EVP (revenu seul) :", evp_type_rs, "\n")
summary(evp_final_rs)


# =======================================================================
# purification -- revenu seul
# =======================================================================

# 1. Modele OLS de base
purif_ols_rs <- lm(purification ~ revenu, data = df)
summary(purif_ols_rs)

# 2. Heteroscedasticite (test de Breusch-Pagan)
bptest(purif_ols_rs)

# 3. Normalite des residus (test de Shapiro-Wilk)
shapiro.test(residuals(purif_ols_rs))

# 4. Autocorrelation spatiale des residus (I de Moran)
purif_moran_rs <- lm.morantest(purif_ols_rs, listw)
purif_moran_rs

# 5. Tests du multiplicateur de Lagrange (pour choisir OLS, SLAG ou SEM)
purif_lm_rs <- lm.LMtests(purif_ols_rs, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
purif_lm_rs

p_lmerr  <- extraire_p(purif_lm_rs, "LMerr", "RSerr")
p_lmlag  <- extraire_p(purif_lm_rs, "LMlag", "RSlag")
p_rlmerr <- extraire_p(purif_lm_rs, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(purif_lm_rs, "RLMlag", "adjRSlag")

# 6. Choix et ajustement du modele spatial final
if (purif_moran_rs$p.value >= 0.05) {
  purif_type_rs <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  purif_type_rs <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  purif_type_rs <- "SLAG"
} else if (p_lmerr < 0.05) {
  purif_type_rs <- "SEM"
} else {
  purif_type_rs <- "OLS"
}

if (purif_type_rs == "SLAG") {
  purif_final_rs <- lagsarlm(purification ~ revenu, data = df, listw = listw)
} else if (purif_type_rs == "SEM") {
  purif_final_rs <- errorsarlm(purification ~ revenu, data = df, listw = listw)
} else {
  purif_final_rs <- purif_ols_rs
}

cat("\n>>> Modele retenu pour purification (revenu seul) :", purif_type_rs, "\n")
summary(purif_final_rs)


# =======================================================================
# temp_regulation -- revenu seul
# =======================================================================

# 1. Modele OLS de base
temp_ols_rs <- lm(temp_regulation ~ revenu, data = df)
summary(temp_ols_rs)

# 2. Heteroscedasticite (test de Breusch-Pagan)
bptest(temp_ols_rs)

# 3. Normalite des residus (test de Shapiro-Wilk)
shapiro.test(residuals(temp_ols_rs))

# 4. Autocorrelation spatiale des residus (I de Moran)
temp_moran_rs <- lm.morantest(temp_ols_rs, listw)
temp_moran_rs

# 5. Tests du multiplicateur de Lagrange (pour choisir OLS, SLAG ou SEM)
temp_lm_rs <- lm.LMtests(temp_ols_rs, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
temp_lm_rs

p_lmerr  <- extraire_p(temp_lm_rs, "LMerr", "RSerr")
p_lmlag  <- extraire_p(temp_lm_rs, "LMlag", "RSlag")
p_rlmerr <- extraire_p(temp_lm_rs, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(temp_lm_rs, "RLMlag", "adjRSlag")

# 6. Choix et ajustement du modele spatial final
if (temp_moran_rs$p.value >= 0.05) {
  temp_type_rs <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  temp_type_rs <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  temp_type_rs <- "SLAG"
} else if (p_lmerr < 0.05) {
  temp_type_rs <- "SEM"
} else {
  temp_type_rs <- "OLS"
}

if (temp_type_rs == "SLAG") {
  temp_final_rs <- lagsarlm(temp_regulation ~ revenu, data = df, listw = listw)
} else if (temp_type_rs == "SEM") {
  temp_final_rs <- errorsarlm(temp_regulation ~ revenu, data = df, listw = listw)
} else {
  temp_final_rs <- temp_ols_rs
}

cat("\n>>> Modele retenu pour temp_regulation (revenu seul) :", temp_type_rs, "\n")
summary(temp_final_rs)


# =======================================================================
# pm10 -- revenu seul
# =======================================================================

# 1. Modele OLS de base
pm10_ols_rs <- lm(pm10 ~ revenu, data = df)
summary(pm10_ols_rs)

# 2. Heteroscedasticite (test de Breusch-Pagan)
bptest(pm10_ols_rs)

# 3. Normalite des residus (test de Shapiro-Wilk)
shapiro.test(residuals(pm10_ols_rs))

# 4. Autocorrelation spatiale des residus (I de Moran)
pm10_moran_rs <- lm.morantest(pm10_ols_rs, listw)
pm10_moran_rs

# 5. Tests du multiplicateur de Lagrange (pour choisir OLS, SLAG ou SEM)
pm10_lm_rs <- lm.LMtests(pm10_ols_rs, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
pm10_lm_rs

p_lmerr  <- extraire_p(pm10_lm_rs, "LMerr", "RSerr")
p_lmlag  <- extraire_p(pm10_lm_rs, "LMlag", "RSlag")
p_rlmerr <- extraire_p(pm10_lm_rs, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(pm10_lm_rs, "RLMlag", "adjRSlag")

# 6. Choix et ajustement du modele spatial final
if (pm10_moran_rs$p.value >= 0.05) {
  pm10_type_rs <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  pm10_type_rs <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  pm10_type_rs <- "SLAG"
} else if (p_lmerr < 0.05) {
  pm10_type_rs <- "SEM"
} else {
  pm10_type_rs <- "OLS"
}

if (pm10_type_rs == "SLAG") {
  pm10_final_rs <- lagsarlm(pm10 ~ revenu, data = df, listw = listw)
} else if (pm10_type_rs == "SEM") {
  pm10_final_rs <- errorsarlm(pm10 ~ revenu, data = df, listw = listw)
} else {
  pm10_final_rs <- pm10_ols_rs
}

cat("\n>>> Modele retenu pour pm10 (revenu seul) :", pm10_type_rs, "\n")
summary(pm10_final_rs)


# =======================================================================
# heatrisk -- revenu seul
# =======================================================================

# 1. Modele OLS de base
heat_ols_rs <- lm(heatrisk ~ revenu, data = df)
summary(heat_ols_rs)

# 2. Heteroscedasticite (test de Breusch-Pagan)
bptest(heat_ols_rs)

# 3. Normalite des residus (test de Shapiro-Wilk)
shapiro.test(residuals(heat_ols_rs))

# 4. Autocorrelation spatiale des residus (I de Moran)
heat_moran_rs <- lm.morantest(heat_ols_rs, listw)
heat_moran_rs

# 5. Tests du multiplicateur de Lagrange (pour choisir OLS, SLAG ou SEM)
heat_lm_rs <- lm.LMtests(heat_ols_rs, listw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag"))
heat_lm_rs

p_lmerr  <- extraire_p(heat_lm_rs, "LMerr", "RSerr")
p_lmlag  <- extraire_p(heat_lm_rs, "LMlag", "RSlag")
p_rlmerr <- extraire_p(heat_lm_rs, "RLMerr", "adjRSerr")
p_rlmlag <- extraire_p(heat_lm_rs, "RLMlag", "adjRSlag")

# 6. Choix et ajustement du modele spatial final
if (heat_moran_rs$p.value >= 0.05) {
  heat_type_rs <- "OLS"
} else if (p_lmlag < 0.05 & p_lmerr < 0.05) {
  heat_type_rs <- if (p_rlmlag < p_rlmerr) "SLAG" else "SEM"
} else if (p_lmlag < 0.05) {
  heat_type_rs <- "SLAG"
} else if (p_lmerr < 0.05) {
  heat_type_rs <- "SEM"
} else {
  heat_type_rs <- "OLS"
}

if (heat_type_rs == "SLAG") {
  heat_final_rs <- lagsarlm(heatrisk ~ revenu, data = df, listw = listw)
} else if (heat_type_rs == "SEM") {
  heat_final_rs <- errorsarlm(heatrisk ~ revenu, data = df, listw = listw)
} else {
  heat_final_rs <- heat_ols_rs
}

cat("\n>>> Modele retenu pour heatrisk (revenu seul) :", heat_type_rs, "\n")
summary(heat_final_rs)


# =======================================================================
# TABLEAU REVENU SEUL
# =======================================================================

extraire_ligne_revenu_seul <- function(nom_var, prefix_type, modele) {
  if (class(modele)[1] == "lm") {
    coefs <- summary(modele)$coefficients
    r2 <- summary(modele)$r.squared
    lag_coeff <- NA
  } else if (prefix_type == "SLAG") {
    coefs <- summary(modele)$Coef
    r2 <- cor(df[[nom_var]], fitted(modele))^2
    lag_coeff <- round(modele$rho, 3)
  } else {
    coefs <- summary(modele)$Coef
    r2 <- cor(df[[nom_var]], fitted(modele))^2
    lag_coeff <- round(modele$lambda, 3)
  }
  
  get_val <- function(terme, colonne) {
    if (terme %in% rownames(coefs)) coefs[terme, colonne] else NA
  }
  
  avec_etoiles <- function(estimation, p_val) {
    etoiles <- dplyr::case_when(
      is.na(p_val) ~ "",
      p_val < 0.001 ~ "***",
      p_val < 0.01  ~ "**",
      p_val < 0.05  ~ "*",
      TRUE ~ ""
    )
    paste0(round(estimation, 3), etoiles)
  }
  
  data.frame(
    variable = nom_var,
    modele = prefix_type,
    intercept = round(get_val("(Intercept)", 1), 3),
    p_intercept = signif(get_val("(Intercept)", 4), 3),
    coef_revenu = round(get_val("revenu", 1), 3),
    p_revenu = signif(get_val("revenu", 4), 3),
    revenu_signif = avec_etoiles(get_val("revenu", 1), get_val("revenu", 4)),
    R2 = round(r2, 3),
    lag_coeff = lag_coeff,
    AIC = round(AIC(modele), 1),
    logLik = round(as.numeric(logLik(modele)), 1)
  )
}

tableau_revenu_seul <- rbind(
  extraire_ligne_revenu_seul("vegetation", veg_type_rs, veg_final_rs),
  extraire_ligne_revenu_seul("canopee", canopee_type_rs, canopee_final_rs),
  extraire_ligne_revenu_seul("EVP", evp_type_rs, evp_final_rs),
  extraire_ligne_revenu_seul("purification", purif_type_rs, purif_final_rs),
  extraire_ligne_revenu_seul("temp_regulation", temp_type_rs, temp_final_rs),
  extraire_ligne_revenu_seul("pm10", pm10_type_rs, pm10_final_rs),
  extraire_ligne_revenu_seul("heatrisk", heat_type_rs, heat_final_rs)
)

cat("\n\n========== TABLEAU REVENU SEUL ==========\n")
print(tableau_revenu_seul, row.names = FALSE)

write.xlsx(tableau_revenu_seul, "tableau_revenu_seul_normalise.xlsx", overwrite = TRUE)

# ---------------------------------------------------------------
### Score de vulnérabilité avec données socio-économiques #### (n = 66 écoles)

# Inversion 
df <- df %>%
  mutate(
    inv_revenu    = 1 - revenu,
    inv_education = 1 - education
  )

# Score envi + socioeco 
df_classement_socioeco <- df %>%
  mutate(across(all_of(INDICATEURS_FAVORABLES), ~ 1 - ., .names = "inv_{.col}")) %>%
  mutate(score_vulnerabilite_socioeco = rowMeans(
    select(., starts_with("inv_"), all_of(INDICATEURS_DEFAVORABLES)),
    na.rm = TRUE
  )) %>%
  arrange(desc(score_vulnerabilite_socioeco))

# histogramme
mediane_score_socioeco <- median(df_classement_socioeco$score_vulnerabilite_socioeco, na.rm = TRUE)

hist_score_socioeco <- ggplot(df_classement_socioeco, aes(x = score_vulnerabilite_socioeco)) +
  geom_histogram(binwidth = 0.05, fill = "#4472C4", color = "white") +
  geom_vline(xintercept = mediane_score_socioeco, linetype = "dashed", color = "red") +
  annotate("text", x = mediane_score_socioeco, y = Inf,
           label = paste0("Médiane = ", round(mediane_score_socioeco, 2)),
           color = "red", vjust = 1.5, hjust = -0.1, size = 3.5) +
  scale_x_continuous(breaks = seq(0, 1, by = 0.2)) +
  scale_y_continuous(labels = scales::label_number(accuracy = 1)) +
  labs(
    x = "Score de vulnérabilité",
    y = "Nombre d'écoles"
  ) +
  theme_minimal()

hist_score_socioeco
ggsave("histogramme_score_vulnerabilite.png", hist_score_socioeco, width = 7, height = 5, dpi = 300)

# Autocorrelation spatiale du score integre (reutilise le listw k=3 deja construit)
moran_score_socioeco <- moran.test(df_classement_socioeco$score_vulnerabilite_socioeco, listw)
moran_score_socioeco

#monte carlo ?
# version Monte Carlo (plus robuste, ne suppose pas la normalite des donnees)
moran_socio_mc <- moran.mc(df_classement_socioeco$score_vulnerabilite_socioeco, listw, nsim = 999)
moran_socio_mc

# top 10 du classement integre
top10_socioeco <- df_classement_socioeco %>%
  mutate(rang = row_number()) %>%
  select(rang, NOM, COMMUNE, type_ecole, LANGUE, score_vulnerabilite_socioeco,
         vegetation, canopee, EVP, purification, temp_regulation, pm10, heatrisk,
         revenu, education) %>%
  head(10)

print(top10_socioeco, width = Inf)

comparaison_classements <- df_classement %>%
  select(NOM, COMMUNE, score_vulnerabilite) %>%
  mutate(rang_env = rank(-score_vulnerabilite)) %>%
  inner_join(
    df_classement_socioeco %>%
      select(NOM, COMMUNE, score_vulnerabilite_socioeco) %>%
      mutate(rang_socioeco = rank(-score_vulnerabilite_socioeco)),
    by = c("NOM", "COMMUNE")
  ) %>%
  mutate(delta_rang = rang_env - rang_socioeco) %>%
  arrange(rang_socioeco)

print(comparaison_classements, n = Inf, width = Inf)

# Correlation de Spearman entre les deux classements
cor.test(comparaison_classements$score_vulnerabilite,
         comparaison_classements$score_vulnerabilite_socioeco,
         method = "spearman")

# export
write.csv(df_classement_socioeco %>%
            select(NOM, COMMUNE, type_ecole, LANGUE, score_vulnerabilite_socioeco),
          "classement_complet_ecoles_socioeco.csv", row.names = FALSE)
write.csv(top10_socioeco, "top10_ecoles_prioritaires_socioeco.csv", row.names = FALSE)
write.csv(comparaison_classements, "comparaison_classements_env_vs_socioeco.csv", row.names = FALSE)