library(here)

csv_path <- here("data", "donnees_standr.csv")
immo_df <- read.csv2(csv_path)

#Pour visualiser les données
library(DT)

# datatable(head(immo_df, 10))

library(sf)


shp_path <- here("data", "EPCI.shp")
epci_sf <- st_read(shp_path)

# et la table attributaire correspondante, avec le package DT
 # datatable(head(epci_sf, 10))

library(mapsf)

# pour voir les données géographiques
mf_map(x = epci_sf)

library(ggplot2)

plot1 <- ggplot(immo_df, aes(x = med_niveau_vis, y = prix_med)) + 
                         geom_point() + 
                         labs(title = "Relation prix median au m² vs niveau de vie median par EPCI") + 
                         xlab("niveau de vie median par EPCI (€)") + ylab("Prix median au m² (€)")+
                         theme(plot.title = element_text(face = "bold"))

plot1 + geom_point(col="blue") + geom_smooth(method = "lm", color="red", se=TRUE) + theme_grey()

library(car)

car::scatterplot(prix_med ~ med_niveau_vis, data = immo_df, id = list(n=10)) 

# Estimation du modèle avec la fonction lm (pour linear model)
lm1 <- lm(data = immo_df, formula = prix_med ~ med_niveau_vis)

summary(lm1)

library(gtsummary)

# En tableau
gtsummary::tbl_regression(lm1)

mod.lm <- lm(formula = prix_med ~ perc_log_vac + perc_maison + 
                      perc_tiny_log + dens_pop + med_niveau_vis + 
                      part_log_suroccup + part_agri_nb_emploi +
                      part_cadre_profintellec_nbemploi, data = immo_df)
          

gtsummary::tbl_regression(mod.lm)

GGally::ggcoef_model(mod.lm)

# Distribution de la variable dépendante :
hist(immo_df$prix_med, freq=FALSE, breaks = 35, col = "darkblue", border = "white", 
     main = "Prix median au m² par EPCI")

library(car)

scatterplotMatrix(~prix_med + perc_log_vac + perc_maison + 
                    perc_tiny_log + dens_pop + med_niveau_vis + 
                    part_log_suroccup + part_agri_nb_emploi +
                    part_cadre_profintellec_nbemploi,
                   data=immo_df, smoother=loessLine,
                  main="Matrice de graphiques cartésiens")

library(correlation)

data_cor <- immo_df %>% select(-SIREN)
immo_cor <- correlation(data_cor, redundant = TRUE)
summary(immo_cor)

library(see)
immo_cor %>%
  summary(redundant = FALSE) %>%
  plot(type="tile", show_labels =TRUE, show_p = TRUE, digits = 1, size_text=3)

library(car)

vif(mod.lm)

# On peut aussi directement l'ajouter au résumé des coefficient obtenu avec gtsummary

library(gtsummary)

mod.lm %>%
  tbl_regression(intercept = TRUE) %>% add_vif()

par(mfrow=c(1,3))
plot(rstudent(mod.lm))
hist(rstudent(mod.lm), breaks = 20, col="darkblue", border="white", main="Analyse visuelle des résidus")
qqPlot(mod.lm)

# Pour étudier la normalité on peut utiliser le test de Shapiro-Wilk
shapiro.test(mod.lm$residuals)

# Pour évaluer l'homoscédasticité on peut utiliser le test de Breusch-Pagan. Le package car propose une fonction pour le réaliser
ncvTest(mod.lm)

immo_df[c(36,266),]

mod.lm2 <- update(mod.lm, subset=-266)
summary(mod.lm2)

par(mfrow=c(1,2))
plot(rstudent(mod.lm2)) 
qqPlot(mod.lm2)

car::compareCoefs(mod.lm, mod.lm2, pvals = TRUE)

# l'option all.x = TRUE permet de garder toutes les lignes de epci_sf,
# même celles qui n'ont pas de correspondance dans immo_df
data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN", all.x = TRUE)


nrow(immo_df)

nrow(epci_sf)

nrow(data_immo)

mf_map(x = data_immo)
mf_map(x = data_immo[is.na(data_immo$prix_med),], col = 'red', add = TRUE)

