10 Régressions multiniveaux

Dans le précédent chapitre, nous avons abordé les modèles à effets mixtes qui permettent d’introduire à la fois des effets fixes et aléatoires (GLMM). Dans ce chapitre, nous poursuivons sur cette voie avec une nouvelle extension des modèles GLM : les modèles multiniveaux. Ces modèles sont simplement une extension des modèles à effets mixtes et permettent de modéliser un phénomène avec une structure hiérarchique des données, tel que décrit dans le chapitre précédent.

Rappel de la structure hiérarchique des données

Exemple à deux niveaux : il s’agit de modéliser un phénomène \(y_{ij}\), soit une variable dépendante Y pour un individu i (niveau 1) niché dans un groupe j (niveau 2). Par exemple, modéliser l’indice de masse corporelle (IMC) de 5000 individus résidant dans 100 quartiers différents.

Exemple à trois niveaux : il s’agit de modéliser un phénomène \(y_{ijk}\), soit une variable dépendante Y pour un individu i (niveau 1), niché dans un groupe j (niveau 2) appartenant à un groupe k (niveau 3). Par exemple, modéliser les notes à un examen de mathématiques d’élèves (niveau 1) nichés dans des classes (niveau 2) nichées dans des écoles (niveau 3).

Nous avons largement décrit précédemment trois principaux types de modèles d’effets mixtes (GLMM) :

Les GLMM avec constantes aléatoires qui permettent d’avoir une constante différente pour chacun des groupes (niveau 2).
Les GLMM avec pentes aléatoires qui permettent de faire varier une variable indépendante au niveau 1 (coefficient) en fonction des groupes au niveau 2.
Les GLMM avec constantes et pentes aléatoires.

Les modèles multiniveaux se différencient des modèles à effets mixtes puisqu’ils permettent d’introduire des variables indépendantes mesurées aux niveaux supérieurs (2 et 3).

Liste des packages utilisés dans ce chapitre

lme4 pour en œuvre des modèles multiniveaux avec une variable dépendante continue.
performance pour obtenir le coefficient intraclasse (ICC).
MuMIn pour obtenir les pseudos R².

10.1 Modèles multiniveaux : deux intérêts majeurs

Les modèles multiniveaux ont deux principaux avantages : analyser la répartition de la variance entre les différents niveaux et introduire des variables explicatives aux différents niveaux du modèle.

10.1.1 Répartition de la variance entre les différents niveaux

Les modèles multiniveaux permettent d’estimer comment se répartit la variance entre les différents niveaux du jeu de données. Dans les deux exemples de l’encadré précédent, ils permettraient de répondre aux questions suivantes :

Quel niveau explique le plus l’IMC, le niveau individuel (niveau 1) ou le niveau contextuel (niveau 2)?
Comment se répartit la variance des notes à l’examen de mathématiques entre les trois niveaux? A-t-on plus de variance pour les individus (niveau 1) ou au sein des classes (niveau 2) ou entre les différentes écoles (niveau 3)?

10.1.2 Estimation des coefficients aux différents niveaux

Les modèles multiniveaux permettent d’estimer simultanément les coefficients de plusieurs variables indépendantes introduites à chacun des niveaux du modèle. Autrement dit, de voir comment les variables indépendantes introduites aux différents niveaux influencent la variable dépendante (Y) mesurée au niveau 1. Si nous reprenons l’exemple à trois niveaux (élèves/classes/écoles), plusieurs facteurs peuvent influencer la réussite ou la performance scolaire des élèves aux différents niveaux :

Variables indépendantes au niveau 1 (élève) : âge, sexe, statut socioéconomique, langue maternelle autre que la langue d’enseignement…
Variables indépendantes au niveau 2 (classe) : nombre d’élèves par classe, programme spécialisé ou pas…
Variables indépendantes au niveau 3 (école) : indice de défavorisation de l’école, école publique ou privée, qualité des infrastructures de l’école (bâtiment, gymnase, cour d’école)…

