Objectif:
Étudier la dépendance entre des variables hydrologiques liées
aux inondations urbaines à Manilles (aux Philippines) à l’aide
des copules, en mettant en évidence des dépendances non
linéaires et/ou de queue, non captées par les corrélations
classiques
1- Chargement du jeu de données
df_flood = read.csv("C://Users//kouak//Downloads//Flood_Prediction_NCR_Philippines.csv")
head(df_flood)
## Date Location Rainfall_mm WaterLevel_m SoilMoisture_pct Elevation_m
## 1 2016-01-01 Quezon City 12.0 0.5 15.3 43
## 2 2016-01-01 Marikina 10.6 1.8 23.2 15
## 3 2016-01-01 Manila 5.7 0.5 15.6 5
## 4 2016-01-01 Pasig 3.7 0.5 5.0 5
## 5 2016-01-02 Quezon City 3.4 0.5 13.3 43
## 6 2016-01-02 Marikina 0.7 0.5 10.7 15
## FloodOccurrence
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0
Description des variables:
Rainfall_mm: la précipitation journalière
(mm)
WaterLevel: le niveau du cours d’eau à la
station la plus proche (m)
SoilMoisture: le niveau d’humidité du sol
(%)
Elevation_m: élévation statique au dessus du
niveau de la mer
Dans la suite du projet nous allons nous concentrer uniquement sur la
ville de Manilles
df_Manila = df_flood[df_flood$Location =="Manila",]
head(df_Manila)
## Date Location Rainfall_mm WaterLevel_m SoilMoisture_pct Elevation_m
## 3 2016-01-01 Manila 5.7 0.5 15.6 5
## 7 2016-01-02 Manila 5.3 0.5 13.2 5
## 11 2016-01-03 Manila 2.4 0.5 12.2 5
## 15 2016-01-04 Manila 4.7 1.1 15.4 5
## 19 2016-01-05 Manila 7.9 1.8 19.0 5
## 23 2016-01-06 Manila 22.0 2.0 18.4 5
## FloodOccurrence
## 3 0
## 7 0
## 11 0
## 15 0
## 19 0
## 23 0
2- Statistiques descriptives
Nous nous concentrons sur les variables
df_var <- df_Manila[, c("Rainfall_mm", "WaterLevel_m", "SoilMoisture_pct")]
summary(df_var)
## Rainfall_mm WaterLevel_m SoilMoisture_pct
## Min. : 0.200 Min. :0.500 Min. : 5.00
## 1st Qu.: 4.900 1st Qu.:0.500 1st Qu.:10.85
## Median : 8.400 Median :1.000 Median :14.80
## Mean : 9.906 Mean :1.273 Mean :15.16
## 3rd Qu.:13.400 3rd Qu.:1.800 3rd Qu.:19.00
## Max. :42.600 Max. :5.600 Max. :36.00
Distributions des variables
par(mfrow = c(1, 3))
# Précipitations
hist(df_var$Rainfall_mm,
prob = TRUE,
main = "Distribution des précipitations",
xlab = "Pluie (mm)",
col = "lightblue",
border = "white")
lines(density(df_var$Rainfall_mm, na.rm = TRUE),
col = "blue", lwd = 2)
#⃣ Niveau d'eau
hist(df_var$WaterLevel,
prob = TRUE,
main = "Distribution du niveau d'eau",
xlab = "Niveau d'eau (m)",
col = "lightgreen",
border = "white")
lines(density(df_var$WaterLevel, na.rm = TRUE),
col = "darkgreen", lwd = 2)
# Humidité du sol
hist(df_var$SoilMoisture,
prob = TRUE,
main = "Distribution de l'humidité du sol",
xlab = "Humidité du sol (%)",
col = "lightcoral",
border = "white")
lines(density(df_var$SoilMoisture, na.rm = TRUE),
col = "red", lwd = 2)
Nuages de points
par=(mfrow=c(1,3))
plot(df_var$Rainfall_mm, df_var$SoilMoisture_pct,
xlab = "Précipitations (mm)",
ylab = "Humidité du sol (%)",
main = "Nuage de points pluie – humidité du sol")
plot(df_var$Rainfall_mm, df_var$WaterLevel_m,
xlab = "Précipitations (mm)",
ylab = "Niveau de l'eau (m)",
main = "Nuage de points pluie – Niveau de l'eau")
plot(df_var$WaterLevel_m, df_var$SoilMoisture_pct,
xlab = "Niveau de l'eau (m)",
ylab = "Humidité du sol (%)",
main = "Niveau de l'eau – humidité du sol")
Sur le premier graphique, l’on observe une dépendance monotone
entre les précipitations et le taux d’humidité du sol avec une forte
dépendance dans la queue de distribution en haut à droite
Sur le second graphique, l’on constate que les faibles niveaux de
pluie (-20mm) n’ont pas d’effet systématique sur le niveau de cours
d’eau. Mais au delà de 20 mm de pluie, l’on constate une dépendance
monotone entre les précipitations et le niveau du cours d’eau avec une
forte dépendance dans la queue en haut à droite
Sur le troisième graphique, l’on ne constate ne soupçonne pas de
lien apparent dépendance entre les deux variables. Nous sommes proche de
la quasi-indépendance
Dans la suite du travail, nous prenons la peine de vérifier ces
observations avec des approches beaucoup plus rigoureuses vues en
cours
3- Pseudo-observations
n <- nrow(df_var)
U1 <- rank(df_var$Rainfall_mm) / (n + 1)
U2 <- rank(df_var$SoilMoisture_pct) / (n + 1)
U3 <- rank(df_var$WaterLevel_m) / (n + 1)
U_mat1 <- cbind(U1, U2)
colnames(U_mat1) <- c("U1", "U2")
U_mat2 <- cbind(U1, U3)
colnames(U_mat2) <- c("U1", "U3")
U_mat3 <- cbind(U2, U3)
colnames(U_mat3) <- c("U2", "U3")
par=(mfrow=c(1,3))
plot(U1, U2,
xlab = "Précipitations (mm)",
ylab = "Humidité du sol",
main = "pseudo-obs pluie – humidité du sol")
plot(U1, U3,
xlab = "Précipitations (mm)",
ylab = "Niveau de l'eau",
main = "Pseudo-obs pluie – Niveau de l'eau")
plot(U2, U3,
xlab = "Niveau de l'eau",
ylab = "Humidité du sol",
main = "pseudo-obs Niveau de l'eau – Humidité du sol")
Les pseudo-observations nous permettent de confirmer nos intuitions
au niveau des dépendances dans les queues. Elles font penser à la copule
de Gumbel ou à la copule de survie de Clayton.
Place aux tests d’indépendance pour confirmer définitivement nos
observations. Ce sera surtout l’occasion de confirmer notre intuition
pour le troisième graphique.
4-Tests d’indépendance
Test de Kendall
cor.test(x=U_mat1[,1], y=U_mat1[,2], method="kendall")
##
## Kendall's rank correlation tau
##
## data: U_mat1[, 1] and U_mat1[, 2]
## z = 21.915, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.3440394
cor.test(x=U_mat2[,1], y=U_mat2[,2], method="kendall")
##
## Kendall's rank correlation tau
##
## data: U_mat2[, 1] and U_mat2[, 2]
## z = 21.702, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.3585618
cor.test(x=U_mat3[,1], y=U_mat3[,2], method="kendall")
##
## Kendall's rank correlation tau
##
## data: U_mat3[, 1] and U_mat3[, 2]
## z = 12.112, p-value < 2.2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.2002208
D’après les tests, il existe une relation de dépendance entre les
deux variables de chaque paire. L’on remarque que ce sont des
dépendances relativement faibles au regards des taux de Kendall
respectifs estimés.
En particulier, pour le troisième couple, contrairement à notre
intuition, l’on a bel et bien une relation de dépendance faible mais
significative ! Mais le problème est que la forme du nuage de points ne
nous permet pas d’identifier aisément le type de copule à laquelle l’on
a affaire! On verra ça.
x-plot
Le test de Kendall nous permet d’apprècier l’existence d’une relation
de dépendance générale entre les variables. Mais nous avons aussi besoin
de confirmer nos observations au niveau des queues. Pour cela, nous
allons réaliser un x-plot ( lire chi-plot ) à partir
des pseudo-observations. En plus, le chi-plot peut permettre aussi de
trancher au niveau de la dépendance globale.
#install.packages("VineCopula")
library(VineCopula)
BiCopChiPlot(u1 = U_mat1[,1], u2 = U_mat1[,2], PLOT = T)
Les valeurs de chi sont toutes positives. Cela confirme une
dépendance globale positive.
Au centre ( autour de lamda=0), l’on a des valeurs de chi
élevées. l’on a donc une dépendance forte au centre. C’est-à-dire les
valeurs typiques tendent à évoluer de manière concordantes
L’on remarque également une dépendance forte au niveau supérieur
droit car pour lambda proche de 1, l’on a une forte concentration et
l’on a les plus hautes valeurs de chi
Toutes ces conditions nous font penser à une copule de Gumbel ou une
copule de survie de Clayton
BiCopChiPlot(u1 = U_mat2[,1], u2 = U_mat2[,2], PLOT = T)
De même l’on a une dépendance positive mais faible en générale
car les valeurs de chi bien qu’étant positive, sont comprises, pour la
plupart entre 0 et 0.5 avec une dépendance relativement forte au
centre.
Elle semble laisser paraître une forte dépendance de queue à
droite
Ce graphe a une structure assez capricieuse. Nous sommes réticents
mais au regard des pseudo observations, nous allons miser sur une copule
de Gumbel ou une copule de survie de Clayton
BiCopChiPlot(u1 = U_mat3[,1], u2 = U_mat3[,2], PLOT = T)
De façon générale nous avons des valeurs de chi positive et alignées
en bande le long de l’axe chi=0. Alors la dépendance est positive mais
elle est trop faible. De plus, en raison de la structure horizontale le
long de l’axe des abscisses, l’on se rapproche ici d’une situation
d’indépendance . Donc nous n’allons pas estimer de copules pour ce
couple de variable dans la suite de notre travail.
5-Estimation empirique des coefficients de dépendance de queue
Fonctions de
calcul des coefficients de queue
# Coefficient de queue de droite
lambda_U_hat <- function(U, q ) {
u <- U[,1]
v <- U[,2]
num <- sum(u > 1-q & v > 1-q)
den <- sum(u > 1-q)
if (den == 0) return(NA_real_)
num / den
}
# coefficient de queue de gauche
lambda_L_hat <- function(U, q ) {
u <- U[,1]
v <- U[,2]
num <- sum(u <= q & v <= q)
den <- sum(u <= q)
if (den == 0) return(NA_real_)
num / den
}
q_grid <- seq(0.01, 0.99, by = 0.01)
Pour le couple
précipitations - humidité du sol
par=(mfrow=c(1,2))
# estimation en queue supérieure
lambda_val1 <- sapply(q_grid, function(q) lambda_U_hat(U_mat1, q))
plot(q_grid, lambda_val1,
type = "l",
lwd = 2,
xlab = "q",
ylab = expression(hat(lambda)[U](q)),
main = "Estimateur empirique de la dépendance en queue supérieure")
abline(h = 0, lty = 2, col = "grey")
# Estimateur en queue inférieure
lambda_val2 <- sapply(q_grid, function(q) lambda_L_hat(U_mat1, q))
plot(q_grid, lambda_val2,
type = "l",
lwd = 2,
xlab = "q",
ylab = expression(hat(lambda)[L](q)),
main = "Estimateur empirique de la dépendance en queue inférieure")
abline(h = 0, lty = 2, col = "grey")
à gauche λL = 0
à droite λU = 1
pour le couple
précipitations - niveau de l’eau
par=(mfrow=c(1,2))
# estimation en queue supérieure
lambda_val3 <- sapply(q_grid, function(q) lambda_U_hat(U_mat2, q))
plot(q_grid, lambda_val3,
type = "l",
lwd = 2,
xlab = "q",
ylab = expression(hat(lambda)[U](q)),
main = "Estimateur empirique de la dépendance en queue supérieure")
abline(h = 0, lty = 2, col = "grey")
# Estimateur en queue inférieure
lambda_val4 <- sapply(q_grid, function(q) lambda_L_hat(U_mat2, q))
plot(q_grid, lambda_val4,
type = "l",
lwd = 2,
xlab = "q",
ylab = expression(hat(lambda)[L](q)),
main = "Estimateur empirique de la dépendance en queue inférieure")
abline(h = 0, lty = 2, col = "grey")
à gauche λL= 0
à droite λU = 1
Nous allons à présent estimer le paramètre de dépendance θ pour
chaque couple. En fait, le paramètre de dépendance contrôle la
force et la structure de dépendance
entre les variables. Pour faire simple, il détermine ” à quel
point les variables sont liées”
6- Estimation par la méthode des moments
##
## Attachement du package : 'copula'
## L'objet suivant est masqué depuis 'package:VineCopula':
##
## pobs
Pour le
couple précipitations - humidité du sol
#estimation de tau pour Gumbel
tau_emp1 <- cor(U_mat1[,1], U_mat1[,2], method = "kendall")
#estimation de tau pour la copule de survie de Clayton
U_surv1 <- 1 - U_mat1
tau_surv1 <- cor(U_surv1[,1], U_surv1[,2], method = "kendall")
theta_gumb_mm1 <- iTau(gumbelCopula(), tau_emp1)
theta_sclay_mm1 <- iTau(claytonCopula(), tau_surv1)
theta_gumb_mm1; theta_sclay_mm1
## [1] 1.524482
## [1] 1.048963
Pour le couple
précipitations - niveau de l’eau
tau_emp2 <- cor(U_mat2[,1], U_mat2[,2], method = "kendall") #estimation de tau pour Gumbel
U_surv2 <- 1 - U_mat2
tau_surv2 <- cor(U_surv2[,1], U_surv2[,2], method = "kendall") #estimation de tau pour la copule de survie de Clayton
theta_gumb_mm2 <- iTau(gumbelCopula(), tau_emp2)
theta_sclay_mm2 <- iTau(claytonCopula(), tau_surv2)
theta_gumb_mm2; theta_sclay_mm2
## [1] 1.558997
## [1] 1.117993
7-Estimation par maximum de Vraissemblance
Pour le
couple précipitations - humidité du sol
# Gumbel
fitGumb1 <- fitCopula(
gumbelCopula(dim = 2),
U_mat1,
method = "ml"
)
# Clayton de survie
U_surv1 <- 1 - U_mat1
fitSClay1 <- fitCopula(
claytonCopula(dim = 2),
U_surv1,
method = "ml"
)
## Call: fitCopula(gumbelCopula(dim = 2), U_mat1, method = "ml")
## Fit based on "maximum likelihood" and 1827 2-dimensional observations.
## Gumbel copula, dim. d = 2
## Estimate Std. Error
## alpha 1.486 0.027
## The maximized loglikelihood is 284.1
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient
## 6 6
## Call: fitCopula(claytonCopula(dim = 2), U_surv1, method = "ml")
## Fit based on "maximum likelihood" and 1827 2-dimensional observations.
## Clayton copula, dim. d = 2
## Estimate Std. Error
## alpha 1.049 0.046
## The maximized loglikelihood is 266.3
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient
## 3 3
Pour le couple
précipitations - niveau de l’eau
# Gumbel
fitGumb2 <- fitCopula(
gumbelCopula(dim = 2),
U_mat2,
method = "ml"
)
# Clayton de survie
U_surv2 <- 1 - U_mat2
fitSClay2 <- fitCopula(
claytonCopula(dim = 2),
U_surv2,
method = "ml"
)
## Call: fitCopula(gumbelCopula(dim = 2), U_mat2, method = "ml")
## Fit based on "maximum likelihood" and 1827 2-dimensional observations.
## Gumbel copula, dim. d = 2
## Estimate Std. Error
## alpha 1.494 0.028
## The maximized loglikelihood is 277.1
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient
## 7 7
## Call: fitCopula(claytonCopula(dim = 2), U_surv2, method = "ml")
## Fit based on "maximum likelihood" and 1827 2-dimensional observations.
## Clayton copula, dim. d = 2
## Estimate Std. Error
## alpha 1.118 0.047
## The maximized loglikelihood is 251.1
## Optimization converged
## Number of loglikelihood evaluations:
## function gradient
## 3 3
Conclusion partielle:
Pour chacun des couples, les deux méthodes donnent des estimations
proches pour les deux copules. On voit que pour les deux copules la
valeur de l’estimée de θ montre bien qu’il y a une dépendence.
Nous allons maintenant comparer les deux modèles à l’aide de
l’AIC/BIC (possible pour l’estimation par maximum de vraisemblance)
8-Comparaison des copules
Les deux modèles ont le même nombre de paramètres, il suffit donc de
regarder la log-vraisemblance que l’on cherche à maximiser
Pour le
couple précipitations - humidité du sol
logLik(fitGumb1); logLik(fitSClay1)
## 'log Lik.' 284.0747 (df=1)
## 'log Lik.' 266.3464 (df=1)
Pour le couple
précipitations - niveau de l’eau
logLik(fitGumb2); logLik(fitSClay2)
## 'log Lik.' 277.0811 (df=1)
## 'log Lik.' 251.132 (df=1)
Il ressort que la copule de Gumbel est le meilleur modèle dans les
deux cas
9- Test de qualité d’ajustement
A cette étape, l’idée est de vérifier si nous avons eu raison de
choisir la copule de Gumbel pour modéliser nos différentes
dépendance
a- Pour le
couple précipitations - humidité du sol
Test de Cramér-Von Mises
Comparaison des lignes de
niveau:
Copule de Gumbel
library(lattice)
library(latticeExtra)
# Calculs
fitGumb1 <- fitCopula(gumbelCopula(), as.matrix(U_mat1), method = "mpl")
emp_cop1 <- empCopula(as.matrix(U_mat1))
# Graphique
contour(fitGumb1@copula, pCopula,
main = "Gumbel : Théorique (noir) vs Empirique (rouge)",
col = "black")
# Superposition
contour(emp_cop1, pCopula,
col = "red",
lty = 2,
add = TRUE)
# Légende (Correction du symbole inattendu ici)
legend("bottomright", legend = c("Théorique", "Empirique"),
col = c("black", "red"), lty = c(1, 2), lwd = 1.5)
Copule de Clayton
# 1. Préparation de la copule théorique (Clayton)
# On récupère l'objet copule ajusté
clay_cop <- fitSClay1@copula
# 2. Création du graphique de base avec la copule théorique
# On utilise contour() qui est la méthode native et stable
contour(clay_cop, pCopula,
main = "Clayton : Théorique (noir) vs Empirique (rouge)",
col = "black",
nlevels = 20)
# 3. Calcul de la copule empirique sur une grille
u <- seq(0, 1, length.out = 30) # Plus de points pour un contour plus lisse
grid <- as.matrix(expand.grid(u1 = u, u2 = u))
# Calcul des valeurs C.n (Copule empirique)
# C.n est la fonction de répartition empirique du package copula
val_emp <- C.n(grid, X = as.matrix(U_mat1))
# Transformation en matrice pour la fonction contour() de base
z_matrix <- matrix(val_emp, nrow = length(u))
# 4. Superposition de la copule empirique
# add = TRUE permet de dessiner par-dessus le premier graphique
contour(u, u, z_matrix,
col = "red",
lty = 2,
add = TRUE,
nlevels = 20,
drawlabels = FALSE)
# 5. Légende
legend("bottomright", legend = c("Clayton Théorique", "Empirique"),
col = c("black", "red"), lty = c(1, 2), bty = "n")
b- Pour le
couple précipitations - niveau de l’eau
Test de Cramér-Von Mises
Comparaison des lignes de
niveau:
Copule de Gumbel
emp_cop2 <- empCopula(as.matrix(U_mat2))
# 1. Copule théorique
contour(fitGumb2@copula, pCopula,
main = "Gumbel : Théorique (noir) vs Empirique (rouge)",
col = "black",
lwd = 1.5)
# 2. Copule empirique (superposition)
contour(emp_cop2, pCopula,
col = "red",
lty = 2,
add = TRUE)
# 3. Légende
legend("bottomright",
legend = c("Théorique", "Empirique"),
col = c("black", "red"),
lty = c(1, 2),
lwd = 1.5)
Copule de Clayton
# 1. Préparation de la copule théorique (Clayton)
# On récupère l'objet copule ajusté
clay_cop <- fitSClay2@copula
# 2. Création du graphique de base avec la copule théorique
# On utilise contour() qui est la méthode native et stable
contour(clay_cop, pCopula,
main = "Clayton : Théorique (noir) vs Empirique (rouge)",
col = "black",
nlevels = 20)
# 3. Calcul de la copule empirique sur une grille
u <- seq(0, 1, length.out = 30) # Plus de points pour un contour plus lisse
grid <- as.matrix(expand.grid(u1 = u, u2 = u))
# Calcul des valeurs C.n (Copule empirique)
# C.n est la fonction de répartition empirique du package copula
val_emp <- C.n(grid, X = as.matrix(U_mat2))
# Transformation en matrice pour la fonction contour() de base
z_matrix <- matrix(val_emp, nrow = length(u))
# 4. Superposition de la copule empirique
# add = TRUE permet de dessiner par-dessus le premier graphique
contour(u, u, z_matrix,
col = "red",
lty = 2,
add = TRUE,
nlevels = 20,
drawlabels = FALSE)
# 5. Légende
legend("bottomright", legend = c("Clayton Théorique", "Empirique"),
col = c("black", "red"), lty = c(1, 2), bty = "n")
Conclusion
Nous conservons pour:
le couple précipitations - humidité du sol, la
copule de Gumbel de paramètre
1.49
le couple précipitations - niveau de l’eau, la
copule de Gumbel de paramètre
1.5