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

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.

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ę!

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