Quelle méthode de comparaison multiple utiliser pour un modèle lmer: lsmeans ou glht?

15

J'analyse un ensemble de données à l'aide d'un modèle à effets mixtes avec un effet fixe (condition) et deux effets aléatoires (participant en raison de la conception et de la paire du sujet). Le modèle a été généré avec le lme4package: exp.model<-lmer(outcome~condition+(1|participant)+(1|pair),data=exp).

Ensuite, j'ai effectué un test de rapport de vraisemblance de ce modèle par rapport au modèle sans l'effet fixe (condition) et j'ai une différence significative. Il y a 3 conditions dans mon ensemble de données, donc je veux faire une comparaison multiple mais je ne sais pas quelle méthode utiliser . J'ai trouvé un certain nombre de questions similaires sur CrossValidated et d'autres forums, mais je suis encore assez confus.

D'après ce que j'ai vu, les gens ont suggéré d'utiliser

1. Le lsmeanspackage - lsmeans(exp.model,pairwise~condition)qui me donne la sortie suivante:

condition     lsmean         SE    df  lower.CL  upper.CL
 Condition1 0.6538060 0.03272705 47.98 0.5880030 0.7196089
 Condition2 0.7027413 0.03272705 47.98 0.6369384 0.7685443
 Condition3 0.7580522 0.03272705 47.98 0.6922493 0.8238552

Confidence level used: 0.95 

$contrasts
 contrast                   estimate         SE    df t.ratio p.value
 Condition1 - Condition2 -0.04893538 0.03813262 62.07  -1.283  0.4099
 Condition1 - Condition3 -0.10424628 0.03813262 62.07  -2.734  0.0219
 Condition2 - Condition3 -0.05531090 0.03813262 62.07  -1.450  0.3217

P value adjustment: tukey method for comparing a family of 3 estimates 

2. Le multcomppackage de deux manières différentes - l'utilisation mcp glht(exp.model,mcp(condition="Tukey"))résultant en

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = outcome ~ condition + (1 | participant) + (1 | pair), 
    data = exp, REML = FALSE)

Linear Hypotheses:
                             Estimate Std. Error z value Pr(>|z|)  
Condition2 - Condition1 == 0  0.04894    0.03749   1.305    0.392  
Condition3 - Condition1 == 0  0.10425    0.03749   2.781    0.015 *
Condition3 - Condition2 == 0  0.05531    0.03749   1.475    0.303  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1
(Adjusted p values reported -- single-step method)

et l'utilisation lsm glht(exp.model,lsm(pairwise~condition))résultant en

Note: df set to 62

     Simultaneous Tests for General Linear Hypotheses

Fit: lmer(formula = outcome ~ condition + (1 | participant) + (1 | pair), 
    data = exp, REML = FALSE)

Linear Hypotheses:
                             Estimate Std. Error t value Pr(>|t|)  
Condition1 - Condition2 == 0 -0.04894    0.03749  -1.305   0.3977  
Condition1 - Condition3 == 0 -0.10425    0.03749  -2.781   0.0195 *
Condition2 - Condition3 == 0 -0.05531    0.03749  -1.475   0.3098  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1
(Adjusted p values reported -- single-step method)

Comme vous pouvez le voir, les méthodes donnent des résultats différents. C'est la première fois que je travaille avec R et les statistiques, donc quelque chose pourrait mal tourner mais je ne sais pas où. Mes questions sont:

Quelles sont les différences entre les méthodes présentées? J'ai lu dans une réponse à une question connexe qu'il s'agit des degrés de liberté ( lsmeansvs. glht). Existe-t-il des règles ou des recommandations quant à l'utilisation de la méthode 1, c'est-à-dire la bonne pour ce type d'ensemble / modèle de données, etc.? Quel résultat dois-je signaler? Sans savoir mieux, j'irais probablement rapporter la valeur de p la plus élevée que j'ai pu jouer en toute sécurité, mais ce serait bien d'avoir une meilleure raison. Merci

