Próbuję użyć, lme4::glmer()
aby dopasować dwumianowy uogólniony model mieszany (GLMM) ze zmienną zależną, która nie jest binarna, ale zmienna ciągła od zera do jednego. Można myśleć o tej zmiennej jako o prawdopodobieństwie; w rzeczywistości jest to prawdopodobieństwo zgłaszane przez ludzi (w eksperymencie, który pomagam analizować). Tj. Nie jest to ułamek „dyskretny”, ale zmienna ciągła.
Moje glmer()
połączenie nie działa zgodnie z oczekiwaniami (patrz poniżej). Dlaczego? Co mogę zrobić?
Późniejsza edycja: moja odpowiedź poniżej jest bardziej ogólna niż oryginalna wersja tego pytania, więc zmodyfikowałem pytanie, aby było bardziej ogólne.
Więcej szczegółów
Najwyraźniej można zastosować regresję logistyczną nie tylko dla DV binarnych, ale także dla ciągłego DV od zera do jednego. Rzeczywiście, kiedy biegnę
glm(reportedProbability ~ a + b + c, myData, family="binomial")
Dostaję komunikat ostrzegawczy
Warning message:
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
ale bardzo rozsądne dopasowanie (wszystkie czynniki są kategoryczne, więc mogę łatwo sprawdzić, czy prognozy modelu są zbliżone do średnich między podmiotami i są).
Jednak tak naprawdę chcę użyć
glmer(reportedProbability ~ a + b + c + (1 | subject), myData, family="binomial")
Daje mi identyczne ostrzeżenie, zwraca model, ale ten model jest wyraźnie bardzo zły; oszacowania ustalonych efektów są bardzo dalekie od glm()
tych i od środków obejmujących wiele przedmiotów. (I muszę dołączyć glmerControl(optimizer="bobyqa")
do glmer
połączenia, w przeciwnym razie w ogóle się nie zbiegnie).
glmmadmb(reportedProbability ~ a + b + c + (1 | subject), myData, family="beta")
, otrzymuję prawidłowe dopasowanie i rozsądne przedziały ufności, ale ostrzeżenie o zbieżności nie powiodło się : - / Próbuję wymyślić, jak zwiększyć liczbę iteracji. Beta może dla mnie działać, ponieważ nie mam przypadków DV = 0 lub DV = 1.+ (1 | rowid)
do mojego glitter call, co daje stabilne oszacowania i stabilne przedziały ufności, niezależnie od mojego wyboru wagi (próbowałem 100 i 500). Próbowałem także uruchomić lmer na logit (reportProbability) i otrzymuję prawie dokładnie to samo. Oba rozwiązania wydają się działać dobrze! Beta MM z glmmadmb daje również bardzo bliskie wyniki, ale z jakiegoś powodu nie zbiega się całkowicie i trwa bez końca. Zastanów się nad opublikowaniem odpowiedzi zawierającej listę tych opcji i wyjaśnienie nieco różnic i zalet / wad! (Wszystkie podane przeze mnie przedziały ufności to Wald).Odpowiedzi:
Warto zacząć od prostszego przypadku bez przypadkowych efektów.
Istnieją cztery sposoby radzenia sobie z ciągłą zmienną odpowiedzi od zera do jednego, która zachowuje się jak ułamek lub prawdopodobieństwo ( jest to nasz najbardziej kanoniczny / oceniany / przeglądany wątek na ten temat, ale niestety nie wszystkie cztery opcje są tam omówione):
n
Logit przekształca odpowiedź i stosuje regresję liniową. Zazwyczaj nie jest to zalecane.
Dopasuj model dwumianowy, a następnie oblicz standardowe błędy z uwzględnieniem nadmiernej dyspersji. Standardowe błędy można obliczyć na różne sposoby:
(a) skalowane błędy standardowe za pomocą oszacowania nadmiernej dyspersji ( jeden , dwa ). Nazywa się to „quasi-dwumianowym” GLM.
(b) solidne błędy standardowe za pomocą estymatora wielowarstwowego ( jeden , dwa , trzy , cztery ). W ekonometrii jest to nazywane „logarytmem ułamkowym”.
(A) i (b) nie są identyczne (patrz ten komentarz oraz sekcje 3.4.1 i 3.4.2 w tej książce , i ten post SO, a także ten i ten ), ale zwykle dają podobne wyniki. Opcja (a) jest realizowana w
glm
następujący sposób:Te same cztery sposoby są dostępne z efektami losowymi.
Za pomocą
weights
argumentu ( jeden , dwa ):Zgodnie z drugim linkiem powyżej dobrym pomysłem może być modelowanie naddyspersji, patrz tam (a także # 4 poniżej).
Używając mieszanego modelu beta:
lub
Jeśli w danych odpowiedzi są dokładne zera lub jedynki, wówczas można użyć modelu beta z zerowym / jednym zawyżeniem
glmmTMB
.Za pomocą transformacji logit odpowiedzi:
Uwzględnianie nadmiernej dyspersji w modelu dwumianowym. Wykorzystuje to inną sztuczkę: dodanie losowego efektu dla każdego punktu danych:
Z jakiegoś powodu nie działa to poprawnie, ponieważ
glmer()
narzeka na liczbę całkowitąp
i daje bzdury. Rozwiązaniem, które wymyśliłem, jest użycie fałszywej stałejweights=k
i upewnienie się, żep*k
zawsze jest liczbą całkowitą. Wymaga to zaokrąglenia,p
ale wybranie odpowiedniok
dużego rozmiaru nie powinno mieć większego znaczenia. Wyniki nie wydają się zależeć od wartościk
.Późniejsza aktualizacja (styczeń 2018 r.): Może to być nieprawidłowe podejście. Zobacz dyskusję tutaj . Muszę to zbadać bardziej.
W moim konkretnym przypadku opcja nr 1 nie jest dostępna.
Opcja nr 2 jest bardzo wolna i ma problemy z konwergencją:Aktualizacja: Próbowałemglmmadmb
uruchomienie zajmuje pięć do dziesięciu minut (i nadal narzeka, że się nie zbiegło!), Natomiastlmer
działa w ułamku sekundy iglmer
zajmuje kilka sekund.glmmTMB
zgodnie z sugestiami @BenBolker i działa prawie tak szybko, jakglmer
bez problemów z konwergencją. Więc tego będę używać.Opcje 3 i 4 dają bardzo podobne szacunki i bardzo podobne przedziały ufności Walda (uzyskane z
confint
). Nie jestem wielkim fanem nr 3, ponieważ to rodzaj oszustwa. A # 4 wydaje się nieco zrzędliwy.Ogromne podziękowania dla @Aarona, który w swoim komentarzu wskazał mi na # 3 i # 4.
źródło
devtools::install_github("glmmTMB/glmmTMB",sub="glmmTMB")
użyciuglmmTMB(p ~ a+b+c + (1|subject), myData, family=list(family="beta",link="logit"))
powinno działać ...glmmTMB
jest szybszy i bardziej stabilny niżglmmADMB
i (nieco) bardziej aktywny rozwój, chociaż nie tak dojrzały.