Dans la même veine, afin d’illustrer l’apport des modèles multiniveaux dans le champ de la géographie de la santé, Philibert et Apparicio (2007, 129) signalent que « pour un modèle à deux niveaux, il s’agit de modéliser \(y_{ij}\), par exemple l’IMC d’un individu i (niveau 1) résidant dans un quartier j (niveau 2). Il est alors possible de mettre des variables explicatives tant au niveau 1 (âge, sexe, revenu, niveau d’éducation, etc.) qu’au niveau 2 (niveau de défavorisation sociale du quartier, offre de services et d’équipements sportifs et récréatifs, caractéristiques de l’environnement urbain, etc.). Dans cet exemple, nous pouvons voir comment la modélisation multiniveau permet d’estimer simultanément les effets environnementaux et individuels de manière à distinguer la contribution de chacun des niveaux (ex. : l’effet du revenu des individus et celui de la défavorisation du quartier) dans l’explication des variations géographiques observées ».

Évaluer les effets de milieu avec des analyses multiniveaux

En santé des populations et en études urbaines, les modèles multiniveaux sont largement mobilisés pour évaluer les effets de milieu (neighbourhoods effects ou area effects en anglais).

Atkinson et Kintrea (2001, 2278) définissent les « effets de milieu comme le changement net dans les potentialités de l’existence (life chances) attribuable au fait de vivre dans un quartier (ou une zone) plutôt qu’un autre » [traduction libre]. Les effets de milieu peuvent être positifs ou négatifs et concerner aussi bien les enfants que les adultes.

Les analyses multiniveaux sont particulièrement adaptées à l’évaluation des effets de milieu. En effet, plusieurs phénomènes — état de santé, comportement ou choix individuels — peuvent être influencés à la fois par des caractéristiques individuelles (âge, sexe, niveau de revenu, niveau d’éducation, etc.) et par des caractéristiques contextuelles (caractéristiques du quartier).

Avec un modèle multiniveau, une fois contrôlées les caractéristiques individuelles (variables indépendantes mesurées au niveau 1), il est alors possible d’évaluer l’effet des caractéristiques du quartier (variables indépendantes mesurées au niveau 2) sur un phénomène \(y_{ij}\) mesuré pour un individu i résidant dans un quartier j.

10.2 Différents types de modèles multiniveaux

10.2.1 Description du jeu de données utilisé

Dans le cadre de cette section, nous présentons uniquement les modèles à deux niveaux, soit celui pour modéliser un phénomène \(y_{ij}\). Pour ce faire, nous utilisons des données tirées d’une étude de Pham et al. (2017). Dans cet article, les auteurs souhaitent évaluer les effets des caractéristiques de la forme urbaine et des caractéristiques socioéconomiques sur la couverture des arbres de rue, et ce, à partir d’un modèle multiniveau. Ils disposent ainsi d’une structure hiérarchique de données avec deux niveaux : les tronçons de rue (niveau 1, n = 10 814) inclus dans un et un seul secteur de recensement (niveau 2, n = 312). La variable dépendante (\(y_{ij}\)) est le pourcentage de la superficie du tronçon de rue qui est couverte par des arbres, calculé à partir d’images satellites à haute résolution (Quickbird, 60 cm, septembre 2008). L’ensemble des variables utilisées pour les modèles sont reportées au tableau 10.1.

Sept variables indépendantes relatives à la forme urbaine sont mesurées pour les tronçons de rue, soit la largeur et la longueur de la rue, l’âge médian des bâtiments (introduit également au carré pour vérifier l’existence d’un effet curvilinéaire; voir la section 7.5.1.1), les pourcentages de bâtiments résidentiels, de duplex et de triplex, le nombre de bâtiments et finalement la distance moyenne entre le bâtiment et la rue. Les variables indépendantes pour les 312 secteurs de recensement (niveau 2) sont extraites du recensement canadien de 2006 (tableau 10.1).

