Questions diverses sur la fonction lmer

[Fabien Vinckier]

Bonjour à tous,

A la demande d'un reviewer, je suis sensé faire une analyse en modèle mixte (auxquels, je dois bien l'avouer, je ne comprends pas grand chose, malgré quelques heures passées sur ce précieux forum et sur le draft de D. Bates)... J'aurais donc plusieurs questions à ce sujet et vous remercie d'avance pour votre aide!

1/ une première question sans doute un peu révélatrice de ma non-compréhension du sujet...

J'utilise un modèle de type

aga1 <- lmer(lrt ~A + B + C + (1 | subj) + (1 | word) , data=u23)

, et j'obtiens ceci :

Code : Tout sélectionner

 Data: u23right
   AIC   BIC logLik deviance REMLdev
 -5279 -5237   2647    -5335   -5293
Random effects:
 Groups   Name        Variance   Std.Dev.
 word     (Intercept) 0.00043719 0.020909
 subj     (Intercept) 0.00480310 0.069304
 Residual             0.00829540 0.091079
Number of obs: 2840, groups: word, 146; subj, 30

Fixed effects:
                    Estimate Std. Error t value
(Intercept)        2.8384507  0.0129920  218.48
A                     0.0052163  0.0020258    2.57
B                    -0.0066632  0.0019706   -3.38
C                    0.0068212  0.0005237   13.02

Correlation of Fixed Effects:
            (Intr) thrsh. length
A          0.081             
B          0.005  0.445       
C         -0.111 -0.556 -0.263


Ce qui me va assez bien, mais lorsque je regarde la table de variance j'obtiens :

Code : Tout sélectionner

    Df  Sum Sq Mean Sq  F value
A  1 1.34008 1.34008 161.5445
B  1 0.00002 0.00002   0.0022
C  1 1.40726 1.40726 169.6436


Comment dois-je interpréter le fait que la F-value soit si basse pour l'une de mes 3 variables? Cela invalide-t-il le résultat?

2/ Par ailleurs, même si j'ai compris que le calcul des p était sujet à caution, le reviewer semble y tenir... J'ai lu quelques méthodes proposées dans des messages précédents, mais je suis aussi tombé sur cette fonction :

http://blog.lib.umn.edu/moor0554/canoem ... s_lrt.html

Cette dernière méthode est-elle standard?

3/ Dans le draft, il est proposé d'utiliser la fonction profile afin de calculer les intervalles de confiance... Peut-être n'ai-je pas la bonne version de cette fonction mais j'obtiens, en suivant l'exemple (Chap1, p16) :

Code : Tout sélectionner

pr1 <- profile(fm1ML)
Error in UseMethod("profile") :
no applicable method for 'profile' applied to an object of class "mer"


Y-a-t-il une autre fonction à charger?

4/ Une remarque plus qu'une question, mais comme j'ai passé pas mal de temps à ruminer là-dessus, je me dis que cela pourra être utile à un autre débutant, par la suite...

Dans le choix du modèle, j'ai voulu utiliser la fonction anova:

Lorsque je fais :

Code : Tout sélectionner

aga2 <- lmer(lrt ~ A + B*C  + (1 | subj) + (1 | word)  ,
        data=u,
            subset= conditions_diverses...)
aga3 <- lmer(lrt ~ A  + (1 | subj) + (1 | word)  ,
        data=u,
            subset= conditions_diverses...)
anova(aga2, aga3)


Cela plante, alors que lorsque je fais :

Code : Tout sélectionner

u14= subset(u, conditions _diverses)
aga2 <- lmer(lrt ~ A + B*C  + (1 | subj) + (1 | word)  ,
        data=u14)
aga3 <- lmer(lrt ~ A  + (1 | subj) + (1 | word)  ,
        data=u14)
anova(aga2, aga3)


Cela marche très bien... Y-a-t-il une explication logique?

5/ Enfin, une dernière question sans doute un peu hors du champ de ce forum mais à tout hasard: quels sont les valeurs à rapporter pour un modèle mixte dans un papier? Estimate et t-value pour chaque variable suffisent-elles, faut-il préciser la f-value? D'autres paramètres du modèle?

En vous remerciant d'avance pour votre aide,

Fabien Vinckier

[Nicolas Péru]

Bonjour,

