8.3 Modèles GLM pour des variables de comptage
Dans cette section, nous présentons les principaux modèles utilisés pour modéliser des variables de comptage. Il peut s’agir de variables comme le nombre d’accidents à une intersection, le nombre de cafés par quartier, le nombre de cas d’une maladie donnée par entité géographique, etc.
8.3.1 Modèle de Poisson
Le modèle GLM de base pour modéliser une variable de comptage est le modèle de Poisson. Pour rappel, la distribution de Poisson a un seul paramètre: \(\lambda\). Il représente le nombre moyen d’évènements observés sur l’intervalle de temps d’une observation, ainsi que la dispersion de la distribution. En conséquence, \(\lambda\) doit être un nombre strictement positif; autrement dit, nous ne pouvons pas observer un nombre négatif d’évènements. Il est donc nécessaire d’utiliser une fonction de lien pour contraindre l’équation de régression sur l’intervalle \([0 ,+\infty]\). La fonction la plus utilisée est le logarithme naturel (log) dont la réciproque est la fonction exponentielle (exp).
Type de variable dépendante | Variable de comptage |
Distribution utilisée | Poisson |
Formulation |
\(Y \sim Poisson(\lambda)\) \(g(\lambda) = \beta_0 + \beta X\) \(g(x) =log(x)\) |
Fonction de lien | log |
Paramètre modélisé | \(\lambda\) |
Paramètres à estimer | \(\beta_0\), \(\beta\) |
Conditions d’application | Absence d’excès de zéros, absence de sur-dispersion ou de sous-dispersion |
8.3.1.1 Interprétation des paramètres
Les coefficients du modèle expriment l’effet du changement d’une unité des variables X sur \(\lambda\) (le nombre de cas) dans l’échelle logarithmique (log-scale). Pour rappel, l’échelle logarithmique est multiplicative : si nous convertissons les coefficients dans leur échelle originale avec la fonction exponentielle, leur effet n’est plus additif, mais multiplicatif. Prenons un exemple concret, admettons que nous avons ajusté un modèle de Poisson à une variable de comptage Y avec deux variables \(X_1\) et \(X_2\) et que nous avons obtenus les coefficients suivants :
\(\beta_0 = 1,8 \text{; } \beta_1 = 0,5 \text{; } \beta_2 = -1,5\)
L’interprétation basique (sur l’échelle logarithmique) est la suivante : une augmentation d’une unité de la variable \(X_1\) est associée avec une augmentation de 0,5 du logarithme du nombre de cas attendus. Une augmentation d’une unité de la variable \(X_2\) est associée avec une réduction de 1,5 unité du logarithme du nombre de cas attendus. Avec la conversion avec la fonction exponentielle, nous obtenons alors :
- \(exp\mbox{(0,5) = 1,649}\), soit une multiplication par 1,649 du nombre de cas attendu (aussi appelé taux d’incident) à chaque augmentation d’une unité de \(X_1\);
- \(exp\mbox{(-1,5) = 0,223}\), soit une division par 4,54 (1/0,223) du nombre de cas attendu (aussi appelé taux d’incident) à chaque augmentation d’une unité de \(X_2\).
Utilisons maintenant notre équation pour effectuer une prédiction si \(X_1 = 1\) et \(X_2 = 1\).
\(\lambda = exp(\mbox{1,8} + (\mbox{0,5}\times 1) + (\mbox{-1,5}\times1)) = \mbox{2,225}\)
Si nous augmentons \(X_1\) d’une unité, nous obtenons alors :
\(\lambda = exp(1,8 + (\mbox{0,5}\times 2) + (\mbox{-1,5}\times1)) = \mbox{3,670}\)
En ayant augmenté d’une unité \(X_1\), nous avons multiplié notre résultat par 1,649 (\(\mbox{2,225} \times \mbox{1,649 = 3,670}\)).
Notez que ces effets se multiplient entre eux. Si nous augmentons à la fois \(X_1\) et \(X_2\) d’une unité chacune, nous obtenons : \(\lambda = exp(1,8 + (\mbox{0,5} \times 2) + (\mbox{-1,5} \times 2)) = \mbox{0,818}\), ce qui correspond bien à \(\mbox{2,225} \times \mbox{1,649} \text{ (effet de }X_1\text{)} \times \mbox{0,223} \text{ (effet de }X_2\text{)} = \mbox{0,818}\).
Il existe des fonctions dans R qui calculent ces prédictions à partir des équations des modèles. Il est cependant essentiel de bien saisir ce comportement multiplicatif induit par la fonction de lien log.
8.3.1.1.1 Conditions d’application
Puisque la distribution de Poisson n’a qu’un seul paramètre, le modèle GLM de Poisson est exposé au même problème potentiel de sur-dispersion que les modèles binomiaux de la section précédente. Référez-vous à la section 8.2.1.2 pour davantage de détails sur le problème posé par la sur-dispersion. Pour détecter une potentielle sur-dispersion dans un modèle de Poisson, il est possible, dans un premier temps, de calculer le ratio entre la déviance du modèle et son nombre de degrés de liberté (SAS Institute Inc 2020b). Ce ratio doit être proche de 1, s’il est plus grand, le modèle souffre de sur-dispersion.
\[\begin{equation} \hat{\phi} = \frac{D(modele)}{N-p} \tag{8.19} \end{equation}\]
avec \(N\) et \(p\) étant respectivement les nombres d’observations et de paramètres. Il est également possible de tester formellement si la sur-dispersion est significative avec un test de dispersion.
La question de l’excès de zéros a été abordée dans la section 2.4.3.7 du chapitre 2. Il s’agit d’une situation où un plus grand nombre de zéros sont présents dans les données que ce que suppose la distribution de Poisson. Dans ce cas, il convient d’utiliser la distribution de Poisson avec excès de zéros.
8.3.1.2 Exemple appliqué dans R
Pour cet exemple, nous reproduisons l’analyse effectuée dans l’article de Cloutier et al. (2014). L’enjeu de cette étude était de modéliser le nombre de piétons blessés autour de plus de 500 carrefours dans les quartiers centraux de Montréal. Pour cela, trois types de variables étaient utilisées : des variables décrivant l’intersection, des variables décrivant les activités humaines dans un rayon d’un kilomètre autour de l’intersection et des variables représentant le trafic routier autour de l’intersection. Un effet direct de ce type d’étude est bien évidemment l’établissement de meilleures pratiques d’aménagement réduisant les risques encourus par les piétons lors de leurs déplacements en ville. Le tableau 8.21 présente l’ensemble des variables utilisées dans l’analyse.
Nom de la variable | Signification | Type de variable | Mesure |
---|---|---|---|
Feux_auto | Présence de feux de circulation | Variable binaire | 0 = absence; 1 = présence |
Feux_piet | Présence de feux de traversée pour les piétons | Variable binaire | 0 = absence; 1 = présence |
Pass_piet | Présence d’un passage piéton | Variable binaire | 0 = absence; 1 = présence |
Terreplein | Présence d’un terre-plein central | Variable binaire | 0 = absence; 1 = présence |
Apaisement | Présence de mesure d’apaisement de la circulation | Variable binaire | 0 = absence; 1 = présence |
LogEmploi | Logarithme du nombre d’emplois dans un rayon d’un kilomètre | Variable continue | Logarithme du nombre d’emploi. Utilisation du logarithme, car la variable est fortement asymétrique |
Densite_pop | Densité de population dans un rayon d’un kilomètre | Variable continue | Habitants par hectare |
Entropie | Diversité des occupations du sol dans un rayon d’un kilomètre (indice d’entropie) | Variable continue | Mesure de 0 à 1; 0 = spécialisation parfaite; 1 = diversité parfaite |
DensiteInter | Densité d’intersections dans un rayon d’un kilomètre (connexité) | Variable continue | Nombre d’intersection par km2 |
Artere | Présence d’une artère à l’intersection | Variable binaire | 0 = absence; 1 = présence |
Long_arterePS | Longueur d’artère dans un rayon d’un kilomètre | Variable continue | Exprimée en mètres |
NB_voies5 | Présence d’une cinq voies à l’intersection | Variable binaire | 0 = absence; 1 = présence |
La distribution originale de la variable est décrite à la figure 8.23. Les barres grises représentent la distribution du nombre d’accidents et les barres rouges une distribution de Poisson ajustée sans variable indépendante (modèle nul). Ce premier graphique peut laisser penser qu’un modèle de Poisson n’est pas nécessairement le plus adapté considérant le grand nombre d’intersections pour lesquelles nous n’avons aucun accident. Cependant, rappelons que la variable Y n’a pas besoin de suivre une distribution de Poisson. Dans un modèle GLM, l’hypothèse que nous formulons est que la variable dépendante (Y) conditionnée par les variables indépendantes (X) suit une certaine distribution (ici Poisson).
# Chargement des données
<- read.csv("data/glm/accident_pietons.csv", sep = ";")
data_accidents # Ajustement d'une distribution de Poisson sans variable indépendante
library(fitdistrplus)
<- fitdist(data_accidents$Nbr_acci,distr = "pois")
model_poisson # Création d'un graphique pour comparer les deux distributions
<- data.frame(x=c(0:19),
dfpoisson y=dpois(0:19, model_poisson$estimate)
)<- data.frame(table(data_accidents$Nbr_acci))
counts names(counts) <- c("nb_accident",'frequence')
$nb_accident <- as.numeric(as.character(counts$nb_accident))
counts$prop <- counts$frequence / sum(counts$frequence)
countsggplot() +
geom_bar(aes(x=nb_accident, weight = prop, fill = "real"), width = 0.6, data = counts)+
geom_bar(aes(x=x, weight = y, fill = "adj"), width = 0.15, data = dfpoisson)+
scale_x_continuous(limits = c(-0.5,7), breaks = c(0:7))+
scale_fill_manual(name = "",
breaks = c("real","adj"),
labels = c("distribution originale", "distribution de Poisson"),
values = c("real" = rgb(0.4,0.4,0.4), "adj" = "red"))+
labs(subtitle = "",
x = "nombre d'accidents",
y = "")
Vérification des conditions d’application
Comme précédemment, notre première étape est de vérifier l’absence de multicolinéarité excessive avec la fonction vif
du package car
.
vif(glm(Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet + Terreplein + Apaisement +
+ Densite_pop + Entropie + DensiteInter +
LogEmploi + Artere + NB_voies5,
Long_arterePS family = poisson(link="log"),
data = data_accidents))
## Feux_auto Feux_piet Pass_piet Terreplein Apaisement LogEmploi Densite_pop
## 2.861708 1.518668 2.321213 1.221683 1.059722 4.763692 1.153096
## Entropie DensiteInter Long_arterePS Artere NB_voies5
## 1.770904 2.040457 4.467841 1.887728 1.520514
Toutes les valeurs de VIF sont inférieures à 5. Notons tout de même que le logarithme de l’emploi et la longueur d’artère dans un rayon d’un kilomètre ont des valeurs de VIF proches de 5. La seconde étape du diagnostic consiste à calculer et à visualiser les distances de Cook.
# Ajustement d'une première version du modèle
<- glm(Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet + Terreplein + Apaisement +
modele + Densite_pop + Entropie + DensiteInter +
LogEmploi + Artere + NB_voies5,
Long_arterePS family = poisson(link="log"),
data = data_accidents)
# Calcul des distances de Cook
<- cooks.distance(modele)
cooksd <- data.frame(
df cook = cooksd,
oid = 1:length(cooksd)
)ggplot(data = df)+
geom_point(aes(x = oid, y = cook), size = 0.5, alpha = 0.5)+
labs(x = "",
y = "Distance de Cook")
La figure 8.24 signale la présence de trois observations avec des valeurs extrêmement fortes de distance de Cook. Nous les isolons dans un premier temps pour les analyser.
<- subset(data_accidents, cooksd>0.1)
cas_etrange print(cas_etrange)
## Nbr_acci Feux_auto Feux_piet Pass_piet Terreplein Apaisement EmpTotBuffer Densite_pop Entropie
## 1 19 1 1 1 1 0 7208.538 5980.923 0.8073926
## 26 7 0 0 1 0 0 1342.625 2751.012 0.0000000
## 428 0 1 1 1 0 1 13122.560 14148.827 0.6643891
## DensiteInter Long_arterePS Artere NB_voies5 log_acci catego_acci catego_acci2 Arret VAG sum_app
## 1 42.41597 6955.00 1 1 2.995732 1 1 0 1 4
## 26 73.35344 2849.66 0 0 2.079442 1 1 1 1 3
## 428 109.25066 4634.20 0 0 0.000000 0 0 0 0 3
## LogEmploi AccOrdinal PopHa
## 1 8.883021 2 5.980923
## 26 7.202382 2 2.751012
## 428 9.482088 0 14.148827
Les deux premiers cas sont des intersections avec de nombreux accidents (respectivement 19 et 7) qui risquent de perturber les estimations du modèle. Le troisième cas ne comprend en revanche aucun accident. Puisqu’il ne s’agit que de trois observations et que leurs distances de Cook sont très nettement supérieures aux autres, nous les retirons du modèle.
<- subset(data_accidents, cooksd<0.1)
data2 # Ajustement d'une première version du modèle
<- glm(Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet + Terreplein + Apaisement +
modele + Densite_pop + Entropie + DensiteInter +
LogEmploi + Artere + NB_voies5,
Long_arterePS family = poisson(link="log"),
data = data2)
# Calcul des distances de Cook
<- cooks.distance(modele)
cooksd <- data.frame(
df cook = cooksd,
oid = 1:length(cooksd)
)ggplot(data = df)+
geom_point(aes(x = oid, y = cook), size = 0.5, alpha = 0.5)+
labs(x = "",
y = "distance de Cook")
La figure 8.25 montre que nous n’avons plus d’observations fortement influentes dans le modèle après avoir retiré les trois observations identifiées précédemment. Nous devons à présent vérifier l’absence de sur-dispersion.
# Calcul du rapport entre déviance et nombre de degrés de liberté du modèle
deviance(modele)/(nrow(data2) - modele$rank)
## [1] 1.674691
Le rapport entre la déviance et le nombre de degrés de liberté du modèle est nettement supérieur à 1, indiquant une sur-dispersion excessive. Nous pouvons confirmer ce résultat avec la fonction dispersiontest
du package AER
.
library(AER)
# Test de sur-dispersion
dispersiontest(modele)
##
## Overdispersion test
##
## data: modele
## z = 5.2737, p-value = 6.686e-08
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion
## 1.891565
Contrairement à la forme classique d’un modèle de Poisson pour laquelle la dispersion attendue est de 1, le test nous indique qu’une dispersion de 1,89 serait mieux ajustée aux données.
Il est également possible d’illuster cet écart à l’aide d’un graphique représentant les valeurs réelles, les valeurs prédites, ainsi que la variance (sous forme de barres d’erreurs) attendue par le modèle (figure 8.26). Nous constatons ainsi que les valeurs réelles (en rouge) ont largement tendance à dépasser la variance attendue par le modèle, surtout pour les valeurs les plus faibles de la distribution.
# Extraction des prédictions du modèle
<- predict(modele, type = "response")
lambdas # Création d'un DataFrame pour contenir la prédiction et les vraies valeurs
<- data.frame(
df1 lambdas = lambdas,
reals = data2$Nbr_acci
)# Calcul de l'intervalle de confiance à 95 % selon la distribution de Poisson
# et stockage dans un second DataFrame
<- seq(0,round(max(lambdas)),1)
seqa <- data.frame(
df2 lambdas = seqa,
lower = qpois(p = 0.025, lambda = seqa),
upper = qpois(p = 0.975, lambda = seqa)
)# Affichage des valeurs réelles et prédites (en rouge)
# et de leur variance selon le modèle (en noir)
ggplot() +
geom_errorbar(data = df2,
mapping = aes(x = lambdas, ymin = lower, ymax = upper),
width = 0.2, color = rgb(0.4,0.4,0.4)) +
geom_point(data = df1,
mapping = aes(x = lambdas, y = reals),
color ="red", size = 0.5) +
labs(x = 'valeurs prédites',
y = "valeurs réelles")
Pour tenir compte de cette particularité des données, nous modifions légèrement le modèle pour utiliser une distribution de quasi-Poisson, intégrant spécifiquement un paramètre de dispersion. Cet ajustement ne modifie pas l’estimation des coefficients du modèle, mais modifie le calcul des erreurs standards et par extension les valeurs de p pour les rendre moins sensibles au problème de sur-dispersion. Une autre approche aurait été de calculer une version robuste des erreurs standards avec le package sandwich
comme nous l’avons fait dans la section 8.2.1 sur le modèle binomial. Après réajustement du modèle, le nouveau paramètre de dispersion estimé est de 1,92.
Les quasi-distributions
Les distributions binomiale et de Poisson ne disposent chacune que d’un paramètre décrivant à la fois leur dispersion et leur espérance. Elles manquent donc de flexibilité et échouent parfois à représenter fidèlement des données avec une forte variance. Il existe donc d’autres distributions, respectivement les distributions quasi-binomiale et quasi-Poisson comprenant chacune un paramètre supplémentaire pour contrôler la dispersion. Bien que cette solution soit attrayante, il ne faut pas perdre de vue que la sur ou la sous dispersion peuvent être causées par l’absence de certaines variables explicatives, la sur-représentation de zéros, ou encore une séparation parfaite de la variable dépendante causée par une variable indépendante.
<- glm(Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet + Terreplein + Apaisement +
modele2 + Densite_pop + Entropie + DensiteInter +
LogEmploi + Artere + NB_voies5,
Long_arterePS family = quasipoisson(link="log"),
data = data2)
Nous pouvons à présent comparer la distribution originale des données et les simulations issues du modèle. Notez que contrairement à la distribution de Poisson simple, il n’existe pas dans R de fonction pour simuler des valeurs issues d’une distribution de quasi-Poisson. Il est cependant possible d’exploiter sa proximité théorique avec la distribution binomiale négative pour définir notre propre fonction de simulation. La figure 8.27 permet de comparer la distribution originale (en gris) et l’intervalle de confiance à 95 % des simulations (en rouge). Nous remarquons que le modèle semble capturer efficacement la forme générale de la distribution originale. À titre de comparaison, nous pouvons effectuer le même exercice avec la distribution de Poisson classique (le code n’est pas montré pour éviter les répétitions). La figure 8.28 montre qu’un simple modèle de Poisson est très éloigné de la distribution originale de Y.
# Définition d'une fonction pour simuler des données quasi-Poisson
<- function(n, lambda, disp) {
rqpois rnbinom(n = n, mu = lambda, size = lambda/(disp-1))
}# Extraction des valeurs prédites par le modèle
<- predict(modele2, type="response")
preds # Génération de 1000 simulations pour chaque prédiction
<- summary(modele2)$dispersion
disp <- 1000
nsim <- lapply(1:length(preds),function(i){
cols <- preds[[i]]
lambda <- round(rqpois(n = nsim, lambda = lambda, disp = disp))
sims return(sims)
})<- do.call(rbind, cols)
mat_sims # Préparation des données pour le graphique (valeurs réelles)
<- data.frame(table(data2$Nbr_acci))
counts names(counts) <- c("nb_accident",'frequence')
$nb_accident <- as.numeric(as.character(counts$nb_accident))
counts$prop <- counts$frequence / sum(counts$frequence)
counts# Préparation des données pour le graphique (valeurs simulées)
<- data.frame(count = 0:25)
df1 <- lapply(1:nsim, function(i){
count_sims <- mat_sims[,i]
sim <- data.frame(table(sim))
cnt <- merge(df1,cnt, by.x="count", by.y = "sim", all.x = T, all.y=F)
df2 $Freq <- ifelse(is.na(df2$Freq),0,df2$Freq)
df2return(df2$Freq)
})<- do.call(cbind,count_sims)
count_sims <- data.frame(
df_sims val = 0:25,
med = apply(count_sims, MARGIN = 1, median),
lower = apply(count_sims, MARGIN = 1, quantile, probs = 0.025),
upper = apply(count_sims, MARGIN = 1, quantile, probs = 0.975)
)
ggplot() +
geom_bar(aes(x=nb_accident, weight = frequence), width = 0.6, data = counts)+
geom_errorbar(aes(x = val, ymin = lower, ymax = upper),
data = df_sims, color = "red", width = 0.6)+
geom_point(aes(x = val, y = med), color = "red", size = 1.3, data = df_sims)+
scale_x_continuous(limits = c(-0.5,7), breaks = c(0:7))+
xlim(-1,12)+
labs(subtitle = "",
x = "nombre d'accidents",
y = "nombre d'occurrences")
La prochaine étape du diagnostic est l’analyse des résidus simulés. La figure 8.29 indique que les résidus du modèle suivent bien une distribution uniforme et qu’aucune valeur aberrante n’est observable.
# Génération de 1000 simulations pour chaque prédiction
<- 1.918757 # trouvable dans le summary(modele2)
disp <- 1000
nsim <- lapply(1:length(preds),function(i){
cols <- preds[[i]]
lambda <- rqpois(n = nsim, lambda = lambda, disp = disp)
sims return(sims)
})<- do.call(rbind, cols)
mat_sims <- createDHARMa(simulatedResponse = mat_sims,
sim_res observedResponse = data2$Nbr_acci,
fittedPredictedResponse = modele2$fitted.values,
integerResponse = T)
plot(sim_res)
Pour affiner notre diagnostic, nous pouvons également comparer les résidus simulés et chaque variable indépendante. La figure 8.30 n’indique aucune relation problématique entre nos variables indépendantes et les résidus.
par(mfrow=c(3,4))
<- c("Feux_auto", "Feux_piet", "Pass_piet", "Terreplein", "Apaisement",
vars "LogEmploi", "Densite_pop", "Entropie", "DensiteInter",
"Long_arterePS", "Artere", "NB_voies5")
for(v in vars){
plotResiduals(sim_res, data2[[v]], xlab = v, main = "", ylab = "résidus")
}
Maintenant que l’ensemble des diagnostics a été effectué, nous pouvons passer à la vérification de la qualité d’ajustement.
Vérification de la qualité d’ajustement
Pour le calcul des pseudo-R2, notez qu’il n’existe pas à proprement parler de loglikelihood pour les quasi-distributions. Pour contourner ce problème, il est possible d’utiliser le loglikelihood d’un simple modèle de Poisson (puisque les coefficients ne changent pas), mais il est important de garder à l’esprit que ces pseudo-R2 seront d’autant plus faibles que la sur-dispersion originale était forte.
<- glm(Nbr_acci ~ 1,
modelnull family = poisson(link="log"),
data = data2)
rsqs(loglike.full = logLik(modele),
loglike.null = logLik(modelnull),
full.deviance = deviance(modele),
null.deviance = deviance(modelnull),
nb.params = modele$rank,
n = nrow(data2))
## $`deviance expliquee`
## [1] 0.4805998
##
## $`McFadden ajuste`
## 'log Lik.' 0.3258375 (df=13)
##
## $`Cox and Snell`
## 'log Lik.' 0.7789704 (df=13)
##
## $Nagelkerke
## 'log Lik.' 0.787958 (df=13)
Le modèle parvient ainsi à expliquer 48 % de la déviance totale. Il obtient un R2 ajusté de McFadden de 0,33 et un R2 de Cox et Snell de 0,78.
# Calcul du RMSE
sqrt(mean((predict(modele2, type = "response") - data2$Nbr_acci)**2))
## [1] 1.838026
# Nombre moyen d'accidents
mean(data2$Nbr_acci)
## [1] 1.503984
L’erreur quadratique moyenne du modèle est de 1,84, ce que signifie qu’en moyenne le modèle se trompe d’environ deux accidents pour chaque intersection. Cette valeur est relativement élevée si nous la comparons avec le nombre moyen d’accidents, soit 1,5. Cela s’explique certainement par le grand nombre de zéros dans la variable Y qui tendent à tirer les prédictions vers le bas.
Interprétation des résultats
L’ensemble des coefficients du modèle sont accessibles via la fonction summary
. Puisque la fonction de lien du modèle est la fonction log
, il est pertinent de convertir les coefficients avec la fonction exp
afin de pouvoir les interpréter sur l’échelle originale (nombre d'accidents
) plutôt que l’échelle logarithmique (log(nombre d'accidents)
). N’oubliez pas que ces effets sont multiplicatifs une fois transformés avec la fonction exp
. Nous pouvons également utiliser les erreurs standards pour calculer des intervalles de confiance à 95 % des exponentiels des coefficients. Le tableau 8.22 présente l’ensemble des informations pertinentes pour l’interprétation des résultats.
# Calcul des coefficients en exponentiel et des intervalles de confiance
<- summary(modele2)$coefficients
tableau <- tableau[,1]
coeffs <- tableau[,2]
err.std <- exp(coeffs)
expcoeff .5 <- exp(coeffs - 1.96*err.std)
exp2<- exp(coeffs - 1.96*err.std)
exp975 <- tableau[,4]
pvals <- cbind(coeffs,err.std,expcoeff,exp2.5,exp975,pvals)
tableauComplet # print(tableauComplet)
Variable | Coeff. | exp(Coeff.) | Val.p | IC 2,5 % exp(Coeff.) | IC 97,5 % exp(Coeff.) | Sign. |
---|---|---|---|---|---|---|
Constante | -3,680 | 0,030 | 0,000 | 0,010 | 0,090 | *** |
Feux_auto | 1,100 | 3,000 | 0,000 | 1,970 | 4,660 | *** |
Feux_piet | 0,330 | 1,390 | 0,009 | 1,090 | 1,790 | ** |
Pass_piet | 0,340 | 1,400 | 0,149 | 0,880 | 2,200 | |
Terreplein | -0,360 | 0,700 | 0,099 | 0,440 | 1,050 | . |
Apaisement | 0,290 | 1,330 | 0,157 | 0,880 | 1,950 | |
LogEmploi | 0,230 | 1,260 | 0,017 | 1,040 | 1,520 |
|
Densite_pop | 0,000 | 1,000 | 0,000 | 1,000 | 1,000 | *** |
Entropie | -0,420 | 0,660 | 0,271 | 0,320 | 1,390 | |
DensiteInter | 0,000 | 1,000 | 0,410 | 1,000 | 1,010 | |
Long_arterePS | 0,000 | 1,000 | 0,684 | 1,000 | 1,000 | |
Artere | 0,030 | 1,030 | 0,842 | 0,780 | 1,360 | |
NB_voies5 | 0,640 | 1,890 | 0,000 | 1,460 | 2,440 | *** |
Parmi les variables décrivant les aménagements de l’intersection, nous constatons que les présences d’un feu de circulation et d’un feu de traversée pour les piétons multiplient le nombre attendu d’accidents à une intersection par 3,0 et 1,39. Par contre, les présences d’un passage piéton, d’un terre-plein ou de mesures d’apaisement n’ont pas d’effets significatifs (valeurs de p > 0,05). Concernant les variables décrivant l’environnement à proximité des intersections, nous observons que la concentration d’emplois et la densité de population contribuent toutes les deux à augmenter le nombre d’accidents à une intersection, bien que leurs effets soient limités. Enfin, la présence d’une rue à cinq voies à l’intersection augmente le nombre d’accidents attendu à l’intersection de 89 %. Nous ne détaillons pas plus les résultats, car nous utilisons le même jeu de données dans les prochaines sections.
8.3.2 Modèle binomial négatif
Dans le cas où une variable de comptage est marquée par une sur ou sous-dispersion, la distribution de Poisson n’est pas en mesure de capturer efficacement sa variance. Pour contourner ce problème, il est possible d’utiliser la distribution binomiale négative plutôt que la distribution de Poisson (ou quasi-Poisson). Cette distribution peut être décrite comme une généralisation de la distribution de Poisson : elle inclut un second paramètre \(\theta\) contrôlant la dispersion. L’intérêt premier de ce changement de distribution est que l’interprétation des paramètres est la même pour les deux modèles, tout en contrôlant directement l’effet d’une potentielle sur ou sous-dispersion.
Type de variable dépendante | Variable de comptage |
Distribution utilisée | Négative binomiale |
Formulation |
\(Y \sim NB(\mu,\theta)\) \(g(\mu) = \beta_0 + \beta X\) \(g(x) =log(x)\) |
Fonction de lien | log |
Paramètre modélisé | \(\mu\) |
Paramètres à estimer | \(\beta_0\), \(\beta\) et \(\theta\) |
Conditions d’application | Absence d’excès de zéros, respect du lien variance-moyenne |
8.3.2.1 Conditions d’application
Les conditions d’application d’un modèle binomial négatif sont presque les mêmes que celles d’un modèle de Poisson. La seule différence est que la condition d’absence de sur ou sous-dispersion est remplacée par une condition de respect du lien espérance-variance. En effet, dans un modèle binomial négatif, le paramètre de dispersion \(\theta\) est combiné avec \(\mu\) (l’espérance) pour exprimer la dispersion de la distribution. Dans le package mgcv
que nous utilisons dans l’exemple, le lien entre \(\mu\), \(\theta\) et la variance est le suivant :
\[\begin{equation} variance = \mu + \mu^{\frac{2}{\theta}} \tag{8.20} \end{equation}\]
Il s’agit donc d’un modèle hétéroscédastique : sa variance n’est pas fixe, mais varie en fonction de sa propre espérance. Si celle-ci augmente, la variance augmente (comme pour un modèle de Poisson), et l’intensité de cette augmentation est contrôlée par le paramètre \(\theta\). Si cette condition n’est pas respectée, l’analyse des résidus simulés révélera un problème de dispersion.
8.3.2.2 Exemple appliqué dans R
Dans l’exemple précédent avec le modèle de Poisson, nous avons observé une certaine sur-dispersion que nous avons contournée en utilisant un modèle de quasi-Poisson. Dans l’article original, les auteurs ont opté pour un modèle binomial négatif, ce que nous reproduisons ici. Les variables utilisées sont les mêmes que pour le modèle de Poisson. Nous utilisons le package mgcv
et sa fonction gam
pour ajuster le modèle.
Vérification des conditions d’application
Nous avons vu précédemment que nos variables indépendantes ne sont pas marquées par une multicolinéarité forte. Il n’est pas nécessaire de recalculer les valeurs de VIF puisque nous utilisons les mêmes données. La première étape du diagnostic est donc de calculer les distances de Cook.
library(mgcv)
# Chargement des données
<- read.csv("data/glm/accident_pietons.csv", sep = ";")
data_accidents # Ajustement d'une première version du modèle
<- gam(Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet +
modelnb + Apaisement +
Terreplein + Densite_pop + Entropie + DensiteInter +
LogEmploi + Artere + NB_voies5,
Long_arterePS family = nb(link="log"),
data = data_accidents)
# Calcul et affichage des distances de Cook
<- cooks.distance(modelnb)
cooksd <- data.frame(
df cook = cooksd,
oid = 1:length(cooksd)
)ggplot(data = df)+
geom_point(aes(x = oid, y = cook), size = 0.5, color = rgb(0.4,0.4,0.4)) +
labs(x = "", y = "distance de Cook")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
Nous observons, dans la figure 8.31, que quatre observations se distinguent très nettement des autres.
<- subset(data_accidents, cooksd > 0.03)
cas_etrange print(cas_etrange)
## Nbr_acci Feux_auto Feux_piet Pass_piet Terreplein Apaisement EmpTotBuffer Densite_pop Entropie
## 1 19 1 1 1 1 0 7208.538 5980.923 0.8073926
## 5 12 1 0 1 0 0 8585.350 8655.430 0.7607852
## 26 7 0 0 1 0 0 1342.625 2751.012 0.0000000
## 34 6 0 0 1 0 0 12516.410 8950.942 0.4300549
## DensiteInter Long_arterePS Artere NB_voies5 log_acci catego_acci catego_acci2 Arret VAG sum_app
## 1 42.41597 6955.00 1 1 2.995732 1 1 0 1 4
## 5 89.11495 6412.27 0 0 2.564949 1 1 0 1 4
## 26 73.35344 2849.66 0 0 2.079442 1 1 1 1 3
## 34 74.91879 8443.01 1 0 1.945910 1 1 1 0 3
## LogEmploi AccOrdinal PopHa
## 1 8.883021 2 5.980923
## 5 9.057813 2 8.655430
## 26 7.202382 2 2.751012
## 34 9.434796 2 8.950942
Il s’agit à nouveau de quatre observations avec un grand nombre d’accidents. Nous décidons de les retirer du jeu de données pour ne pas fausser les résultats concernant l’ensemble des autres intersections. Dans une analyse plus détaillée, il serait judicieux de chercher à comprendre pourquoi ces quatre observations sont particulièrement accidentogènes.
<- subset(data_accidents, cooksd < 0.03)
data2 # Ajustement d'une première version du modèle
<- gam(Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet +
modelnb + Apaisement +
Terreplein + Densite_pop + Entropie + DensiteInter +
LogEmploi + Artere + NB_voies5,
Long_arterePS family = nb(link="log"),
data = data2)
# Calcul et affichage des distances de Cook
<- cooks.distance(modelnb)
cooksd <- data.frame(
df cook = cooksd,
oid = 1:length(cooksd)
)ggplot(data = df)+
geom_point(aes(x = oid, y = cook), size = 0.5, color = rgb(0.4,0.4,0.4)) +
labs(x = "", y = "distance de Cook")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
Après avoir retiré ces quatre observations, les distances de Cook (figure 8.32) ne révèlent plus d’observations fortement influentes dans le modèle. La prochaine étape du diagnostic est donc d’analyser les résidus simulés.
# Extraction de la valeur de theta
<- modelnb$family$getTheta(T)
theta <- 1000
nsim # Extraction des valeurs prédites par le modèle
<- predict(modelnb, type = "response")
mus # Calcul des simulations
<- lapply(1:length(mus),function(i){
cols <- mus[[i]]
mu <- rnbinom(n = nsim, mu = mu, size = theta)
sims return(sims)
})<- do.call(rbind, cols)
mat_sims # Calcul des résidus simulés
<- createDHARMa(simulatedResponse = mat_sims,
sim_res observedResponse = data2$Nbr_acci,
fittedPredictedResponse = mus,
integerResponse = T)
# Affichage du diagnostic
plot(sim_res)
La figure 8.33 présentant le diagnostic des résidus simulés, montre que ces derniers suivent bien une distribution uniforme et aucun problème de dispersion ni de valeurs aberrantes. La figure 8.34 permet de comparer la distribution originale de la variable Y et les simulations issues du modèle (intervalles de confiance représentés en bleu). Nous constatons que le modèle parvient bien à reproduire la distribution originale, et ce, même pour les valeurs les plus extrèmes de la distribution.
# Extraction des valeurs prédites par le modèle
<- predict(modelnb, type="response")
mus # Génération de 1000 simulations pour chaque prédiction
<- modelnb$family$getTheta(T)
theta <- 1000
nsim <- lapply(1:length(mus),function(i){
cols <- mus[[i]]
mu <- round(rnbinom(n = nsim, mu = mu, size = theta))
sims return(sims)
})<- do.call(rbind, cols)
mat_sims # Préparation des données pour le graphique (valeurs réelles)
<- data.frame(table(data2$Nbr_acci))
counts names(counts) <- c("nb_accident",'frequence')
$nb_accident <- as.numeric(as.character(counts$nb_accident))
counts$prop <- counts$frequence / sum(counts$frequence)
counts# Préparation des données pour le graphique (valeurs simulées)
<- data.frame(count = 0:25)
df1 <- lapply(1:nsim, function(i){
count_sims <- mat_sims[,i]
sim <- data.frame(table(sim))
cnt <- merge(df1,cnt, by.x="count", by.y = "sim", all.x = T, all.y=F)
df2 $Freq <- ifelse(is.na(df2$Freq),0,df2$Freq)
df2return(df2$Freq)
})<- do.call(cbind,count_sims)
count_sims <- data.frame(
df_sims val = 0:25,
med = apply(count_sims, MARGIN = 1, median),
lower = apply(count_sims, MARGIN = 1, quantile, probs = 0.025),
upper = apply(count_sims, MARGIN = 1, quantile, probs = 0.975)
)# Affichage du graphique
ggplot() +
geom_bar(aes(x=nb_accident, weight = frequence), width = 0.6, data = counts)+
geom_errorbar(aes(x = val, ymin = lower, ymax = upper),
data = df_sims, color = "blue", width = 0.6)+
geom_point(aes(x = val, y = med), color = "blue", size = 1.3, data = df_sims)+
scale_x_continuous(limits = c(-0.5,7), breaks = c(0:7))+
xlim(-1,12)+
labs(subtitle = "",
x = "nombre d'accidents",
y = "nombre d'occurrences")
À titre de comparaison, nous pouvons à nouveau réaliser le graphique permettant de visualiser si la variance attendue par le modèle est proche de celle effectivement observée dans les données. Nous avons constaté avec ce graphique, lorsque nous ajustions un modèle de Poisson, que la variance des données était trop grande comparativement à celle attendue par le modèle (figure 8.26).
# Extraction des prédictions du modèle
<- predict(modelnb, type = "response")
mus # Création d'un DataFrame pour contenir la prédiction et les vraies valeurs
<- data.frame(
df1 mus = mus,
reals = data2$Nbr_acci
)# Calcul de l'intervalle de confiance à 95 % selon la distribution de Poisson
# et stockage dans un second DataFrame
<- seq(0,round(max(mus)),1)
seqa <- data.frame(
df2 mus = seqa,
lower = qnbinom(p = 0.025, mu = seqa, size = theta),
upper = qnbinom(p = 0.975, mu = seqa, size = theta)
)# Affichage des valeurs réelles et prédites (en rouge)
# et de leur variance selon le modèle (en noir)
ggplot() +
geom_errorbar(data = df2,
mapping = aes(x = mus, ymin = lower, ymax = upper),
width = 0.2, color = rgb(0.4,0.4,0.4)) +
geom_point(data = df1,
mapping = aes(x = mus, y = reals),
color ="red", size = 0.5) +
labs(x = 'valeurs prédites',
y = "valeurs réelles")
Nous pouvons ainsi constater à la figure 8.35 que le modèle binomial négatif autorise une variance bien plus large que le modèle de Poisson et est ainsi mieux ajusté aux données.
Vérification de la qualité d’ajustement
# Calcul des pseudo R2
rsqs(loglike.full = logLik(modelnb),
loglike.null = logLik(modelnull),
full.deviance = deviance(modelnb),
null.deviance = modelnb$null.deviance,
nb.params = modelnb$rank,
n = nrow(data2))
## $`deviance expliquee`
## [1] 0.458052
##
## $`McFadden ajuste`
## 'log Lik.' 0.384353 (df=14)
##
## $`Cox and Snell`
## 'log Lik.' 0.8304773 (df=14)
##
## $Nagelkerke
## 'log Lik.' 0.8399731 (df=14)
# Calcul du RMSE
sqrt(mean((predict(modelnb, type = "response") - data2$Nbr_acci)**2))
## [1] 1.825278
Le modèle parvient à expliquer 45 % de la déviance. Il obtient un R2 ajusté de McFadden de 0,14 et un R2 de Nagelkerke de 0,42. L’erreur moyenne quadratique de la prédiction est de 1,82, ce qui est identique au modèle de Poisson ajusté précédemment.
Interprétation des résultats
Il est possible d’accéder à l’ensemble des coefficients du modèle via la fonction summary
. À nouveau, les coefficients doivent être convertis avec la fonction exponentielle (du fait de la fonction de lien log) et interprétés comme des effets multiplicatifs. Le tableau 8.24 présente les coefficients estimés par le modèle. Les résultats sont très similaires à ceux du modèle de quasi-Poisson original. Nous notons cependant que la variable représentant la présence d’un feu pour piéton n’est plus significative au seuil de 0,05.
Variable | Coeff. | exp(Coeff.) | Val.p | IC 2,5 % exp(Coeff.) | IC 97,5 % exp(Coeff.) | Sign. |
---|---|---|---|---|---|---|
Constante | -3,880 | 0,020 | 0,000 | 0,010 | 0,080 | *** |
Feux_auto | 1,130 | 3,100 | 0,000 | 2,030 | 4,710 | *** |
Feux_piet | 0,350 | 1,420 | 0,016 | 1,060 | 1,900 |
|
Pass_piet | 0,220 | 1,240 | 0,300 | 0,830 | 1,880 | |
Terreplein | -0,340 | 0,710 | 0,155 | 0,440 | 1,140 | |
Apaisement | 0,240 | 1,270 | 0,315 | 0,790 | 2,030 | |
LogEmploi | 0,230 | 1,260 | 0,025 | 1,030 | 1,550 |
|
Densite_pop | 0,000 | 1,000 | 0,000 | 1,000 | 1,000 | *** |
Entropie | -0,170 | 0,840 | 0,669 | 0,380 | 1,860 | |
DensiteInter | 0,000 | 1,000 | 0,925 | 0,990 | 1,010 | |
Long_arterePS | 0,000 | 1,000 | 0,587 | 1,000 | 1,000 | |
Artere | 0,110 | 1,110 | 0,497 | 0,820 | 1,510 | |
NB_voies5 | 0,700 | 2,010 | 0,000 | 1,480 | 2,750 | *** |
8.3.3 Modèle de Poisson avec excès fixe de zéros
Dans le cas où la variable Y comprendrait significativement plus de zéros que ce que suppose une distribution de Poisson, il est possible d’utiliser la distribution de Poisson avec excès de zéros. Pour rappel, cette distribution ajoute un paramètre p contrôlant pour la proportion de zéros dans la distribution. Du point de vue conceptuel, cela revient à formuler l’hypothèse suivante : dans les données que nous avons observées, deux processus distincts sont à l’œuvre. Le premier est issu d’une distribution de Poisson et l’autre produit des zéros qui s’ajoutent aux données. Les zéros produits par la distribution de Poisson sont appelés les vrais zéros, alors que ceux produits par le second phénomène sont appelés les faux zéros.
Type de variable dépendante | Variable de comptage |
Distribution utilisée | Poisson avec excès de zéros |
Formulation |
\(Y \sim ZIP(\mu,\theta)\) \(g(\lambda) = \beta_0 + \beta X\) \(g(x) =log(x)\) |
Fonction de lien | log |
Paramètre modélisé | \(\lambda\) |
Paramètres à estimer | \(\beta_0\), \(\beta\) et p |
Conditions d’application | Absence de sur-dispersion |
Dans cette formulation, p est fixé. Nous n’avons donc aucune information sur ce qui produit les zéros supplémentaires, mais seulement leur proportion totale dans le jeu de données.
8.3.3.1 Interprétation des paramètres
L’interprétation des paramètres est identique à celle d’un modèle de Poisson. Le paramètre p représente la proportion de faux zéros dans la variable Y une fois que les variables indépendantes sont contrôlées.
8.3.3.2 Exemple appliqué dans R
La variable de comptage des accidents des piétons que nous avons utilisée dans les deux exemples précédents semble être une bonne candidate pour une distribution de Poisson avec excès de zéros. En effet, nous avons pu constater une sur-dispersion dans le modèle de Poisson original, ainsi qu’un nombre important d’intersections sans accident. Tentons donc d’améliorer notre modèle en ajustant un excès fixe de zéros. Nous utilisons la fonction gamlss
du package gamlss
.
library(gamlss)
<- gamlss(formula = Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet +
modelzi + Apaisement + LogEmploi + Densite_pop +
Terreplein + DensiteInter + Long_arterePS + Artere + NB_voies5,
Entropie sigma.formula = ~1,
family = ZIP(mu.link = "log", sigma.link="logit"),
data = data_accidents)
## GAMLSS-RS iteration 1: Global Deviance = 1516.53
## GAMLSS-RS iteration 2: Global Deviance = 1514.656
## GAMLSS-RS iteration 3: Global Deviance = 1514.52
## GAMLSS-RS iteration 4: Global Deviance = 1514.508
## GAMLSS-RS iteration 5: Global Deviance = 1514.506
## GAMLSS-RS iteration 6: Global Deviance = 1514.506
<- glm(formula = Nbr_acci ~ 1,
modelnull family = poisson(link="log"),
data = data_accidents)
# Constante pour p
<- modelzi$sigma.coefficients
coeff_p cat("Coefficient pour p =", round(coeff_p,4))
## Coefficient pour p = -1.4376
# Calcul de la déviance expliquée
1 - deviance(modelzi) / deviance(modelnull)
## [1] 0.08267513
# Calcul de la probabilité de base p d'être un faux 0
# en appliquant l'inverse de la fonction logistique
exp(-coeff_p) / (1+exp(-coeff_p))
## (Intercept)
## 0.80809
Nous constatons immédiatement que le modèle avec excès fixe de zéros est peu ajusté aux données. Cette version du modèle ne parvient à capter que 8 % de la déviance, ce qui s’explique facilement, car nous n’avons donné aucune variable au modèle pour distinguer les vrais et les faux zéros. Pour cela, nous devons passer au prochain modèle : Poisson avec excès ajusté de zéros. Notons tout de même que d’après ce modèle, 81 % des observations seraient des faux zéros.
8.3.4 Modèle de Poisson avec excès ajusté de zéros
Nous avons vu dans le modèle précédent, que l’excès de zéro était conceptualisé comme la combinaison de deux phénomènes, l’un issu d’une distribution de Poisson que nous souhaitons modéliser, et l’autre générant des zéros supplémentaires. Il est possible d’aller plus loin que de simplement contrôler la proportion de zéros supplémentaires en modélisant explicitement ce second processus en ajoutant une deuxième équation au modèle. Cette deuxième équation a pour enjeu de modéliser p (la proportion de 0) à l’aide de variables indépendantes, ces dernières pouvant se retrouver dans les deux parties du modèle. L’idée étant que, pour chaque observation, le modèle évalue sa probabilité d’être un faux zéro (partie binomiale), et le nombre attendu d’accidents.
Type de variable dépendante | Variable de comptage |
Distribution utilisée | Poisson avec excès de zéros |
Formulation |
\(Y \sim ZIP(\mu,\theta)\) \(g(\mu) = \beta_0 + \beta X\) \(s(p) = \alpha_0 + \alpha_X\) \(g(x) =log(x)\) \(s(x) = log(\frac{x}{1-x})\) |
Fonction de lien | log et logistique |
Paramètre modélisé | \(\mu\) et p |
Paramètres à estimer | \(\beta_0\), \(\beta\), \(\alpha_0\) et \(\alpha\) |
Conditions d’application | Absence de sur-dispersion |
8.3.4.1 Interprétation des paramètres
L’interprétation des paramètres \(\beta_0\) et \(\beta\) est identique à celle d’un modèle de Poisson. Les paramètres \(\alpha_0\) et \(\alpha\) sont identiques à ceux d’un modèle binomial. Plus spécifiquement, ces derniers paramètres modélisent la probabilité d’observer des valeurs supérieures à zéro.
8.3.4.2 Exemple appliqué
Nous avons vu, dans l’exemple précédent, que l’utilisation du modèle avec excès fixe de zéros pour les données d’accident des piétons aux intersections ne donnait pas de résultats satisfaisants. Nous tentons ici d’améliorer le modèle en ajoutant les variables indépendantes significatives du modèle Poisson dans la seconde équation de régression destinée à détecter les faux zéros.
Vérification des conditions d’application
Pour un modèle de Poisson avec excès de zéros, il n’est pas possible de calculer les distances de Cook. Nous devons donc directement passer à l’analyse des résidus simulés.
# Ajuster une première version du modèle
<- gamlss(formula = Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet +
modelza + Apaisement + LogEmploi + Densite_pop +
Terreplein + DensiteInter + Long_arterePS + Artere + NB_voies5,
Entropie sigma.formula = ~1 + Feux_auto + Feux_piet + Densite_pop + NB_voies5,
family = ZIP(mu.link = "log", sigma.link="logit"),
data = data_accidents)
## GAMLSS-RS iteration 1: Global Deviance = 1505.155
## GAMLSS-RS iteration 2: Global Deviance = 1488.658
## GAMLSS-RS iteration 3: Global Deviance = 1483.304
## GAMLSS-RS iteration 4: Global Deviance = 1482.085
## GAMLSS-RS iteration 5: Global Deviance = 1481.868
## GAMLSS-RS iteration 6: Global Deviance = 1481.832
## GAMLSS-RS iteration 7: Global Deviance = 1481.827
## GAMLSS-RS iteration 8: Global Deviance = 1481.825
## GAMLSS-RS iteration 9: Global Deviance = 1481.825
# Extraire la prédiction des valeurs lambda
<- predict(modelza, type = "response", what = "mu")
lambdas # Extraire la prédiction des valeurs p
<- predict(modelza, type = "response", what = "sigma")
ps # Calculer la combinaison de ces deux éléments
<- lambdas * ps
preds # Effectuer les 1000 simulations
<- 1000
nsim <- lapply(1:length(lambdas),function(i){
cols <- lambdas[[i]]
lambda <- ps[[i]]
p <- rZIP(n = nsim, mu = lambda, sigma = p)
sims return(sims)
})<- do.call(rbind, cols)
mat_sims # Calculer les résidus simulés
<- createDHARMa(simulatedResponse = mat_sims,
sim_res observedResponse = data_accidents$Nbr_acci,
fittedPredictedResponse = preds,
integerResponse = T)
plot(sim_res)
La figure 8.36 indique deux problèmes importants dans le modèle : la présence de valeurs aberrantes ainsi qu’un potentiel problème de dispersion. Nous commençons donc par identifier ces valeurs aberrantes.
# Identification des outliers
<- outliers(sim_res,return = "logical",lowerQuantile = 0.001,
isOutlier upperQuantile = 0.999)
<- subset(data_accidents, isOutlier)
cas_etrange print(cas_etrange)
## Nbr_acci Feux_auto Feux_piet Pass_piet Terreplein Apaisement EmpTotBuffer Densite_pop Entropie
## 1 19 1 1 1 1 0 7208.538 5980.923 0.8073926
## 2 13 1 1 1 0 0 15198.958 11696.805 0.6432018
## 5 12 1 0 1 0 0 8585.350 8655.430 0.7607852
## 26 7 0 0 1 0 0 1342.625 2751.012 0.0000000
## 483 0 0 0 1 0 0 2665.846 7208.356 0.6251996
## DensiteInter Long_arterePS Artere NB_voies5 log_acci catego_acci catego_acci2 Arret VAG sum_app
## 1 42.41597 6955.00 1 1 2.995732 1 1 0 1 4
## 2 86.50445 9581.51 1 0 2.639057 1 1 0 1 4
## 5 89.11495 6412.27 0 0 2.564949 1 1 0 1 4
## 26 73.35344 2849.66 0 0 2.079442 1 1 1 1 3
## 483 44.49561 1684.07 0 0 0.000000 0 0 1 1 4
## LogEmploi AccOrdinal PopHa
## 1 8.883021 2 5.980923
## 2 9.628982 2 11.696805
## 5 9.057813 2 8.655430
## 26 7.202382 2 2.751012
## 483 7.888277 0 7.208356
Nous retirons des données les quelques observations pouvant avoir une trop forte influence sur le modèle. Après réajustement, la figure 8.37 nous informe que nous n’avons plus de valeurs aberrantes restantes ni de fort problème de dispersion. En revanche, le premier quantile des résidus tant à être plus faible que ce que nous aurions pu nous attendre d’une distribution uniforme. Ce constat laisse penser que le modèle a du mal à bien identifier les faux zéros. Ce résultat n’est pas étonnant, car aucune variable n’avait été identifiée à cette fin dans l’article original (Cloutier et al. 2014) qui utilisait un modèle binomial négatif.
<- subset(data_accidents, isOutlier==FALSE)
data2 # Ajuster une première version du modèle
<- gamlss(formula = Nbr_acci ~ Feux_auto + Feux_piet + Pass_piet +
modelza + Apaisement + LogEmploi + Densite_pop +
Terreplein + DensiteInter + Long_arterePS + Artere + NB_voies5,
Entropie sigma.formula = ~1 + Feux_auto + Feux_piet + Densite_pop + NB_voies5,
family = ZIP(mu.link = "log", sigma.link="logit"),
data = data2)
## GAMLSS-RS iteration 1: Global Deviance = 1397.751
## GAMLSS-RS iteration 2: Global Deviance = 1386.151
## GAMLSS-RS iteration 3: Global Deviance = 1382.5
## GAMLSS-RS iteration 4: Global Deviance = 1381.285
## GAMLSS-RS iteration 5: Global Deviance = 1380.893
## GAMLSS-RS iteration 6: Global Deviance = 1380.777
## GAMLSS-RS iteration 7: Global Deviance = 1380.745
## GAMLSS-RS iteration 8: Global Deviance = 1380.736
## GAMLSS-RS iteration 9: Global Deviance = 1380.733
## GAMLSS-RS iteration 10: Global Deviance = 1380.732
# Extraire la prédiction des valeurs lambda
<- predict(modelza, type = "response", what = "mu")
lambdas # Extraire la prédiction des valeurs p
<- predict(modelza, type = "response", what = "sigma")
ps # Calculer la combinaison de ces deux éléments
<- lambdas * ps
preds # Effectuer les 1000 simulations
<- 1000
nsim <- lapply(1:length(lambdas),function(i){
cols <- lambdas[[i]]
lambda <- ps[[i]]
p <- rZIP(n = nsim, mu = lambda, sigma = p)
sims return(sims)
})<- do.call(rbind, cols)
mat_sims # Calculer les résidus simulés
<- createDHARMa(simulatedResponse = mat_sims,
sim_res observedResponse = data2$Nbr_acci,
fittedPredictedResponse = preds,
integerResponse = T)
plot(sim_res)
Nous pouvons une fois encore comparer des simulations issues du modèle et de la distribution originale de la variable Y. La figure 8.38 montre clairement que les simulations du modèle (en bleu) sont très éloignées dans la distribution originale (en gris), ce qui remet directement en question la pertinence de ce modèle.
# Extraire la prédiction des valeurs lambda
<- predict(modelza, type = "response", what = "mu")
lambdas # Extraire la prédiction des valeurs p
<- predict(modelza, type = "response", what = "sigma")
ps # Génération de 1000 simulations pour chaque prédiction
<- 1000
nsim <- lapply(1:length(lambdas),function(i){
cols <- lambdas[[i]]
lambda <- ps[[1]]
p <- round(rZIP(nsim,mu=lambda, sigma = p))
sims return(sims)
})<- do.call(rbind, cols)
mat_sims # Préparation des données pour le graphique (valeurs réelles)
<- data.frame(table(data2$Nbr_acci))
counts names(counts) <- c("nb_accident",'frequence')
$nb_accident <- as.numeric(as.character(counts$nb_accident))
counts$prop <- counts$frequence / sum(counts$frequence)
counts# Préparation des données pour le graphique (valeurs simulées)
<- data.frame(count = 0:25)
df1 <- lapply(1:nsim, function(i){
count_sims <- mat_sims[,i]
sim <- data.frame(table(sim))
cnt <- merge(df1,cnt, by.x="count", by.y = "sim", all.x = T, all.y=F)
df2 $Freq <- ifelse(is.na(df2$Freq),0,df2$Freq)
df2return(df2$Freq)
})<- do.call(cbind,count_sims)
count_sims <- data.frame(
df_sims val = 0:25,
med = apply(count_sims, MARGIN = 1, median),
lower = apply(count_sims, MARGIN = 1, quantile, probs = 0.025),
upper = apply(count_sims, MARGIN = 1, quantile, probs = 0.975)
)# Affichage du graphique
ggplot() +
geom_bar(aes(x=nb_accident, weight = frequence), width = 0.6, data = counts)+
geom_errorbar(aes(x = val, ymin = lower, ymax = upper),
data = df_sims, color = "blue", width = 0.6)+
geom_point(aes(x = val, y = med), color = "blue", size = 1.3, data = df_sims)+
scale_x_continuous(limits = c(-0.5,7), breaks = c(0:7))+
xlim(-1,12)+
labs(subtitle = "",
x = "nombre d'accidents",
y = "nombre d'occurrences")
Vérification la qualité d’ajustement
<- glm(Nbr_acci ~ 1,
modelenull family = poisson(link="log"),
data = data2)
# Calcul des R2
rsqs(loglike.full = logLik(modelza),
loglike.null = logLik(modelenull),
full.deviance = deviance(modelza),
null.deviance = deviance(modelenull),
nb.params = modelza$sigma.df + modelza$mu.df,
n = nrow(data2)
)
## $`deviance expliquee`
## [1] 0.08807561
##
## $`McFadden ajuste`
## 'log Lik.' 0.3484515 (df=18)
##
## $`Cox and Snell`
## 'log Lik.' 0.7955335 (df=18)
##
## $Nagelkerke
## 'log Lik.' 0.805948 (df=18)
# Calcul du RMSE
sqrt(mean((preds - data2$Nbr_acci)**2))
## [1] 2.576983
Le modèle avec excès de zéro ajusté ne parvient à expliquer que 11 % de la déviance totale. Il obtient toutefois des valeurs de R2 assez hautes (McFadden ajusté : 0,36, Nagerlkerke : 0,82). Son RMSE est très élevé (2,6), comparativement à celui que nous avons obtenu avec le modèle binomial négatif (1,9). Considérant ces éléments, ce modèle est nettement moins informatif que le modèle binomial négatif et ne devrait pas être retenu. Nous montrons tout de même ici comment interpréter ces résultats.
Interprétation des résultats
L’ensemble des coefficients du modèle sont accessibles avec la fonction summary
. Les coefficients dédiés à la partie Poisson (appelée Mu dans le résumé) doivent être analysés et interprétés de la même manière que s’ils provenaient d’un modèle de Poisson. Les coefficients appartenant à la partie logistique (appelé Sigma dans le résumé) doivent être analysés et interprétés de la même manière que s’ils provenaient d’un modèle logistique.
# Extraction des résultats
<- summary(modelza) base_table
## ******************************************************************
## Family: c("ZIP", "Poisson Zero Inflated")
##
## Call: gamlss(formula = Nbr_acci ~ Feux_auto + Feux_piet +
## Pass_piet + Terreplein + Apaisement + LogEmploi +
## Densite_pop + Entropie + DensiteInter + Long_arterePS +
## Artere + NB_voies5, sigma.formula = ~1 + Feux_auto +
## Feux_piet + Densite_pop + NB_voies5, family = ZIP(mu.link = "log",
## sigma.link = "logit"), data = data2)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: log
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.472e+00 5.483e-01 -4.508 8.21e-06 ***
## Feux_auto 5.628e-01 1.909e-01 2.948 0.003359 **
## Feux_piet 3.905e-01 1.102e-01 3.545 0.000431 ***
## Pass_piet 3.389e-01 1.926e-01 1.760 0.079087 .
## Terreplein -3.318e-01 1.667e-01 -1.991 0.047078 *
## Apaisement 2.164e-01 1.534e-01 1.411 0.158933
## LogEmploi 1.638e-01 7.644e-02 2.142 0.032668 *
## Densite_pop 7.841e-05 1.533e-05 5.116 4.51e-07 ***
## Entropie -1.308e-01 3.094e-01 -0.423 0.672799
## DensiteInter 3.677e-03 2.082e-03 1.766 0.078108 .
## Long_arterePS 8.103e-06 2.054e-05 0.395 0.693332
## Artere 3.142e-02 1.156e-01 0.272 0.785871
## NB_voies5 4.761e-01 1.055e-01 4.511 8.10e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: logit
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.1597263 0.5259636 2.205 0.02793 *
## Feux_auto -1.6591393 0.6038115 -2.748 0.00622 **
## Feux_piet 0.0653804 0.7645555 0.086 0.93189
## Densite_pop -0.0001252 0.0000532 -2.353 0.01903 *
## NB_voies5 -1.4651605 0.8694733 -1.685 0.09261 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 500
## Degrees of Freedom for the fit: 18
## Residual Deg. of Freedom: 482
## at cycle: 10
##
## Global Deviance: 1380.732
## AIC: 1416.732
## SBC: 1492.595
## ******************************************************************
# Multiplication par 1000 des coefficients de population
# (effet pour 1000 habitants)
8,1] <- 1000 * base_table[8,1]
base_table[8,2] <- 1000 * base_table[8,2]
base_table[17,1] <- 1000 * base_table[17,1]
base_table[17,2] <- 1000 * base_table[17,2]
base_table[# Multiplication par 1000 des coefficients de longueur artère
# (effet pour 1 km)
11,1] <- 1000 * base_table[11,1]
base_table[11,2] <- 1000 * base_table[11,2]
base_table[# Calcul des exponentiels des variables indépendantes
# et des intervalles de confiance
<- exp(base_table[,1])
expcoeff .5 <- exp(base_table[,1] - 1.96 * base_table[,2])
expcoeff2.5 <- exp(base_table[,1] + 1.96 * base_table[,2])
expcoeff97<- cbind(base_table, expcoeff, expcoeff2.5,expcoeff97.5)
base_table # Calculer une colonne indiquant le niveau de significativité
<- case_when(
sign 4] < 0.001 ~ "***",
base_table[,4] >= 0.001 & base_table[,4]<0.01 ~ "**",
base_table[,4] >= 0.01 & base_table[,4]<0.05 ~ "*",
base_table[,4] >= 0.05 & base_table[,4]<0.1 ~ ".",
base_table[,TRUE ~ ""
)# Arrondir à trois décimales
<- round(base_table,3)
base_table # Enlever les colonnes de valeurs de t et d'erreur standard
<- base_table[,c(1,4,5,6,7)]
base_table <- cbind(base_table, sign)
base_table # Remplacer les 0 dans la colonne pval
2] <- ifelse(base_table[,2]=="0","<0.001",base_table[,2])
base_table[,# Séparer le tout en deux tableaux
<- base_table[1:13,]
part_poiss <- base_table[14:18,]
part_logit # Mettre les bons noms de colonnes
colnames(part_poiss) <- c("Coeff.","Val.p","Exp(Coeff.)",
"IC 2,5 % exp(Coeff.)","IC 97,5 % exp(Coeff.)", "Sign.")
colnames(part_logit) <- c("Coeff.","Val.p","RC","IC 2,5 % RC","IC 97,5 % RC", "Sign.")
Nous rapportons les résultats de ce modèle de Poisson avec excès de zéro ajusté dans les tableaux 8.27 et 8.28.
Variable | Coeff. | Val.p | Exp(Coeff.) | IC 2,5 % exp(Coeff.) | IC 97,5 % exp(Coeff.) | Sign. |
---|---|---|---|---|---|---|
(Intercept) | -2.472 | <0.001 | 0.084 | 0.029 | 0.247 | *** |
Feux_auto | 0.563 | 0.003 | 1.756 | 1.208 | 2.552 | ** |
Feux_piet | 0.391 | <0.001 | 1.478 | 1.191 | 1.834 | *** |
Pass_piet | 0.339 | 0.079 | 1.403 | 0.962 | 2.047 | . |
Terreplein | -0.332 | 0.047 | 0.718 | 0.518 | 0.995 |
|
Apaisement | 0.216 | 0.159 | 1.242 | 0.919 | 1.677 | |
LogEmploi | 0.164 | 0.033 | 1.178 | 1.014 | 1.368 |
|
Densite_pop | 0.078 | <0.001 | 1.082 | 1.05 | 1.115 | *** |
Entropie | -0.131 | 0.673 | 0.877 | 0.478 | 1.609 | |
DensiteInter | 0.004 | 0.078 | 1.004 | 1 | 1.008 | . |
Long_arterePS | 0.008 | 0.693 | 1.008 | 0.968 | 1.05 | |
Artere | 0.031 | 0.786 | 1.032 | 0.823 | 1.294 | |
NB_voies5 | 0.476 | <0.001 | 1.61 | 1.309 | 1.98 | *** |
Variable | Coeff. | Val.p | RC | IC 2,5 % RC | IC 97,5 % RC | Sign. |
---|---|---|---|---|---|---|
(Intercept) | 1.16 | 0.028 | 3.189 | 1.138 | 8.941 |
|
Feux_auto | -1.659 | 0.006 | 0.19 | 0.058 | 0.621 | ** |
Feux_piet | 0.065 | 0.932 | 1.068 | 0.239 | 4.777 | |
Densite_pop | -0.125 | 0.019 | 0.882 | 0.795 | 0.979 |
|
NB_voies5 | -1.465 | 0.093 | 0.231 | 0.042 | 1.27 | . |
Nous observons ainsi que la présence d’un feu de circulation divise par 5 les chances de rentrer dans la catégorie d’intersection où des accidents peuvent se produire. De même, la densité de population réduit les chances de passer dans cette catégorie de 11 %.
Concernant les coefficients pour la partie Poisson du modèle, nous observons que les présences d’un feu de circulation et d’un feu pour piéton contribuent à multiplier respectivement par 2 et 1,5 le nombre attendu d’accidents à une intersection. De même, la présence d’un axe de circulation à cinq voies augmente de 57 % le nombre d’accidents. Enfin, la densité de population est aussi associée à une augmentation du nombre d’accidents : pour 1 000 habitants supplémentaires autour de l’intersection, nous augmentons le nombre d’accidents attendu de 9 %.
8.3.5 Conclusion sur les modèles destinés à des variables de comptage
Dans cette section, nous avons vu que modéliser une variable de comptage ne doit pas toujours être réalisé avec une simple distribution de Poisson. Il est nécessaire de tenir compte de la sur ou sous-dispersion potentielle ainsi que de l’excès de zéros. Nous n’avons cependant pas couvert tous les cas. Il est en effet possible d’ajuster des modèles avec une distribution binomiale négative avec excès de zéros (avec le package gamlss
), ainsi que des modèles de Hurdle. Ces derniers ont une approche différente de celle proposée par les distributions ajustées pour tenir compte de l’excès de zéro que nous détaillons dans l’encadré « pour aller plus loin » ci-dessous. Le processus de sélection du modèle peut être résumé avec la figure 8.39. Notez que même en suivant cette procédure, rien ne garantit que votre modèle final reflète bien les données que vous étudiez. L’analyse approfondie des résidus et des prédictions du modèle est la seule façon de déterminer si oui ou non le modèle est fiable.
Modèle de Hurdle versus modèle avec excès de zéro
Les modèles de Hurdle sont une autre catégorie de modèles GLM. Ils peuvent être décrits avec la formulation suivante :
\[\begin{equation}\left\{\begin{array}{c} Y \sim \text {Binomial}(p) \text { si } y=0 \\ \text { logit }(p)=\beta_{a} X_{a} \\ Y \sim \text { TrPoisson}(\lambda) \text { si } y>0 \\ \log (\lambda)=\beta_{b} X_{b} \end{array}\right. \tag{8.21} \end{equation}\]
Nous constatons qu’un modèle de Hurdle utilise deux distributions, la première est une distribution binomiale dont l’objectif est de prédire si les observations sont à 0 ou au-dessus de 0. La seconde est une distribution strictement positive (supérieure à 0), il peut s’agir d’une distribution tronquée de Poisson, tronquée binomiale négative, Gamma, log-normale ou autre, dépendamment du phénomène modélisé. Puisque le modèle fait appel à deux distributions, deux équations de régression sont utilisées, l’une pour prédire p (la probabilité d’observer une valeur au-dessus de 0) et l’autre l’espérance (moyenne) de la seconde distribution.
En d’autres termes, un modèle de Hurdle modélise les données à zéro et les données au-delà de 0 comme deux processus différents (chacun avec sa propre distribution). Cette approche se distingue des modèles avec excès de zéros qui utilisent une seule distribution pour décrire l’ensemble des données. D’après un modèle avec excès de zéro, il existe de vrais et de faux zéros que l’on tente de distinguer. Dans un modèle de Hurdle, l’idée est que les zéros constituent une limite. Nous modélisons la probabilité de dépasser cette limite et ensuite la magnitude du dépassement de cette limite.
Prenons un exemple pour rendre la distinction plus concrète. Admettons que nous utilisons un capteur capable de mesurer la concentration de particules fines dans l’air. D’après les spécifications du fabricant, le capteur est capable de mesurer des taux de concentration à partir de 0,001 µg/m3. Dans une ville avec des niveaux de concentration très faibles, il est très fréquent que le capteur enregistre des valeurs à zéro. Considérant ce phénomène, il serait judicieux de modéliser le processus avec un modèle de Hurdle Gamma puisque les 0 représentent une limite qui n’a pas été franchie : le seuil de détection du capteur. Nous traitons donc différemment les secteurs au-dessous et au-dessus de ce seuil. Si nous reprenons notre exemple sur les accidents des piétons à des intersections, il est plus judicieux, dans ce cas, de modéliser le phénomène avec un modèle avec excès de zéro puisque nous pouvons observer zéro accident à une intersection dangereuse (vrai zéro) et zéro accident à une intersection sur laquelle aucun piéton ne traverse jamais (faux zéro).