Tableau 10.1: Statistiques descriptives pour les variables des modèles multiniveaux
Nom	Intitulé	Type	Niveau	Moy.	Écart type
PCTArb	Arbres sur le tronçon de rue (%)	VD	1	7,2	10,7
Width	Largeur des rues	VI	1	16,0	7,3
Length	Longueur de rues	VI	1	136,0	87,8
AgeMed	Âge médian des bâtiments	VI	1	1 952,7	28,3
ResiPCT	Bâtiments résidentiels (%)	VI	1	83,5	28,0
DuTriPct	Duplex ou triplex (%)	VI	1	41,8	39,3
NoLog	Nombre de bâtiments	VI	1	14,0	14,4
Setback	Distance entre le bâtiment et la rue	VI	1	7,2	4,3
ValLog	Valeur moyenne des logements (milliers de dollars)	VI	2	267,6	80,0
UDipPCT	Diplômés universitaires (%)	VI	2	16,9	9,6
PCTFRAVI	Personnes à faible revenu (%)	VI	2	30,3	11,5
PCTIMGRE	Immigrants récents (%)	VI	2	10,0	7,3
AvecEnf	Ménages avec enfants (%)	VI	2	34,8	12,6
FranPCT	Langue maternelle française (%)	VI	2	66,9	24,1

10.2.2 Démarche classique pour les modèles multiniveaux

La démarche habituelle en analyse multiniveau est de réaliser plusieurs modèles, allant du plus simple au plus complexe. Cette stratégie permet habituellement de bien cerner la répartition de la variance entre les différents niveaux et l’apport des variables explicatives introduites aux différents niveaux. De la sorte, cinq types de modèles peuvent être construits :

Le modèle vide (appelé aussi modèle inconditionnel) qui comprend des constantes aléatoires au niveau 2, mais aucune variable explicative.
Le modèle avec uniquement les variables indépendantes au niveau 1 et des constantes aléatoires au niveau 2.
Le modèle complet avec les variables indépendantes aux deux niveaux et des constantes aléatoires.
Le modèle complet avec les variables indépendantes aux deux niveaux, incluant une interaction entre une variable indépendante mesurée au niveau 1 et une autre mesurée au niveau 2.
Le modèle avec les variables indépendantes aux deux niveaux et des constantes et pentes aléatoires.

Dans les sous-sections suivantes, nous détaillons les quatre premiers modèles en prenant soin de montrer les similitudes qu’ils partagent avec les modèles à effets mixtes vus précédemment. Notez d’emblée que les trois premiers modèles sont les plus fréquemment utilisés.

10.2.2.1 Modèle vide

Comme son nom l’indique, le modèle vide ne comprend aucune variable explicative. Il consiste simplement à faire varier la constante du niveau 1 avec des effets aléatoires au niveau 2, ce qui explique qu’il est souvent comparé à une ANOVA avec des effets aléatoires. En d’autres termes, ce modèle correspond à un GLMM avec constantes aléatoires dans lequel aucune variable indépendante n’est incluse au niveau 1. D’ailleurs, si vous comparez l’équation 10.1 avec l’équation 9.2 au chapitre précédent, vous constaterez que seul le paramètre \(\beta_1 x_1\) a été ôté et qu’il comprend aussi deux variances : l’une fixe au niveau 1 (\(\sigma_e\)) et l’autre aléatoire (stochastique) au niveau 2 (\(\sigma_{\upsilon}\)).

\[ \begin{aligned} &Y \sim Normal(\mu,\sigma_e)\\ &g(\mu) = \beta_0 + \upsilon \\ &\upsilon \sim Normal(0, \sigma_{\upsilon}) \\ &g(x) = x \end{aligned} \tag{10.1}\]

Quel est alors l’intérêt de réaliser un modèle si simple? À partir des deux variances, il est possible de calculer le coefficient de corrélation intraclasse (intraclass-correlation (ICC) en anglais) qui est le rapport entre la variance aléatoire et la somme des variances des deux niveaux, soit fixe et aléatoire (équation 10.2). Ce coefficient varie ainsi de 0 à 1 et indique la proportion de la variance de la variable dépendante qui est imputable au niveau 2. Tous(tes) les auteur(e)s s’entendent sur le fait qu’il est impératif de commencer une analyse de multiniveau en calculant ce modèle vide qui nous informe de la répartition de la variance entre les deux niveaux (Raudenbush et Bryk 2002; Gelman et Hill 2006; Tabachnick, Fidell et Ullman 2007; Bressoux 2010). Nous pourrons ensuite analyser l’évolution de ce coefficient dans les modèles subséquents.

