Mam pytania dotyczące specyfikacji i interpretacji GLMM. 3 pytania są zdecydowanie statystyczne, a 2 dotyczą bardziej R. Piszę tutaj, ponieważ ostatecznie myślę, że problemem jest interpretacja wyników GLMM.
Obecnie próbuję dopasować GLMM. Korzystam z danych amerykańskiego spisu ludności z bazy danych podłużnych dróg . Moje obserwacje są traktatami spisowymi. Moją zmienną zależną jest liczba wolnych lokali mieszkalnych i interesuje mnie związek między wakatem a zmiennymi społeczno-ekonomicznymi. Przykład jest prosty, wykorzystując tylko dwa ustalone efekty: procent populacji niebiałej (rasa) i mediana dochodu gospodarstwa domowego (klasa) plus ich interakcja. Chciałbym uwzględnić dwa zagnieżdżone efekty losowe: traktaty w ciągu dziesięcioleci i dziesięcioleci, tj. (Dekada / traktat). Rozważam te losowe w celu kontroli autokorelacji przestrzennej (tj. Między traktami) i czasowej (tj. Między dekadami). Jednak interesuje mnie dekada jako stały efekt, więc włączam ją również jako stały czynnik.
Ponieważ moja zmienna niezależna jest nieujemną zmienną zliczającą liczby całkowite, próbowałem dopasować poissony i ujemne dwumianowe GLMM. Korzystam z dziennika wszystkich jednostek mieszkaniowych jako przesunięcia. Oznacza to, że współczynniki są interpretowane jako wpływ na współczynnik pustostanów, a nie całkowitą liczbę wolnych domów.
Obecnie mam wyniki dla Poissona i ujemnego dwumianowego GLMM oszacowanego za pomocą glmer i glmer.nb z lme4 . Interpretacja współczynników ma dla mnie sens na podstawie mojej wiedzy na temat danych i obszaru badań.
Jeśli chcesz dane i skrypt , są one na moim Githubie . Skrypt zawiera więcej opisowych badań, które przeprowadziłem przed zbudowaniem modeli.
Oto moje wyniki:
Model Poissona
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson ( log )
Formula: R_VAC ~ decade + P_NONWHT + a_hinc + P_NONWHT * a_hinc + offset(HU_ln) + (1 | decade/TRTID10)
Data: scaled.mydata
AIC BIC logLik deviance df.resid
34520.1 34580.6 -17250.1 34500.1 3132
Scaled residuals:
Min 1Q Median 3Q Max
-2.24211 -0.10799 -0.00722 0.06898 0.68129
Random effects:
Groups Name Variance Std.Dev.
TRTID10:decade (Intercept) 0.4635 0.6808
decade (Intercept) 0.0000 0.0000
Number of obs: 3142, groups: TRTID10:decade, 3142; decade, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.612242 0.028904 -124.98 < 2e-16 ***
decade1980 0.302868 0.040351 7.51 6.1e-14 ***
decade1990 1.088176 0.039931 27.25 < 2e-16 ***
decade2000 1.036382 0.039846 26.01 < 2e-16 ***
decade2010 1.345184 0.039485 34.07 < 2e-16 ***
P_NONWHT 0.175207 0.012982 13.50 < 2e-16 ***
a_hinc -0.235266 0.013291 -17.70 < 2e-16 ***
P_NONWHT:a_hinc 0.093417 0.009876 9.46 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) dc1980 dc1990 dc2000 dc2010 P_NONWHT a_hinc
decade1980 -0.693
decade1990 -0.727 0.501
decade2000 -0.728 0.502 0.530
decade2010 -0.714 0.511 0.517 0.518
P_NONWHT 0.016 0.007 -0.016 -0.015 0.006
a_hinc -0.023 -0.011 0.023 0.022 -0.009 0.221
P_NONWHT:_h 0.155 0.035 -0.134 -0.129 0.003 0.155 -0.233
convergence code: 0
Model failed to converge with max|grad| = 0.00181132 (tol = 0.001, component 1)
Negatywny model dwumianowy
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(25181.5) ( log )
Formula: R_VAC ~ decade + P_NONWHT + a_hinc + P_NONWHT * a_hinc + offset(HU_ln) + (1 | decade/TRTID10)
Data: scaled.mydata
AIC BIC logLik deviance df.resid
34522.1 34588.7 -17250.1 34500.1 3131
Scaled residuals:
Min 1Q Median 3Q Max
-2.24213 -0.10816 -0.00724 0.06928 0.68145
Random effects:
Groups Name Variance Std.Dev.
TRTID10:decade (Intercept) 4.635e-01 6.808e-01
decade (Intercept) 1.532e-11 3.914e-06
Number of obs: 3142, groups: TRTID10:decade, 3142; decade, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.612279 0.028946 -124.79 < 2e-16 ***
decade1980 0.302897 0.040392 7.50 6.43e-14 ***
decade1990 1.088211 0.039963 27.23 < 2e-16 ***
decade2000 1.036437 0.039884 25.99 < 2e-16 ***
decade2010 1.345227 0.039518 34.04 < 2e-16 ***
P_NONWHT 0.175216 0.012985 13.49 < 2e-16 ***
a_hinc -0.235274 0.013298 -17.69 < 2e-16 ***
P_NONWHT:a_hinc 0.093417 0.009879 9.46 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) dc1980 dc1990 dc2000 dc2010 P_NONWHT a_hinc
decade1980 -0.693
decade1990 -0.728 0.501
decade2000 -0.728 0.502 0.530
decade2010 -0.715 0.512 0.517 0.518
P_NONWHT 0.016 0.007 -0.016 -0.015 0.006
a_hinc -0.023 -0.011 0.023 0.022 -0.009 0.221
P_NONWHT:_h 0.154 0.035 -0.134 -0.129 0.003 0.155 -0.233
Testy Poissona DHARMa
One-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.044451, p-value = 8.104e-06
alternative hypothesis: two-sided
DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
data: simulationOutput
ratioObsExp = 1.3666, p-value = 0.159
alternative hypothesis: more
Negatywne dwumianowe testy DHARMa
One-sample Kolmogorov-Smirnov test
data: simulationOutput$scaledResiduals
D = 0.04263, p-value = 2.195e-05
alternative hypothesis: two-sided
DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model
data: simulationOutput2
ratioObsExp = 1.376, p-value = 0.174
alternative hypothesis: more
Działki DHARMa
Poissona
Ujemny dwumianowy
Pytania statystyczne
Ponieważ wciąż zastanawiam się nad GLMM, nie mam pewności co do specyfikacji i interpretacji. Mam parę pytań:
Wygląda na to, że moje dane nie obsługują przy użyciu modelu Poissona, dlatego lepiej jest mi z dwumianem ujemnym. Jednak konsekwentnie otrzymuję ostrzeżenia, że moje ujemne modele dwumianowe osiągają limit iteracji, nawet gdy zwiększę limit maksymalny. „W pliku theta.ml (Y, mu, wagi = obiekt @ lub $ wagi, limit = limit, osiągnięty limit iteracji.” Dzieje się tak przy użyciu kilku różnych specyfikacji (tj. Modeli minimalnych i maksymalnych dla efektów stałych i losowych). Próbowałem również usunąć wartości odstające w mojej grupie zależnej (brutto, wiem!), Ponieważ najwyższy 1% wartości to bardzo duże wartości odstające (dolny zakres 99% od 0-1012, górny 1% od 1013-5213). nie mają żadnego wpływu na iteracje i mają bardzo niewielki wpływ na współczynniki. Nie zamieszczam tutaj tych szczegółów. Współczynniki między Poissonem a dwumianem ujemnym są również bardzo podobne. Czy ten brak konwergencji stanowi problem? Czy ujemny model dwumianowy jest dobrze dopasowany? Uruchomiłem również ujemny model dwumianowyAllFit i nie wszystkie optymalizatory generują to ostrzeżenie (bobyqa, Nelder Mead i nlminbw nie zrobiły tego).
Wariancja dla mojego stałego efektu dekady jest konsekwentnie bardzo niska lub wynosi 0. Rozumiem, że to może oznaczać, że model jest przeregulowany. Usunięcie dekady ze stałych efektów zwiększa dziesięcioletnią wariancję efektu losowego do 0,2620 i nie ma większego wpływu na współczynniki efektu stałego. Czy jest coś złego w pozostawieniu tego? Świetnie interpretuję to jako po prostu niepotrzebne do wyjaśnienia wariancji obserwacji.
Czy te wyniki wskazują, że powinienem wypróbować modele z zerowym napełnieniem? DHARMa wydaje się sugerować, że inflacja zerowa może nie być problemem. Jeśli uważasz, że powinienem spróbować mimo to, patrz poniżej.
Pytania R.
Byłbym skłonny wypróbować modele z zerowym napełnieniem, ale nie jestem pewien, który pakiet sugeruje zagnieżdżone losowe efekty dla z napompowanego zerowo Poissona i ujemnych dwumianowych GLMM. Użyłbym glmmADMB do porównania AIC z modelami z napompowaniem zerowym, ale jest ograniczony do pojedynczego efektu losowego, więc nie działa dla tego modelu. Mógłbym spróbować MCMCglmm, ale nie znam statystyki bayesowskiej, więc to też nie jest atrakcyjne. Jakieś inne opcje?
Czy mogę wyświetlić współczynniki potęgowane w podsumowaniu (modelu), czy też muszę to zrobić poza podsumowaniem, tak jak to tutaj zrobiłem?
źródło
decade
zarówno ustalonych, jak i losowych nie ma sensu. Albo ustaw go jako stały i uwzględniaj tylko(1 | decade:TRTID10)
losowo (co jest równoznaczne z(1 | TRTID10)
założeniem, że twójTRTID10
nie ma takich samych poziomów przez różne dziesięciolecia), albo usuń go ze stałych efektów. Mając tylko 4 poziomy, lepiej być naprawionym: zwykle zaleca się dopasowanie efektów losowych, jeśli jeden ma 5 poziomów lub więcej.bobyqa
optymalizator i nie wygenerował on żadnego ostrzeżenia. W czym zatem problem? Po prostu użyjbobyqa
.bobyqa
wynika, że lepiej łączy się z domyślnym optymalizatorem (i myślę, że gdzieś przeczytałem, że stanie się domyślny w przyszłych wersjachlme4
). Nie sądzę, że musisz się martwić o brak konwergencji z domyślnym optymalizatorem, jeśli się on zbiegabobyqa
.Odpowiedzi:
Uważam, że istnieje kilka ważnych problemów, którymi należy się zająć przy szacowaniu.
Z tego, co zebrałem, badając twoje dane, twoje jednostki nie są pogrupowane geograficznie, tj. Obszary spisowe w obrębie hrabstw. Zatem użycie traktów jako czynnika grupującego nie jest odpowiednie do uchwycenia przestrzennej heterogeniczności, ponieważ oznacza to, że masz taką samą liczbę osobników jak grupy (lub inaczej mówiąc, wszystkie twoje grupy mają tylko jedną obserwację). Korzystanie ze strategii modelowania wielopoziomowego pozwala nam oszacować wariancję na poziomie indywidualnym, jednocześnie kontrolując wariancję między grupami. Ponieważ w każdej grupie jest tylko jedna osoba, wariancja między grupami jest taka sama, jak wariancja na poziomie indywidualnym, co przeczy celowi podejścia wielopoziomowego.
Z drugiej strony współczynnik grupowania może reprezentować powtarzające się pomiary w czasie. Na przykład, w przypadku badań podłużnych, wyniki matematyczne danej osoby mogą być odtwarzane corocznie, dlatego mielibyśmy roczną wartość dla każdego ucznia przez n lat (w tym przypadku czynnikiem grupującym jest uczeń, ponieważ mamy n liczba obserwacji „zagnieżdżonych” wśród studentów). W twoim przypadku powtórzyłeś pomiary dla każdego obszaru spisu ludności przez
decade
. Zatem możesz użyć swojejTRTID10
zmiennej jako czynnika grupującego, aby uchwycić „wariancję między dekadami”. Prowadzi to do 3142 obserwacji zagnieżdżonych w 635 traktach, z około 4 i 5 obserwacjami na jeden spis ludności.Jak wspomniano w komentarzu, użycie
decade
jako czynnika grupującego nie jest zbyt odpowiednie, ponieważ masz tylko około 5 dekad dla każdego obszaru spisu, a ich efekt można lepiej uchwycić, wprowadzającdecade
jako zmienną towarzyszącą.Po drugie, aby ustalić, czy dane powinny być modelowane przy użyciu modelu dwumianowego Poissona, czy ujemnego (lub podejścia o zawyżonym zeru). Zastanów się nad nadmierną dyspersją danych. Podstawową cechą rozkładu Poissona jest równomierność, co oznacza, że średnia jest równa wariancji rozkładu. Patrząc na twoje dane, jest całkiem jasne, że istnieje duża naddyspersja. Rozbieżności są znacznie większe niż średnie.
Niemniej jednak, aby ustalić, czy ujemny dwumian jest bardziej odpowiedni statystycznie, standardową metodą jest wykonanie testu współczynnika prawdopodobieństwa między modelem Poissona i ujemnym dwumianowym modelem, co sugeruje, że negbin jest lepiej dopasowany.
Po ustaleniu tego, kolejny test mógłby rozważyć, czy podejście wielopoziomowe (model mieszany) jest uzasadnione przy użyciu podobnego podejścia, co sugeruje, że wersja wielopoziomowa zapewnia lepsze dopasowanie. (Podobny test można zastosować do porównania dopasowania połysku przy założeniu rozkładu Poissona do obiektu glmer.nb, o ile modele są takie same.)
Jeśli chodzi o szacunki modeli Poissona i NB, to w rzeczywistości oczekuje się, że będą bardzo do siebie podobne, przy czym głównym rozróżnieniem są błędy standardowe, tj. Jeżeli występuje nadmierna dyspersja, model Poissona ma tendencję do dostarczania stronniczych błędów standardowych. Biorąc twoje dane jako przykład:
Zauważ, że wszystkie oszacowania współczynników są bardzo podobne, a główną różnicą jest tylko znaczenie jednej z Twoich zmiennych towarzyszących, a także różnica wariancji efektów losowych, co sugeruje, że wariancja na poziomie jednostki została przechwycona przez parametr nadmiernej dyspersji w nb model (
theta
wartość w obiekcie glmer.nb) przechwytuje niektóre wariancje między traktami uchwycone przez efekty losowe.W odniesieniu do współczynników wykładniczych (i powiązanych przedziałów ufności) można użyć następujących opcji:
Ostatnie przemyślenia dotyczące zerowej inflacji. Nie ma wielopoziomowej implementacji (przynajmniej jestem tego świadomy) modelu poissona lub negbina nadmuchanego zera, który pozwala określić równanie dla nadmuchanego zera składnika mieszaniny.
glmmADMB
modelu pozwala na oszacowanie parametru inflacji stała zerowy. Alternatywnym podejściem byłoby użyciezeroinfl
funkcji wpscl
pakiecie, chociaż nie obsługuje ona modeli wielopoziomowych. W ten sposób można porównać dopasowanie dwumianu ujemnego jednopoziomowego do dwumianu napełnionego ujemnie jednopoziomowego. Istnieje prawdopodobieństwo, że jeśli zerowa inflacja nie będzie znacząca dla modeli jednopoziomowych, jest prawdopodobne, że nie będzie znacząca dla specyfikacji wielopoziomowej.Uzupełnienie
Jeśli obawiasz się autokorelacji przestrzennej, możesz to kontrolować za pomocą jakiejś formy regresji ważonej geograficznie (chociaż uważam, że wykorzystuje to dane punktowe, a nie obszary). Alternatywnie możesz pogrupować traktaty spisowe według dodatkowego czynnika grupującego (stany, powiaty) i uwzględnić to jako efekt losowy. Wreszcie i nie jestem pewien, czy jest to w pełni wykonalne, może być możliwe włączenie zależności przestrzennej, na przykład wykorzystując średnią liczbę
R_VAC
sąsiadów pierwszego rzędu jako zmienną towarzyszącą. W każdym razie przed takimi podejściami rozsądne byłoby ustalenie, czy rzeczywiście występuje autokorelacja przestrzenna (przy użyciu I, testów LISA i podobnych metod Global Morana).źródło
brms
może pasować do zerowych pompowanych modeli dwumianowych z efektami losowymi.brms
i porównać go z modelem glmer.nb, jak opisano powyżej. Spróbuję również uwzględnić miejsce zdefiniowane w spisie (zasadniczo gmina, ~ 170 grup) jako czynnik grupujący dla efektów losowych (tylko 5 hrabstw w danych, więc nie użyję tego). Będę również testować autokorelację przestrzenną reszt za pomocą Global Morana I. Opiszę, kiedy to zrobię.brms
i porównywanie ich.