Différences entre PROC Mixed et lme / lmer en R - degrés de liberté

Remarque: cette question est une rediffusion, car ma question précédente a dû être supprimée pour des raisons juridiques.

En comparant PROC MIXED de SAS avec la fonction lmedu nlmepackage dans R, je suis tombé sur des différences assez confuses. Plus précisément, les degrés de liberté dans les différents tests diffèrent entre PROC MIXEDet lme, et je me suis demandé pourquoi.

Commencez à partir de l'ensemble de données suivant (code R ci-dessous):

• ind: facteur indiquant l'individu où la mesure est prise
• fac: organe où la mesure est prise
• trt: facteur indiquant le traitement
• y: une variable de réponse continue

L'idée est de construire les modèles simples suivants:

y ~ trt + (ind): indcomme facteur aléatoire y ~ trt + (fac(ind)): facimbriqué indcomme facteur aléatoire

Notez que le dernier modèle devrait provoquer des singularités, car il n'y a qu'une seule valeur de ypour chaque combinaison de indet fac.

Premier modèle

En SAS, je construis le modèle suivant:

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM ind /s;
run;

Selon les tutoriels, le même modèle dans R utilisant nlmedevrait être:

> require(nlme)
> options(contrasts=c(factor="contr.SAS",ordered="contr.poly"))
> m2<-lme(y~trt,random=~1|ind,data=Data)

Les deux modèles donnent les mêmes estimations pour les coefficients et leur SE, mais lorsqu'ils effectuent un test F pour l'effet de trt, ils utilisent une quantité différente de degrés de liberté:

SAS : 
Type 3 Tests of Fixed Effects 
Effect Num DF Den DF     F  Value Pr > F 
trt         1      8  0.89        0.3724 

R : 
> anova(m2)
            numDF denDF  F-value p-value
(Intercept)     1     8 70.96836  <.0001
trt             1     6  0.89272  0.3812

Question 1: Quelle est la différence entre les deux tests? Les deux sont équipés de REML et utilisent les mêmes contrastes.

REMARQUE: j'ai essayé différentes valeurs pour l'option DDFM = (y compris BETWITHIN, qui devrait théoriquement donner les mêmes résultats que lme)

Deuxième modèle

En SAS:

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM fac(ind) /s;
run;

Le modèle équivalent dans R devrait être:

> m4<-lme(y~trt,random=~1|ind/fac,data=Data)

Dans ce cas, il existe des différences très étranges:

• R rentre sans se plaindre, alors que SAS note que la toile de jute finale n'est pas définie positive (ce qui ne me surprend pas un peu, voir ci-dessus)
• Le SE sur les coefficients diffère (est plus petit en SAS)
• Encore une fois, le test F a utilisé une quantité différente de DF (en fait, en SAS, cette quantité = 0)

Sortie SAS:

Effect     trt Estimate Std Error  DF t Value Pr > |t| 
Intercept        0.8863    0.1192  14    7.43 <.0001 
trt       Cont  -0.1788    0.1686   0   -1.06 . 

Sortie R:

> summary(m4)
...
Fixed effects: y ~ trt 
               Value Std.Error DF   t-value p-value
(Intercept)  0.88625 0.1337743  8  6.624963  0.0002
trtCont     -0.17875 0.1891855  6 -0.944840  0.3812
...

(Notez que dans ce cas, les tests F et T sont équivalents et utilisent le même DF.)

Fait intéressant, lors de l'utilisation lme4dans R, le modèle ne convient même pas:

> require(lme4)
> m4r <- lmer(y~trt+(1|ind/fac),data=Data)
Error in function (fr, FL, start, REML, verbose)  : 
  Number of levels of a grouping factor for the random effects
must be less than the number of observations

Question 2 : Quelle est la différence entre ces modèles avec des facteurs imbriqués? Sont-ils spécifiés correctement et si oui, comment les résultats sont-ils si différents?

Données simulées dans R:

