Dlaczego SAS PROC GLIMMIX daje mi BARDZO różne losowe nachylenia niż glmer (lme4) dla dwumianowego glmm

12

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.

Wyjście R. Wyjście SAS 1 Wyjście SAS 2 Wyjście SAS 3


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.

R vs SAS dla PS

Nova
źródło
Zauważam, że IDto nie jest czynnik w R; sprawdź i zobacz, czy to coś zmieni.
Aaron opuścił Stack Overflow
Widzę, że dla obu tych wartości korzystasz z aproksymacji Laplace'a. Jakie są ich wyniki wiarygodności dziennika?
usεr11852
1
Czy sprawdziłeś, że modelujesz zmienną zależną w tym samym kierunku?
Peter Flom - Przywróć Monikę
1
Nawiasem mówiąc, Peter zmierza do tego, że domyślnie z danymi dwumianowymi oznaczonymi jako 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.
Jake Westfall,
1
@EricaN Jedną rzecz, o której mi właśnie przypomniałeś, jest to, że powinieneś porównywać losowe efekty z R do tych w SAS, używając tej 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ć.
Jake Westfall,

Odpowiedzi:

11

Wydaje się, że nie powinienem oczekiwać, że losowe nachylenia będą podobne między pakietami, zgodnie z Zhang i in. 2011. W swoim artykule na temat dopasowania uogólnionych liniowych modeli mieszanych efektów dla odpowiedzi binarnych przy użyciu różnych pakietów statystycznych opisują:

Abstrakcyjny:

Uogólniony liniowy model efektów mieszanych (GLMM) jest popularnym paradygmatem rozszerzania modeli danych przekrojowych na ustawienie wzdłużne. W przypadku modelowania odpowiedzi binarnych różne pakiety oprogramowania, a nawet różne procedury w pakiecie mogą dawać zupełnie inne wyniki. W tym raporcie opisujemy podejścia statystyczne, które leżą u podstaw tych różnych procedur, i omawiamy ich mocne i słabe strony, gdy są stosowane w celu dopasowania skorelowanych odpowiedzi binarnych. Następnie ilustrujemy te rozważania, stosując procedury zaimplementowane w niektórych popularnych pakietach oprogramowania do danych symulowanych i rzeczywistych badań. Nasze wyniki symulacji wskazują na brak niezawodności większości rozważanych procedur, co niesie ze sobą znaczące implikacje dla zastosowania takich popularnych pakietów oprogramowania w praktyce.

Mam nadzieję, że @BenBolker wraz z zespołem rozważy moje pytanie jako głos za tym, by R zawierał empiryczne błędy standardowe i zdolność kwadratury Gaussa-Hermity do modeli z kilkoma losowymi wartościami nachylenia, aby migotać, ponieważ wolę interfejs R i chciałbym móc zastosować kilka dalszych analiz w tym programie. Na szczęście, mimo że R i SAS nie mają porównywalnych wartości dla losowych stoków, ogólne trendy są takie same. Dziękuję wszystkim za Twój wkład. Naprawdę doceniam poświęcony czas i uwagę!

Nova
źródło
przepraszam: co to jest „standardowy błąd standardowy”? Masz na myśli standardowe błędy składników wariancji? Czy masz na myśli standardowe błędy kanapkowe?
Ben Bolker
przepraszam ... oznaczało SE empiryczne / kanapkowe. Zredagowałem swoją odpowiedź.
Nova,
@BenBolker Czy to kiedykolwiek zostało włączone?
Lepidopterist
Nie. Wciąż próbuję wymyślić, w jaki sposób będę wspierać taki rozwój, ponieważ technicznie nie jest to część mojego programu badawczego ...
Ben Bolker
4

Mieszanka odpowiedzi i komentarza / więcej pytań:

