Et avec R ?

La syntaxe de lme4

Pour ajuster les modèles étudiés avec les données de l’étude sur la privation de sommeil, nous avons utilisé lmer(). Cette fonction du package lme4 (Bates et al. 2015) reprend la syntaxe de lm() et permet d’ajouter des effets aléatoires entre parenthèses.

La forme générale est (termes | groupe) : à gauche de la barre, on indique ce qui varie ; à droite, selon quel groupe. Voici quelques écritures utiles :

Formule Modèle
y ~ x + (1 | g) Intercept aléatoire par groupe g
y ~ x + (x | g) Intercept et pente de x aléatoires, corrélés
y ~ x + (x || g) Intercept et pente aléatoires, non corrélés
y ~ x + (0 + x | g) Pente aléatoire seule

Un aide-mémoire plus complet à la fin de cette section (aide-mémoire). On utilise aussi le package lmerTest (Kuznetsova et al. 2017), qui remplace lmer() par une version identique MAIS qui ajoute des p-values pour les effets fixes (pratique).

Inférence statistique

Tests statistiques

Note

Par défaut, lmer() estime les paramètres par maximum de vraisemblance restreinte (REML). C’est le choix recommandé lorsque l’on souhaite estimer et rapporter un modèle. Si l’on veut comparer des modèles qui diffèrent par leurs effets fixes, il faut utiliser le maximum de vraisemblance (ML). Avec anova(m1, m2), les modèles sont automatiquement réajustés par ML. On peut empêcher ce réajustement avec refit = FALSE, ce qui est légitime lorsque seuls les effets aléatoires diffèrent.

Une fois notre modèle estimé, nous pouvons vouloir tester un effet fixe. Avec un test de Wald, on compare la statistique \(t = \hat\beta / \mathrm{SE}(\hat\beta)\) à une loi de Student. Encore faut-il connaître le nombre de degrés de liberté à utiliser !

Dans un modèle mixte, l’information vient à la fois des observations et des groupes. Le “bon” nombre de degrés de liberté n’est connu exactement que dans des plans équilibrés simples. C’est pourquoi les auteurs de lme4 ont choisi de ne pas donner de p-values pour ces tests. Le package lmerTest, lui, propose deux approximations :

  • Satterthwaite, utilisée par défaut : elle est rapide et satisfaisante dans la plupart des cas ;
  • Kenward-Roger : elle corrige aussi l’erreur standard et peut être plus précise en petit échantillon, mais elle coûte davantage en calcul et nécessite le package pbkrtest.

Comparons-les pour notre modèle :

m_randInterceptSlope <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy) 

summary(m_randInterceptSlope)$coefficients # Satterthwaite
             Estimate Std. Error       df   t value     Pr(>|t|)
(Intercept) 251.40510   6.824597 16.99973 36.838090 1.171558e-17
Days         10.46729   1.545790 16.99998  6.771481 3.263824e-06
summary(m_randInterceptSlope, ddf = "Kenward-Roger")$coefficients # Kenward-Roger
             Estimate Std. Error df   t value     Pr(>|t|)
(Intercept) 251.40510   6.824597 17 36.838090 1.171003e-17
Days         10.46729   1.545790 17  6.771481 3.263808e-06

Si l’on souhaite tester un facteur à plusieurs modalités ou des interactions, anova(m_randInterceptSlope) produit une table de tests F. Avec lmerTest, ces tests sont de type III par défaut.

ImportantDéfinition

Le test du rapport de vraisemblance (LRT) compare deux modèles emboîtés : le modèle réduit correspond à l’hypothèse nulle, tandis que le modèle complet correspond à une hypothèse plus générale.

Pour tester Days, nous comparons le modèle complet à un modèle sans Days. Comme leurs effets fixes diffèrent, anova() les réajuste par ML (section 5.2) :

