Jestem użytkownikiem bardziej zaznajomionym z R. Próbowałem oszacować losowe zbocza (współczynniki selekcji) dla około 35 osobników w ciągu 5 lat dla czterech zmiennych siedlisk. Zmienna odpowiedzi określa, czy lokalizacja była siedliskiem „używanym” (1), czy „dostępnym” (0) („używaj” poniżej).
Korzystam z 64-bitowego komputera z systemem Windows.
W wersji R 3.1.0 używam poniższych danych i wyrażeń. PS, TH, RS i HW są stałymi efektami (znormalizowana, mierzona odległość od typów siedlisk). lme4 V 1.1-7.
str(dat)
'data.frame': 359756 obs. of 7 variables:
$ use : num 1 1 1 1 1 1 1 1 1 1 ...
$ Year : Factor w/ 5 levels "1","2","3","4",..: 4 4 4 4 4 4 4 4 3 4 ...
$ ID : num 306 306 306 306 306 306 306 306 162 306 ...
$ PS: num -0.32 -0.317 -0.317 -0.318 -0.317 ...
$ TH: num -0.211 -0.211 -0.211 -0.213 -0.22 ...
$ RS: num -0.337 -0.337 -0.337 -0.337 -0.337 ...
$ HW: num -0.0258 -0.19 -0.19 -0.19 -0.4561 ...
glmer(use ~ PS + TH + RS + HW +
(1 + PS + TH + RS + HW |ID/Year),
family = binomial, data = dat, control=glmerControl(optimizer="bobyqa"))
glmer daje mi oszacowania parametrów dla ustalonych efektów, które mają dla mnie sens, a losowe zbocza (które interpretuję jako współczynniki selekcji dla każdego typu siedliska) również mają sens, gdy badam jakościowo dane. Prawdopodobieństwo dziennika dla modelu wynosi -3050,8.
Jednak większość badań w dziedzinie ekologii zwierząt nie wykorzystuje R, ponieważ w przypadku danych o lokalizacji zwierząt autokorelacja przestrzenna może powodować, że standardowe błędy są podatne na błąd typu I. Podczas gdy R używa standardowych błędów opartych na modelu, preferowane są błędy standardowe empiryczne (również białe Hubera lub sandwich).
Podczas gdy R obecnie nie oferuje tej opcji (według mojej wiedzy - PROSZĘ, popraw mnie, jeśli się mylę), SAS ma - chociaż nie mam dostępu do SAS, kolega zgodził się pozwolić mi pożyczyć swój komputer, aby ustalić, czy standardowe błędy zmieniają się znacząco, gdy stosowana jest metoda empiryczna.
Po pierwsze, chcieliśmy upewnić się, że przy użyciu standardowych błędów opartych na modelu, SAS wygeneruje oszacowania podobne do R - aby mieć pewność, że model jest określony w ten sam sposób w obu programach. Nie obchodzi mnie, czy są dokładnie takie same - po prostu podobne. Próbowałem (SAS V 9.2):
proc glimmix data=dat method=laplace;
class year id;
model use = PS TH RS HW / dist=bin solution ddfm=betwithin;
random intercept PS TH RS HW / subject = year(id) solution type=UN;
run;title;
Próbowałem także różnych innych form, takich jak dodawanie linii
random intercept / subject = year(id) solution type=UN;
random intercept PS TH RS HW / subject = id solution type=UN;
Próbowałem bez podawania
solution type = UN,lub komentowanie
ddfm=betwithin;Bez względu na to, jak określimy model (i próbowaliśmy na wiele sposobów), nie mogę uzyskać losowych nachyleń w SAS, aby zdalnie przypominały te wyjściowe z R - nawet jeśli ustalone efekty są wystarczająco podobne. A kiedy mam na myśli inny, to znaczy, że nawet znaki nie są takie same. Prawdopodobieństwo dziennika -2 w SAS wyniosło 71344,94.
Nie mogę załadować mojego pełnego zestawu danych; więc stworzyłem zestaw danych z zabawkami, zawierający tylko rekordy od trzech osób. SAS daje mi dane wyjściowe za kilka minut; w R zajmuje to ponad godzinę. Dziwne. Dzięki temu zestawowi danych zabawek otrzymuję teraz różne oszacowania ustalonych efektów.
Moje pytanie: czy ktoś może rzucić światło na to, dlaczego szacunkowe losowe nachylenia mogą być tak różne dla R i SAS? Czy mogę coś zrobić w R lub SAS, aby zmodyfikować mój kod, aby wywołania dawały podobne wyniki? Wolę zmienić kod w SAS, ponieważ „wierzę”, że moje R szacuje więcej.
Naprawdę martwię się tymi różnicami i chcę dotrzeć do sedna tego problemu!
Moje dane wyjściowe z zabawkowego zestawu danych, który wykorzystuje tylko trzy z 35 osób w pełnym zbiorze danych dla R i SAS, są uwzględnione jako pliki JPEG.
EDYCJA I AKTUALIZACJA:
Jak @JakeWestfall pomógł odkryć, stoki w SAS nie zawierają ustalonych efektów. Kiedy dodam stałe efekty, oto wynik - porównanie nachyleń R ze spadkami SAS dla jednego ustalonego efektu, „PS”, między programami: (Współczynnik wyboru = losowe nachylenie). Zwróć uwagę na zwiększoną zmienność w SAS.
0si i 1s Rmodeluje prawdopodobieństwo odpowiedzi „1”, podczas gdy SAS modeluje prawdopodobieństwo odpowiedzi „0”. Aby model SAS był prawdopodobieństwem „1”, musisz zapisać zmienną odpowiedzi jako use(event='1'). Oczywiście, nawet bez tego, uważam, że nadal powinniśmy oczekiwać takich samych oszacowań wariancji efektu losowego, jak również tych samych oszacowań efektu stałego, choć z odwróconymi znakami.
ranef()funkcji zamiast coef(). Pierwszy daje rzeczywiste efekty losowe, a drugi daje losowe efekty plus wektor efektów stałych. To wyjaśnia, dlaczego liczby przedstawione w twoim poście różnią się, ale nadal istnieje znaczna rozbieżność, której nie potrafię całkowicie wyjaśnić.
IDto nie jest czynnik w R; sprawdź i zobacz, czy to coś zmieni.