\[ \rho = \frac{\sigma_{\upsilon}}{\sigma_{\upsilon} + \sigma_{e}} \tag{10.2}\]

Les résultats du modèle vide (inconditionnel) à partir des données de Pham et al. (2017) sont présentés au tableau 10.2. La variance du niveau 1 est de 92,93 contre 19,82 au niveau 2. Le coefficient de corrélation intraclasse vaut alors : \(\mbox{19,82} / \mbox{(19,82} + \mbox{92,93}) = \mbox{0,1758}\). Cela signifie que 18 % de la variance de la variable dépendante sont imputables au niveau 2 (des secteurs de recensement) et 82 % au niveau 1 (des tronçons). Nous verrons comment évolue ce coefficient dans les modèles subséquents.

Tableau 10.2: Résultats du modèle vide (modèle 1)
Paramètre	Coefficient	Erreur type	Valeur de T
Effets fixes (niveau 1)
Constante	7,337	0,277	26,466
Répartition de la variance
Variance (niveau 1)	19,818
Variance (niveau 2)	92,925
Coefficient de corrélation intraclasse	0,176
Qualité d’ajustement du modèle
AIC	80 305,219
R² marginal	0,000
R² conditionnel	0,176

10.2.2.2 Modèle avec les variables indépendantes du niveau 1

Dans ce second modèle, nous introduisons uniquement les variables explicatives au niveau 1. Par conséquent, ce modèle est tout simplement un modèle à effets mixtes (GLMM) avec des constantes aléatoires largement décrit à la section 9.2.1). Si vous comparez l’équation du modèle vide (équation 10.1) avec l’équation de ce modèle (équation 10.3), vous constaterez que le paramètre \(\beta_1 x_1\) a été ajouté. Il correspond au coefficient pour la variable indépendante \(X_1\) mesurée au niveau 1 (effet fixe). Nous pourrions alors ajouter d’autres paramètres pour les autres variables indépendantes du modèle, soit \(\beta_1 x_1 + \beta_2 x_2+ \ldots + \beta_k x_k\) (\(k\) étant le nombre de variables explicatives mesurées au niveau 1, effets fixes).

\[ \begin{aligned} &Y \sim Normal(\mu,\sigma_e)\\ &g(\mu) = \beta_0 + \beta_1 x_1 + \upsilon \\ &\upsilon \sim Normal(0, \sigma_{\upsilon})) \\ &g(x) = x \end{aligned} \tag{10.3}\]

Les résultats du second modèle sont présentés au tableau 10.3.

La répartition de la variance entre les deux niveaux. La variance du niveau 1 est désormais de 15,263 contre 80,317 au niveau 2, ce qui permet d’obtenir un coefficient de corrélation intraclasse de 0,1597. Cela signifie que de 16 % de la variance de la variable dépendante sont imputables au niveau 2 (des secteurs de recensement), une fois contrôlées les caractéristiques des tronçons.

La qualité d’ajustement du modèle. Dans le chapitre précédent sur les GLMM, nous avons largement décrit plusieurs mesures de la qualité d’ajustement du modèle, notamment l’AIC et les R² marginal et conditionnel. À titre de rappel, voici comment interpréter ces mesures :

Plus la valeur de l’AIC est faible, mieux le modèle est ajusté. En comparant les valeurs d’AIC du modèle vide et du modèle avec les variables explicatives du niveau 1 (80 305 versus 78 785), nous constatons, sans surprise, que ce dernier modèle est plus performant.
Le R² marginal indique la proportion de la variance expliquée uniquement si les effets fixes sont pris en compte (ici, 0,129). Quant au R² conditionnel, il indique la proportion de la variance expliquée à la fois par les effets fixes et aléatoires (ici, 0,268). L’écart important entre les deux R² signale que les secteurs de recensement (effets aléatoires, niveau 2) jouent un rôle important dans le modèle.