m_rs_0 <- lmer(Reaction ~ 1 + (Days | Subject), data = sleepstudy)
anova(m_rs_0, m_randInterceptSlope)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
m_rs_0 5 1785.476 1801.441 -887.7379 1775.476 NA NA NA
m_randInterceptSlope 6 1763.939 1783.097 -875.9697 1751.939 23.53654 1 1.2e-06

Intervalles de confiance

Pour accompagner les estimations, nous pouvons aussi calculer des intervalles de confiance. La fonction confint() propose trois méthodes :

  • method = "profile" est la méthode par défaut. Elle utilise le profil de vraisemblance, est fiable et fournit aussi des intervalles pour les écarts-types et les corrélations des effets aléatoires ;
  • method = "Wald" utilise \(\hat\beta \pm 1{,}96\,\mathrm{SE}\). Elle est rapide, mais fournit des intervalles uniquement pour les effets fixes ;
  • method = "boot" utilise un bootstrap paramétrique avec nsim simulations. Elle est souvent la plus fiable en petit échantillon si le modèle est bien spécifié, mais elle est lente.

Calculons les intervalles par profil pour m_randInterceptSlope :

ic_sleep <- confint(m_randInterceptSlope, oldNames = FALSE)
ic_sleep
                                   2.5 %      97.5 %
sd_(Intercept)|Subject        14.3813850  37.7158318
sd_Days|Subject                3.8011641   8.7533825
cor_Days.(Intercept)|Subject  -0.4815008   0.6849863
sigma                         22.8982669  28.8579965
(Intercept)                  237.6806955 265.1295147
Days                           7.3586533  13.5759188

Choisir les effets aléatoires

Revenons à l’étude sur la privation de sommeil : faut-il inclure une pente aléatoire ? Nous pouvons comparer m_randIntercept, avec seulement un intercept aléatoire, à m_randInterceptSlope, qui comporte aussi une pente aléatoire. Les deux modèles ont les mêmes effets fixes, donc un LRT par REML convient ; nous indiquons refit = FALSE pour garder cet ajustement :

m_randIntercept <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy) 
lrt_pente <- anova(m_randIntercept, m_randInterceptSlope, refit = FALSE)
lrt_pente
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
m_randIntercept 4 1794.465 1807.237 -893.2325 1786.465 NA NA NA
m_randInterceptSlope 6 1755.628 1774.786 -871.8141 1743.628 42.83681 2 0

Le modèle m_randInterceptSlope comporte deux paramètres supplémentaires : \(\tau_1\) pour l’écart-type de la pente et \(\rho\) pour sa corrélation avec l’intercept. C’est pourquoi la sortie indique 2 degrés de liberté. Mais, pour un test de variance, il y a un piège.

AvertissementTester une variance : le problème de la frontière

Sous l’hypothèse nulle \(H_0 : \tau_1^2 = 0\), la variance est sur la frontière de son espace : elle ne peut pas être négative. La statistique du LRT ne suit alors plus une loi du \(\chi^2\) à 2 degrés de liberté. Sous les conditions asymptotiques considérées ici, sa loi est un mélange à parts égales de \(\chi^2_1\) et \(\chi^2_2\) (Self et Liang 1987; Stram et Lee 1994). La p-value “naïve” est donc conservatrice, c’est-à-dire trop grande.

Pour tester un intercept aléatoire seul, le mélange asymptotique est \(\frac12\chi^2_0 + \frac12\chi^2_1\) : la p-value correspondante est la moitié de la p-value naïve.

Calculons les deux p-values pour voir ce que cela change ici :

lr <- lrt_pente$Chisq[2]
p_melange <- 0.5 * pchisq(lr, df = 1, lower.tail = FALSE) + 0.5 * pchisq(lr, df = 2, lower.tail = FALSE)
c(statistique = lr, p_melange = p_melange)
 statistique    p_melange 
4.283681e+01 2.792530e-10

Ici, \(p<0,05\) donc nous rejetons l’hypothèse nulle et concluons que la pente aléatoire semble nécessaire.