data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN")
nrow(data_immo)

mf_map(x = data_immo)

# Pour exporter 

# st_write(data_immo, "data_immo.gpkg", driver = "GPKG")

# La palette de couleurs :
cols_v1 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272",  "#fb6a4a", "#de2d26")

# Carte du prix médian 
mf_map(x = data_immo, 
       var = "prix_med", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       breaks=quantile(data_immo$prix_med,seq(0,1, by=1/11)), 
       pal=cols_v1, 
       leg_title = "Prix médian", 
       leg_val_rnd = 1)
mf_title("Prix Médian du logement par EPCI France Métropolitaine") #titre

data_immo$res_reg <- mod.lm$residuals

# Définition d'une palette de couleur
cols_v1 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272",  "#fb6a4a", "#de2d26")

# Réalisation de la carte
mf_map(x = data_immo, 
       var = "res_reg", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       breaks=quantile(data_immo$res_reg,seq(0,1, by=1/11)), 
       pal=cols_v1, 
       leg_title = "Résidus de régression\nlinéaire 'classique'", 
       leg_val_rnd = 1)
mf_title("Résidus modèle lm") #titre

library(spdep)
# Création de la liste des voisins : avec l'option queen = TRUE, 
# sont considérés comme voisins 2 polygones possédant au moins 1 sommet commun
#help(poly2nb)
neighbours_epci <- poly2nb(data_immo, queen = TRUE) 

# cette opération se fait aussi avec un objet sp 
#shape_nbq <- poly2nb(shape, queen = TRUE)

# Obtention des coordonnées des centroïdes
#coord<-coordinates(shape)
coord <- st_coordinates(st_centroid(data_immo))

# Faire un graphe de voisinage
plot(neighbours_epci, coord)

# si on prend le 1er élément de neighbours_epci, on voit qu'il a pour voisins les polygones 62, 74 etc.
neighbours_epci[[1]]

# ce qu'on peut vérifier sur la carte :
neighbours1 <- data_immo[c(1,62,74,338,1135,1136,1137,1140),]
neighbours1$index <- rownames(neighbours1)
mf_map(x = neighbours1)
mf_label(x = neighbours1, var = "index")

# la fonction nb2listw attribue des poids à chaque voisin
# par ex. si un polygone a 4 voisins, ils auront chacun un poids de 1/4 = 0.25
#help("nb2listw")
neighbours_epci_w <- nb2listw(neighbours_epci)
# les poids sont stockés dans le 3ème élément de neighbours_epci_w
# par ex. si on veut connaître les poids des voisins du 1er élément :
neighbours_epci_w[[3]][1]

# La palette de couleurs :
cols_v1 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272",  "#fb6a4a", "#de2d26")

# Carte du prix médian 
mf_map(x = data_immo, 
       var = "prix_med", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       breaks=quantile(data_immo$prix_med,seq(0,1, by=1/11)), 
       pal=cols_v1, 
       leg_title = "Prix médian", 
       leg_val_rnd = 1)
mf_title("Prix Médian du logement par EPCI France Métropolitaine") #titre

library(spdep)

# Pour l'occasion on va standardiser notre prix médian. 
# Cela permettra par la suite de le comparer à d'autres variables si d'autres analyses d'autocorrélation spatiale sont réalisées
data_immo$prix_med_z <- (data_immo$prix_med-mean(data_immo$prix_med))/sd(data_immo$prix_med)


# Test de Moran :

# On indique dans un premier temps la variable que l'on souhaite analyser
# Puis la matrice de voisinage
# L'argument zero.policy=TRUE permet de préciser que l'on souhaite intégrer à l'analyse les entités spatiales qui n'auraient pas de voisins
# L'argument randomisation=FALSE transmet à la fonction que nous supposons que la distribution est normale. Dans le cas contraire on devrait partir sur une solution de type Monte-Carlo qui repose sur la randomisation
moran.test(data_immo$prix_med_z, neighbours_epci_w, zero.policy = TRUE, randomisation = FALSE)

moran.plot(data_immo$prix_med_z,neighbours_epci_w, 
           labels = TRUE, 
           xlab="prix medians par epci", 
           ylab="moyenne du prix médians par epci des voisins")

