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

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

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

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

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

5. Czy mogę wyświetlić współczynniki potęgowane w podsumowaniu (modelu), czy też muszę to zrobić poza podsumowaniem, tak jak to tutaj zrobiłem?

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ć swojej TRTID10zmiennej 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 decadejako 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ąc decadejako 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.

library(dplyr)    
 dispersionstats <- scaled.mydata %>%
 + group_by(decade) %>%
 + summarise(
 + means = mean(R_VAC),
 + variances = var(R_VAC),
 + ratio = variances/means)

##   dispersionstats
##   # A tibble: 5 x 5
##   decade     means variances     ratio 
##    <int>     <dbl>     <dbl>     <dbl> 
## 1   1970  45.43513   4110.89  90.47822 
## 2   1980 103.52365  17323.34 167.33707 
## 3   1990 177.68038  62129.65 349.67087 
## 4   2000 190.23150  91059.60 478.67784 
## 5   2010 247.68246 126265.60 509.78821 

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.

library(MASS)
library(lmtest)

modelformula <- formula(R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln))

poismodel <- glm(modelformula, data = scaled.mydata, family = "poisson")   
nbmodel <- glm.nb(modelformula, data = scaled.mydata)

lrtest(poismodel, nbmodel)

## Likelihood ratio test

##  Model 1: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)  
## Model 2: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   8 -154269
## 2   9  -17452  1 273634  < 2.2e-16 ***
##  ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

library(lme4)

glmmformula <- update(modelformula, . ~ . + (1|TRTID10))

nbglmm <- glmer.nb(glmmformula, data = scaled.mydata)

lrtest(nbmodel, nbglmm)

## Model 1: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)
## Model 2: R_VAC ~ factor(decade) + P_NONWHT + a_hinc + (1 | TRTID10) +
##     P_NONWHT:a_hinc + offset(HU_ln)
##   #Df LogLik Df Chisq Pr(>Chisq)
## 1   9 -17452
## 2  10 -17332  1 239.3  < 2.2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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:

poissonglmm <- glmer(glmmformula, data = scaled.mydata)
summary(poissonglmm)

## Random effects:
##  Groups  Name        Variance Std.Dev.
## TRTID10 (Intercept) 0.2001   0.4473
## Number of obs: 3142, groups:  TRTID10, 635

## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -2.876013   0.020602 -139.60   <2e-16 ***
## factor(decade)1980  0.092597   0.007602   12.18   <2e-16 ***
## factor(decade)1990  0.903543   0.007045  128.26   <2e-16 ***
## factor(decade)2000  0.854821   0.006913  123.65   <2e-16 ***
## factor(decade)2010  0.986126   0.006723  146.67   <2e-16 ***
## P_NONWHT           -0.125500   0.014007   -8.96   <2e-16 ***
## a_hinc             -0.107335   0.001480  -72.52   <2e-16 ***
## P_NONWHT:a_hinc     0.160937   0.003117   51.64   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(nbglmm)
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  TRTID10 (Intercept) 0.09073  0.3012
## Number of obs: 3142, groups:  TRTID10, 635

## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -2.797861   0.056214  -49.77  < 2e-16 ***
## factor(decade)1980  0.118588   0.039589    3.00  0.00274 **
## factor(decade)1990  0.903440   0.038255   23.62  < 2e-16 ***
## factor(decade)2000  0.843949   0.038172   22.11  < 2e-16 ***
## factor(decade)2010  1.068025   0.037376   28.58  < 2e-16 ***
## P_NONWHT            0.020012   0.089224    0.22  0.82253
## a_hinc             -0.129094   0.008109  -15.92  < 2e-16 ***
## P_NONWHT:a_hinc     0.149223   0.018967    7.87 3.61e-15 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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 ( thetawartość 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:

fixed <- fixef(nbglmm)
confnitfixed <- confint(nbglmm, parm = "beta_", method = "Wald") # Beware: The Wald method is less accurate but much, much faster.

# The exponentiated coefficients are also known as Incidence Rate Ratios (IRR)
IRR <- exp(cbind(fixed, confintfixed)
IRR
##                         fixed      2.5 %     97.5 %
## (Intercept)        0.06094028 0.05458271 0.06803835
## factor(decade)1980 1.12590641 1.04184825 1.21674652
## factor(decade)1990 2.46807856 2.28979339 2.66024515
## factor(decade)2000 2.32553168 2.15789585 2.50619029
## factor(decade)2010 2.90962703 2.70410073 3.13077444
## P_NONWHT           1.02021383 0.85653208 1.21517487
## a_hinc             0.87889172 0.86503341 0.89297205
## P_NONWHT:a_hinc    1.16093170 1.11856742 1.20490048

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. glmmADMBmodelu pozwala na oszacowanie parametru inflacji stała zerowy. Alternatywnym podejściem byłoby użycie zeroinflfunkcji w psclpakiecie, 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_VACsą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).