Analyses multiniveaux avec R
Autoformation
S’autoformer en statistique via des tutoriels et des capsules vidéo
Analyses multiniveaux avec R
Ce tutoriel a pour objectif de vous initier aux analyses multiniveaux (ou MLM pour « multilevel modeling »), soit une méthode statistique qui permet de tester les relations entre des variables lorsque nos données ont une structure hiérarchique, comme c’est le cas par exemple avec des données longitudinales.
Les analyses de ce tutoriel seront exécutées sur R. Pour vous pratiquer, il suffit de suivre le code ci-dessous !
Qu’est-ce qu’un modèle multiniveau ?
Plusieurs types de données, notamment en sciences sociales, peuvent être structurées en fonction de différents niveaux (c.-à-d., groupes et/ou structures hiérarchiques). Par exemple (voir Figure 1), pensons à la satisfaction au travail des employés dans une organisation, variable qui semble être grandement influencée par les gens avec qui un employé travaille quotidiennement (collègues, gestionnaire, etc.). Partageant un microclimat organisationnel similaire, il serait plausible de penser que deux employés dans la même équipe de travail ou dans le même département seraient davantage similaires que deux employés choisis aléatoirement dans l’organisation. En effet, le fait d’appartenir au même groupe (« cluster » en anglais) créé une certaine dépendance entre les employés d’une même équipe de travail. Les analyses multiniveaux permettent justement de contrôler pour cette dépendance en nichant les observations dans leur groupe. Les modèles multiniveaux sont utilisés pour modéliser les liens entre les variables dépendantes et indépendantes en tenant compte des corrélations entre les observations d’un même groupe.
Figure 1. Exemple de modèle multiniveau avec différents groupes d’individus.
Les modèles multiniveaux sont également communs pour analyser des données longitudinales. Par exemple, si on cherche à étudier la performance des employés d’une organisation sur 3 ans, en prenant des mesures mi-annuelle – les différents scores de performance d’un même individu au cours des trois années sont corrélés entre eux. Une analyse multiniveau permettrait de nicher ces données intra-individuelle pour chaque employé respectif, et contrôlerait ainsi pour la dépendance de ces observations. L’analyse multiniveau permet également de distinguer la variance provenant des différents niveaux d’analyses, soit le temps (niveau 1) et l’individu (niveau 2). La variance comprend donc une composante intra-individuelle (variance des scores d’un même individu d’un temps à l’autre) et une composante inter-individuelle (variance des scores moyens d’un individu à l’autre).
Figure 2. Exemple de modèle multiniveau avec données longitudinales.
Effets fixes et effets aléatoires
Une particularité des analyses multiniveaux est qu’elles permettent de distinguer les effets fixes des effets aléatoires (variances d’effets). De façon générale, on assume que les prédicteurs sont fixes, ce qui signifie qu’ils représentent l’effet moyen pour tous les individus. Toutefois, les analyses peuvent permettre aux intercepts et aux pentes des prédicteurs au niveau 1 d’être aléatoires – c’est à quoi on réfère quand on parle d’effets aléatoires. Les effets aléatoires permettent d’observer la variance dans les effets.
Par exemple, reprenons notre exemple examinant la performance des employés sur une période de 3 ans avec des mesures mi-annuelles. Ici, l’effet fixe serait de dire qu’en moyenne, à la première mesure de performance, les employés ont une cote de performance de 70 et qu’à chaque évaluation (demi-année), la cote augmente de 5 points. On peut l’illustrer sous forme de régression : cote de performance = 5*demi-année + 70. Un effet aléatoire dans ce cas-ci indiquerait deux choses. D’abord, il peut avoir une variabilité dans la cote de performance de base de chaque employé (ici, on réfère à l’effet aléatoire de l’intercept, donc du score quand tous les prédicteurs sont à 0). Certains peuvent avoir une cote de base plus haute, d’autres peuvent en avoir une plus basse. L’effet aléatoire de l’intercept vient capter cette variabilité au baseline. Deuxièment, la progression de la performance des employés peuvent varier d’une personne à l’autre et cela sera également identifier par les effets aléatoires. Bien que certains peuvent augmenter leur cote de 5 points chaque 6 mois, d’autres peuvent l’augmenter de 8 points et certains peuvent la descendre de 6 points. L’effet aléatoire de la pente vient capter la variabilité dans la progression de chaque individu au cours des temps de mesure.
Quand utiliser les modèles multiniveaux ?
Les modèles multiniveaux sont les analyses de choix lorsque nos observations étudiées ne sont pas indépendantes les unes des autres. Donc, si un ou une chercheur(e) décide d’observer les comportements d’employés au sein de différents départements ou de les étudier longitudinalement, le devis de recherche comprend une dépendance entre les données (via le département ou via l’individu). Pour des raisons méthodologiques, il est donc plus adéquat d’utiliser des analyses multiniveaux afin de contrôler pour cette dépendance. De plus, les modèles multiniveaux sont appropriés lorsqu’un ou une chercheur(e) désire modéliser la variation à plusieurs niveaux dans un échantillon, comme pour déterminer quel pourcentage de la variance des scores est due au niveau 1 ou au niveau 2.
Comment compléter un modèle multiniveau étape par étape sur R ?
Comme pour toutes analyses complétées sur R, il faut d’abord télécharger et appeler le package comprenant les fonctions nécessaires pour réaliser les analyses. Dans notre cas, nous appelons le package lme4. Par le fait même, pour les fins de ce tutoriel, nous utiliserons un jeu de données inclus dans ce package. Le jeu de données « sleepstudy » observe le temps de réaction moyen (en millisecondes) par jour de participants manquant de sommeil. Dans ce jeu de données, le temps de réaction (Reaction) est la variable dépendante et est observée au niveau 1. Chaque participant a des entrées quotidiennes ce qui laisse comprendre que les jours (Days) sont également au niveau 1. Finalement, les participants (Subject) eux-mêmes (c.-à-d., les différences entre individus) sont au niveau 2.
Tel qu’illustré dans la Figure 3, le jeu de données « sleepstudy » a un format long, ce qui signifie qu’il y a une ligne (ou rangée) par observation (c.-à-d., temps de mesure) et plusieurs lignes par personne, soit notre variable de niveau 2.
Figure 3. Jeu de données “sleepstudy”.
#### Téléchargement du package et des données ####
install.packages("lme4")
library(lme4)
df<-sleepstudy
View(df)
1) Corrélation intraclasse: décomposition de la variance
La prochaine étape est de décomposer la variance de notre variable dépendante, le temps de réaction, afin d’examiner quelle partie de la variance est intra-individuelle (et donc varie selon les jours de mesures) et quelle partie de la variance est inter-individuelles (et donc due aux différences individuelles). Cette étape permet de mieux comprendre ce que notre modèle doit expliquer. Pour ce faire, nous testerons le modèle nul qui comprend un terme fixe, soit la moyenne de toutes les observations du temps de réaction. Le modèle nul nous permettra également de calculer la variance à chaque niveau.
#### Modèle nul: moyenne de VD quand tous les facteurs sont à 0, décomposition de la variance ####
#### Temps de réaction divisé en intra vs inter individuelle
# null model
formule_null <- Reaction ~ 1 + (1|Subject)
model_null <- lmer(formule_null, data=df, subset=Days>=2)
summary(model_null)
Nous avons d’abord indiqué la formule pour le modèle nul. Comme pour toute formule d’analyse multiniveau, la première partie (celle qui n’est pas entre parenthèses) est semblable à une formule de régression et représente les effets fixes de notre modèle. Dans le cas d’un modèle nul, aucun prédicteur n’est inclus puisqu’on cherche simplement à décomposer la variance de notre variable dépendante. Ainsi, nous écrivons Reaction ~ 1, où le 1 représente l’intercept (et donc, la moyenne) du temps de réaction. La partie entre parenthèses, (1|Subject) représente l’effet aléatoire de l’intercept, qui varie pour chaque « Subject ». Dans le cas de notre modèle nul, on indique donc que le temps de réaction varie par individu.
La deuxième étape est d’exécuter l’analyse multiniveau avec la fonction lmer, dans laquelle nous devons indiquer la formule établie à l’étape précédente (formule_null), le jeu de données utilisé (data=df), et, au besoin, un sous-échantillon à analyser. Dans notre exemple, puisque les jours 0 et 1 étaient des jours de pratiques, nous échantillonnerons seulement les jours 2 et plus pour les analyses (subset=Days>=2).
Suite à ce premier modèle, les résultats dans la section des effets fixes suggèrent que la moyenne du temps de réaction parmi toutes les observations est M = 307.99 ms. Au niveau des effets aléatoires, les résultats indiquent une variance inter-individuelle (niveau 2) de 1650 (SD = 40.62) et une variance intra-individuelle (niveau 1) de 1690 (SD = 41.11). À la lumière de ces résultats, nous pouvons maintenant calculer l’ICC (pour « interclass correlation »), coefficient qui indique la corrélation entre les observations d’un même individu (ou groupe). Plus cette corrélation est élevée, plus les données d’une même personne sont corrélées entre elles.
# icc = var inter / (var inter + var intra)
icc <- 1650/(1650+1690)
icc
Dans notre cas, l’ICC est de .494 – cela indique 49% de la variance dans le temps de réaction est due à l’individu (c.-à-d., expliquée par les différences individuelles).
2) Intercept aléatoire, pente fixe
Nous pouvons maintenant ajouter un prédicteur, soit les différents jours de collecte (Days) tout en maintenant une pente (« slope ») fixe.
#### Ajouter les jours comme prédicteurs, pente fixe ####
formule_fixeds <- Reaction ~ 1 + Days + (1|Subject)
model_fixeds <- lmer(formule_fixeds, data=df, subset=Days>=2)
summary(model_fixeds)
Pour clarifier la formule, la section Reaction ~ 1 + Days est similaire à la formule d’une régression normale et établit la pente fixe des jours. La seconde partie (1|Subject) représente les effets aléatoires (« random effects »), ici un intercept qui varie entre personnes. Dans le cas suivant, on indique que l’intercept (représenté par le 1) est aléatoire. La partie à droite de la barre verticale | représente la variable de niveau 2, soit les individus, à travers lesquels les effets varient. Pour faciliter la compréhension, voici un exemple graphique illustrant une pente fixe et des intercepts aléatoires : tous les individus ont la même pente (i.e., pente fixe), mais chaque individu a son propre intercept (i.e., intercept aléatoire). La ligne pointillée indique les effets fixes, soit la trajectoire moyenne.
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1746.9 41.80
Residual 913.1 30.22
Number of obs: 144, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 245.097 11.829 20.72
Days 11.435 1.099 10.40
Regardons d’abord les résultats des effets fixes. En lien avec l’intercept, celui-ci indique qu’au Jour 0, le temps de réaction moyen est de 245.1 ms. Également, à chaque augmentation d’une journée, le temps de réaction augmente de 11.44 ms. On peut l’imaginer comme ceci : Reaction = 11.44(Days) + 245.1.
Au niveau des effets aléatoires, les résultats indique qu’au Jour 0, le temps de réaction a une variance inter-individuelle de 1746.9 (SD = 41.80). Cette variance représente la variabilité de temps de réaction qu’ont les participants au Jour 0.
3) Intercept aléatoire, pente aléatoire
Nous pouvons maintenant tester le modèle avec l’intercept et la pente aléatoire.
#### Jours comme prédicteur, pente aléatoire ####
formule_randoms <- Reaction ~ 1 + Days + (1+Days|Subject)
model_randoms <- lmer(formule_randoms, data=df, subset=Days>=2)
summary(model_randoms)
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 992.69 31.507
Days 45.77 6.766 -0.25
Residual 651.59 25.526
Number of obs: 144, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 245.097 9.260 26.468
Days 11.435 1.845 6.197
L’interprétation des résultats au niveau des effets fixes est la même que celle mentionnée plus haut (Reaction = 11.44*Days + 245.1). Les effets aléatoires indiquent qu’au Jour 0, le temps de réaction a une variance inter-individuelle de 992.69 (SD = 31.51). De plus, la progression du temps de réaction à travers les jours (c.-à-d., la pente de Days) a une variance de 45.77 (SD = 6.77). La corrélation de r = -.25 entre l’intercept et la pente indique que plus les participants ont un temps de réaction élevé au Jour 0, moins leur temps de réaction tend à augmenter au fil des jours.
4) Comparaison des modèles à pente fixe vs aléatoire
Il nous reste maintenant à voir si notre modèle avec la pente aléatoire est mieux ajusté avec les données en comparaison avec celui où la pente est fixe. Pour ce faire, nous utiliserons la fonction anova().
#### Comparaison des deux modèles: la pente aléatoire est-elle plus adéquate? ####
anova(model_randoms, model_fixeds)
model_fixeds: Reaction ~ 1 + Days + (1 | Subject)
model_randoms: Reaction ~ 1 + Days + (1 + Days | Subject)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
model_fixeds 4 1446.5 1458.4 -719.25 1438.5
model_randoms 6 1425.2 1443.0 -706.58 1413.2 25.332 2 3.156e-06 ***
Le test significatif nous indique que de modéliser une pente aléatoire mène effectivement à un meilleur fit que s’il est fixe. Au niveau de l’interprétation, on en comprend donc que les gens varient non seulement au Jour 0, mais que la progression de leur temps de réaction au fil des jours est également significativement différente.
Syntaxe complète
#### Téléchargement du package et des données ####
install.packages("lme4")
library(lme4)
df<-sleepstudy
View(df)
#### Modèle nul: moyenne de VD quand tous les facteurs sont à 0, décomposition de la variance ####
#### Temps de réaction divisé en intra vs inter individuelle
# null model
formule_null <- Reaction ~ 1 + (1|Subject)
model_null <- lmer(formule_null, data=df, subset=Days>=2)
summary(model_null)
# icc = var inter / (var inter + var intra)
icc <- 1650/(1650+1690)
icc
#### Ajouter les jours comme prédicteurs, pente fixe ####
formule_fixeds <- Reaction ~ 1 + Days + (1|Subject)
model_fixeds <- lmer(formule_fixeds, data=df, subset=Days>=2)
summary(model_fixeds)
#### Jours comme prédicteur, pente aléatoire ####
formule_randoms <- Reaction ~ 1 + Days + (1+Days|Subject)
model_randoms <- lmer(formule_randoms, data=df, subset=Days>=2)
summary(model_randoms)
#### Comparaison des deux modèles: la pente aléatoire est-elle plus adéquate? ####
anova(model_randoms, model_fixeds)
PARTAGEZ VOS IMPRESSIONS & POSEZ VOS QUESTIONS ICI !