Różnice między PROC Mixed i lme / lmer w R - stopnie swobody

Uwaga: to pytanie jest repost, ponieważ moje poprzednie pytanie musiało zostać usunięte ze względów prawnych.

Porównując PROC MIXED z SAS z funkcją lmez nlmepakietu w R, natknąłem się na pewne dość mylące różnice. Mówiąc dokładniej, stopnie swobody w różnych testach różnią się między PROC MIXEDi lmezastanawiałem się, dlaczego.

Zacznij od następującego zestawu danych (kod R podany poniżej):

• ind: współczynnik wskazujący osobę, w której wykonywany jest pomiar
• fac: organ, na którym dokonywany jest pomiar
• trt: czynnik wskazujący na leczenie
• y: jakaś zmienna ciągłej odpowiedzi

Chodzi o zbudowanie następujących prostych modeli:

y ~ trt + (ind): indjako czynnik losowy y ~ trt + (fac(ind)): faczagnieżdżony indjako czynnik losowy

Zauważ, że ostatni model powinien powodować osobliwości, ponieważ istnieje tylko 1 wartość ydla każdej kombinacji indi fac.

Pierwszy model

W SAS buduję następujący model:

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

Zgodnie z samouczkami ten sam model w użyciu R nlmepowinien być:

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

Oba modele dają takie same oszacowania dla współczynników i ich SE, ale przeprowadzając test F dla efektu trt, używają innej ilości stopni swobody:

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

Pytanie 1: Jaka jest różnica między obydwoma testami? Oba są dopasowane za pomocą REML i używają tych samych kontrastów.

UWAGA: Próbowałem różnych wartości dla opcji DDFM = (w tym BETWITHIN, które teoretycznie powinny dać takie same wyniki jak lme)

Drugi model

W SAS:

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

Odpowiednikiem modelu w R powinno być:

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

W tym przypadku istnieją bardzo dziwne różnice:

• R pasuje bez narzekania, podczas gdy SAS zauważa, że ​​ostateczny hessian nie jest zdecydowanie pozytywny (co mnie nie dziwi, patrz wyżej)
• SE na współczynniki różnią się (jest mniejszy w SAS)
• Ponownie, test F wykorzystał inną ilość DF (w rzeczywistości w SAS ta ilość = 0)

Wyjście 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 . 

R Wyjście:

> 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
...

(Należy pamiętać, że w tym przypadku test F i T są równoważne i używają tego samego DF.)

Co ciekawe, przy użyciu lme4w R model nawet nie pasuje:

> 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

Pytanie 2 : Jaka jest różnica między tymi modelami z zagnieżdżonymi czynnikami? Czy są określone poprawnie, a jeśli tak, to dlaczego wyniki są tak różne?

Dane symulowane w 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")

Dane symulowane:

   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

W przypadku pierwszego pytania domyślna metoda znalezienia pliku df w SAS nie jest zbyt inteligentna; szuka warunków w efekcie losowym, które syntaktycznie obejmują ustalony efekt, i wykorzystuje to. W tym przypadku, ponieważ trtnie można go znaleźć w ind, nie robi to dobrze. Nigdy nie próbowałem BETWITHINi nie znam szczegółów, ale albo opcja Satterthwaite ( satterth), albo użycie ind*trtjako efektu losowego daje prawidłowe wyniki.

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;

Co do drugiego pytania, kod SAS nie do końca pasuje do kodu R; ma tylko termin fac*ind, podczas gdy kod R ma termin zarówno dla, jak indi fac*ind. (Zobacz dane wyjściowe Variance Components, aby to zobaczyć.) Dodanie tego daje taki sam SE dla trtwszystkich modeli w Q1 i Q2 (0.1892).

Jak zauważasz, jest to dziwny model do dopasowania, ponieważ fac*indtermin ma jedną obserwację dla każdego poziomu, więc jest równoważny z błędem. Znajduje to odzwierciedlenie w danych wyjściowych SAS, gdzie fac*indtermin ma zerową wariancję. To właśnie mówi komunikat o błędzie z lme4; przyczyną błędu jest to, że najprawdopodobniej coś źle sprecyzowałeś, ponieważ termin błędu dołączasz do modelu na dwa różne sposoby. Co ciekawe, istnieje jedna niewielka różnica w modelu NLM; w jakiś sposób znajduje warunek wariancji dla fac*indterminu oprócz warunku błędu, ale zauważysz, że suma tych dwóch wariancji jest równa warunkowi błędu zarówno SAS, jak i nlme bez tego fac*indterminu. Jednak SE dla trtpozostaje taka sama (0,1892), jak trtjest zagnieżdżonaind, więc te warunki o niższej wariancji nie mają na to wpływu.

Wreszcie, ogólna uwaga na temat stopni swobody w tych modelach: są one obliczane po dopasowaniu modelu, a zatem różnice w stopniach swobody między różnymi programami lub opcjami programu niekoniecznie oznaczają, że model jest inaczej dopasowany. W tym celu należy spojrzeć na oszacowania parametrów, zarówno parametry efektu stałego, jak i parametry kowariancji.

Również stosowanie aproksymacji t i F z określoną liczbą stopni swobody jest dość kontrowersyjne. Istnieje nie tylko kilka sposobów przybliżenia df, ale niektórzy uważają, że praktyka robienia tego nie jest dobrym pomysłem. Kilka słów porady:

1. Jeśli wszystko jest zrównoważone, porównaj wyniki tradycyjną metodą najmniejszych kwadratów, tak jak powinny się zgodzić. Jeśli jest bliski wyważeniu, oblicz je sam (zakładając równowagę), aby mieć pewność, że te, których używasz, znajdują się we właściwym polu.

2. Jeśli masz dużą próbkę, stopnie swobody nie mają większego znaczenia, ponieważ rozkłady zbliżają się do normalnej i chi-kwadrat.

3. Sprawdź metody wnioskowania Douga Batesa. Jego starsza metoda oparta jest na symulacji MCMC; jego nowsza metoda opiera się na profilowaniu prawdopodobieństwa.