10.4 Mise en œuvre dans R
Pour mettre en œuvre des modèles multiniveaux avec une variable dépendante continue, nous utilisons la fonction lmer du package lme4. Pour d’autres distributions, nous pouvons utiliser la fonction glmer implémentant différentes familles de modèles GLM, notamment binomiale (modèle multiniveau logistique), gaussien, Gamma, inverse gaussien, Poisson, Quasi-poisson, etc. Comme pour les modèles GLMM, lorsque d’autres distributions sont nécessaires, il est possible d’utiliser le package gamlss.
10.4.1 Le modèle vide
Dans le code R ci-dessous, la syntaxe lmer(PCTArb ~ 1 + (1| SRNOM), data = Multiniveau) permet de construire le modèle vide avec la variable indépendante PCTArb et SRNOM comme variable définissant les groupes au niveau 2, soit les 312 secteurs de recensement. À titre de rappel, le modèle vide ne comprend aucune variable indépendante.
library("lme4")
library("MuMIn")
# chargement du jeu de données
load("data/multiniveau/dataArbres.RData")
# MODÈLE 1 : modèle vide (sans prédicteurs)
#------------------------------------------------------
# Écrire Y ~ 1 signifie que le modèle est vide
# 1| SRNOM : signifie que l'on fait varier la constante avec la variable SRNOM
Modele1 <- lmer(PCTArb ~ 1 + (1| SRNOM), data = Multiniveau)
# Nombre de groupes
cat("nombre de groupes =", length(unique(Multiniveau$SRNOM)))## nombre de groupes = 312La fonction summary(Modele1) permet d’afficher les résultats du modèle. Dans la section intitulée Random effects, la variance pour le niveau 2 (SRNOM (Intercept)) est de 19,82 contre 92,93 pour le niveau 1 (Residual). Le coefficient de corrélation intraclasse (ICC) est donc égal à \(\mbox{19,82 / (19,82+92,93)} \times \mbox{100 = 17,58}\)%.
# Résultats du modèle
summary(Modele1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: PCTArb ~ 1 + (1 | SRNOM)
## Data: Multiniveau
##
## REML criterion at convergence: 80299.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9413 -0.5295 -0.2235 0.2175 8.4695
##
## Random effects:
## Groups Name Variance Std.Dev.
## SRNOM (Intercept) 19.82 4.452
## Residual 92.93 9.640
## Number of obs: 10814, groups: SRNOM, 312
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.3373 0.2772 314.9183 26.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Notez qu’il est possible d’obtenir directement la valeur de l’ICC avec la fonction icc(Modele1) du package performance et les statistiques d’ajustement du modèle avec les fonctions logLik, AIC et BIC.
# Calcul de l'ICC (coefficient intraclasse)
performance::icc(Modele1)## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.176
## Unadjusted ICC: 0.176ICC1 <- performance::icc(Modele1)
cat("Part de la variance de la variable dépendante imputable au niveau 2 : ",
round(ICC1$ICC_adjusted*100,2), "%", sep="")## Part de la variance de la variable dépendante imputable au niveau 2 : 17.58%# Qualité d'ajustement du modèle
cat("Statistiques d'ajustement du modèle :",
"\n-2 Log V = ", -2*logLik(Modele1),
"\nAIC =", AIC(Modele1),
"\nBIC =", BIC(Modele1))## Statistiques d'ajustement du modèle :
## -2 Log V = 80299.22
## AIC = 80305.22
## BIC = 80327.0810.4.2 Modèle avec les variables indépendantes du niveau 1
Le second modèle consiste à introduire les variables indépendantes mesurées pour les tronçons de rue (niveau 1). Notez comment sont centrées préalablement les variables explicatives.
# Centrage des variables indépendantes
VINiv1 <- c("Width","Length","AgeMed","AgeMed2","ResiPCT","DuTriPct","NoLog","Setback")
for (e in VINiv1){
e.c <- paste(e, ".c", sep="")
Multiniveau[[e.c]] <- Multiniveau[[e]] - mean(Multiniveau[[e]])
}
# MODÈLE 2 : modèle avec les prédicteurs au niveau 1 (rues)
# ------------------------------------------------------
Modele2 <- lmer(PCTArb ~
# Variables indépendantes au niveau 1
Width.c+Length.c+AgeMed.c+AgeMed2.c+ResiPCT.c+DuTriPct.c+NoLog.c+Setback.c+
(1| SRNOM), data = Multiniveau)
# Résultats du modèle
summary(Modele2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: PCTArb ~ Width.c + Length.c + AgeMed.c + AgeMed2.c + ResiPCT.c + DuTriPct.c + NoLog.c + Setback.c + (1 | SRNOM)
## Data: Multiniveau
##
## REML criterion at convergence: 78763.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9065 -0.5536 -0.1941 0.2569 9.4205
##
## Random effects:
## Groups Name Variance Std.Dev.
## SRNOM (Intercept) 15.26 3.907
## Residual 80.32 8.962
## Number of obs: 10814, groups: SRNOM, 312
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.228e+00 2.479e-01 3.183e+02 29.151 < 2e-16 ***
## Width.c -1.292e-01 1.272e-02 1.075e+04 -10.160 < 2e-16 ***
## Length.c 1.085e-02 1.717e-03 1.073e+04 6.322 2.69e-10 ***
## AgeMed.c 1.103e+00 1.856e-01 1.080e+04 5.946 2.83e-09 ***
## AgeMed2.c -2.950e-04 4.791e-05 1.080e+04 -6.158 7.62e-10 ***
## ResiPCT.c 4.699e-02 3.466e-03 1.080e+04 13.558 < 2e-16 ***
## DuTriPct.c -1.299e-02 2.683e-03 1.070e+04 -4.842 1.30e-06 ***
## NoLog.c 1.473e-01 1.057e-02 1.080e+04 13.938 < 2e-16 ***
## Setback.c 2.018e-01 2.295e-02 1.080e+04 8.792 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Wdth.c Lngth. AgMd.c AgMd2. RsPCT. DTrPc. NoLg.c
## Width.c -0.003
## Length.c 0.011 -0.216
## AgeMed.c 0.000 0.010 -0.014
## AgeMed2.c 0.002 -0.011 0.013 -1.000
## ResiPCT.c 0.056 0.095 0.208 0.023 -0.024
## DuTriPct.c -0.010 0.022 0.086 -0.074 0.077 0.025
## NoLog.c -0.030 0.156 -0.785 -0.008 0.009 -0.269 -0.127
## Setback.c 0.048 -0.018 -0.146 0.007 -0.008 -0.014 0.035 0.038
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling# Calcul de l'ICC (coefficient intraclasse)
performance::icc(Modele2)## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.160
## Unadjusted ICC: 0.139ICC2 <- performance::icc(Modele2)
cat("Part de la variance de la variable dépendante ",
"\nimputable au niveau 2 : ", round(ICC2$ICC_adjusted*100,2), "%", sep="")## Part de la variance de la variable dépendante
## imputable au niveau 2 : 15.97%# Calcul des R2 conditionnel et marginal avec les fonctions
# r.squaredGLMM ou r2_nakagawa du package performance
r.squaredGLMM(Modele2)## R2m R2c
## [1,] 0.1292872 0.2683329r2_nakagawa(Modele2)## # R2 for Mixed Models
##
## Conditional R2: 0.268
## Marginal R2: 0.129# Qualité d'ajustement du modèle
cat("Statistiques d'ajustement du modèle",
"\n-2 Log L = ", -2*logLik(Modele2),
"\nAIC =", AIC(Modele2),
"\nBIC =", BIC(Modele2))## Statistiques d'ajustement du modèle
## -2 Log L = 78763.83
## AIC = 78785.83
## BIC = 7886610.4.3 Modèle avec les variables indépendantes aux niveaux 1 et 2
Le troisième modèle comprend à la fois les variables indépendantes mesurées aux deux niveaux (tronçons et secteurs de recensement).
# MODÈLE 3 : modèle complet avec les prédicteurs aux niveaux 1 et 2
# ------------------------------------------------------
Modele3 <- lmer(PCTArb ~
# Variables indépendantes au niveau 1
Width.c+Length.c+AgeMed.c+AgeMed2.c+ResiPCT.c+DuTriPct.c+NoLog.c+Setback.c+
# Variables indépendantes au niveau 2
ValLog+UDipPCT+PCTFRAVI+PCTIMGRE+AvecEnf+FranPCT+
(1| SRNOM), data = Multiniveau)
# Résultats du modèle
summary(Modele3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: PCTArb ~ Width.c + Length.c + AgeMed.c + AgeMed2.c + ResiPCT.c +
## DuTriPct.c + NoLog.c + Setback.c + ValLog + UDipPCT + PCTFRAVI + PCTIMGRE + AvecEnf + FranPCT + (1 | SRNOM)
## Data: Multiniveau
##
## REML criterion at convergence: 78742.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0461 -0.5558 -0.1939 0.2622 9.4190
##
## Random effects:
## Groups Name Variance Std.Dev.
## SRNOM (Intercept) 12.12 3.482
## Residual 80.35 8.964
## Number of obs: 10814, groups: SRNOM, 312
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -5.175e-01 3.227e+00 3.136e+02 -0.160 0.87271
## Width.c -1.319e-01 1.271e-02 1.076e+04 -10.375 < 2e-16 ***
## Length.c 1.076e-02 1.717e-03 1.073e+04 6.265 3.87e-10 ***
## AgeMed.c 1.097e+00 1.854e-01 1.078e+04 5.920 3.31e-09 ***
## AgeMed2.c -2.936e-04 4.785e-05 1.078e+04 -6.136 8.75e-10 ***
## ResiPCT.c 4.649e-02 3.482e-03 1.078e+04 13.352 < 2e-16 ***
## DuTriPct.c -1.268e-02 2.677e-03 1.061e+04 -4.737 2.20e-06 ***
## NoLog.c 1.478e-01 1.057e-02 1.079e+04 13.985 < 2e-16 ***
## Setback.c 1.944e-01 2.307e-02 1.079e+04 8.428 < 2e-16 ***
## ValLog 1.591e-02 3.856e-03 3.127e+02 4.126 4.75e-05 ***
## UDipPCT 1.405e-02 3.546e-02 3.292e+02 0.396 0.69221
## PCTFRAVI -8.837e-02 2.958e-02 3.282e+02 -2.988 0.00302 **
## PCTIMGRE 2.367e-01 4.860e-02 3.217e+02 4.870 1.76e-06 ***
## AvecEnf 5.778e-04 3.226e-02 3.141e+02 0.018 0.98572
## FranPCT 5.213e-02 1.638e-02 3.163e+02 3.183 0.00160 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling# Qualité d'ajustement du modèle
cat("Statistiques d'ajustement du modèle",
"\n-2 Log L = ", -2*logLik(Modele3),
"\nAIC =", AIC(Modele3), "\nBIC =", BIC(Modele3))## Statistiques d'ajustement du modèle
## -2 Log L = 78742.85
## AIC = 78776.85
## BIC = 78900.75# Calcul de l'ICC (coefficient intraclasse)
performance::icc(Modele3)## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.131
## Unadjusted ICC: 0.110ICC3 <- performance::icc(Modele3)
cat("Part de la variance de la variable dépendante ",
"\nimputable au niveau 2 : ", round(ICC3$ICC_adjusted*100,2), "%", sep="")## Part de la variance de la variable dépendante
## imputable au niveau 2 : 13.11%# Calcul des R2 conditionnel et marginal avec les fonctions
# r.squaredGLMM ou r2_nakagawa du package performance
r.squaredGLMM(Modele3)## R2m R2c
## [1,] 0.1598477 0.269979r2_nakagawa(Modele3)## # R2 for Mixed Models
##
## Conditional R2: 0.270
## Marginal R2: 0.16010.4.4 Modèle complet avec une interaction
Le quatrième modèle consiste à ajouter au modèle complet une interaction entre deux variables des deux niveaux.
# Variance d'interaction
Multiniveau$PCTFRAVI_Setback <- Multiniveau$PCTFRAVI * Multiniveau$Setback.c
# MODÈLE 4 : interaction aux deux niveaux
# ------------------------------------------------------
Modele4 <- lmer(PCTArb ~
# Variables indépendantes au niveau 1
Width.c+Length.c+AgeMed.c+AgeMed2.c+ResiPCT.c+DuTriPct.c+NoLog.c+Setback.c+
# Variables indépendantes au niveau 2
ValLog+UDipPCT+PCTFRAVI+PCTIMGRE+AvecEnf+FranPCT+
# Variable d'interaction
PCTFRAVI_Setback+
(1| SRNOM), data = Multiniveau)
# Résultats du modèle
summary(Modele4)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: PCTArb ~ Width.c + Length.c + AgeMed.c + AgeMed2.c + ResiPCT.c +
## DuTriPct.c + NoLog.c + Setback.c + ValLog + UDipPCT + PCTFRAVI +
## PCTIMGRE + AvecEnf + FranPCT + PCTFRAVI_Setback + (1 | SRNOM)
## Data: Multiniveau
##
## REML criterion at convergence: 78732.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0148 -0.5568 -0.1922 0.2598 9.4261
##
## Random effects:
## Groups Name Variance Std.Dev.
## SRNOM (Intercept) 11.83 3.439
## Residual 80.24 8.958
## Number of obs: 10814, groups: SRNOM, 312
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -9.484e-03 3.198e+00 3.132e+02 -0.003 0.99764
## Width.c -1.357e-01 1.273e-02 1.076e+04 -10.660 < 2e-16 ***
## Length.c 1.079e-02 1.716e-03 1.073e+04 6.289 3.32e-10 ***
## AgeMed.c 1.092e+00 1.852e-01 1.078e+04 5.894 3.88e-09 ***
## AgeMed2.c -2.923e-04 4.780e-05 1.078e+04 -6.114 1.00e-09 ***
## ResiPCT.c 4.608e-02 3.480e-03 1.078e+04 13.239 < 2e-16 ***
## DuTriPct.c -1.268e-02 2.674e-03 1.059e+04 -4.742 2.14e-06 ***
## NoLog.c 1.454e-01 1.057e-02 1.079e+04 13.747 < 2e-16 ***
## Setback.c 3.443e-03 4.759e-02 1.076e+04 0.072 0.94233
## ValLog 1.582e-02 3.820e-03 3.119e+02 4.141 4.46e-05 ***
## UDipPCT 9.404e-03 3.515e-02 3.288e+02 0.268 0.78920
## PCTFRAVI -7.883e-02 2.937e-02 3.322e+02 -2.684 0.00764 **
## PCTIMGRE 2.194e-01 4.829e-02 3.261e+02 4.543 7.82e-06 ***
## AvecEnf -7.063e-03 3.199e-02 3.136e+02 -0.221 0.82542
## FranPCT 4.987e-02 1.623e-02 3.161e+02 3.072 0.00231 **
## PCTFRAVI_Setback 8.169e-03 1.779e-03 1.047e+04 4.591 4.46e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling# Qualité d'ajustement du modèle
cat("Statistiques d'ajustement du modèle",
"\n-2 Log L = ", -2*logLik(Modele4),
"\nAIC =", AIC(Modele3), "\nBIC =", BIC(Modele4))## Statistiques d'ajustement du modèle
## -2 Log L = 78732.66
## AIC = 78776.85
## BIC = 78899.85# Calcul de l'ICC (coefficient intraclasse)
performance::icc(Modele4)## # Intraclass Correlation Coefficient
##
## Adjusted ICC: 0.128
## Unadjusted ICC: 0.108ICC4 <- performance::icc(Modele4)
cat("Part de la variance de la variable dépendante ",
"\nimputable au niveau 2 : ", round(ICC4$ICC_adjusted*100,2), "%", sep="")## Part de la variance de la variable dépendante
## imputable au niveau 2 : 12.85%# Calcul des R2 conditionnel et marginal avec les fonctions
# r.squaredGLMM ou r2_nakagawa du package performance
r.squaredGLMM(Modele4)## R2m R2c
## [1,] 0.1628372 0.270394r2_nakagawa(Modele4)## # R2 for Mixed Models
##
## Conditional R2: 0.270
## Marginal R2: 0.16310.4.5 Comparaison des quatre modèles
Pour comparer les modèles, nous utilisons habituellement les statistiques d’ajustement du modèle vues plus haut, soit le maximum de vraisemblance (−2 Log-likelihood), l’AIC, l’ICC et les R2 marginal et conditionnel.
c_logLik <- c(logLik(Modele1),logLik(Modele2),logLik(Modele3),logLik(Modele4))
ICC <- c(performance::icc(Modele1)$ICC_adjusted,
performance::icc(Modele2)$ICC_adjusted,
performance::icc(Modele3)$ICC_adjusted,
performance::icc(Modele4)$ICC_adjusted)
R2m <- c(r.squaredGLMM(Modele1)[1],
r.squaredGLMM(Modele2)[1],
r.squaredGLMM(Modele3)[1],
r.squaredGLMM(Modele4)[1])
R2c <- c(r.squaredGLMM(Modele1)[2],
r.squaredGLMM(Modele2)[2],
r.squaredGLMM(Modele3)[2],
r.squaredGLMM(Modele4)[2])
print(data.frame(
Modele = c("Modèle 1 (vide)",
"Modèle 2 (VI : niv. 1)",
"Modèle 3 (VI : niv. 1 et 2)",
"Modèle 4 (interaction niv. 1 et 2"),
dl = AIC(Modele1, Modele2, Modele3, Modele4)$df,
Moins2LogLik = round(-2*c_logLik,0),
AIC = round(AIC(Modele1, Modele2, Modele3, Modele4)$AIC,0),
ICC = round(ICC,4),
R2marg = round(R2m,3),
R2cond = round(R2c,3)
))## Modele dl Moins2LogLik AIC ICC R2marg R2cond
## 1 Modèle 1 (vide) 3 80299 80305 0.1758 0.000 0.176
## 2 Modèle 2 (VI : niv. 1) 11 78764 78786 0.1597 0.129 0.268
## 3 Modèle 3 (VI : niv. 1 et 2) 17 78743 78777 0.1311 0.160 0.270
## 4 Modèle 4 (interaction niv. 1 et 2 18 78733 78769 0.1285 0.163 0.270Vous constaterez ci-dessus que les valeurs d’AIC et de -2 log de vraisemblance diminuent des modèles 1 à 4, signalant une amélioration progressive des modèles. Cela se traduit aussi par une augmentation du R2 conditionnel incluant à la fois les effets fixes et aléatoires. Sans surprise, la valeur du coefficient de corrélation intraclasse diminue du modèle vide au modèle complet : plus nous ajoutons de variables dépendantes, plus la capacité explicative du niveau 2 diminue.
Il est également judicieux de vérifier si un modèle est significativement différent du modèle précédent avec la fonction anova qui compare les différences de leurs déviances. En guise d’exemple, la différence de déviance de 59 (\(\mbox{78 625}-\mbox{78 684}=\mbox{59}\)) entre les modèles 3 et 2 (modèle complet versus modèle GLMM) avec six degrés de liberté – puisque le modèle 3 inclut six variables indépendantes de plus que le précédent (\(\mbox{17}-\mbox{11}=\mbox{6}\)) – est significative (p < 0,001). Cela indique que le modèle 3 est plus performant que le précédent.
anova(Modele1, Modele2, Modele3, Modele4)## Data: Multiniveau
## Models:
## Modele1: PCTArb ~ 1 + (1 | SRNOM)
## Modele2: PCTArb ~ Width.c + Length.c + AgeMed.c + AgeMed2.c + ResiPCT.c + DuTriPct.c + NoLog.c + Setback.c + (1 | SRNOM)
## Modele3: PCTArb ~ Width.c + Length.c + AgeMed.c + AgeMed2.c + ResiPCT.c + DuTriPct.c + NoLog.c + Setback.c + ValLog + UDipPCT + PCTFRAVI + PCTIMGRE + AvecEnf + FranPCT + (1 | SRNOM)
## Modele4: PCTArb ~ Width.c + Length.c + AgeMed.c + AgeMed2.c + ResiPCT.c + DuTriPct.c + NoLog.c + Setback.c + ValLog + UDipPCT + PCTFRAVI + PCTIMGRE + AvecEnf + FranPCT + PCTFRAVI_Setback + (1 | SRNOM)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Modele1 3 80304 80326 -40149 80298
## Modele2 11 78706 78786 -39342 78684 1614.351 8 < 2.2e-16 ***
## Modele3 17 78659 78783 -39313 78625 59.131 6 6.758e-11 ***
## Modele4 18 78640 78771 -39302 78604 21.166 1 4.213e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1