library(dplyr)
library(sf)
library(mapsf)
library(ggplot2)
library(GWmodel)Analyse locale des relations entre géodiversité et biodiversité en Occitanie
Note
Les données et une partie des analyses présentées ci-dessous proviennent de l’article suivant :
- Gakou M., Souhil T., Barragan-Jason G., Feuillet T., Bétard F., Cauchoix M., Unveiling the spatial link between geodiversity and biodiversity: a multi-taxon study in the South of France. Soumis à Landscape ecology.
1 Géodiversité et biodiversité : notions et hypothèses
La géodiversité représente l’hétérogénéité des paysages en termes de caractéristiques hydrologiques, lithologiques, pédologiques et géomorphologique (voir ici une manière de la calculer en R). Comme elle illustre la variété des habitats potentiels au sein d’un territoire, on pose l’hypothèse qu’elle est positivement corrélée à la biodiversité taxonomique (en d’autres termes la richesse spécifique d’un groupe taxonomique donné). Néanmoins, cette relation positive est aussi susceptible de varier selon les contextes spatiaux, par exemple en fonction de la pression anthropique et donc de la proximité aux zones urbanisées, ou encore du fait d’être située ou non dans un espace protégé. On pose ainsi une seconde hypothèse selon laquelle cette relation entre géodiversité et biodiversité est spatialement non stationnaire.
2 Explorer la non-stationnarité spatiale : l’apport de la GWR
La GWR (geographically weighted regression) est une méthode de régression locale à coefficients spatialement variables. Pour chaque observation localisée, on estime un modèle linéaire sur la base du voisinage uniquement. Pour plus de détail, voir cette page.
L’équation de la GWR est la suivante :
\[ y_i = \beta_0(u_i, v_i) + \sum_{k=1}^{p} \beta_k(u_i, v_i) x_{ik} + \varepsilon_i \]
Où :
\((u_i, v_i)\) sont les coordonnées spatiales du point \(i\),
\(\beta_0(u_i, v_i)\) est l’ordonnée à l’origine spécifique à chaque point,
\(\beta_k(u_i, v_i)\) sont les coefficients estimés localement,
\(x_{ik}\) sont les variables explicatives,
\(\varepsilon_i\) est le terme d’erreur.
Cette méthode est parfaitement adaptée pour tester notre hypothèse de non-stationnarité spatiale des relations géodiversité-diversité et nou allons donc l’employer ci-dessous.
3 Les données mobilisées
Toutes les données de cette analyse ont été calculées au sein de mailles régulières de 10km de large, sur l’étendue de la région Occitanie (soit 595 mailles). Nous allons importer la table et décrire les variables une à une, mais activons dans un premier temps les packages nécessaires :
Importons maintenant la table, qui peut être téléchargée ici :
data <- st_read("data/gwrBota.gpkg", quiet = TRUE)3.1 Description des variables
La table contient cinq variables que nous allons décrire.
id_cell: Identifiant de la maillebiodiv: Variable dépendante. Il s’agit d’une estimation non-paramétrique de la richesse en espèces de plantes vasculaires dans chaque maille. Les données d’occurrence proviennent de l’INPN et de PatriNat, principalement à partir d’observations citoyennes et d’initiatives de science participative. Pour contrer le biais de sous- ou sur-obsevations de certaines espèces, nous avons estimé le Jackknife d’ordre 1, qui correspond précisément à notre variable dépendante (voir ici des précisions sur cet estimateur) : \[ \text{biodiv} = S_{\text{obs}} + \frac{n-1}{n} f_1 \]
\(\text{biodiv}\) : Estimation de la richesse spécifique par le Jackknife d’ordre 1,
\(S_{\text{obs}}\) : Nombre d’espèces observées dans la maille,
\(n\) : Nombre total d’observations dans la maille,
\(f_1\) : Nombre d’espèces observées une seule fois dans la maille (singletons).
geodiv: Prédicteur. Il s’agit de la géodiversité, estimée ainsi :
\[ \text{geodiv} = H_{\text{litho}} + H_{\text{pedo}} + H_{\text{hydro}} \]
\(H_{\text{litho}} = -\sum p_i \log p_i\) : Indice de Shannon de diversité lithologique,
\(H_{\text{pedo}} = -\sum p_i \log p_i\) : Indice de Shannon de diversité pédologique,
\(H_{\text{hydro}} = -\sum p_i \log p_i\) : Indice de Shannon de diversité hydrologique.
naturalite: Covariable décrivant le gradient de pression anthropique, estimée dans le cadre du projet CartNat.altiRange: Covariable décrivant l’ampltiude de l’altitude dans chaque maille.
3.2 Cartographie des variables
mf_theme(fg = "#593196", bg = "grey99", tab = FALSE, pos = "center")
par(mfrow = c(1, 2))
mf_map(data, col = NA, border = NA, expandBB = c(.15,0,0,0))
mf_shadow(data, add = TRUE)
mf_map(
x = data,
var = "biodiv",
type = "choro",
breaks = "q6",
pal = "Red-Yellow",
leg_pos = "bottomleft",
leg_horiz = TRUE,
leg_val_rnd = 0,
leg_box_cex = c(.75, 1),
leg_title = "Jackknife 1",
add = TRUE
)
mf_scale(size = 50, pos = "bottomright")
mf_title("Biodiversité (Trachéophytes) en Occitanie")
mf_map(data, col = NA, border = NA, expandBB = c(.15,0,0,0))
mf_shadow(data, add = TRUE)
mf_map(
x = data,
var = "geodiv",
type = "choro",
breaks = "q6",
pal = "Red-Yellow",
leg_horiz = TRUE,
leg_pos = "bottomleft",
leg_box_cex = c(.75, 1),
leg_title = "Shannon",
add = TRUE
)
mf_scale(size = 50, pos = "bottomright")
mf_title("Géoodiversité en Occitanie")
4 Régression globale
Commençons par estimer un modèle linéaire classique (global, sans considération de l’information spatiale) pour modéliser la relation entre nos deux variables d’intérêt. On transforme la variable Y en log car sa distribution est un peu tirée vers la droite :
library(ggplot2)
ggplot(data, aes(x = geodiv, y = log(biodiv))) +
geom_point()+
geom_smooth(method = "lm", col = "red", se = TRUE)
summary(lm(log(biodiv) ~ geodiv, data = data))
Call:
lm(formula = log(biodiv) ~ geodiv, data = data)
Residuals:
Min 1Q Median 3Q Max
-1.07140 -0.25191 -0.00624 0.28083 0.85336
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.32189 0.08584 73.648 < 2e-16 ***
geodiv 0.11263 0.01904 5.915 5.62e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3315 on 593 degrees of freedom
Multiple R-squared: 0.05571, Adjusted R-squared: 0.05412
F-statistic: 34.98 on 1 and 593 DF, p-value: 5.619e-09
Résultat
On constate que la relation est bien positive et significative, mais qu’elle n’est pas très forte en moyenne.
5 GWR et MGWR avec GWmodel
5.1 GWR univariée
Commençons par l’estimation d’une GWR univariée, c’est-à-dire avec seule la géodiversité comme prédicteur, sans covariable. Il faut d’abord construire la matrice de distance entre chaque observations (i.e. maille) :
matDist <- gw.dist(dp.locat = st_coordinates(st_centroid(data)))Il convient ensuite de déterminer la bandwidth (soit un nombre de voisins, soit une distance) optimale, c’est-à-dire celle qui minimise un critère d’information défini (comme l’AIC). Dans le cas d’un maillage régulier, une distance ou un nombre de plus proches voisins sont équivalents. Il faut également choisir une pondération du voisinage (gaussienne, bicarrée, etc., qui correspond au noyau). Nous allons sélectionner ici une pondération exponentielle (inversement proportionnelle à la distance). Notez que nous transformons Y en log.
# Exponential
bw <- bw.gwr(data = data,
approach = "AICc", # Critère d'optimisation
kernel = "exponential", # Nature du noyau
adaptive = TRUE, # Nombre de voisins et non distance
dMat = matDist,
formula = log(biodiv) ~ geodiv)Adaptive bandwidth (number of nearest neighbours): 375 AICc value: 249.8773
Adaptive bandwidth (number of nearest neighbours): 240 AICc value: 206.4548
Adaptive bandwidth (number of nearest neighbours): 155 AICc value: 150.6997
Adaptive bandwidth (number of nearest neighbours): 104 AICc value: 91.08254
Adaptive bandwidth (number of nearest neighbours): 71 AICc value: 31.83995
Adaptive bandwidth (number of nearest neighbours): 52 AICc value: -21.24602
Adaptive bandwidth (number of nearest neighbours): 38 AICc value: -68.90637
Adaptive bandwidth (number of nearest neighbours): 32 AICc value: -99.73669
Adaptive bandwidth (number of nearest neighbours): 25 AICc value: -131.986
Adaptive bandwidth (number of nearest neighbours): 24 AICc value: -137.9546
Adaptive bandwidth (number of nearest neighbours): 20 AICc value: -169.0905
Adaptive bandwidth (number of nearest neighbours): 21 AICc value: -161.0213
Adaptive bandwidth (number of nearest neighbours): 23 AICc value: -142.6887
Adaptive bandwidth (number of nearest neighbours): 21 AICc value: -161.0213
Adaptive bandwidth (number of nearest neighbours): 22 AICc value: -146.1872
Adaptive bandwidth (number of nearest neighbours): 21 AICc value: -161.0213
Adaptive bandwidth (number of nearest neighbours): 21 AICc value: -161.0213
Adaptive bandwidth (number of nearest neighbours): 20 AICc value: -169.0905 On peut maintenant estimer la GWR univariée :
gwrUniv <- gwr.basic(formula = log(biodiv) ~ geodiv,
data = data,
bw = bw,
adaptive = TRUE,
dMat = matDist,
kernel = "exponential")
print(gwrUniv) ***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2025-12-12 12:02:40.121893
Call:
gwr.basic(formula = log(biodiv) ~ geodiv, data = data, bw = bw,
kernel = "exponential", adaptive = TRUE, dMat = matDist)
Dependent (y) variable: biodiv
Independent variables: geodiv
Number of data points: 595
***********************************************************************
* Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3Q Max
-1.07140 -0.25191 -0.00624 0.28083 0.85336
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.32189 0.08584 73.648 < 2e-16 ***
geodiv 0.11263 0.01904 5.915 5.62e-09 ***
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3315 on 593 degrees of freedom
Multiple R-squared: 0.05571
Adjusted R-squared: 0.05412
F-statistic: 34.98 on 1 and 593 DF, p-value: 5.619e-09
***Extra Diagnostic information
Residual sum of squares: 65.1808
Sigma(hat): 0.3315374
AIC: 378.7559
AICc: 378.7966
BIC: -183.9127
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: exponential
Adaptive bandwidth: 20 (number of nearest neighbours)
Regression points: the same locations as observations are used.
Distance metric: A distance matrix is specified for this model calibration.
****************Summary of GWR coefficient estimates:******************
Min. 1st Qu. Median 3rd Qu. Max.
Intercept 5.707076 6.143518 6.334424 6.597268 7.0170
geodiv -0.050468 0.077951 0.108946 0.135415 0.2192
************************Diagnostic information*************************
Number of data points: 595
Effective number of parameters (2trace(S) - trace(S'S)): 59.94072
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 535.0593
AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): -169.0905
AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): -212.5251
BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): -611.9107
Residual sum of squares: 22.93096
R-square value: 0.6677931
Adjusted R-square value: 0.6305076
***********************************************************************
Program stops at: 2025-12-12 12:02:40.194095 On constate (i) que le modèle est de meilleure qualité (R² ~ 0,67) et (ii) que les coefficients associés à la géodiverité varient de -0,04 à 0,16. Il faut spatialiser ces coefficients pour bien visualiser les structures de non-stationnarité spatiale, en tenant compte également de la significativité locale des coefficients. Il y a différentes manières de cartographier ces betas locaux avec leur t-value, en voici une :
data <- data %>%
mutate(
geodiv_b = gwrUniv$SDF$geodiv,
geodiv_t = gwrUniv$SDF$geodiv_TV
)
dataSig <- data %>% filter(abs(geodiv_t) > 1.96)
mf_shadow(data)
mf_map(x = data,
var = "geodiv_b",
type = "choro",
border = "grey80",
pal = "Red-Yellow",
breaks = "q6",
leg_pos = "bottomright",
leg_adj = c(0, 3),
lwd = 0.3, add = TRUE)
mf_map(x = dataSig,
border = "black",
lwd = 0.8,
col = NA,
add = TRUE)
mf_legend(type = "symb", val = "t-value > |1,96|",
pos = "bottomright", title = "",
pal = "black", lwd = .8, pch = 0, cex = 1.9)
mf_scale(size = 50, pos = "bottomleft")
mf_title("Betas locaux (GWR) associés à la géodiversité (Y = biodiversité)")
5.2 MGWR
L’effet de la géodiversité est susceptible d’être affecté par les effets d’autres variables. Nous allons donc maintenant ajouter deux covariables dans le modèle (la naturalité et l’amplitude d’altitude) de manière à contrôler ces autres effets et à raisonner toutes choses égales par ailleurs.
Seulement, ces deux autres effets sont aussi potentiellement non-stationnaires. Mais pas forcément… La meilleure manière est d’estimer une MGWR (multiscale GWR, voir Fotheringham et al., (2017)). Cette méthode optimise une bandwidth différente pour chaque variable indépendante par le biais d’un algorithme de rétro-ajustement (qui consiste à optimiser chaque bw tout en gardant les autres fixes, jusqu’à convergence) :
\[ y_i = \beta_0(u_i, v_i) + \sum_{k=1}^{p} \beta_{bwk}(u_i, v_i) x_{ik} + \varepsilon_i \] où :
- \(y_i\) est la variable dépendante pour l’observation \(i\).
- \((u_i, v_i)\) représente les coordonnées spatiales de l’observation \(i\).
- \(\beta_0(u_i, v_i)\) est l’ordonnée à l’origine, qui peut varier spatialement.
- \(\beta_{bwk}(u_i, v_i)\) est le coefficient de la variable explicative \(x_{ik}\), estimé localement avec une bande passante spécifique \(bw_k\).
- \(x_{ik}\) est la valeur de la \(k\)-ième variable explicative pour l’observation \(i\).
- \(\varepsilon_i\) est le terme d’erreur.
Chaque coefficient \(\beta_{bw_j}(u_i, v_i)\) est estimé avec une bande passante \(bw_j\) qui lui est propre, ce qui permet d’examiner les relations à différentes échelles spatiales, et d’oberserver plus finement les variations spatiales des processus.
Il se trouve que si une des variables est stationnaire, alors sa bande passante \(bw_j\) sera très large, c’est-à-dire globale (un seul beta moyen sera estimé), et on se trouvera alors dans le cas d’une GWR mixte. La MGWR est aujourd’hui considérée comme l’approche à privilégier en GWR, quelles que soient le nombre de variables ou les hypothèses de stationnarité.
Voici comment l’estimer avec GWmodel :
mgwrBiodiv <- gwr.multiscale(formula = log(biodiv) ~ geodiv + naturalite + altiRange,
data = data,
adaptive = TRUE,
bws0 = rep(10,4), # Valeurs des bw d'initialisation pour chaque Xk
dMats = list(matDist,matDist,matDist,matDist), # Liste des matrices de distances
var.dMat.indx = 1:4, # Index des matrices de distances
predictor.centered = TRUE,
kernel = "exponential")print(mgwrBiodiv) ***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2025-12-12 12:02:40.590578
Call:
gwr.multiscale(formula = log(biodiv) ~ geodiv + naturalite +
altiRange, data = data, kernel = "exponential", adaptive = TRUE,
dMats = list(matDist, matDist, matDist, matDist), var.dMat.indx = 1:4,
bws0 = rep(10, 4), predictor.centered = TRUE)
Dependent (y) variable: biodiv
Independent variables: geodiv naturalite altiRange
Number of data points: 595
***********************************************************************
* Multiscale (PSDM) GWR *
***********************************************************************
*********************Model calibration information*********************
Kernel function: exponential
Adaptive bandwidths for each coefficient(number of nearest neighbours):
(Intercept) geodiv naturalite altiRange
Bandwidth 7 41 9 7
*************Summary of multiscale GWR coefficient estimates:***************
Min. 1st Qu. Median 3rd Qu. Max.
Intercept 5.945981 6.229072 6.355433 6.528865 6.8800
geodiv 0.052403 0.080460 0.097406 0.118965 0.1540
naturalite -0.194433 -0.072316 -0.017991 0.050153 0.1213
altiRange -0.179005 0.055473 0.109120 0.169049 0.3470
************************Diagnostic information*************************
Number of data points: 595
Effective number of parameters (2trace(S) - trace(S'S)): 185.7421
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 409.2579
AICc value: -346.8924
AIC value: -523.1202
BIC value: -492.1753
Residual sum of squares: 11.8967
R-square value: 0.8276494
Adjusted R-square value: 0.7492363
***********************************************************************
Program stops at: 2025-12-12 12:03:49.343552 Ajoutons les vecteurs des coefficients estimées à la table initiale :
data <- data %>%
mutate(
m_inter_b = mgwrBiodiv$SDF$Intercept,
m_inter_t = mgwrBiodiv$SDF$Intercept_TV,
m_geodiv_b = mgwrBiodiv$SDF$geodiv,
m_geodiv_t = mgwrBiodiv$SDF$geodiv_TV,
m_natur_b = mgwrBiodiv$SDF$naturalite,
m_natur_t = mgwrBiodiv$SDF$naturalite_TV,
m_alti_b = mgwrBiodiv$SDF$altiRange,
m_alti_t = mgwrBiodiv$SDF$altiRange_TV
)Puis cartographions ces coefficients :
par(mfrow = c(2, 2))
dataSigInter <- data %>% filter(abs(m_inter_t) > 1.96)
mf_map(data, col = NA, border = NA, expandBB = c(0,0,0,0.25))
mf_shadow(data, add = TRUE)
mf_map(x = data,
var = "m_inter_b",
type = "choro",
border = "grey80",
pal = "Red-Yellow",
breaks = "q6",
leg_title = paste("Bandwidth :\n",mgwrBiodiv$GW.arguments$bws[1],"voisins"),
leg_pos = "bottomright",
leg_adj = c(0, 3),
lwd = 0.3, add = TRUE)
mf_map(x = dataSigInter,
border = "black",
lwd = 0.8,
col = NA,
add = TRUE)
mf_legend(type = "symb", val = "t-value > |1,96|",
pos = "bottomright", title = "",
pal = "black", lwd = .8, pch = 0, cex = 1)
mf_scale(size = 50, pos = "bottomleft")
mf_title("Constantes locales (Y = biodiversité)")
dataSigGeod <- data %>% filter(abs(m_geodiv_t) > 1.96)
mf_map(data, col = NA, border = NA, expandBB = c(0,0,0,0.25))
mf_shadow(data, add = TRUE)
mf_map(x = data,
var = "m_geodiv_b",
type = "choro",
border = "grey80",
pal = "Red-Yellow",
breaks = "q6",
leg_title = paste("Bandwidth:\n",mgwrBiodiv$GW.arguments$bws[2],"voisins"),
lwd = 0.3,
leg_pos = "bottomright",
leg_adj = c(0, 3),
add = TRUE)
mf_map(x = dataSigInter,
border = "black",
lwd = 0.8,
col = NA,
add = TRUE)
mf_legend(type = "symb", val = "t-value > |1,96|",
pos = "bottomright", title = "",
pal = "black", lwd = .8, pch = 0, cex = 1)
mf_scale(size = 50, pos = "bottomleft")
mf_title("Betas locaux (GWR) associés à la géodiversité")
dataSigNatur <- data %>% filter(abs(m_natur_t) > 1.96)
mf_map(data, col = NA, border = NA, expandBB = c(0,0,0,0.25))
mf_shadow(data, add = TRUE)
mf_map(x = data,
var = "m_natur_b",
type = "choro",
border = "grey80",
pal = mf_get_pal(c(3,2), c("Teal", "YlOrRd"), neutral = "white"),
breaks = "q6",
leg_title = paste("Bandwidth:\n",mgwrBiodiv$GW.arguments$bws[3],"voisins"),
lwd = 0.3,
leg_pos = "bottomright",
leg_adj = c(0, 3),
add = TRUE)
mf_map(x = dataSigNatur,
border = "black",
lwd = 0.8,
col = NA,
add = TRUE)
mf_legend(type = "symb", val = "t-value > |1,96|",
pos = "bottomright", title = "",
pal = "black", lwd = .8, pch = 0, cex = 1)
mf_scale(size = 50, pos = "bottomleft")
mf_title("Betas locaux (GWR) associés à la naturalité")
dataSigAlti <- data %>% filter(abs(m_alti_t) > 1.96)
mf_map(data, col = NA, border = NA, expandBB = c(0,0,0,0.25))
mf_shadow(data, add = TRUE)
mf_map(x = data,
var = "m_alti_b",
type = "choro",
border = "grey80",
pal = c("lightblue", "grey85", "#E8CE8F", "#CE9758", "#AA592C", "#7D0112"),
breaks = "q6",
leg_title = paste("Bandwidth:\n",mgwrBiodiv$GW.arguments$bws[4],"voisins"),
lwd = 0.3,
leg_pos = "bottomright",
leg_adj = c(0, 3),
add = TRUE)
mf_map(x = dataSigAlti,
border = "black",
lwd = 0.8,
col = NA,
add = TRUE)
mf_legend(type = "symb", val = "t-value > |1,96|",
pos = "bottomright", title = "",
pal = "black", lwd = .8, pch = 0, cex = 1)
mf_scale(size = 50, pos = "bottomleft")
mf_title("Betas locaux (GWR) associés à la variation d'altitude")
Résultat
On constate que la non-stationnarité spatiale de la géodiversité présente des patrons différents en MGWR par rapport à la GWR univariée. C’est normal car elle prend maintenant en compte les effets partiels des autres variables (qui influencent en partie la relation). C’est ce résultat dont il faut tenir compte. Pour les deux autres co-variables, on constate qu’une partie des mailles seulement est significative. En certains lieux, ces variables n’ont pas d’effet significatif.
6 MGWR avec mgwrsar
Le package mgwrsar permet d’estimer des GWR avec autocorrélation spatiale (spatial autoregressive GWR). On peut ainsi coupler la prise en compte au sein d’un même modèle à la fois de l’hérogénéité spatiale des coefficients et de la dépendance spatiale de la variable endogène, avec la possibilité de considérer des échelles différentes entre ces deux phénomènes (Géniaux et Martinetti, 2018). Ce package propose également des méthodes de GWR avec sélection automatique des points focaux qui rend possible l’estimation de GWR sur très gros échantillons (Géniaux, 2024), qu’on abordera pas ici puisque notre échantillon est très petit. Enfin il introduit des algorithmes d’estimation plus performants de la multiscale GWR (Géniaux, 2025) et d’une variante de la MGWR qui permet une adaptation des bandes passantes à la fois entre variables mais également en fonction de la localisation des observations. C’est par ailleurs le seul package R qui propose des méthodes de prédiction adaptées à la Multiscale GWR pour faire de la prédiction hors échantillon.
Nous allons estimer dans un premier temps une GWR multivariée avec et sans autocorrelation spatiale à partir de ce package. Pour cela, on créé une matrice de coordonnées, la variable endogène logbiodiv dans le dataframe et une matrice d’interaction spatiale W de dimension n x n avec des 0 sur la diagonale et des poids positifs pour les 4 premiers voisins.
library(mgwrsar)
### on créé une matrice de coordonnées
coords=as.matrix(st_coordinates(st_centroid(data)))
### on créé un dataframe classique
mydata=as.data.frame(data)
## on créé la variable log(biodiv) et on standardise tous les termes
mydata$logbiodiv=log(mydata$biodiv)
mydata$logbiodiv=scale(mydata$logbiodiv)
mydata$geodiv=scale(mydata$geodiv)
mydata$naturalite=scale(mydata$naturalite)
mydata$altiRange=scale(mydata$altiRange)
myformula=as.formula('logbiodiv ~ geodiv + naturalite + altiRange')
### on créé une matrice d'interaction spatiale des 4 premiers voisins normalisés par ligne et avec une diagonale nulle.
W=kernel_matW(H=4,kernels='rectangle',coords=coords,NN=4,adaptive=TRUE,diagnull=TRUE)On estime alors une GWR multivariée standard qui nous sert ici d’étalon pour comparer aux estimations par GWR avec autocorrélation spatiale.
#### version GWR
res_GWR=golden_search_bandwidth(formula =myformula ,data =mydata,coords=coords,fixed_vars=NULL,kernels=c('bisq'), Model='GWR',control=list(W=W,adaptive=T,criterion='AICc',verbose=F),lower.bound=10, upper.bound=300,tolerance=0.001)
res_GWR$minimum # bandwidth = 47 [1] 47res_GWR$objective # AICc 918.8497[1] 918.8497model_GWR=res_GWR$model
summary(model_GWR)Call:
MGWRSAR(formula = myformula, data = mydata, coords = coords,
fixed_vars = fixed_vars, kernels = c("bisq"), H = c(res,
H2), Model = Model, control = list(W = W, adaptive = T,
criterion = "AICc", verbose = F))
Model: GWR
Kernels function: bisq
Kernels adaptive: YES
Kernels type: GD
Bandwidth: 47
Computation time: 0.307
Use of parallel computing: FALSE ncore = 1
Use of rough kernel: NO
Use of Target Points: NO
Number of data points: 595
Varying parameters: Intercept geodiv naturalite altiRange
Intercept geodiv naturalite altiRange
Min. -3.411803 -0.605674 -2.124409 -2.1929
1st Qu. -0.604011 0.022691 -0.303337 0.0690
Median -0.131958 0.181138 -0.088605 0.3808
Mean 0.025863 0.163826 -0.071417 0.4632
3rd Qu. 0.823034 0.292479 0.246920 0.6738
Max. 2.367201 0.582428 0.821547 4.1782
Effective degrees of freedom: 491.3263
AICc: 918.8497
Residual sum of squares: 106.7035
RMSE: 0.423478 Maintenant nous allons estimer des modèles par GWR avec autocorrélation spatiale en utilisant la matrice W et en considérant les 2 types de modèles suivant :
\(y=\lambda(u_i,v_i) Wy+\beta_v(u_i,v_i)X_v+\,\epsilon_i\;\;\;\;\;\;\;\; (MGWRSAR_1_0_kv)\)
\(y=\lambda Wy+\beta_v(u_i,v_i)X_v+\,\epsilon_i\;\;\;\;\;\;\;\; (MGWRSAR_0_0_kv)\)
Dans le modèle MGWRSAR_1_0_kv, tous les coefficients varient dans l’espace, y compris le coefficient d’autocorrélation spatiale lambda, alors que dans le modèle MGWRSAR_0_0_kv, le coefficient d’autocorrélation spatiale lambda est constant.
#### version GWR-SAR avec lambda variable dans l'espace
res_MGWRSAR_1_0_kv=golden_search_bandwidth(formula = myformula,data = mydata,coords=coords,fixed_vars=NULL,kernels=c('bisq'), Model='MGWRSAR_1_0_kv',control=list(W=W,adaptive=T,criterion='AICc',verbose=F),lower.bound=10, upper.bound=300,tolerance=0.001)
res_MGWRSAR_1_0_kv$minimum # bandwidth = 298 [1] 298res_MGWRSAR_1_0_kv$objective # AICc 977.0346[1] 1012.912#### version GWR-SAR avec lambda constant
res_MGWRSAR_0_0_kv=golden_search_bandwidth(formula = myformula,data = mydata,coords=coords,fixed_vars=NULL,kernels=c('bisq'), Model='MGWRSAR_0_0_kv',control=list(W=W,adaptive=T,criterion='AICc',verbose=F),lower.bound=10, upper.bound=300,tolerance=0.001)
res_MGWRSAR_0_0_kv$minimum # bandwidth = 79 [1] 50res_MGWRSAR_0_0_kv$objective # AICc 941.0762[1] 910.5948### rajout coefficient dans le data original
mycoef<-res_MGWRSAR_0_0_kv$model@Betav
colnames(mycoef)=paste0('Beta_GWRSAR_',colnames(mycoef))
data<-bind_cols(data,mycoef)On peut voir à la valeur des AICc obtenues qu’on privilégierait plutôt le modèle avec lambda constant (MGWRSAR_0_0_kv) entre les deux modèles avec autocorrélation spatiale, même si une confirmation formelle demanderait l’utilisation d’un test par bootstrap via la fonction mgwrsar_bootstrap_test().
Rajout des plots du modèle res_MGWRSAR_0_0_kv avec mapsf ?
On peut voir ici que l’utilisation d’une matrice d’interaction spatiale basées sur les 4 premiers voisins ne permet ni d’améliorer les critères d’informations de type AICc ni les RMSE. A priori c’est une mauvaise piste, même avec les autres matrices d’interaction spatiale que nous avons testées. (voir les vignettes Introduction with House Price Data et GWR and Mixed GWR with spatial autocorrelation pour plus d’informations sur les GWR avec autocorrélation spatiale).
Pour terminer, nous allons estimer les deux types d’estimateurs de type multiscale GWR que propose ce package via des approches dites Top-Down-Scale. Le premier algorithme intitulé tds_mwgr (Top-Down-Scale-MGWR) vise a estimer exactement le même modèle qu’avec la fonction gwr.multiscale du package GWmodel, mais en s’appuyant sur une méthode d’optimisation plus précise et plus rapide.
model_tds_mgwr=tds_mgwr(formula=myformula,Model='tds_mgwr',data=mydata,coords=coords,kernels='bisq',control_tds=list(nns=30,verbose=T,get_AIC=T),control=list(adaptive=T,doMC=T,ncore=10))
i= 1 Starting model
########################################################################
STAGE 1 : find unique bandwidth for each covariate
using Top Down Scale Approach with backfitting ...
########################################################################
1 Intercept
AICc : 1320.628 geodiv
AICc : 1299.708 naturalite
AICc : 1285.634 altiRange
AICc : 1266.494
AICc : 1266.494 stable 0 0 0 0 delta_rmse 0.2390377
H= 437 437 437 437
2 Intercept
AICc : 1201.922 geodiv
AICc : 1198.055 naturalite
AICc : 1192.792 altiRange
AICc : 1182.805
AICc : 1182.805 stable 0 0 0 0 delta_rmse 0.07482683
H= 321 321 321 321
3 Intercept
AICc : 1127.021 geodiv
AICc : 1124.889 naturalite
AICc : 1121.751 altiRange
AICc : 1115.051
AICc : 1115.051 stable 0 0 0 0 delta_rmse 0.06410837
H= 236 236 236 236
4 Intercept
AICc : 1065.338 geodiv
AICc : 1064.001 naturalite
AICc : 1058.618 altiRange
AICc : 1050.466
AICc : 1050.466 stable 0 0 0 0 delta_rmse 0.06472046
H= 173 173 173 173
5 Intercept
AICc : 1010.604 geodiv
AICc : 1009.978 naturalite
AICc : 1000.293 altiRange
AICc : 990.1034
AICc : 990.1034 stable 0 0 0 0 delta_rmse 0.06337666
H= 127 149 127 127
6 Intercept
AICc : 960.8788 geodiv
AICc : 961.0162 naturalite
AICc : 949.0141 altiRange
AICc : 937.5186
AICc : 937.5186 stable 0 1 0 0 delta_rmse 0.05856968
H= 94 149 94 94
7 Intercept
AICc : 912.4747 geodiv
AICc : 912.6234 naturalite
AICc : 900.8824 altiRange
AICc : 890.4559
AICc : 890.4559 stable 0 0 0 0 delta_rmse 0.05768064
H= 69 173 69 69
8 Intercept
AICc : 870.6943 geodiv
AICc : 871.1034 naturalite
AICc : 865.0742 altiRange
AICc : 856.1128
AICc : 856.1128 stable 0 0 0 0 delta_rmse 0.05922573
H= 50 202 50 50
9 Intercept
AICc : 845.0892 geodiv
AICc : 846.1411 naturalite
AICc : 849.4767 altiRange
AICc : 838.0593
AICc : 838.0593 stable 0 1 0 0 delta_rmse 0.03650563
H= 37 202 59 37
10 Intercept
AICc : 822.3523 geodiv
AICc : 823.3104 naturalite
AICc : 823.4576 altiRange
AICc : 813.3705
AICc : 813.3705 stable 0 2 0 0 delta_rmse 0.06371739
H= 27 202 69 27
11 Intercept
AICc : 801.0333 geodiv
AICc : 801.8813 naturalite
AICc : 801.0367 altiRange
AICc : 796.4467
AICc : 796.4467 stable 0 0 0 0 delta_rmse 0.05411151
H= 20 236 80 20
12 Intercept
AICc : 796.1566 geodiv
AICc : 796.5947 naturalite
AICc : 796.0906 altiRange
AICc : 795.4866
AICc : 795.4866 stable 1 0 0 1 delta_rmse -0.003829808
H= 20 275 94 20
13 Intercept
AICc : 794.4368 geodiv
AICc : 794.6749 naturalite
AICc : 795.6953 altiRange
AICc : 796.1856
AICc : 796.1856 stable 2 1 0 2 delta_rmse -0.001745192
H= 20 275 109 20
14 Intercept
AICc : 794.9717 geodiv
AICc : 795.114 naturalite
AICc : 796.0657 altiRange
AICc : 796.6083
AICc : 796.6083 stable 3 2 0 3 delta_rmse -0.001495699
H= 20 275 127 20
15 Intercept
AICc : 795.1931 geodiv
AICc : 795.439 naturalite
AICc : 797.4083 altiRange
AICc : 797.5682
AICc : 797.5682 stable 4 0 0 4 delta_rmse -0.00294258
H= 20 321 149 20
16 Intercept
AICc : 795.8264 geodiv
AICc : 796.0243 naturalite
AICc : 797.9424 altiRange
AICc : 798.3777
AICc : 798.3777 stable 5 1 0 5 delta_rmse -0.001587235
H= 20 321 173 20
17 Intercept
AICc : 796.6705 geodiv
AICc : 796.81 naturalite
AICc : 798.51 altiRange
AICc : 798.8331
AICc : 798.8331 stable 6 2 0 6 delta_rmse -0.00108856
H= 20 321 202 20
18 Intercept
AICc : 797.2142 geodiv
AICc : 797.3336 naturalite
AICc : 797.7367 altiRange
AICc : 798.354
AICc : 798.354 stable 7 3 1 7 delta_rmse 0.0004392464
H= 20 321 202 20
Time = 152.132 model_tds_mgwr@AICc # AICc = 798.354
798.354 summary(model_tds_mgwr)Call:
tds_mgwr(formula = myformula, data = mydata, coords = coords,
Model = "tds_mgwr", kernels = "bisq", control_tds = list(nns = 30,
verbose = T, get_AIC = T), control = list(adaptive = T,
doMC = T, ncore = 10))
Model: tds_mgwr
Method for spatial autocorrelation:
Kernels function: bisq
Kernels adaptive: YES
Kernels type: GD
Bandwidth: 20 321 202 20
Computation time: 152.133
Use of parallel computing: FALSE ncore = 1
Use of rough kernel: NO
Use of Target Points: NO
Number of data points: 595
Varying parameters: Intercept geodiv naturalite altiRange
Intercept geodiv naturalite altiRange
Min. -1.180314 0.109984 -0.356723 -1.1976
1st Qu. -0.501681 0.154186 -0.303352 0.1538
Median -0.189269 0.179235 -0.206409 0.3721
Mean 0.050884 0.194088 -0.140626 0.4205
3rd Qu. 0.692781 0.243694 0.072018 0.7022
Max. 1.528831 0.289369 0.162427 2.1928
Effective degrees of freedom: 484.3587
AICc: 798.354
Residual sum of squares:
RMSE: 0.3760669 ## les différentes bandes passantes candidates utilisées
model_tds_mgwr@V # 595 510 437 375 321 275 236 202 173 149 127 109 94 80 69 59 50 43 37 32 27 23 20 17 15 13 11 9 8 7 6 [1] 595 510 437 375 321 275 236 202 173 149 127 109 94 80 69 59 50 43 37
[20] 32 27 23 20 17 15 13 11 9 8 7 6## les bandes passantes sélectionnées
model_tds_mgwr@H # Intercept 20 geodiv 321 naturalite 202 altiRange 20 Intercept geodiv naturalite altiRange
20 321 202 20 On peut voir que l’algorithme tds_mgwr est 11 fois plus rapide que gwr.multiscale (GWmodel) même sur petit échantillon et permet d’atteindre un AICc notoirement plus bas (798.35 contre 933.78). Le vecteur de bande passante obtenu est très différent puisque les échelles des variables geodiv et naturalite sont au final beaucoup plus larges (321 et 202 avec tds_mgwr contre 41 et 9 avec gwr.multiscale ). Comme le montre geniaux (2025) sur un jeux de simulation très variés, l’estimateur tds_mgwr est plus robuste à la présence de pattern spatiaux complexes et fournit des estimations plus fiables que l’algorithme initial proposé par Fotheringam et al 2017 implémenté dans le package R GWmodel et dans la librairie python mgwr.
On peut voir que le coefficient geodiv agit à une échelle plus globale et on n’observe plus d’inversion de signes dans les effets de cette variable dont l’effet reste positif et varie moins (0.11-0.29) contre (-0.28 - 0.270) pour la GWR et (-0.16 - 0.24) pour la gwr.multiscale (optimum local).
7 atds_mgwr
Le deuxième algorithme intitulé atds_mgwr (Adaptive-Top-Down-Scale-MGWR) combine un algorithme de rétroajustement et un algorithme de gradient boosting qui permet de corriger séquentiellement les estimations en considérant toutes une série de bandes passantes et pas seulement une seule par variable. Cela permet de corriger localement les bandes passantes et de les adapter à des patterns spatiaux complexes. Pour que cette correction de deuxième étape puisse fournir effectivement de meilleures estimations, il faut naturellement être en présence de patterns spatiaux complexes et avoir des échantillons de taille suffisantes (>1000 obs pour un modèle multivarié). Pour cette raison, la deuxième étape atds_mgwr ne produit pas ici de meilleurs résultats et conserve, à un bémol près, l’estimation de première étape fournit par tds_mgwr.
model_atds_mgwr=tds_mgwr(formula=myformula,Model='atds_mgwr',data=mydata,coords=coords,kernels='bisq',control_tds=list(nns=30,verbose=F,get_AIC=T,model_stage1=model_tds_mgwr),control=list(adaptive=T,doMC=T,ncore=10))
model_atds_mgwr@AICc # AICc = 798.5862[1] 798.5862summary(model_atds_mgwr)Call:
tds_mgwr(formula = myformula, data = mydata, coords = coords,
Model = "atds_mgwr", kernels = "bisq", control_tds = list(nns = 30,
verbose = F, get_AIC = T, model_stage1 = model_tds_mgwr),
control = list(adaptive = T, doMC = T, ncore = 10))
Model: atds_mgwr
Method for spatial autocorrelation:
Kernels function: bisq
Kernels adaptive: YES
Kernels type: GD
Bandwidth: 20 437 321 20
Computation time: 403.293
Use of parallel computing: FALSE ncore = 1
Use of rough kernel: NO
Use of Target Points: NO
Number of data points: 595
Varying parameters: Intercept geodiv naturalite altiRange
Intercept geodiv naturalite altiRange
Min. -1.180314 0.111680 -0.352654 -1.1976
1st Qu. -0.501681 0.146550 -0.303614 0.1538
Median -0.189269 0.178127 -0.190293 0.3721
Mean 0.050884 0.194456 -0.128943 0.4205
3rd Qu. 0.692781 0.252398 0.070352 0.7022
Max. 1.528831 0.293660 0.191956 2.1928
Effective degrees of freedom: 484.8096
AICc: 798.5862
Residual sum of squares:
RMSE: 0.3760037 ### rajout coefficient dans le data original de l'estimation tds_mgwr
mycoef<-model_atds_mgwr@Betav
colnames(mycoef)=paste0('Beta_tds_mgwr_',colnames(mycoef))
data<-bind_cols(data,mycoef)PLOT des coefficient tds_mgwr avec mapsf ?
apply(model_GWR@Betav,2,range) Intercept geodiv naturalite altiRange
[1,] -3.411803 -0.6056741 -2.124409 -2.192854
[2,] 2.367201 0.5824284 0.821547 4.178232apply(model_GWR@Betav,2,range) Intercept geodiv naturalite altiRange
[1,] -3.411803 -0.6056741 -2.124409 -2.192854
[2,] 2.367201 0.5824284 0.821547 4.178232