Quelles informations peut-on tirer des coefficients de régression? Les variables indépendantes relatives à la forme urbaine les plus importantes sont : le pourcentage de bâtiments résidentiels (ResiPCT), la largeur de la rue (Width) et le nombre de bâtiments (NoLog). Aussi, la distance entre le bâtiment et la rue (Setback) est associée positivement avec la variable dépendante. En effet, à chaque ajout d’un mètre de la distance moyenne entre les bâtiments et le tronçon de rue, la couverture des arbres sur le tronçon augmente de 0,202 point de pourcentage, toutes choses étant égales par ailleurs.

Tableau 10.3: Résultats du modèle avec les variables indépendantes au niveau 1 (modèle 2)
Paramètre	Coefficient	Erreur type	Valeur de T
Effets fixes (niveau 1)
Constante	-1 028,618	179,736	-5,723
Width	-0,129	0,013	-10,160
Length	0,011	0,002	6,322
AgeMed	1,103	0,186	5,946
AgeMed2	0,000	0,000	-6,158
ResiPCT	0,047	0,003	13,558
DuTriPct	-0,013	0,003	-4,842
NoLog	0,147	0,011	13,938
Setback	0,202	0,023	8,792
Répartition de la variance
Variance (niveau 1)	15,263
Variance (niveau 2)	80,317
Coefficient de corrélation intraclasse	0,160
Qualité d’ajustement du modèle
AIC	78 785,827
R² marginal	0,129
R² conditionnel	0,268

Remarquez la valeur de la constante : −1028,618. À titre de rappel, la constante est la valeur que prend la variable dépendante quand toutes les variables indépendantes sont égales à 0. Or, il est impossible qu’elles soient toutes égales à zéro.

Centrage des variables quantitatives mesurées au niveau 1

En analyse multiniveau, il est très courant et souvent recommandé de centrer les variables explicatives quantitatives au niveau 1. Deux options sont alors possibles :

Pour une variable indépendante donnée, les observations sont centrées sur leur moyenne générale, c’est-à-dire la moyenne de l’ensemble des observations du jeu de données, soit \(X_{ij} - \bar{X}\). Dans ce cas, la constante est donc la valeur que prend la variable \(Y\) quand toutes les variables indépendantes sont égales à leur moyenne respective.
Chaque observation est centrée sur la moyenne de son groupe respectif, soit \(X_{ij} - \bar{X}_{.j}\).

Comme signalé par Bressoux (2010, 328), « dans le premier cas, la variance des pentes sera estimée pour l’individu moyen dans la distribution générale, tandis que dans le second elle est estimée pour l’individu moyen de chaque groupe ».

Autrement dit, comparativement à un modèle sans centrage, les valeurs des coefficients pour les variables indépendantes sont les mêmes dans le premier cas (seule la valeur de la constante va changer) tandis qu’elles sont différentes dans le second cas.

Attention, il ne faut pas appliquer de centrage sur une variable qualitative, qu’elle soit dichotomique, nominale ou ordinale.

Pourquoi la pratique du centrage en analyse multiniveau est si courante?

Dans la plupart des livres sur les régressions multiniveaux, le centrage est recommandé, notamment dans l’ouvrage classique de Raudenbush et Bryk (2002). Rappelons que ces modèles sont largement utilisés en éducation avec une structure hiérarchique classique élève/école/classe. Nous nous intéressons alors à l’individu moyen (l’élève), ce qui explique que le centrage est habituellement appliqué. Par exemple, ne pas centrer l’âge des élèves fait que la constante qui est obtenue est peu interprétable : difficile d’évaluer la note moyenne à un examen quand la variable âge de l’élève a la valeur de 0, tout comme les autres variables explicatives quantitatives relatives à l’élève.

Centrage et réduction de l’ensemble des variables du modèle

