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ą lme
z nlme
pakietu 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 MIXED
i lme
zastanawiał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)
: ind
jako czynnik losowy
y ~ trt + (fac(ind))
: fac
zagnieżdżony ind
jako czynnik losowy
Zauważ, że ostatni model powinien powodować osobliwości, ponieważ istnieje tylko 1 wartość y
dla każdej kombinacji ind
i 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 nlme
powinien 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 lme4
w 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
źródło
Odpowiedzi:
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ż
trt
nie można go znaleźć wind
, nie robi to dobrze. Nigdy nie próbowałemBETWITHIN
i nie znam szczegółów, ale albo opcja Satterthwaite (satterth
), albo użycieind*trt
jako efektu losowego daje prawidłowe wyniki.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, jakind
ifac*ind
. (Zobacz dane wyjściowe Variance Components, aby to zobaczyć.) Dodanie tego daje taki sam SE dlatrt
wszystkich modeli w Q1 i Q2 (0.1892).Jak zauważasz, jest to dziwny model do dopasowania, ponieważ
fac*ind
termin ma jedną obserwację dla każdego poziomu, więc jest równoważny z błędem. Znajduje to odzwierciedlenie w danych wyjściowych SAS, gdziefac*ind
termin 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 dlafac*ind
terminu 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 tegofac*ind
terminu. Jednak SE dlatrt
pozostaje taka sama (0,1892), jaktrt
jest 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:
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.
Jeśli masz dużą próbkę, stopnie swobody nie mają większego znaczenia, ponieważ rozkłady zbliżają się do normalnej i chi-kwadrat.
Sprawdź metody wnioskowania Douga Batesa. Jego starsza metoda oparta jest na symulacji MCMC; jego nowsza metoda opiera się na profilowaniu prawdopodobieństwa.
źródło