library(rgeoda)

# Pour utiliser la fonction local_moran proposée par le package rgeoda 2 pré-requis:

# 1- utiliser la fonction queen_weights du package rgeoda pour calculer une matrice de contiguïté de type queen 
queen_w <- queen_weights(data_immo)
# 2- Sortir la variable à étudier dans un vecteur
prix_med_z = data_immo["prix_med_z"]

lisa <- local_moran(queen_w, prix_med_z)

# Pour visualiser les résultats des LISA il faut les stocker dans des objets ou dans les base de données pour les représenter 
lisa_colors <- lisa_colors(lisa)
lisa_labels <- lisa_labels(lisa)
lisa_clusters <- lisa_clusters(lisa)
lisa_value <- lisa_values(lisa)
lisa_pvalue<- lisa_pvalues(lisa)


# Carte de Moran locaux
data_immo$moranlocalvalue <- lisa_values(lisa)

cols_v2 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272",  "#fb6a4a", "#de2d26")

mf_map(x = data_immo, 
       var = "moranlocalvalue", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       pal=cols_v2, 
       leg_title = "Local Moran", 
       leg_val_rnd = 1)
mf_title("Carte des LISA du prix médian du logement")

# Carte des p-value des moran locaux
data_immo$moranlocalpvalue<- lisa_pvalues(lisa)

library(dplyr)

# Pour plus de lisibilité dans la carte on va faire des classes des p-value
data_immo<- data_immo %>%  mutate(lisapvalue_fac = case_when(moranlocalpvalue <= 0.002 ~ "[0.001;0.002[",
                           moranlocalpvalue <= 0.01 ~ "[0.002;0.01[",
                           moranlocalpvalue <= 0.05 ~ "[0.01;0.05[",
                           TRUE ~ "[0.05;0.5]")) %>%
  mutate(lisapvalue_fac = factor(lisapvalue_fac,
                        levels = c("[0.001;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.5]")))

mypal <- mf_get_pal(n = 4, palette = "Reds")

mf_map(x = data_immo, 
       var = "lisapvalue_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       pal=mypal, 
       leg_title = "P-value Local Moran")
mf_title("Carte des significativité des LISA")

# Carte des clusters lisa avec rgeoda

plot(st_geometry(data_immo), 
     col=sapply(lisa_clusters, function(x){return(lisa_colors[[x+1]])}), 
     border = "#333333", lwd=0.2)
title(main = "Cluster LISA")
legend('bottomleft', legend = lisa_labels, fill = lisa_colors, border = "#eeeeee")

# En utilisant le package mapsf

data_immo$testmoran <- sapply(lisa_clusters, function(x){return(lisa_labels[[x+1]])})

colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")

mf_map(x = data_immo, 
       var = "testmoran", 
       type = "typo", 
       border = "black", 
       lwd = 0.1, 
       pal= colors,
       val_order = c("Not significant","Low-Low","Low-High","High-Low","High-High"),
       leg_title = "Lisa cluster")
mf_title("LISA clusters")

class(data_immo)

sp_data_immo <- as(data_immo, "Spatial")
class(sp_data_immo)

library(GWmodel)

# Pour construire la matrice de distances entre centroïdes des EPCI :
dm.calib <- gw.dist(dp.locat = coordinates(sp_data_immo))