schvaba986
la source

Réponses:

17

Pas une réponse complète ...

La différence entre glht(myfit, mcp(myfactor="Tukey"))les deux autres méthodes est que cette méthode utilise une statistique "z" (distribution normale), tandis que les autres utilisent une statistique "t" (distribution de Student). La statistique "z" est identique à une statistique "t" avec un degré de liberté infini. Cette méthode est asymptotique et fournit des valeurs de p plus petites et des intervalles de confiance plus courts que les autres. Les valeurs de p peuvent être trop petites et les intervalles de confiance peuvent être trop courts si l'ensemble de données est petit.

Lorsque j'exécute, lsmeans(myfit, pairwise~myfactor)le message suivant apparaît:

Loading required namespace: pbkrtest

Cela signifie que lsmeans(pour un lmermodèle) utilise le pbkrtestpackage qui implémente la méthode de Kenward & Rogers pour les degrés de liberté de la statistique "t". Cette méthode vise à fournir de meilleures valeurs de p et intervalles de confiance que l'asymptotique (il n'y a pas de différence lorsque le degré de liberté est grand).

Maintenant, à propos de la différence entre lsmeans(myfit, pairwise~myfactor)$contrastset glht(myfit, lsm(pairwise~factor), je viens de faire quelques tests et mes observations sont les suivantes:

  • lsmest une interface entre le lsmeanspackage et le multcomppackage (voir ?lsm)

  • pour un design équilibré il n'y a pas de différence entre les résultats

  • pour une conception déséquilibrée, j'ai observé de petites différences entre les résultats (les erreurs standard et le rapport t)

Malheureusement, je ne sais pas quelle est la cause de ces différences. Il ressemble à des lsmappels lsmeansuniquement pour obtenir la matrice d'hypothèses linéaires et les degrés de liberté, mais lsmeansutilise une méthode différente pour calculer les erreurs standard.

Stéphane Laurent
la source
Merci pour la réponse détaillée! J'ai complètement raté la différence de statistique de test ... Vous mentionnez que les valeurs peuvent être trop petites et les CI trop étroits pour la méthode asymptotique. Mon ensemble de données comprend environ 30 participants, donc je suppose que je m'en tiendrai à la statistique t. Lorsque vous dites que la méthode de Kenward & Rogers conduit à de meilleures valeurs de p, voulez-vous dire plus précises ou plus petites? Donc, les différences sont dues à des différences dans les méthodes de calcul df et SE et non à l'utilisation incorrecte de l'une d'entre elles avec mon modèle, si je vous ai bien compris. Existe-t-il un moyen de choisir la "meilleure" méthode ici?
schvaba986
11
(Je suis le développeur du package lsmeans ) lsmeansutilise le package pbkrtest, qui prévoit (1) les calculs df de Kenward-Rogers et (2) une matrice de covariance ajustée avec un biais réduit dans les estimations. Si vous définissez d'abord lsm.options(disable.pbkrtest=TRUE), alors l' lsmeansappel avec adjust="mvt"donnera les mêmes résultats que glht, à l'exception de légères différences dues à l'algorithme randomisé utilisé par les deux packages pour la distribution t multivariée.
Russ Lenth
3
Cependant, je suggère l'ajustement "mvt" sans désactiver pbkrtest, en raison de l'ajustement du biais et du fait que sans df, les valeurs asymptotiques (z) supposent essentiellement un df infini, produisant ainsi des valeurs de P irréalistes.
Russ Lenth
3
Soit dit en passant, la summaryméthode pour glhtpermet diverses méthodes de test de réduction en plus de l'ajustement de multiplicité en une étape (CI simultanés) par défaut. Sur un point complètement différent, si vous avez plus d'un facteur, vous lsmpouvez créer les types de comparaisons habituels assez facilement, tout en mcpne pouvant pas le faire du tout.
Russ Lenth