Jak dopasować model mieszany ze zmienną odpowiedzi od 0 do 1?

15

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 glmerpołączenia, w przeciwnym razie w ogóle się nie zbiegnie).

ameba mówi Przywróć Monikę
źródło
1
A może najpierw przekształcisz prawdopodobieństwa? Czy możesz uzyskać coś, co jest bliżej normalnie dystrybuowane, powiedzmy, transformacja logitów? A może arcsin-sqrt? To byłoby moje preferencje zamiast używania blasku. Lub w swoim rozwiązaniu do hakowania możesz również spróbować dodać losowy efekt dla każdej obserwacji, aby uwzględnić niedostateczną dyspersję ze względu na wybór wag.
Aaron opuścił Stack Overflow
Dzięki. Tak, mogę zalogować DV, a następnie użyć mieszanego modelu Gaussa (lmer), ale jest to również rodzaj hackowania i przeczytałem, że nie jest to zalecane. Spróbuję losowego efektu dla każdej obserwacji! W tej chwili próbuję mieszanego modelu beta; lme4 nie może sobie z tym poradzić, ale glmmadmb potrafi. Kiedy biegam 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.
ameba mówi Przywróć Monikę
Nie wiem dla glmera, ale dla glm może to pomóc: stats.stackexchange.com/questions/164120/… :
1
@Aaron: Próbowałem dodać + (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).
Amoeba mówi: Przywróć Monikę
1
I są absolutnie pewni swojej wartości, takiej jak 0,9, czy też mają jakiś „margines błędu”? Czy możesz założyć, że zaufanie zgłaszane przez różne podmioty jest równie dokładne?

Odpowiedzi:

21

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):

  1. p=m/nnnN.

    glm(p ~ a+b+c, myData, family="binomial", weights=n)
  2. pp01

    betareg(p ~ a+b+c, myData)
  3. Logit przekształca odpowiedź i stosuje regresję liniową. Zazwyczaj nie jest to zalecane.

    lm(log(p/(1-p)) ~ a+b+c, myData)
  4. 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 glmnastępujący sposób:

    glm(p ~ a+b+c, myData, family="quasibinomial")

Te same cztery sposoby są dostępne z efektami losowymi.

  1. Za pomocą weightsargumentu ( jeden , dwa ):

    glmer(p ~ a+b+c + (1|subject), myData, family="binomial", weights=n)

    Zgodnie z drugim linkiem powyżej dobrym pomysłem może być modelowanie naddyspersji, patrz tam (a także # 4 poniżej).

  2. Używając mieszanego modelu beta:

    glmmadmb(p ~ a+b+c + (1|subject), myData, family="beta")

    lub

    glmmTMB(p ~ a+b+c + (1|subject), myData, 
            family=list(family="beta",link="logit"))

    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.

  3. Za pomocą transformacji logit odpowiedzi:

    lmer(log(p/(1-p)) ~ a+b+c + (1|subject), myData)
  4. Uwzględnianie nadmiernej dyspersji w modelu dwumianowym. Wykorzystuje to inną sztuczkę: dodanie losowego efektu dla każdego punktu danych:

    myData$rowid = as.factor(1:nrow(myData))
    glmer(p ~ a+b+c + (1|subject) + (1|rowid), myData, family="binomial",
          glmerControl(optimizer="bobyqa"))

    Z jakiegoś powodu nie działa to poprawnie, ponieważ glmer()narzeka na liczbę całkowitą pi daje bzdury. Rozwiązaniem, które wymyśliłem, jest użycie fałszywej stałej weights=ki upewnienie się, że p*kzawsze jest liczbą całkowitą. Wymaga to zaokrąglenia, pale wybranie odpowiednio kdużego rozmiaru nie powinno mieć większego znaczenia. Wyniki nie wydają się zależeć od wartości k.

    k = 100
    glmer(round(p*k)/k ~ a+b+c + (1|subject) + (1|rowid), myData, 
          family="binomial", weights=rowid*0+k, glmerControl(optimizer="bobyqa"))

    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ą: glmmadmburuchomienie zajmuje pięć do dziesięciu minut (i nadal narzeka, że ​​się nie zbiegło!), Natomiast lmerdziała w ułamku sekundy i glmerzajmuje kilka sekund. Aktualizacja: Próbowałem glmmTMBzgodnie z sugestiami @BenBolker i działa prawie tak szybko, jak glmerbez 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.

ameba mówi Przywróć Monikę
źródło
1
Ładna odpowiedź, dobrze wyjaśniona i połączona z modelami bez efektów losowych. Nie nazwałbym jednak oszustwem # 3 (transformacja), tego rodzaju transformacje są bardzo powszechne w tego typu analizach. Powiedziałbym zamiast tego, że zarówno nr 3, jak i 4 przyjmują założenia dotyczące związku między rozkładem danych, a więc także związku między średnią a wariancją, i tylko dlatego, że nr 4 modeluje w skali, w której dane zebrane w dniu nie oznacza, że ​​te założenia będą lepsze.
Aaron opuścił Stack Overflow
1
# 3 zakłada, że ​​logarytm prawdopodobieństw jest normalny ze stałą wariancją, a # 4 zakłada, że ​​wariancja jest proporcjonalna do p (1-p). Z twojego opisu dopasowania wydają się być na tyle podobne, że nie mają większego znaczenia. A # 3 jest prawie na pewno bardziej standardowy (w zależności od odbiorców), więc jeśli diagnostyka jest rozsądna, to właśnie ja bym wolał.
Aaron opuścił Stack Overflow
1
inną możliwością jest użycie glmmTMB ; po instalacji przy devtools::install_github("glmmTMB/glmmTMB",sub="glmmTMB")użyciu glmmTMB(p ~ a+b+c + (1|subject), myData, family=list(family="beta",link="logit"))powinno działać ...
Ben Bolker
@BenBolker Thanks! Czy istnieje powód, aby preferować glmmTMB zamiast glmmADMB (dla modeli beta) lub odwrotnie? Czy jeden z tych pakietów jest nowszy czy bardziej aktywnie rozwijany? Poza tym, czy mogę zapytać, jakie podejście spośród wymienionych w tej odpowiedzi - gluss gaussowski po transformacji logit, beta glmm lub dwumianowy glmm z terminem (1 | rowid) - czy uważasz, że ogólnie jest bardziej odpowiedni?
ameba mówi Przywróć Monikę
1
Wolę wersję beta GLMM, jeśli jest to wykonalne - to model statystyczny, który ma mierzyć zmiany proporcji między zmiennymi towarzyszącymi / grupami. glmmTMBjest szybszy i bardziej stabilny niż glmmADMBi (nieco) bardziej aktywny rozwój, chociaż nie tak dojrzały.
Ben Bolker,