Dlaczego otrzymuję bardzo różne wyniki dla poly (raw = T) vs. poly ()?

10

Chcę modelować dwie różne zmienne czasowe, z których niektóre są silnie współliniowe w moich danych (wiek + kohorta = okres). Robiąc to, miałem problemy z lmerinterakcjami i poly(), ale prawdopodobnie nie jest to ograniczone lmer, otrzymałem takie same wyniki z nlmeIIRC.

Oczywiście brakuje mi zrozumienia tego, co robi funkcja poli (). Rozumiem, co poly(x,d,raw=T)robi i pomyślałem, że bez raw=Tniego powstają wielomiany ortogonalne (nie mogę powiedzieć, że naprawdę rozumiem, co to oznacza), co ułatwia dopasowanie, ale nie pozwala na bezpośrednią interpretację współczynników.
I czytać , bo używam funkcji przewidywania, prognozy powinny być takie same.

Ale tak nie jest, nawet jeśli modele zbiegają się normalnie. Używam zmiennych wyśrodkowanych i najpierw pomyślałem, że może wielomian ortogonalny prowadzi do wyższej korelacji efektu stałego z kolinearnym terminem interakcji, ale wydaje się, że jest porównywalny. Wkleiłem tutaj dwa podsumowania modeli .

Te wykresy mają nadzieję ilustrują zakres różnicy. Użyłem funkcji przewidywania, która jest dostępna tylko w dev. wersja lme4 (słyszałem o tym tutaj ), ale ustalone efekty są takie same w wersji CRAN (i same również wydają się wyłączone, np. ~ 5 dla interakcji, gdy mój DV ma zakres 0-4).

Lmer był

cohort2_age =lmer(churchattendance ~ 
poly(cohort_c,2,raw=T) * age_c + 
ctd_c + dropoutalive + obs_c + (1+ age_c |PERSNR), data=long.kg)

Prognozę poprawiono tylko w przypadku efektów, na fałszywych danych (wszystkie inne predyktory = 0), gdzie zaznaczyłem zakres obecny w oryginalnych danych jako ekstrapolacja = F.

predict(cohort2_age,REform=NA,newdata=cohort.moderates.age)

Mogę podać więcej kontekstu, jeśli zajdzie taka potrzeba (nie udało mi się łatwo stworzyć powtarzalnego przykładu, ale oczywiście mogę się bardziej postarać), ale myślę, że jest to bardziej podstawowy zarzut: wyjaśnij mi tę poly()funkcję, proszę bardzo.

Surowe wielomiany

Surowe wielomiany

Wielomiany ortogonalne (obcięte, nielipcowane w Imgur )

Wielomiany ortogonalne

Ruben
źródło

Odpowiedzi:

10

Myślę, że jest to błąd w funkcji przewidywania (i stąd moja wina), którego w rzeczywistości nlme nie udostępnia. ( Edycja : należy naprawić w najnowszej wersji R-forge lme4.) Przykład poniżej ...

Myślę, że twoje zrozumienie ortogonalnych wielomianów jest prawdopodobnie w porządku. Trudną rzeczą, którą musisz o nich wiedzieć, jeśli próbujesz napisać metodę prognozowania dla klasy modeli, jest to, że podstawa dla wielomianów ortogonalnych jest definiowana na podstawie danego zestawu danych, więc jeśli naiwnie (tak jak ja!) ) użyj, model.matrixaby spróbować wygenerować macierz projektową dla nowego zestawu danych, otrzymasz nową podstawę - co nie ma już sensu w przypadku starych parametrów. Dopóki nie naprawię tego, może być konieczne umieszczenie pułapki, która mówi ludziom, że predictnie działa z ortogonalnymi podstawami wielomianowymi (lub podstawami splajnu, które mają tę samą właściwość).

d <- expand.grid(x=seq(0,1,length=50),f=LETTERS[1:10])
set.seed(1001)
u.int <- rnorm(10,sd=0.5)
u.slope <- rnorm(10,sd=0.2)
u.quad <- rnorm(10,sd=0.1)
d <- transform(d,
               ypred = (1+u.int[f])+
               (2+u.slope[f])*x-
               (1+u.quad[f])*x^2)
d$y <- rnorm(nrow(d),mean=d$ypred,sd=0.2)
ggplot(d,aes(x=x,y=y,colour=f))+geom_line()+
    geom_line(aes(y=ypred),linetype=2)

library(lme4)
fm1 <- lmer(y~poly(x,2,raw=TRUE)+(1|f)+(0+x|f)+(0+I(x^2)|f),
            data=d)


fm2 <- lmer(y~poly(x,2)+(1|f)+(0+x|f)+(0+I(x^2)|f),
            data=d)
newdat <- data.frame(x=unique(d$x))
plot(predict(fm1,newdata=newdat,REform=NA))
lines(predict(fm2,newdata=newdat,REform=NA),col=2)
detach("package:lme4")

library(nlme)
fm3 <- lme(y~poly(x,2,raw=TRUE),
           random=list(~1|f,~0+x|f,~0+I(x^2)|f),
            data=d)
VarCorr(fm3)

fm4 <- lme(y~poly(x,2),
           random=list(~1|f,~0+x|f,~0+I(x^2)|f),
            data=d)

newdat <- data.frame(x=unique(d$x))
lines(predict(fm3,newdata=newdat,level=0),col=4)
lines(predict(fm4,newdata=newdat,level=0),col=5)
Ben Bolker
źródło
Dziękuję, to uspokajające. Powtórzmy: przeczytałem, że nie można przyjmować ustalonych efektów wielomianu ortogonalnego według wartości nominalnej, ale czasami wydają się one niesamowicie duże. Na przykład, jeśli uruchomię interakcję dwóch wielomianów sześciennych, otrzymam ustalone efekty dla wielomianów i ich interakcji w zakresie od -22 do -127400. Wydaje mi się to dalekie, szczególnie biorąc pod uwagę, że wszystkie ustalone efekty są negatywne. Czy zmieniona funkcja przewidywania miałaby sens z tymi stałymi efektami, czy też modele fałszywie się zbiegały, czy w końcu ich coś jest nie tak?
Ruben
Ponownie podejrzewam (ale oczywiście nie wiem na pewno), że wszystko jest w porządku. Orth. wielomiany są dobre do testowania stabilności liczbowej i testowania hipotez, ale (jak się domyślasz) rzeczywiste wartości parametrów mogą być trudniejsze do interpretacji. Obecna wersja lme4-devel (właśnie opublikowałem wersję, która powinna przejść testy, odbudowanie na r-forge może zająć około 24 godzin, chyba że możesz sam zbudować z SVN) powinno dać ci pasujące prognozy między wielomianami surowymi / orto. Alternatywą jest wyśrodkowanie i skalowanie ciągłych predyktorów à la Schielzeth 2010 Methods in Ecology & Evolution ...
Ben Bolker
Tak, te dwa wielomiany są teraz całkowicie zgodne. Wielkie dzięki! Przeskalowałem i wyśrodkowałem moje predyktory, ale niektóre modele nie pasowały do ​​surowych wielomianów.
Ruben