Il est à noter que certains auteurs centrent et réduisent l’ensemble des variables du modèle. À titre de rappel, le centrage consiste à soustraire à chaque valeur la moyenne de la variable; la réduction, à la diviser par l’écart-type de la variable (section 2.5.5.2). Pour chaque variable, la moyenne est alors égale à 0 et l’écart-type à 1. Les coefficients s’interprètent alors en termes d’augmentation d’une unité d’écart-type tant pour la VI que la VD. Ils correspondent alors à des coefficients de régression standardisés (abordés dans la section 7.4.2). Ce processus de centrage et de réduction des variables peut être motivé par des problèmes de convergence du modèle (lorsque l’algorithme d’optimisation n’arrive pas à trouver une solution pour produire les coefficients).

Par conséquent, nous vous proposons de centrer les variables du niveau 1 de notre jeu de données. Si vous comparez les tableaux 10.3 et 10.4, vous constaterez que les valeurs relatives aux coefficients, aux mesures de la répartition de la variance et à la qualité d’ajustement du modèle sont les mêmes. Seule la valeur de la constante change : elle passe de −1028,618 à 7,228. Elle s’interprète désormais de la façon suivante : si toutes les variables explicatives sont égales à leurs moyennes respectives, alors le pourcentage de la superficie du tronçon couverte par des arbres est égal à 7,228 %.

Tableau 10.4: Résultats du modèle avec les variables indépendantes centrées au niveau 1 (modèle 2)
Paramètre	Coefficient	Erreur type	Valeur de T
Effets fixes (niveau 1)
Constante	7,228	0,248	29,151
Width.c	-0,129	0,013	-10,160
Length.c	0,011	0,002	6,322
AgeMed.c	1,103	0,186	5,946
AgeMed2.c	0,000	0,000	-6,158
ResiPCT.c	0,047	0,003	13,558
DuTriPct.c	-0,013	0,003	-4,842
NoLog.c	0,147	0,011	13,938
Setback.c	0,202	0,023	8,792
Répartition de la variance
Variance (niveau 1)	15,263
Variance (niveau 2)	80,317
Coefficient de corrélation intraclasse	0,160
Qualité d’ajustement du modèle
AIC	78 785,827
R² marginal	0,129
R² conditionnel	0,268

10.2.2.3 Modèle complet avec les variables indépendantes aux niveaux 1 et 2

Le troisième type de modèle consiste à introduire à la fois les variables explicatives mesurées au niveau 1 et au niveau 2 (équation 10.4). Il est communément appelé le modèle complet. Les variables explicatives du niveau 2 sont aussi considérées comme des effets fixes.

\[ \begin{aligned} &Y \sim Normal(\mu,\sigma_e)\\ &g(\mu) = \underbrace{\beta_0 + \beta_1 x_1}_{\mbox{effets fixes (niveau 1)}}+ \underbrace{\beta_2 z_2}_{\mbox{effets fixes (niveau 2)}}+\epsilon \\ &\upsilon \sim Normal(0, \sigma_{\upsilon}) \\ &g(x) = x \end{aligned} \tag{10.4}\]

Les résultats du troisième modèle sont présentés au tableau 10.5. Ce modèle permet d’évaluer les effets des caractéristiques socioéconomiques (mesurés au niveau des secteurs de recensement) sur la couverture des arbres des îlots, une fois contrôlées les caractéristiques de la forme urbaine des tronçons. Rappelons que dans ce modèle, les constantes sont aléatoires et les variables indépendantes au niveau 1 sont centrées.

Quelles informations peut-on tirer des coefficients de régression du niveau 2? D’emblée, deux caractéristiques n’ont pas d’effet significatif sur la variable dépendante, soit les pourcentages de diplômés universitaires et de ménages avec enfants. Par contre, toutes choses étant égales par ailleurs, la valeur moyenne des logements et le pourcentage d’immigrants récents sont associés à une augmentation de la couverture végétale. À l’inverse, le pourcentage de personnes à faible revenu est associé à une diminution de la couverture végétale.

