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 = 312

La 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 ' ' 1

Notez 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
##   Conditional ICC: 0.176
ICC1 <- 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.08

10.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
##   Conditional ICC: 0.139
ICC2 <- 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.2683329
r2_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 = 78866

10.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
##   Conditional ICC: 0.110
ICC3 <- 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.269979
r2_nakagawa(Modele3)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.270
##      Marginal R2: 0.160

10.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
##   Conditional ICC: 0.108
ICC4 <- 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.270394
r2_nakagawa(Modele4)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.270
##      Marginal R2: 0.163

10.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.270

Vous 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 + 
## Modele2:     DuTriPct.c + NoLog.c + Setback.c + (1 | SRNOM)
## Modele3: PCTArb ~ Width.c + Length.c + AgeMed.c + AgeMed2.c + ResiPCT.c + 
## Modele3:     DuTriPct.c + NoLog.c + Setback.c + ValLog + UDipPCT + PCTFRAVI + 
## Modele3:     PCTIMGRE + AvecEnf + FranPCT + (1 | SRNOM)
## Modele4: PCTArb ~ Width.c + Length.c + AgeMed.c + AgeMed2.c + ResiPCT.c + 
## Modele4:     DuTriPct.c + NoLog.c + Setback.c + ValLog + UDipPCT + PCTFRAVI + 
## Modele4:     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