Si le modèle comporte plusieurs termes aléatoires, la fonction ranova() de lmerTest automatise leur comparaison en les retirant tour à tour :

ranova(m_randInterceptSlope)
npar logLik AIC LRT Df Pr(>Chisq)
6 -871.8141 1755.628 NA NA NA
Days in (Days | Subject) 4 -893.2325 1794.465 42.83681 2 0

La ligne Days in (Days | Subject) teste la suppression de la pente aléatoire et de sa corrélation avec l’intercept. On retrouve donc le LRT à 2 degrés de liberté ci-dessus. Avec le modèle non corrélé m_randIntSlope_uncorrelated, ranova() teste séparément chaque variance, avec 1 degré de liberté pour chacune.

Aide-mémoire lme4, lmerTest et nlme

Syntaxe des effets aléatoires dans lme4

Syntaxe Signification
(1 | g) Intercept aléatoire par groupe g
(x | g) ou (1 + x | g) Intercept et pente de x aléatoires, corrélés
(x || g) Intercept et pente aléatoires, non corrélés (équivaut à (1 | g) + (0 + x | g))
(0 + x | g) Pente aléatoire seule, sans intercept aléatoire

Fonctions utiles

Fonction Rôle
fixef(m) Effets fixes \(\hat{\boldsymbol\beta}\)
ranef(m) Effets aléatoires prédits
coef(m) Coefficients propres à chaque groupe (\(\hat\beta + \hat u_i\))
VarCorr(m), sigma(m) Composantes de variance, écart-type résiduel
predict(m, re.form = NA) Prédictions marginales (effets fixes seuls)
predict(m, newdata, allow.new.levels = TRUE) Prédictions conditionnelles ; un groupe nouveau reçoit un effet aléatoire nul
confint(m, method = "profile") IC par profil (aussi "Wald", "boot")
anova(m1, m2) test du rapport de vraisemblance (réajuste par ML si nécessaire ; refit = FALSE pour garder REML)
anova(m, type = 3) (lmerTest) tests F des effets fixes
ranova(m) (lmerTest) tests du rapport de vraisemblance des termes aléatoires

nlme : effets aléatoires et structure résiduelle

nlme::lme(
    fixed = y ~ temps, random = ~ temps | sujet, data = donnees,
    correlation = nlme::corAR1(form = ~ temps | sujet), # ou corCAR1, corCompSymm, corSymm...
    weights = nlme::varIdent(form = ~ 1 | temps)
) # une variance par temps

nlme::gls(y ~ temps,
    data = donnees, # sans effet aléatoire
    correlation = nlme::corSymm(form = ~ 1 | sujet)
)

nlme::getVarCov(m, type = "marginal") # matrice de covariance marginale V_i
nlme::ACF(m, resType = "normalized") # autocorrélation des résidus

Les références

Bates, Douglas, Martin Mächler, Ben Bolker, et Steve Walker. 2015. « Fitting Linear Mixed-Effects Models Using lme4 ». Journal of Statistical Software 67 (1): 1‑48. https://doi.org/10.18637/jss.v067.i01.
Kuznetsova, Alexandra, Per B. Brockhoff, et Rune H. B. Christensen. 2017. « lmerTest Package: Tests in Linear Mixed Effects Models ». Journal of Statistical Software 82 (13): 1‑26. https://doi.org/10.18637/jss.v082.i13.
Self, Steven G., et Kung-Yee Liang. 1987. « Asymptotic Properties of Maximum Likelihood Estimators and Likelihood Ratio Tests Under Nonstandard Conditions ». Journal of the American Statistical Association 82 (398): 605‑10. https://doi.org/10.1080/01621459.1987.10478472.
Stram, Daniel O., et Jae Won Lee. 1994. « Variance Components Testing in the Longitudinal Mixed Effects Model ». Biometrics 50 (4): 1171‑77. https://doi.org/10.2307/2533455.