Tableau 10.5: Résultats du modèle avec les variables indépendantes aux niveaux 1 et 2 (modèle 3)
Paramètre	Coefficient	Erreur type	Valeur de T
Effets fixes (niveau 1 : tronçons)
Constante	-0,518	3,227	-0,160
Width.c	-0,132	0,013	-10,375
Length.c	0,011	0,002	6,265
AgeMed.c	1,097	0,185	5,920
AgeMed2.c	0,000	0,000	-6,136
ResiPCT.c	0,046	0,003	13,352
DuTriPct.c	-0,013	0,003	-4,737
NoLog.c	0,148	0,011	13,985
Setback.c	0,194	0,023	8,428
Effets fixes (niveau 2 : secteurs de recensement)
ValLog	0,016	0,004	4,126
UDipPCT	0,014	0,035	0,396
PCTFRAVI	-0,088	0,030	-2,988
PCTIMGRE	0,237	0,049	4,870
AvecEnf	0,001	0,032	0,018
FranPCT	0,052	0,016	3,183
Répartition de la variance
Variance (niveau 1)	12,121
Variance (niveau 2)	80,347
Coefficient de corrélation intraclasse	0,131
Qualité d’ajustement du modèle
AIC	78 776,845
R² marginal	0,160
R² conditionnel	0,270

10.2.2.4 Modèle avec une interaction entre deux niveaux

Dans la section 7.5.4, nous avons vu comment introduire des variables d’interaction dans une régression linéaire multiple, soit entre deux variables continues (section 7.5.4.1), soit entre une variable continue et une variable dichotomique (section 7.5.4.2), soit entre deux variables dichotomiques (section 7.5.4.3). En analyse multiniveau, il peut être pertinent d’introduire une interaction entre une variable mesurée au niveau 1 et une autre mesurée au niveau 2 (équation 10.5).

\[ \begin{aligned} &Y \sim Normal(\mu,\sigma_e)\\ &g(\mu) = \underbrace{\beta_0 + \beta_1 x_1}_{\mbox{effets fixes (niv. 1)}}+ \underbrace{\beta_2 z_2}_{\mbox{effets fixes (niv. 2)}}+ \underbrace{\beta_3 (x_1 \times z_2)}_{\mbox{interaction (niv. 1 et 2)}}+ \epsilon \\ &\upsilon \sim Normal(0, \sigma_{\upsilon})) \\ &g(x) = x \end{aligned} \tag{10.5}\]

Dans le tableau 10.6, nous introduisons une variable d’interaction entre la distance entre le bâtiment et la rue (Setback.c) et le pourcentage de personnes à faible revenu (PCTFRAVI). On constate alors que PCTFRAVI est associé négativement avec la variable dépendante (\(\beta =\) −0,079, t = −2,684). Toutefois, lorsqu’elle est mise en interaction avec la variable Setback.c, cette variable est significative et positive (\(\beta =\) 0,008, t = 4,591).

Tableau 10.6: Résultats du modèle avec une variable d’interaction entre les deux niveaux 1 et 2 (modèle 4)
Paramètre	Coefficient	Erreur type	Valeur de T
Effets fixes (niveau 1 : tronçons)
Constante	-0,009	3,198	-0,003
Width.c	-0,136	0,013	-10,660
Length.c	0,011	0,002	6,289
AgeMed.c	1,092	0,185	5,894
AgeMed2.c	0,000	0,000	-6,114
ResiPCT.c	0,046	0,003	13,239
DuTriPct.c	-0,013	0,003	-4,742
NoLog.c	0,145	0,011	13,747
Setback.c	0,003	0,048	0,072
Effets fixes (niveau 2 : secteurs de recensement)
ValLog	0,016	0,004	4,141
UDipPCT	0,009	0,035	0,268
PCTFRAVI	-0,079	0,029	-2,684
PCTIMGRE	0,219	0,048	4,543
AvecEnf	-0,007	0,032	-0,221
FranPCT	0,050	0,016	3,072
Variable d’interaction (niv. 1 et 2)
Setback X PCTFRAVI	0,008	0,002	4,591
Répartition de la variance
Variance (niveau 1)	11,829
Variance (niveau 2)	80,239
Coefficient de corrélation intraclasse	0,128
Qualité d’ajustement du modèle
AIC	78 768,659
R² marginal	0,163
R² conditionnel	0,270

10.3 Conditions d’application des régressions multiniveaux

