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

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

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/n$$n$n$N$

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

2. $p$$p$$0$$1$

   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.