je ne vais pas être en mesure de répondre sur tes questions directement liées aux GLMM. Mais une chose me dérange quant à ta question sur les p-value. Effectivement, j'ai également vu apparaitre ici et ailleurs que les p-values peuvent être remises en question pour les modèles mixtes (pour les autres aussi mais pour des raisons différentes). Si c'est le cas, alors il existe certainement des références à ce sujet (Pinheiro et Bates ? Zuur et al 2009?).
Ce n'est pas parce que le reviewer y tient qu'il faut le mettre. Il faut simplement pouvoir argumenter ses choix avec des références à l'appuis sinon tu vas potentiellement perpétuer une erreur, ce qui n'est pas le but d'un article.

Après bien sûr tu es libre de tes choix :D

Nicolas

[Renaud Lancelot]

J'utilise un modèle de type

aga1 <- lmer(lrt ~A + B + C + (1 | subj) + (1 | word) , data=u23)

, et j'obtiens ceci :

Code : Tout sélectionner

 Data: u23right
   AIC   BIC logLik deviance REMLdev
 -5279 -5237   2647    -5335   -5293
Random effects:
 Groups   Name        Variance   Std.Dev.
 word     (Intercept) 0.00043719 0.020909
 subj     (Intercept) 0.00480310 0.069304
 Residual             0.00829540 0.091079
Number of obs: 2840, groups: word, 146; subj, 30

Fixed effects:
                    Estimate Std. Error t value
(Intercept)        2.8384507  0.0129920  218.48
A                     0.0052163  0.0020258    2.57
B                    -0.0066632  0.0019706   -3.38
C                    0.0068212  0.0005237   13.02

Correlation of Fixed Effects:
            (Intr) thrsh. length
A          0.081             
B          0.005  0.445       
C         -0.111 -0.556 -0.263


Ce qui me va assez bien, mais lorsque je regarde la table de variance j'obtiens :

Code : Tout sélectionner

    Df  Sum Sq Mean Sq  F value
A  1 1.34008 1.34008 161.5445
B  1 0.00002 0.00002   0.0022
C  1 1.40726 1.40726 169.6436


Comment dois-je interpréter le fait que la F-value soit si basse pour l'une de mes 3 variables? Cela invalide-t-il le résultat?

L'analyse de variance est séquentielle. Avec un plan d'observation déséquilibré, cela lui retire une grande partie de son intérêt sinon tout car les SS dépendent de l'ordre d'introduction des variables. Il vaut mieux utiliser un test du rapport des vraisemblances (avec modèles ajustés avec REML = FALSE):

Code : Tout sélectionner

aga0 <- lmer(lrt ~ A + B + C + (1 | subj) + (1 | word) , REML = FALSE, data = u23)
aga1 <- lmer(lrt ~ A + B + (1 | subj) + (1 | word), REML = FALSE, data = u23)
aga2 <- lmer(lrt ~ A + C + (1 | subj) + (1 | word), REML = FALSE, data = u23)
aga3 <- lmer(lrt ~ B + C + (1 | subj) + (1 | word), REML = FALSE, data = u23)
anova(aga0, aga1)
anova(aga0, aga2)
anova(aga0, aga3)


2/ Par ailleurs, même si j'ai compris que le calcul des p était sujet à caution, le reviewer semble y tenir... J'ai lu quelques méthodes proposées dans des messages précédents, mais je suis aussi tombé sur cette fonction :

http://blog.lib.umn.edu/moor0554/canoem ... s_lrt.html

Cette dernière méthode est-elle standard?

Utiliser le test du rapport des vraisemblances qui est classique (même si pas parfait).

3/ Dans le draft, il est proposé d'utiliser la fonction profile afin de calculer les intervalles de confiance... Peut-être n'ai-je pas la bonne version de cette fonction mais j'obtiens, en suivant l'exemple (Chap1, p16) :

Code : Tout sélectionner

pr1 <- profile(fm1ML)
Error in UseMethod("profile") :
no applicable method for 'profile' applied to an object of class "mer"


Y-a-t-il une autre fonction à charger?

C'est dans le packafe lme4a qui est sur R-Forge. Je viens de vérifier et apparemment D Bates est en train de tout réorganiser dans ses packages, avec un nouveau package ModelMatrix faisant les calculs des matrices des plans d'observation des effets fixes et aléatoires. Les fonctions de lme4a ont été mises dans lme4b (!). Voici un exemple traité