Dopasowałem zestaw danych „zabawki” do trzech różnych opcji optymalizacji. (* Uwaga 1: Prawdopodobnie bardziej użyteczne byłoby dla celów porównawczych wykonanie małego zestawu danych poprzez podpróbkowanie od roku każdego roku i identyfikatora, niż przez podpróbkowanie zmiennych grupujących. W tej chwili wiemy, że GLMM nie będzie działał szczególnie dobrze przy tak małej liczbie zmiennych grupujących. Możesz to zrobić za pomocą:

library(plyr)
subdata <- ddply(fulldata,c("year","id"),
    function(x) x[sample(nrow(x),size=round(nrow(x)*0.1)),])

Kod dopasowania partii:

Ntoy <- readRDS("Newton_toy.RDS")
library(lme4)
fitfun <- function(opt) {
    tt <- system.time(fit1 <- glmer(use ~  ps + th + rs + hw +
                                    (1 + ps + th + rs + hw |id/year),
                                    family = binomial, data = Ntoy,
                                    control=glmerControl(optimizer=opt),
                                    verbose=100))
    return(list(time=tt,fit=fit1))
}

opts <- c("nloptwrap","nlminbwrap","bobyqa")
## use for() instead of lapply so we can checkpoint more easily
res <- setNames(vector("list",length(opts)),opts)
for (i in opts) {
    res[[i]] <- fitfun(i)
    save("res",file="Newton_batch.RData")
}

Następnie przeczytałem wyniki w nowej sesji:

load("Newton_batch.RData")
library(lme4)

Upływający czas i dewiacja:

cbind(time=unname(sapply(res,function(x) x$time["elapsed"])),
          dev=sapply(res,function(x) deviance(x$fit)))
##                time      dev
## nloptwrap  1001.824 6067.706
## nlminbwrap 3495.671 6068.730
## bobyqa     4945.332 6068.731

Odchylenia te są znacznie niższe niż odchylenie zgłoszone przez OP z R (6101.7), i nieco poniżej odchyleń zgłoszonych przez OP z SAS (6078.9), chociaż porównywanie odchyleń między pakietami nie zawsze jest rozsądne.

Byłem naprawdę zaskoczony, że SAS zebrał tylko około 100 ocen funkcji!

Czasy wahają się od 17 minut ( nloptwrap) do 80 minut ( bobyqa) na Macbooku Pro, zgodnie z doświadczeniem OP. Odchylenie jest nieco lepsze nloptwrap.

round(cbind(sapply(res,function(x) fixef(x$fit))),3)
##             nloptwrap nlminbwrap bobyqa
## (Intercept)    -5.815     -5.322 -5.322
## ps             -0.989      0.171  0.171
## th             -0.033     -1.342 -1.341
## rs              1.361     -0.140 -0.139
## hw             -2.100     -2.082 -2.082

Odpowiedzi wyglądają zupełnie inaczej nloptwrap- chociaż standardowe błędy są dość duże ...

round(coef(summary(res[[1]]$fit)),3)
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -5.815      0.750  -7.750    0.000
## ps            -0.989      1.275  -0.776    0.438
## th            -0.033      2.482  -0.013    0.989
## rs             1.361      2.799   0.486    0.627
## hw            -2.100      0.490  -4.283    0.000

(kod tutaj zawiera ostrzeżenia o year:idtym, że powinienem wyśledzić)

Ciąg dalszy nastąpi ... ?

Ben Bolker
źródło
czy byłoby bardziej pomocne, gdybym wysłał ci pełny zestaw danych? Jedynym problemem jest to, że konwergencja zajmuje około 9 godzin z pełnym zestawem danych, więc twoja sugestia dotycząca próbkowania jest dobra. Próbowałem przekształcić dane za pomocą transformacji dziennika, ale skumulowany wykres resztkowy jest nadal brzydki - czy uważasz, że wykres resztkowy wyjaśnia część problemu związanego z tymi danymi? Wreszcie - czy twoje wyniki w SAS były podobne do R?
Nova