Data <- structure(list(y = c(1.05, 0.86, 1.02, 1.14, 0.68, 1.05, 0.22, 
1.07, 0.46, 0.65, 0.41, 0.82, 0.6, 0.49, 0.68, 1.55), ind = structure(c(1L, 
2L, 3L, 1L, 3L, 4L, 4L, 2L, 5L, 6L, 7L, 8L, 6L, 5L, 7L, 8L), .Label = c("1", 
"2", "3", "4", "5", "6", "7", "8"), class = "factor"), fac = structure(c(1L, 
1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L), .Label = c("l", 
"r"), class = "factor"), trt = structure(c(2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Cont", 
"Treat"), class = "factor")), .Names = c("y", "ind", "fac", "trt"
), row.names = c(NA, -16L), class = "data.frame")

Données simulées:

   y ind fac   trt
1.05   1   l Treat
0.86   2   l Treat
1.02   3   l Treat
1.14   1   r Treat
0.68   3   r Treat
1.05   4   l Treat
0.22   4   r Treat
1.07   2   r Treat
0.46   5   r  Cont
0.65   6   l  Cont
0.41   7   l  Cont
0.82   8   l  Cont
0.60   6   r  Cont
0.49   5   l  Cont
0.68   7   r  Cont
1.55   8   r  Cont

Pour la première question, la méthode par défaut dans SAS pour trouver le df n'est pas très intelligente; il recherche des termes dans l'effet aléatoire qui incluent syntaxiquement l'effet fixe, et l'utilise. Dans ce cas, puisque trtne se trouve pas dans ind, il ne fait pas la bonne chose. Je n'ai jamais essayé BETWITHINet je ne connais pas les détails, mais soit l'option Satterthwaite ( satterth) ou l'utilisation ind*trtcomme effet aléatoire donne des résultats corrects.

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s ddfm=satterth;
    RANDOM ind /s;
run;

PROC MIXED data=Data;
    CLASS ind fac trt;
    MODEL y = trt /s;
    RANDOM ind*trt /s;
run;

Quant à la deuxième question, votre code SAS ne correspond pas tout à fait à votre code R; il n'a qu'un terme pour fac*ind, tandis que le code R a un terme pour indet fac*ind. (Voir la sortie Composants de variance pour voir ceci.) L'ajout de cela donne le même SE pour trttous les modèles à la fois au T1 et au T2 (0,1889).

Comme vous le constatez, il s'agit d'un modèle étrange à adapter car le fac*indterme a une observation pour chaque niveau, il est donc équivalent au terme d'erreur. Cela se reflète dans la sortie SAS, où le fac*indterme a une variance nulle. C'est aussi ce que le message d'erreur de lme4 vous dit; la raison de l'erreur est que vous avez probablement mal spécifié quelque chose car vous incluez le terme d'erreur dans le modèle de deux manières différentes. Fait intéressant, il existe une légère différence dans le modèle nlme; il s'agit en quelque sorte de trouver un terme de variance pour le fac*indterme en plus du terme d'erreur, mais vous remarquerez que la somme de ces deux variances est égale au terme d'erreur de SAS et de nlme sans le fac*indterme. Cependant, la SE pour trtreste la même (0,1892) que celle trtimbriquée dansind, de sorte que ces termes de variance inférieurs ne l'affectent pas.

Enfin, une note générale sur les degrés de liberté dans ces modèles: ils sont calculés une fois que le modèle est adapté, et donc les différences dans les degrés de liberté entre les différents programmes ou options d'un programme ne signifient pas nécessairement que le modèle est ajusté différemment. Pour cela, il faut regarder les estimations des paramètres, à la fois paramètres à effet fixe et paramètres de covariance.

De plus, l'utilisation des approximations t et F avec un nombre donné de degrés de liberté est assez controversée. Non seulement il existe plusieurs façons d'approximer le df, certains pensent que la pratique de le faire n'est pas une bonne idée de toute façon. Quelques conseils:

1. Si tout est équilibré, comparez les résultats avec la méthode traditionnelle des moindres carrés, comme ils devraient être d'accord. Si elle est proche de l'équilibre, calculez-les vous-même (en supposant un équilibre) afin de vous assurer que ceux que vous utilisez sont dans le bon stade.

2. Si vous avez un échantillon de grande taille, les degrés de liberté importent peu car les distributions se rapprochent d'une normale et d'un chi carré.

3. Découvrez les méthodes de Doug Bates pour l'inférence. Son ancienne méthode est basée sur la simulation MCMC; sa nouvelle méthode est basée sur le profilage de la probabilité.