Code : Tout sélectionner

> install.packages(c("lme4b", "MatrixModels"),
+                  repos = "http://R-Forge.R-project.org")
Avis dans install.packages(c("lme4b", "MatrixModels"), repos = "http://R-Forge.R-project.org") :
  l'argument 'lib' manque : 'C:/R/RLIBS' est utilisé
essai de l'URL 'http://R-Forge.R-project.org/bin/windows/contrib/2.11/lme4b_0.999375-60.zip'
Content type 'application/zip' length 1066887 bytes (1.0 Mb)
URL ouverte
downloaded 1.0 Mb

essai de l'URL 'http://R-Forge.R-project.org/bin/windows/contrib/2.11/MatrixModels_0.2-1.zip'
Content type 'application/zip' length 150865 bytes (147 Kb)
URL ouverte
downloaded 147 Kb

le package 'lme4b' a été décompressé et les sommes MD5 ont été vérifiées avec succés
le package 'MatrixModels' a été décompressé et les sommes MD5 ont été vérifiées avec succés

Les packages téléchargés sont dans
        C:\Temp\RtmpdtHyDw\downloaded_packages
> library(lme4b)
Le chargement a nécessité le package : Matrix
Le chargement a nécessité le package : lattice

Attachement du package : 'Matrix'

The following object(s) are masked from 'package:base':

    det

Le chargement a nécessité le package : minqa
Le chargement a nécessité le package : Rcpp

Attachement du package : 'lme4b'

The following object(s) are masked from 'package:stats':

    AIC

Message d'avis :
les méthodes S3 ‘plot.coef.mer’, ‘dotplot.coef.mer’, ‘plot.ranef.mer’, ‘qqmath.ranef.mer’, ‘dotplot.ranef.mer’ sont déclarées dans NAMESPACE mais sont introuvables
> library(MatrixModels)

Attachement du package : 'MatrixModels'

The following object(s) are masked from 'package:Matrix':

    model.Matrix

> data(Orthodont, package = "nlme")
> fm1 <- lmer1(distance ~ age + Sex + (1 | Subject), data = Orthodont)
> prof1 <- profile(fm1)
> confint(prof1)
                 2.5 %     97.5 %
.sig01       1.2745452  2.4042453
.lsig        0.2061015  0.5148922
(Intercept) 16.0848469 19.3285790
age          0.5387505  0.7816198
SexFemale   -3.8096622 -0.8323833
>

4/ Une remarque plus qu'une question, mais comme j'ai passé pas mal de temps à ruminer là-dessus, je me dis que cela pourra être utile à un autre débutant, par la suite...

Dans le choix du modèle, j'ai voulu utiliser la fonction anova:

Lorsque je fais :

Code : Tout sélectionner

aga2 <- lmer(lrt ~ A + B*C  + (1 | subj) + (1 | word)  ,
        data=u,
            subset= conditions_diverses...)
aga3 <- lmer(lrt ~ A  + (1 | subj) + (1 | word)  ,
        data=u,
            subset= conditions_diverses...)
anova(aga2, aga3)


Cela plante, alors que lorsque je fais :

Code : Tout sélectionner

u14= subset(u, conditions _diverses)
aga2 <- lmer(lrt ~ A + B*C  + (1 | subj) + (1 | word)  ,
        data=u14)
aga3 <- lmer(lrt ~ A  + (1 | subj) + (1 | word)  ,
        data=u14)
anova(aga2, aga3)


Cela marche très bien... Y-a-t-il une explication logique?

L'argument subset est probablement peu fiable. Faire explicitement les jeux de données comme dans le second cas.

5/ Enfin, une dernière question sans doute un peu hors du champ de ce forum mais à tout hasard: quels sont les valeurs à rapporter pour un modèle mixte dans un papier? Estimate et t-value pour chaque variable suffisent-elles, faut-il préciser la f-value? D'autres paramètres du modèle?

Les coefficients et leur intervalle de confiance, par exemple. Si c'est pour juger de l'effet global d'une variable qualitative à plusieurs modalités, test du rapport des vraisemblances comme ci-dessus.