# Définition de la bande passante (bandwidth en anglais)
# Cette fonction n'est pas compatible avec le format sf il faut un Spatial DataFrame (contrairement à spgwr qui peut travailler avec un format sf)
bw_g <- bw.gwr(data = sp_data_immo, 
              approach = "AICc", 
              kernel = "gaussian", 
              adaptive = TRUE,
              dMat = dm.calib,
              formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

mod.gwr_g <- gwr.robust(data = sp_data_immo, 
                   dMat = dm.calib,
                   bw = bw_g,
                   kernel = "gaussian",
                   filtered = FALSE, # un des problèmes de la GWR est de gérer des individus "aberrants" au niveau local. 2 méthodes ont été définies pour gérer cela. 
                                    # Méthode 1 (argument TRUE) on filtre en fonction des individus standardisés. L'objectif est de détecter les individus dont les résidus sont très élevés et de les exclure.
                                    # Methode 2 (argument FALSE) on diminue le poids des observations aux résidus élevés.
                   adaptive = TRUE,
                   formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

bw_tri <- bw.gwr(data = sp_data_immo, 
              approach = "AICc", 
              kernel = "tricube", 
              adaptive = TRUE, 
              dMat = dm.calib,
              formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

mod.gwr_tri <- gwr.robust(data = sp_data_immo, 
                   dMat = dm.calib,
                   bw = bw_tri,
                   kernel = "gaussian",
                   filtered = FALSE,
                   adaptive = TRUE,
                   formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

Best_gwr <- cbind(
  rbind(bw_g, bw_tri),
  rbind(mod.gwr_g$GW.diagnostic$gw.R2,mod.gwr_tri$GW.diagnostic$gw.R2),
  rbind(mod.gwr_g$GW.diagnostic$AIC,mod.gwr_tri$GW.diagnostic$AIC)) %>% 
  `colnames<-`(c("Nb Voisins","R2","AIC")) %>% 
  `rownames<-`(c("GAUSSIAN","TRICUBE"))

Best_gwr

# Pour voir les différents éléments qui composent notre modèle de GWR
summary(mod.gwr_g)

mod.gwr_g

# Pour visualiser ce fichier dans R
#View(mod.gwr_g$SDF@data)

#Pour voir à quoi il ressemble de manière interactive
datatable(mod.gwr_g$SDF@data)

# Pour voir les variables qui le constituent
names(mod.gwr_g$SDF@data)

res_gwr <- mod.gwr_g$SDF$Stud_residual
data_immo$res_gwr <- res_gwr

# carto résidus v2
mf_map(x = data_immo, 
       var = "res_gwr", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       breaks = quantile(data_immo$res_gwr, seq(0,1, by=1/11)), 
       pal = cols_v2, 
       leg_title = "Résidus GWR", 
       leg_val_rnd = 1)
mf_title("Résidus GWR")

# On ajoute à data_immo les coefficients

data_immo$agri.coef=mod.gwr_g$SDF$part_agri_nb_emploi
data_immo$perc_maison.coef <- mod.gwr_g$SDF$perc_maison
data_immo$dens_pop.coef=mod.gwr_g$SDF$dens_pop
data_immo$med_vie.coef=mod.gwr_g$SDF$med_niveau_vis
data_immo$logvac.coef=mod.gwr_g$SDF$perc_log_vac
data_immo$tinylog.coef=mod.gwr_g$SDF$perc_tiny_log
data_immo$suroccup.coef=mod.gwr_g$SDF$part_log_suroccup
data_immo$cadre.coef=mod.gwr_g$SDF$part_cadre_profintellec_nbemploi

# Réaliser la collection des cartes

par(mfrow = c(2, 4))

mf_map(x = data_immo, var = "agri.coef", type = "choro", pal= "Earth")
mf_title("Coefficients Agriculteurs")

mf_map(x = data_immo, var = "perc_maison.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Maison")

mf_map(x = data_immo, var = "dens_pop.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de dens pop")

mf_map(x = data_immo, var = "med_vie.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Médianne niveau de vie")

mf_map(x = data_immo, var = "logvac.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Logements vacants")

mf_map(x = data_immo, var = "tinylog.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Petits logements")

mf_map(x = data_immo, var = "suroccup.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de logement suroocupés")

mf_map(x = data_immo, var = "cadre.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Cadre")

data_immo$agri.t <- mod.gwr_g$SDF$part_agri_nb_emploi_TV
data_immo$maison.t <- mod.gwr_g$SDF$perc_maison_TV
data_immo$dens.t <- mod.gwr_g$SDF$dens_pop_TV
data_immo$medvie.t <- mod.gwr_g$SDF$med_niveau_vis_TV
data_immo$logvac.t <- mod.gwr_g$SDF$perc_log_vac_TV
data_immo$tinylog.t <- mod.gwr_g$SDF$perc_tiny_log_TV
data_immo$suroccup.t <- mod.gwr_g$SDF$part_log_suroccup_TV
data_immo$cadre.t <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi_TV     

#Définir contrib max
df<- as.data.frame(data_immo)
# On passe les t-values en valeurs absolues pour voir la plus grande contribution dans un sens sens ou dans l'autre
data_immo$contribmax<- colnames(df[, c(30:37)])[max.col(abs(df[, c(30:37)]),ties.method="first")]


par(mfrow = c(1, 1))

# Carte
mf_map(x = data_immo, var = "contribmax", type = "typo", pal= "Zissou 1")
mf_title("Carte des variables contribuant le plus par epci")

data_immo$r2local=mod.gwr_g$SDF$Local_R2

mf_map(x = data_immo, var = "r2local", type = "choro", pal= "Reds")
mf_title("R2 Locaux")

# Pour rappel si on a plus de 200 individus et le t-value > |1.96| on pourra considérer le coefficient comme significatif au seuil de 0.05 (95% chances que ce ne soit pas dû au hasard)

data_immo$nbsignif_t<- rowSums(abs(df[, c(30:37)]) > 1.96)

mf_map(x = data_immo, var = "nbsignif_t", type = "typo", pal= "Reds")
mf_title("Nombre des Betas significatif par EPCI (t-value)")

#  Les p-value ne sont pas fournit dans le modèle de la GWR, on pourrait les calculer à partir de t-value et de l'erreur standard mais le package GWmodel propose une fonction pour les obttenir

pvalue<-gwr.t.adjust(mod.gwr_g)

# On ajoute les p-value à notre fichier
data_immo$agri.p <- pvalue$SDF$part_agri_nb_emploi_p 
data_immo$maison.p <- pvalue$SDF$perc_maison_p
data_immo$dens.p <- pvalue$SDF$dens_pop_p
data_immo$medvie.p <- pvalue$SDF$med_niveau_vis_p
data_immo$logvac.p <- pvalue$SDF$perc_log_vac_p
data_immo$tinylog.p <- pvalue$SDF$perc_tiny_log_p
data_immo$suroccup.p <- pvalue$SDF$part_log_suroccup_p
data_immo$cadre.p <- pvalue$SDF$part_cadre_profintellec_nbemploi_p

df<- as.data.frame(data_immo)
data_immo$nbsignif_p<- rowSums(df[, c(41:48)] < 0.05)

mf_map(x = data_immo, var = "nbsignif_p", type = "typo", pal= "Reds")

# Ici nous représenterons les p-value avec un decoupage par classe de significativité et seulement les p-value de 2 VI

par(mfrow = c(1, 2))

# Par exemple les p-value des coefficients de la variable part de l'emploi agriculteur
data_immo<- data_immo %>%  mutate(agri.p_fac = case_when(agri.p <= 0.002 ~ "[0;0.002[",
                           agri.p <= 0.01 ~ "[0.002;0.01[",
                           agri.p <= 0.05 ~ "[0.01;0.05[",
                           agri.p <= 0.1 ~ "[0.05;0.1[",
                           TRUE ~ "[0.1;1]")) %>%
  mutate(agri.p_fac = factor(agri.p_fac,
                        levels = c("[0;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.1[", "[0.1;1]")))

mypal2 <- mf_get_pal(n = 5, palette = "SunsetDark")

mf_map(x = data_immo, 
       var = "agri.p_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       pal=mypal2, 
       leg_title = "Classe P-value")
mf_title("P-value du coefficient de la part d'emploi agriculteurs")

# Pour la densité de population
data_immo<- data_immo %>%  mutate(dens.p_fac = case_when(dens.p<= 0.002 ~ "[0;0.002[",
                           dens.p <= 0.01 ~ "[0.002;0.01[",
                           dens.p <= 0.05 ~ "[0.01;0.05[",
                           dens.p <= 0.1 ~ "[0.05;0.1[",
                           TRUE ~ "[0.1;1]")) %>%
  mutate(dens.p_fac = factor(dens.p_fac,
                        levels = c("[0;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.1[", "[0.1;1]")))

mf_map(x = data_immo, 
       var = "dens.p_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       pal=mypal2, 
       leg_title = "Classe P-value")
mf_title("P-value du coefficient de la densité de population")