load("data/Lyon.Rdata")
library(spdep)
## Matrice de contiguïté selon le partage d'un nœud (Queen)
<- poly2nb(LyonIris, queen = TRUE)
nb_queen <- nb2listw(nb_queen, zero.policy = TRUE, style = "W")
w_queen
# I de Moran sur la variable lden elon l’hypothèse de la normalité
moran.test(LyonIris$Lden, listw=w_queen, zero.policy = TRUE, randomisation = FALSE)
# I de Moran sur la variable lden selon la normalité selon l’hypothèse de la randomisation
moran.test(LyonIris$Lden, listw=w_queen, zero.policy = TRUE, randomisation = TRUE)
# I de Moran sur la variable lden selon des permutations Monte-Carlo
moran.mc(LyonIris$Lden, listw=w_queen, zero.policy = TRUE, nsim = 999)
11 Correction des exercices
11.1 Exercices du chapitre 1
11.1.1 Exercice 1
11.1.2 Exercice 2
library(spdep)
load("data/Lyon.Rdata")
# Modèle MCO
<- "Lden ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed"
formule <- lm(formule, data = LyonIris)
modele_mco
# Résultats du modèle
summary(modele_mco)
## Matrice de contiguïté (Rook)
<- poly2nb(LyonIris, queen = FALSE)
nb_rook <- nb2listw(nb_rook, zero.policy = TRUE, style = "W")
w_rook
# I de Moran sur les résidus du modèle global (MCO)
lm.morantest(modele_mco, w_rook)
# Cartographie des résidus
$MCO.Residus <- modele_mco$residuals
LyonIristmap_mode("plot")
<- list(text.separator = "à",
legende_parametres text.less.than = "Moins de",
text.or.more = "et plus",
decimal.mark = ",",
big.mark = " ")
tm_shape(LyonIris)+
tm_polygons(
col="gray25", lwd=.5,
fill = "MCO.Residus",
fill.scale = tm_scale_intervals(
n = 6,
midpoint = 0,
style = "pretty",
values = "-brewer.rd_bu",
label.format = legende_parametres
),fill.legend = tm_legend(
title = "MCO",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center")
)+
) tm_layout(frame = FALSE)
11.2 Exercices du chapitre 2
11.2.1 Exercice 1
library(sf, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(tmap, quietly = TRUE)
library(dplyr, quietly = TRUE)
library(ggpubr, quietly = TRUE)
library(cowplot, quietly = TRUE)
library(spdep, quietly = TRUE)
library(DHARMa, quietly = TRUE)
library(nlme, quietly = TRUE)
library(car, quietly = TRUE)
library(sf, quietly = TRUE)
load("data/Lyon.Rdata")
<- st_centroid(LyonIris)
LyonIris <- st_coordinates(LyonIris)
XY $X <- XY[,1]
LyonIris$Y <- XY[,2]
LyonIris
## Construction du modèle MCO
<- lm(NO2 ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
modele_MCO data = LyonIris)
summary(modele_MCO)
## Autocorrélation spatiale des résidus du modèle MCO
$residual_lm <- residuals(modele_MCO)
LyonIris<- knearneigh(LyonIris, k = 10)
knn_list <- knn2nb(knn_list, sym = TRUE)
nb_list <- nb2listwdist(nb_list, LyonIris, style = 'W', type = 'idw', alpha = 2)
listw moran.mc(LyonIris$residual_lm, listw, nsim = 999)
# Exponentielle
<- gls(NO2 ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
modelExp data = LyonIris,
correlation = corExp(form = ~X+Y, nugget = TRUE))
# Gaussienne
<- gls(NO2 ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
modelGau data = LyonIris,
correlation = corGaus(form = ~X+Y, nugget = TRUE))
# Sphérique
<- gls(NO2 ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
modelSph data = LyonIris,
correlation = corSpher(form=~X+Y, nugget = TRUE))
# Rationelle quadratique
<- gls(NO2 ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
modelRatio data = LyonIris,
correlation = corRatio(form = ~X+Y, nugget = TRUE))
# Comparaison des modèles avec l'AIC
AIC(modelExp, modelGau, modelSph, modelRatio)
<- modelExp$modelStruct$corStruct
struc print(struc)
## Structure exponentielle
<- function(n, r, d) {
corExpo return((1 - n) * exp(-r / d))
}
<- seq(0, 10000, 50)
distances
<- data.frame(
df_corr distance = distances,
correlation = corExpo(n = 0.01718,
d = 3279.694,
r = distances
)
)
ggplot(df_corr) +
geom_line(aes(x = distance, y = correlation)) +
theme_bw() +
labs(title = 'Structure de corrélation du modèle GLS (exponentielle)')
# Résidus du modèle GLS exponentielle
$residual_gls <- residuals(modelExp, type = "normalized")
LyonIrismoran.mc(LyonIris$residual_gls, listw, nsim = 999)
qqPlot(LyonIris$residual_gls)
# subset(LyonIris, LyonIris$residual_gls > 2)
# Comparaison des modèles MCO et GLS
library(knitr, quietly = TRUE)
library(kableExtra, quietly = TRUE)
<- summary(modele_MCO)
s1 <- summary(modelGau)
s2
<- data.frame(s1$coefficients[,c(1,2,4)])
df1 <- data.frame(s2$tTable[,c(1,2,4)])
df2
kable(cbind(df1, df2), digits = 3,
col.names = c('coeff.', 'err. std.', 'val. p',
'coeff.', 'err. std.', 'val. p')) %>%
kable_classic() %>%
add_header_above(c(' ' = 1, "Modèle linéaire" = 3, "Modèle GLS (exp)" = 3))
11.2.2 Exercice 2
load("data/Lyon.Rdata")
library(spdep)
library(spatialreg)
library(tmap)
## Construction du modèle MCO
<- lm(NO2 ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
modele_MCO data = LyonIris)
# Construire un matrice de contiguïté (Queen)
<- poly2nb(LyonIris, queen = TRUE)
nbs <- nb2listw(nbs, style = 'B', zero.policy = TRUE)
listw
# I de Moran sur les résidus du modèles MCO
moran.mc(residuals(modele_MCO), listw = listw, nsim = 999)
# Construire le modèle CAR
<- spautolm(
modele_car ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
NO2 data = LyonIris,
listw = listw,
family = 'CAR',
zero.policy = TRUE
)
summary(modele_car)
# I de Moran sur les résidus du modèle CAR
moran.mc(residuals(modele_car), listw = listw, nsim = 999)
# Cartographie de l'effet spatial
$car_term <- modele_car$fit$signal_stochastic
LyonIris
<- list(text.separator = "à",
legende_parametres decimal.mark = ",",
big.mark = " ")
tm_shape(LyonIris) +
tm_polygons(col = NULL,
fill ="car_term",
fill.scale = tm_scale_intervals(
n = 7, style = "fisher", midpoint = 1,
label.format = legende_parametres,
values = "-brewer.rd_bu"),
fill.legend = tm_legend(
title = "Terme spatial \ndu modèle CAR",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center"))
+
) tm_layout(frame=FALSE) +
tm_scalebar(breaks = c(0,10),
position = tm_pos_out("right", "bottom"))
11.2.3 Exercice 3
library(ggplot2, quietly = TRUE)
library(tmap, quietly = TRUE)
library(dplyr, quietly = TRUE)
library(spaMM, quietly = TRUE)
library(spdep, quietly = TRUE)
library(DHARMa, quietly = TRUE)
load("data/Lyon.Rdata")
## Construction du modèle GLM gaussien
<- glm(NO2 ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
modele_glm family = gaussian,
data = LyonIris)
# Construire un matrice de contiguïté (Queen)
<- poly2nb(LyonIris, queen = TRUE)
nbs <- nb2listw(nbs, style = 'B', zero.policy = TRUE)
listw
# I de Moran sur les résidus du modèles MCO
moran.mc(residuals(modele_glm), listw = listw, nsim = 999)
# Construire le modèle GLMM de type CAR
$oid <- 1:nrow(LyonIris)
LyonIris<- poly2nb(LyonIris, queen = TRUE)
queen_nb <- matrix(0, ncol = length(queen_nb), nrow = length(queen_nb))
wmat for(i in 1:length(queen_nb)){
<- 1
wmat[i, queen_nb[[i]]]
}
<- fitme(
model_car ~ Pct0_14 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed +
NO2 adjacency(1|oid),
data = st_drop_geometry(LyonIris),
adjMatrix = wmat,
family = gaussian()
)
summary(modele_car)
# I de Moran sur les résidus du modèle CAR
moran.mc(residuals(modele_car), listw = listw, nsim = 999)
# Cartographie de l'effet spatial
<- LyonIris
pred_df <- c('Pct0_14', 'Pct_65', 'Pct_Img', 'Pct_brevet')
vars for(var in vars){
<- 0
pred_df[[var]]
}
$pred <- as.numeric(predict(model_car,
pred_dfnewdata = pred_df,
type = 'link'))
$pred <- pred_df$pred - mean(pred_df$pred)
pred_df
tm_shape(pred_df) +
tm_polygons(col = NULL,
fill ="pred",
fill.scale = tm_scale_intervals(
n = 7, style = "fisher",
label.format = legende_parametres,
values = "-brewer.rd_bu"),
fill.legend = tm_legend(
title = "Effet aléatoire \nspatial",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center"))
+
) tm_layout(frame=FALSE)
11.3 Exercices du chapitre 3
11.3.1 Exercice 1
library(sf)
library(spatialreg)
load("data/Lyon.Rdata")
# Matrice de contiguïté selon le partage d'un segment
<- poly2nb(LyonIris, queen=FALSE)
nb_rook <- nb2listw(nb_rook, zero.policy = TRUE, style = "W")
w_rook
# Modèles
<- "PM25 ~ Pct0_14+Pct_65+Pct_Img+Pct_brevet+NivVieMed"
formule
# Modèles MOC, SLX, SDM, SDEM, SAR, SEM, Manski
<- lm(formule, data = LyonIris)
MCO <- lmSLX(formule, listw = w_rook, data = LyonIris)
SLX <- sacsarlm(formule, listw=w_rook, data = LyonIris, type="sacmixed")
Manski <- lagsarlm(formule, listw = w_rook, data = LyonIris, type = "mixed")
SDM <- errorsarlm(formule, listw=w_rook, data = LyonIris, etype = 'emixed')
SDEM <- lagsarlm(formule,listw=w_rook, data = LyonIris, type = 'lag')
SAR <- errorsarlm(formule, listw=w_rook, data = LyonIris) SEM
11.3.2 Exercice 2
summary(lm.LMtests(model = MCO,
listw = w_rook,
test = c("LMlag","LMerr","RLMlag","RLMerr")))
11.3.3 Exercice 3
## Valeurs d'AIC et de BIC
<- AIC(MCO, SLX, SDM, SDEM, SAR, SEM)
AICs <- BIC(MCO, SLX, SDM, SDEM, SAR, SEM)
BICs
## Autocorrélation spatiale des résidus
<- moran.mc(resid(MCO), w_rook, nsim = 999)
IMoran_MCO <- moran.mc(resid(SLX), w_rook, nsim = 999)
IMoran_SLX <- moran.mc(resid(SAR), w_rook, nsim = 999)
IMoran_SAR <- moran.mc(resid(SEM), w_rook, nsim = 999)
IMoran_SEM <- moran.mc(resid(SDM), w_rook, nsim = 999)
IMoran_SDM <- moran.mc(resid(SDEM), w_rook, nsim = 999)
IMoran_SDEM
<- c(IMoran_MCO$statistic, IMoran_SLX$statistic,
MoranI_s $statistic, IMoran_SEM$statistic,
IMoran_SAR$statistic, IMoran_SDEM$statistic)
IMoran_SDM
<- c(IMoran_MCO$p.value, IMoran_SLX$p.value,
MoranI_p $p.value, IMoran_SEM$p.value,
IMoran_SAR$p.value, IMoran_SDEM$p.value)
IMoran_SDM
## Tableau
<- data.frame(Modele = c("MCO", "SLX", "SAR", "SEM", "SDM", "SDEM"),
Comparaison AIC = round(AICs$AIC,2),
BIC = round(BICs$BIC,2),
dl = AICs$df,
MoranI = round(MoranI_s,3),
MoranIp = round(MoranI_p,3))
Comparaison
11.4 Exercices du chapitre 4
11.4.1 Exercice 1
11.4.2 Exercice 2
11.4.3 Exercice 3
11.5 Exercices du chapitre 5
11.5.1 Exercice 1
11.5.2 Exercice 2
11.5.3 Exercice 3
11.6 Exercices du chapitre 6
11.6.1 Exercice 1
library(sf)
library(mgcv)
library(car)
library(DHARMa)
library(mgcViz)
# Chargement du jeu de données sur l'agglomération Lyonnaise
load("data/Lyon.Rdata")
# Vérification de la multicolinéarité (VIF)
vif(lm(PM25 ~ Lden+Pct_65+Pct_Img+Pct_brevet+NivVieMed, data = LyonIris))
# Construction du modèle GLM gaussien
<- gam(PM25 ~ Pct0_14+Pct_65+Pct_Img+Pct_brevet+NivVieMed,
modele_gam1 family = gaussian, data = LyonIris)
# Résidus simulés du modèle gaussien
<- simulateResiduals(modele_gam1, plot = FALSE)
gam1_res plot(gam1_res)
# Construction du modèle GLM avec une distribution de Student
<- gam(PM25 ~ Pct0_14+Pct_65+Pct_Img+Pct_brevet+NivVieMed,
modele_gam2 family = scat, data = LyonIris)
# Résidus simulés avec une distribution de Student
<- simulateResiduals(modele_gam2, plot = FALSE)
gam2_res plot(gam2_res)
# Comparaison des deux modèles
AIC(modele_gam1, modele_gam2)
# Résultats du modèle gaussien
summary(modele_gam2)
# Matrice de contiguïté selon le partage d'un nœud
<- poly2nb(LyonIris, queen = TRUE)
queen_nb <- nb2listw(queen_nb, style = 'W', zero.policy = TRUE)
queen_w
# I de Moran sur les résidus gaussiens
moran.mc(residuals(modele_gam2,
type = 'pearson'),
listw = queen_w, nsim = 999,
zero.policy = TRUE)
# I de Moran sur les résidus simulés
moran.mc(residuals(gam2_res,
type = 'pearson'),
listw = queen_w, nsim = 999,
zero.policy = TRUE)
11.6.2 Exercice 2
# Chargement du jeu de données sur l'agglomération Lyonnaise
load("data/Lyon.Rdata")
# Ajout des coordonnées x et y dans la couche sf LyonIris
<- st_coordinates(st_centroid(LyonIris))
XY $X <- XY[,1]
LyonIris$Y <- XY[,2]
LyonIris
# Modèle GAM avec la distribution de student et une spline
# sur les coordonnées x et y avec la distribution de Student
<- gam(NO2 ~ Pct0_14+Pct_65+Pct_Img+Pct_brevet+NivVieMed+
modele_gam_spline s(X,Y, bs = 'gp', m = 3),
data = LyonIris, family = scat)
summary(modele_gam_spline)
11.7 Exercices du chapitre 7
11.7.1 Exercice 1
11.7.2 Exercice 2
11.7.3 Exercice 3
11.8 Exercices du chapitre 8
11.8.1 Exercice 1
library(sf)
library(spgwr)
load("data/Lyon.Rdata")
# Ajout des coordonnées x et y
<- st_coordinates(st_centroid(LyonIris))
xy $X <- xy[,1]
LyonIris$Y <- xy[,2]
LyonIris
# Optimisation du nombre de voisins avec le CV
<- "PM25 ~ Pct0_14+Pct_65+Pct_Img+Pct_brevet+NivVieMed"
formule <- gwr.sel(formule,
bwaCV_voisins data = LyonIris,
method = "cv",
gweight=gwr.bisquare,
adapt=TRUE,
verbose = FALSE,
RMSE = TRUE,
longlat = FALSE,
coords=cbind(LyonIris$X,LyonIris$Y))
# Optimisation du nombre de voisins avec l'AIC
<- "PM25 ~ Pct0_14+Pct_65+Pct_Img+Pct_brevet+NivVieMed"
formule <- gwr.sel(formule,
bwaCV_voisins data = LyonIris,
method = "CV",
gweight=gwr.bisquare,
adapt=TRUE,
verbose = FALSE,
RMSE = TRUE,
longlat = FALSE,
coords=cbind(LyonIris$X,LyonIris$Y))
# Réalisation de la GWR
<- gwr(formule,
modele_gwr data = LyonIris,
adapt=bwaCV_voisins,
gweight=gwr.bisquare,
hatmatrix=TRUE,
se.fit=TRUE,
coords=cbind(LyonIris$X,LyonIris$Y),
longlat=FALSE)
# Affichage des résultats
modele_gwr
11.8.2 Exercice 2
# Modèle MCO
<- lm(formule, data = LyonIris)
modele_global
# Comparaison des R2
<- summary(modele_global)$r.squared
r2_global <- sum((modele_gwr$lm$y - modele_gwr$SDF$pred)^2)
rss <- sum((modele_gwr$lm$y - mean(modele_gwr$SDF$pred))^2)
tss <- 1 - (rss / tss)
R2_GWRquasiglobal cat("R2 global (LM) = ", round(r2_global, 3),
"\nR2 quasi-global (GWR) : ", round(R2_GWRquasiglobal, 3))
# Comparaison des R2
anova(modele_gwr)
LMZ.F1GWR.test(modele_gwr)
LMZ.F2GWR.test(modele_gwr)
# Les coefficients du modèle GWR varient-ils significativement dans l'espace?
LMZ.F3GWR.test(modele_gwr)
11.8.3 Exercice 3
library(tmap)
# Récupération du R carré locaux et valeurs de t locales
$GWR.R2 <- modele_gwr$SDF$localR2
LyonIris<- c("Pct0_14", "Pct_65", "Pct_Img", "Pct_brevet", "NivVieMed")
vars_indep for(e in vars_indep){
# Nom des nouvelles variables
<- paste0("GWR.", "B_", e)
var.coef <- paste0("GWR.", "T_", e)
var.t # Récupération des coefficients pour les variables indépendantes
<- modele_gwr$SDF[[e]]
LyonIris[[var.coef]] # Calcul des valeurs de t pour les variables indépendantes
<- modele_gwr$SDF[[e]] / modele_gwr$SDF[[paste0(e, "_se")]]
LyonIris[[var.t]]
}
# Cartographie des R2 locaux
<- list(text.separator = "à",
legendes_parametres decimal.mark = ",",
big.mark = " ")
tm_shape(LyonIris) +
tm_polygons(
col = "gray25", lwd = .2,
fill = "GWR.R2",
fill.scale = tm_scale_intervals(
style = "quantile", n = 5,
values = "brewer.yl_or_br",
label.format = legende_parametres),
fill.legend = tm_legend(
title = "R2 locaux",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center"))
+
) tm_layout(frame = FALSE) +
tm_scale(breaks = c(0, 5))
# Cartographie des valeurs de t locales
## Paramètres pour la légende
<- list(text.separator = "à",
legende.parametres text.less.than = "Moins de",
text.or.more = "et plus",
decimal.mark = ",",
big.mark = " ")
## Cartographie
= c(-Inf, -3.29, -2.58, -1.96, 1.96, 2.58, 3.29, Inf)
classes_intervalles
<- tm_shape(LyonIris) +
carte1 tm_polygons(
col = "gray25", lwd = .2,
fill = "GWR.T_Pct0_14",
fill.scale = tm_scale_intervals(
breaks = classes_intervalles,
midpoint = 0,
values = "-brewer.rd_bu",
label.format = legende_parametres),
fill.legend = tm_legend(
title = "Moins de 15 ans (%)",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center"))
+
) tm_layout(frame = FALSE)
<- tm_shape(LyonIris) +
carte2 tm_polygons(
col = "gray25", lwd = .2,
fill = "GWR.T_Pct_65",
fill.scale = tm_scale_intervals(
breaks = classes_intervalles,
midpoint = 0,
values = "-brewer.rd_bu",
label.format = legende_parametres),
fill.legend = tm_legend(
title = "65 ans et plus (%)",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center"))
+
) tm_layout(frame = FALSE)
<- tm_shape(LyonIris) +
carte3 tm_polygons(
col = "gray25", lwd = .2,
fill = "GWR.T_Pct_Img",
fill.scale = tm_scale_intervals(
breaks = classes_intervalles,
midpoint = 0,
values = "-brewer.rd_bu",
label.format = legende_parametres),
fill.legend = tm_legend(
title = "Immigrants (%)",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center"))
+
) tm_layout(frame = FALSE)
<- tm_shape(LyonIris) +
carte4 tm_polygons(
col = "gray25", lwd = .2,
fill = "GWR.B_Pct_brevet",
fill.scale = tm_scale_intervals(
breaks = classes_intervalles,
midpoint = 0,
values = "-brewer.rd_bu",
label.format = legende_parametres),
fill.legend = tm_legend(
title = "Faible scolarité (%)",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center"))
+
) tm_layout(frame = FALSE)
<- tm_shape(LyonIris) +
carte5 tm_polygons(
col = "gray25", lwd = .2,
fill = "GWR.T_NivVieMed",
fill.scale = tm_scale_intervals(
breaks = classes_intervalles,
midpoint = 0,
values = "-brewer.rd_bu",
label.format = legende_parametres),
fill.legend = tm_legend(
title = "Niveau de vie (€1000)",
frame = FALSE,
position = tm_pos_out("right", pos.v = "center"))
+
) tm_layout(frame = FALSE)
tmap_arrange(carte1, carte2, carte3, carte4, carte5, ncol = 2, nrow=3)
11.9 Exercices du chapitre 9
11.9.1 Exercice 1
11.9.2 Exercice 2
11.9.3 Exercice 3
11.10 Exercices du chapitre 10
11.10.1 Exercice 1
library(sf, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(tmap, quietly = TRUE)
library(dplyr, quietly = TRUE)
library(ggpubr, quietly = TRUE)
library(mgcv, quietly = TRUE)
library(spdep, quietly = TRUE)
library(collapse, quietly = TRUE)
library(gratia, quietly = TRUE)
load("data/Lyon.Rdata")
<- st_coordinates(st_centroid(LyonIris))
xy $X <- xy[,1]
LyonIris$Y <- xy[,2]
LyonIris
# Modèle non spatial avec une distribution gaussienne
<- gam(NO2 ~
model1 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed,
Pct0_14 data = LyonIris, method = 'REML',
family = gaussian(link = "identity"))
summary(model1)
# Modèle intégrant uniquement une constante variant spatialement
<- gam(NO2 ~
model2 + Pct_65 + Pct_Img + Pct_brevet + NivVieMed +
Pct0_14 s(X,Y, bs = 'gp'),
data = LyonIris, method = 'REML',
family = gaussian(link = "identity"))
summary(model2)
# Modèle intégrant intégrant les coefficients variant spatialement
<- gam(NO2 ~
model3 s(X,Y, bs = 'gp', k = 25, by = Pct0_14) +
s(X,Y, bs = 'gp', k = 25, by = Pct_brevet) +
s(X,Y, bs = 'gp', k = 25, by = Pct_65) +
s(X,Y, bs = 'gp', k = 25, by = NivVieMed),
data = LyonIris, method = 'REML',
family = gaussian(link = "identity"))
summary(model3)
round(concurvity(model3),3)
# Modèle intégrant intégrant avec uniquement le niveau de vie variant spatialement
<- gam(NO2 ~ Pct0_14 + Pct_brevet + Pct_65 +
model3 s(X,Y, bs = 'gp', k = 25, by = NivVieMed),
data = LyonIris, method = 'REML',
family = gaussian(link = "identity"))
summary(model3)
<- AIC(model1, model2, model3)
comp_df $REML <- c(model1$gcv.ubre, model2$gcv.ubre, model3$gcv.ubre)
comp_df$R2 <- c(summary(model1)$r.sq,
comp_dfsummary(model2)$r.sq,
summary(model3)$r.sq)
::kable(comp_df, row.names = TRUE, digits = 2)
knitr
$residual_1 <- residuals(model1, type = 'scaled.pearson')
LyonIris$residual_2 <- residuals(model2, type = 'scaled.pearson')
LyonIris$residual_3 <- residuals(model3, type = 'scaled.pearson')
LyonIris
<- knearneigh(st_centroid(LyonIris), k = 10)
knn_list <- knn2nb(knn_list, sym = TRUE)
nb_list <- nb2listwdist(nb_list, st_centroid(LyonIris), style = 'W', type = 'idw', alpha = 2)
listw
<- moran.mc(LyonIris$residual_1, listw = listw, nsim = 999)
test1 <- moran.mc(LyonIris$residual_2, listw = listw, nsim = 999)
test2 <- moran.mc(LyonIris$residual_3, listw = listw, nsim = 999)
test3
<- data.frame(
tableau variable = c('residus1', 'residus2', 'residus3'),
I = c(test1$statistic, test2$statistic, test3$statistic),
P = c(test1$p.value, test2$p.value, test3$p.value)
)print(tableau)
summary(model3)
# Cartographie des effets spatiaux
library(terra)
library(gratia)
library(matrixStats)
<- LyonIris %>%
mask_data st_union() %>%
vect()
# Préparation d'un ensemble de points de prédictions formant une grille
# régulière de quadrats de 200 mètres de coté
<- st_bbox(mask_data)
bbox <- expand.grid(
coords X = seq(bbox$xmin, bbox$xmax, 100),
Y = seq(bbox$ymin, bbox$ymax, 100)
)
# Création du dataframe pour les prédictions
<- coords
pred_df $Pct0_14 <- 0
pred_df$Pct_brevet <- 0
pred_df$Pct_65 <- 0
pred_df$NivVieMed <- mean(LyonIris$NivVieMed)
pred_df
# Prédictions à partir du modèle
<- predict(model3, newdata = pred_df, type = 'terms')
preds <- data.frame(preds)
preds <- cbind(pred_df[, c("X", "Y")], preds)
preds names(preds) <- c("X", "Y", "Pct0_14", "Pct_brevet", "Pct_65", "spline_NivVieMed")
# Conversion des prédictions en objets raster pour la cartographie
<- c("X", "Y","spline_NivVieMed")
vars <- as.matrix(preds[, vars])
data_rast <- rast(data_rast[,c(1,2,3)],type = 'xyz')
raster_NivVieMed crs(raster_NivVieMed) <- crs(LyonIris)
# Cartographie de l'effet spatial de niveau de vie
<- list(text.separator = "à",
legende_parametres decimal.mark = ",",
big.mark = " "
)tmap_options(frame = FALSE, legend.frame = FALSE)
tm_shape(mask(raster_NivVieMed, mask_data)) +
tm_raster(
col.scale = tm_scale_intervals(
n = 6,
midpoint = 0,
style = "pretty",
values = '-brewer.rd_bu',
label.format = legende_parametres),
col.legend = tm_legend(title = "Niveau de vie",
position = tm_pos_out(pos.h = "left", pos.v = "top")
) )
11.10.2 Exercice 2
load("data/Lyon.Rdata")
library(sf, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(tmap, quietly = TRUE)
library(dplyr, quietly = TRUE)
library(spdep, quietly = TRUE)
library(DHARMa, quietly = TRUE)
library(sdmTMB, quietly = TRUE)
# Récupération en mise à l'échelle des coordonnées spatiales (division par 100)
<- c('Pct0_14', 'Pct_65', 'Pct_Img', 'Pct_brevet', 'NivVieMed')
vars <- st_coordinates(st_centroid(LyonIris))
XY <- LyonIris[, c("NO2", vars)]
LyonIris $X <- XY[,1] / 100
LyonIris$Y <- XY[,2] / 100
LyonIris
# Standardisation des différentes variables quantitatives
for(var in vars){
paste0('cr_', var)]] <-
LyonIris[[scale(LyonIris[[var]], center = TRUE, scale = TRUE)
}
# Réalisation de la mesh
<- make_mesh(st_drop_geometry(LyonIris), c('X', 'Y'), cutoff = 20)
mesh1 <- make_mesh(st_drop_geometry(LyonIris), c('X', 'Y'), cutoff = 1)
mesh2 <- make_mesh(st_drop_geometry(LyonIris), c('X', 'Y'), cutoff = 5)
mesh3
plot(mesh1)
plot(mesh2)
plot(mesh3)
11.10.3 Exercice 3
# Premier modèle avec uniquement la constante variant spatialement
<- sdmTMB(
model1 data = st_drop_geometry(LyonIris),
formula = NO2 ~
+ cr_Pct_65 + cr_Pct_Img + cr_Pct_brevet + cr_NivVieMed,
cr_Pct0_14 mesh = mesh2,
family = gaussian(link = "identity"),
spatial = TRUE
)
# Analyse des résidus simulés du modèle
<- simulate(model1, nsim = 500, type = "mle-mvn")
simulations <- dharma_residuals(simulations, model1, return_DHARMa = TRUE)
resids plot(resids)
# Second modèle qui fait aussi varier spatialement les variables indépendantes
<- sdmTMB(
model2 data = st_drop_geometry(LyonIris),
formula = NO2 ~
+ cr_Pct_65 + cr_Pct_Img + cr_Pct_brevet + cr_NivVieMed,
cr_Pct0_14 mesh = mesh2,
spatial = TRUE,
spatial_varying = ~ 1 + cr_Pct0_14 + cr_Pct_65 +
+ cr_Pct_brevet + cr_NivVieMed,
cr_Pct_Img family = gaussian(link = "identity"),
)
# Identification des variables qui varient spatialement
tidy(model2, effects = 'ran_pars', conf.int = TRUE) %>%
::kable(digits = 3)
knitr
# Troisième modèle avec uniquement les variables à faire
# varier spatialement
<- sdmTMB(
model3 data = st_drop_geometry(LyonIris),
formula = NO2 ~
+ cr_Pct_65 + cr_Pct_Img + cr_Pct_brevet + cr_NivVieMed,
cr_Pct0_14 mesh = mesh2,
spatial = TRUE,
spatial_varying = ~ 1 + cr_Pct_Img + cr_NivVieMed,
family = gaussian(link = "identity"),
)
# analyse des résidus simulés du modèle
<- simulate(model3, nsim = 500, type = "mle-mvn")
simulations <- dharma_residuals(simulations, model3, return_DHARMa = TRUE)
resids plot(resids)
<- cAIC(model1)
caic1 <- cAIC(model2)
caic2 <- cAIC(model3)
caic3 print(c(caic1, caic2, caic3))
# Extraction des coefficients et de leurs intervalles de confiance
<- tidy(model3, effects = 'fixed', conf.int = TRUE)
tab_fixe
# Récupération des écarts types
<- c('Pct0_14', 'Pct_65', 'Pct_Img', 'Pct_brevet', 'NivVieMed')
vars <- sapply(vars, function(x){
sd_values sd(LyonIris[[x]])
})<- data.frame(
sd_df sd = c(1,sd_values),
var = c('(Intercept)', paste0('cr_',vars))
)# Application des écarts types aux effets estimés
<- match(tab_fixe$term, sd_df$var)
ids $estimate_std <- tab_fixe$estimate
tab_fixeprint(tab_fixe)
for(col in c("estimate", "conf.low", "conf.high")){
<- tab_fixe[[col]] / sd_df$sd[ids]
tab_fixe[[col]]
}$std.error <- NULL
tab_fixe
# Extraction des prédictions
<- predict(model3, type = 'link')
preds
# Récupération de la constante spatiale
$rand_intercept <- preds$omega_s
LyonIris
# Extraction de 999 simulations des valeurs locales des coefficients spatiaux
library(matrixStats)
set.seed(888)
<- predict(model3, type = 'link', nsim = 999, sims_var = 'zeta_s')
preds2
# Extraction de 999 simulations des valeurs locales de la constante spatiale
<- predict(model3, type = 'link', nsim = 999, sims_var = 'omega_s')
preds3
# Pour chacun des coefficients spatiaux, déterminer si la variation spatiale est
# significativement différente de 0 au seuil 0.05
<- sapply(preds2, function(tab){
pvals <- rowQuantiles(tab, probs = c(0.025, 0.975))
qtl <- ((qtl[,1] < 0) & (qtl[,2] < 0)) | ((qtl[,1] > 0) & (qtl[,2] > 0))
test return(test)
})<- data.frame(pvals)
pvals
# idem avec la constante spatiale
<- rowQuantiles(preds3, probs = c(0.025, 0.975))
qtl <- ((qtl[,1] < 0) & (qtl[,2] < 0)) | ((qtl[,1] > 0) & (qtl[,2] > 0))
test $rand_intercept_sign <- test
LyonIristable(LyonIris$rand_intercept_sign)
# Récupération des coefficients spatiaux pour nos différentes variables
# avec un effet spatial
<- c('cr_Pct_Img', 'cr_NivVieMed')
sp_vars for(var in sp_vars){
# récupération du coefficient fixe de base
<- tab_fixe$estimate_std[tab_fixe$term == var]
base_coeff # récupération des variations spatiales du coefficient
<- preds[[paste0('zeta_s_', var)]]
sp_coeff # addition et conversion vers l'échelle originale des données
<- (base_coeff + sp_coeff) / sd_df$sd[sd_df$var == var]
loc_coeff paste0('sp_coeff_', var)]] <- loc_coeff
LyonIris[[paste0('sp_sign_', var)]] <- pvals[[var]]
LyonIris[[
}table(LyonIris$sp_sign_cr_Pct_Img)
table(LyonIris$sp_sign_cr_NivVieMed)
<- list(text.separator = "à",
legende_parametres decimal.mark = ",",
big.mark = " ",
digits = 2)
tmap_options(frame = FALSE, legend.frame = FALSE)
tm_shape(LyonIris) +
tm_polygons(
col = NULL,
fill = 'rand_intercept',
fill.scale = tm_scale_intervals(
style = 'fisher', n = 7, midpoint = 0,
values = '-brewer.rd_bu',
label.format = legende_parametres),
fill.legend = tm_legend(
title = "Constante aléatoire spatial",
position = tm_pos_out("center", "bottom"),
orientation = "landscape")
+
) tm_shape(subset(LyonIris, !LyonIris$rand_intercept_sign)) +
tm_polygons(col = NULL, fill = 'grey') +
tm_add_legend(type = 'polygons',
position = tm_pos_out("center", "bottom"),
labels = 'non significatif',
col = NA)
tm_shape(LyonIris) +
tm_polygons(
col = NULL,
fill = 'sp_coeff_cr_Pct_Img',
fill.scale = tm_scale_intervals(
style = 'fisher', n = 7, midpoint = 0,
values = '-brewer.rd_bu',
label.format = legende_parametres),
fill.legend = tm_legend(
title = "Constante aléatoire spatial",
position = tm_pos_out("center", "bottom"),
orientation = "landscape")
+
) tm_shape(subset(LyonIris, !LyonIris$sp_sign_cr_Pct_Img)) +
tm_polygons(col = NULL, fill = 'grey') +
tm_add_legend(type = 'polygons',
position = tm_pos_out("center", "bottom"),
labels = 'non significatif',
col = NA)
tm_shape(LyonIris) +
tm_polygons(
col = NULL,
fill = 'sp_coeff_cr_NivVieMed',
fill.scale = tm_scale_intervals(
style = 'fisher', n = 7, midpoint = 0,
values = '-brewer.rd_bu',
label.format = legende_parametres),
fill.legend = tm_legend(
title = "Constante aléatoire spatial",
position = tm_pos_out("center", "bottom"),
orientation = "landscape")
+
) tm_shape(subset(LyonIris, !LyonIris$sp_sign_cr_NivVieMed)) +
tm_polygons(col = NULL, fill = 'grey') +
tm_add_legend(type = 'polygons',
position = tm_pos_out("center", "bottom"),
labels = 'non significatif',
col = NA)