Puisque les modèles multiniveaux sont une extension des modèles à effets mixtes (GLMM), nous retrouvons globalement les mêmes conditions d’application (voir la section 9.3), dont les principales sont :

l’absence de multicolinéarité excessive,
la normalité des résidus,
l’absence d’observations trop influentes dans le modèle.

Combien de groupes au niveau 2? Dans la section 9.3, nous avons vu que dans un modèle GLMM, plusieurs auteur(e)s, notamment Gelman et Hill (2006), préconisent un minimum de cinq groupes dans un modèle à effets mixtes. Toutefois, dans un modèle complet d’une régression multiniveau, nous introduisons aussi des variables indépendantes au niveau 2. Par conséquent, le nombre de groupes doit être augmenté significativement, et ce, idéalement proportionnellement au nombre des variables indépendantes ajoutées au niveau 2. En ce sens, Bressoux (2010, 325) conseille d’avoir au moins 10 groupes pour chaque variable indépendante mesurée au niveau 2. Toujours selon Bressoux (2010, 325), certains auteurs recommandent même 20 groupes par variable indépendante au niveau 2. En conséquence, bien qu’aucune règle de pouce soit clairement admise, il est clair qu’un modèle complet multiniveau nécessite un nombre de groupes conséquent.

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 ['lmerMod']
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 t value
(Intercept)   7.3373     0.2772   26.47

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
  Unadjusted 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 ['lmerMod']
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 t value
(Intercept)  7.228e+00  2.479e-01  29.151
Width.c     -1.292e-01  1.272e-02 -10.160
Length.c     1.085e-02  1.717e-03   6.322
AgeMed.c     1.103e+00  1.856e-01   5.946
AgeMed2.c   -2.950e-04  4.791e-05  -6.158
ResiPCT.c    4.699e-02  3.466e-03  13.558
DuTriPct.c  -1.299e-02  2.683e-03  -4.842
NoLog.c      1.473e-01  1.057e-02  13.938
Setback.c    2.018e-01  2.295e-02   8.792

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.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 ['lmerMod']
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 t value
(Intercept) -5.175e-01  3.227e+00  -0.160
Width.c     -1.319e-01  1.271e-02 -10.375
Length.c     1.076e-02  1.717e-03   6.265
AgeMed.c     1.097e+00  1.854e-01   5.920
AgeMed2.c   -2.936e-04  4.785e-05  -6.136
ResiPCT.c    4.649e-02  3.482e-03  13.352
DuTriPct.c  -1.268e-02  2.677e-03  -4.737
NoLog.c      1.478e-01  1.057e-02  13.985
Setback.c    1.944e-01  2.307e-02   8.428
ValLog       1.591e-02  3.856e-03   4.126
UDipPCT      1.405e-02  3.546e-02   0.396
PCTFRAVI    -8.837e-02  2.958e-02  -2.988
PCTIMGRE     2.367e-01  4.860e-02   4.870
AvecEnf      5.778e-04  3.226e-02   0.018
FranPCT      5.213e-02  1.638e-02   3.183
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.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 ['lmerMod']
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 t value
(Intercept)      -0.0094840  3.1984732  -0.003
Width.c          -0.1356543  0.0127250 -10.660
Length.c          0.0107911  0.0017159   6.289
AgeMed.c          1.0915688  0.1851983   5.894
AgeMed2.c        -0.0002923  0.0000478  -6.114
ResiPCT.c         0.0460769  0.0034803  13.239
DuTriPct.c       -0.0126785  0.0026735  -4.742
NoLog.c           0.1453667  0.0105748  13.747
Setback.c         0.0034430  0.0475943   0.072
ValLog            0.0158168  0.0038196   4.141
UDipPCT           0.0094041  0.0351462   0.268
PCTFRAVI         -0.0788297  0.0293729  -2.684
PCTIMGRE          0.2193705  0.0482873   4.543
AvecEnf          -0.0070626  0.0319907  -0.221
FranPCT           0.0498673  0.0162308   3.072
PCTFRAVI_Setback  0.0081695  0.0017794   4.591
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.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 R² 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 R² 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

Vrai
Faux

Vrai
Faux

Vrai
Faux