[Fabien Vinckier]

Fantastique, je vais essayer tout ça de ce pas.

Merci beaucoup!

Fabien

[Stéphane Laurent]

Bonjour,
J'ajoute une question : comment récupérer la variance résiduelle quand on a ajusté un modèle avec lmer(). Par exemple :

Code : Tout sélectionner

> summary(fit)
Linear mixed model fit by REML
Formula: values ~ (1 | Lot) + (1 | Lot.Pompe)
   Data: DataR
    AIC    BIC logLik deviance REMLdev
 -19030 -19006   9519   -19050  -19038
Random effects:
 Groups    Name        Variance   Std.Dev.
 Lot.Pompe (Intercept) 2.1877e-05 0.0046773
 Lot       (Intercept) 2.4117e-06 0.0015530
 Residual              1.2552e-04 0.0112034
Number of obs: 3111, groups: Lot.Pompe, 24; Lot, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.637171   0.001339     476

Avec VarCorr(fit) je récupère les variances correspondant à Lot et Lot.Pompe, mais pas la Residual.

[Renaud Lancelot]

C'est bien caché dans la sortie de VarCorr:

Code : Tout sélectionner

> library(lme4)
Le chargement a nécessité le package : Matrix
Le chargement a nécessité le package : lattice

Attachement du package : 'Matrix'

The following object(s) are masked from 'package:base':

    det


Attachement du package : 'lme4'

The following object(s) are masked from 'package:stats':

    AIC, BIC

> fm1 <- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)
> summary(fm1)
Linear mixed model fit by REML
Formula: Reaction ~ Days + (Days | Subject)
   Data: sleepstudy
  AIC  BIC logLik deviance REMLdev
 1756 1775 -871.8     1752    1744
Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Subject  (Intercept) 612.092  24.7405       
          Days         35.072   5.9221  0.066
 Residual             654.941  25.5918       
Number of obs: 180, groups: Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  251.405      6.825   36.84
Days          10.467      1.546    6.77

Correlation of Fixed Effects:
     (Intr)
Days -0.138


Quand on regarde la sortie de VarCorr:

Code : Tout sélectionner

> (vc <- VarCorr(fm1))
$Subject
            (Intercept)      Days
(Intercept)  612.091776  9.603985
Days           9.603985 35.071536
attr(,"stddev")
(Intercept)        Days
  24.740489    5.922123
attr(,"correlation")
            (Intercept)       Days
(Intercept)  1.00000000 0.06554896
Days         0.06554896 1.00000000

attr(,"sc")
[1] 25.59182

Donc:

Code : Tout sélectionner

> attr(vc, "sc")
[1] 25.59182

NB: c'est l'écart-type, et non pas la variance.

[Stéphane Laurent]

Gracias!

[Stéphane Laurent]

Ca me semblerait idéal que la somme des composantes de la variance soit égale à la variance des données. Sauriez-vous pourquoi ce n'est pas le cas ?

[Stéphane Laurent]

Une autre question "diverse" : quand je stocke une sortie de lmer() dans, disons, "fit", la commande model.matrix(fit) semble retourner uniquement la matrice des effets fixes, est-il possible d'obtenir celle des effets aléatoires ?

[Renaud Lancelot]

Il y a une toute nouvelle fonction d'extraction dans la toute nouvelle version du package Matrix (version 1.0): getME(). Elle extrait et éventuellement transforme de manière adéquate les objets stockés dans les slots de l'objet de classe "mer" retourné par lmer and friends. Voir ?mer pour avoir la liste de ces slots. Pour la transposée de la matrice Z du plan d'observation des effets aléatoires:

Code : Tout sélectionner

> fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
> Zt <- getME(fm1, "Zt")
> dim(Zt)
[1]  36 180


Apparemment, on peut directement obtenir Z mais ce n'est pas documenté:

Code : Tout sélectionner

> Z <- getME(fm1, "Z")
> dim(Z)
[1] 180  36
> all.equal(t(Z), Zt)
[1] TRUE


Les auteurs recommandent d'utiliser getME plutôt que d'aller chercher directement les données dans les slots, car (i) les objets stockés peuvent l'être sous une forme transformée et (ii) la fonction sera mise à jour pour la prochaine version majeure de lme4 dans laquelle la liste des slots sera différente.