Saya punya beberapa pertanyaan tentang spesifikasi dan interpretasi GLMM. 3 pertanyaan pasti statistik dan 2 lebih spesifik tentang R. Saya posting di sini karena pada akhirnya saya pikir masalahnya adalah interpretasi hasil GLMM.

Saat ini saya mencoba untuk menyesuaikan GLMM. Saya menggunakan data sensus AS dari Longitudinal Tract Database . Pengamatan saya adalah trus sensus. Variabel dependen saya adalah jumlah unit rumah kosong dan saya tertarik pada hubungan antara lowongan dan variabel sosial-ekonomi. Contoh di sini sederhana, hanya menggunakan dua efek tetap: persen populasi bukan kulit putih (ras) dan pendapatan rumah tangga rata-rata (kelas), ditambah interaksi mereka. Saya ingin memasukkan dua efek acak bersarang: traktat dalam dekade dan dekade, yaitu (dekade / traktat). Saya mempertimbangkan ini secara acak dalam upaya untuk mengendalikan autokorelasi spasial (yaitu antara traktat) dan temporal (yaitu antara dekade). Namun, saya juga tertarik pada dekade sebagai efek tetap, jadi saya memasukkannya juga sebagai faktor tetap.

Karena variabel independen saya adalah variabel hitung bilangan bulat non-negatif, saya sudah mencoba menyesuaikan poisson dan GLMM binomial negatif. Saya menggunakan log total unit perumahan sebagai offset. Ini berarti koefisien ditafsirkan sebagai efek pada tingkat kekosongan, bukan jumlah total rumah kosong.

Saat ini saya memiliki hasil untuk Poisson dan GLMM binomial negatif yang diperkirakan menggunakan glmer dan glmer.nb dari lme4 . Interpretasi koefisien masuk akal bagi saya berdasarkan pengetahuan saya tentang data dan bidang studi.

Jika Anda ingin data dan skrip ada di Github saya . Script mencakup lebih banyak penyelidikan deskriptif yang saya lakukan sebelum membangun model.

Inilah hasil saya:

Model Poisson

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)

Model binomial negatif

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

Tes Poisson 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

Tes DHARMa binomial negatif

    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

Plot DHARMa

Poisson

Binomial negatif

Pertanyaan statistik

Karena saya masih mencari tahu GLMMs saya merasa tidak aman tentang spesifikasi dan interpretasi. Saya punya beberapa pertanyaan:

1. Tampaknya data saya tidak mendukung menggunakan model Poisson dan oleh karena itu saya lebih baik dengan binomial negatif. Namun, saya secara konsisten mendapatkan peringatan bahwa model binomial negatif saya mencapai batas iterasi mereka, bahkan ketika saya meningkatkan batas maksimum. "Dalam theta.ml (Y, mu, bobot = objek @ resp $ bobot, batas = batas,: batas iterasi tercapai." Ini terjadi menggunakan beberapa spesifikasi yang berbeda (yaitu model minimum dan maksimal untuk efek tetap dan acak). Saya juga telah mencoba menghilangkan outlier dalam ketergantungan saya (kotor, saya tahu!), Karena 1% nilai teratas adalah outlier yang sangat banyak (kisaran 99% bawah 0-1012, atas 1% dari 1013-5213). t tidak memiliki efek pada iterasi dan sangat sedikit berpengaruh pada koefisien juga. Saya tidak memasukkan detail itu di sini. Koefisien antara Poisson dan binomial negatif juga sangat mirip. Apakah kurangnya konvergensi ini merupakan masalah? Apakah model binomial negatif cocok? Saya juga menjalankan model binomial negatif menggunakanAllFit dan tidak semua pengoptimal melempar peringatan ini (bobyqa, Nelder Mead, dan nlminbw tidak).

2. Varian untuk efek tetap dasawarsa saya secara konsisten sangat rendah atau 0. Saya mengerti ini bisa berarti modelnya terlalu bagus. Mengambil dekade dari efek tetap memang meningkatkan varians efek acak dekade ke 0,2620 dan tidak memiliki banyak efek pada koefisien efek tetap. Apakah ada yang salah dengan membiarkannya? Saya baik-baik saja menafsirkannya sebagai hanya tidak diperlukan untuk menjelaskan antara perbedaan pengamatan.

3. Apakah hasil ini menunjukkan saya harus mencoba model nol-inflasi? DHARMa tampaknya menyarankan inflasi nol mungkin bukan masalah. Jika Anda pikir saya harus tetap mencoba, lihat di bawah.

Pertanyaan R

4. Saya akan bersedia untuk mencoba model zero-inflated, tetapi saya tidak yakin paket mana yang menerapkan efek acak untuk Poisson zero-inflated dan GLMM binomial negatif. Saya akan menggunakan glmmADMB untuk membandingkan AIC dengan model zero-inflated, tetapi terbatas pada efek acak tunggal sehingga tidak berfungsi untuk model ini. Saya bisa mencoba MCMCglmm, tetapi saya tidak tahu statistik Bayesian sehingga juga tidak menarik. Ada opsi lain?

5. Dapatkah saya menampilkan koefisien eksponensial dalam ringkasan (model), atau apakah saya harus melakukannya di luar ringkasan seperti yang saya lakukan di sini?

Saya percaya ada beberapa masalah penting yang harus diatasi dengan estimasi Anda.

Dari apa yang saya kumpulkan dengan memeriksa data Anda, unit Anda tidak dikelompokkan secara geografis, yaitu trus sensus di dalam kabupaten. Jadi, menggunakan traktat sebagai faktor pengelompokan tidak sesuai untuk menangkap heterogenitas spasial karena ini berarti Anda memiliki jumlah individu yang sama dengan kelompok (atau dengan kata lain, semua grup Anda hanya memiliki satu pengamatan masing-masing). Menggunakan strategi pemodelan bertingkat memungkinkan kita untuk memperkirakan varians tingkat individu, sambil mengendalikan varians antar-kelompok. Karena grup Anda hanya memiliki satu individu, varians antar-grup Anda sama dengan varians level individual Anda, sehingga mengalahkan tujuan pendekatan multilevel.

Di sisi lain, faktor pengelompokan dapat mewakili pengukuran berulang dari waktu ke waktu. Misalnya, dalam kasus studi longitudinal, nilai "matematika" seseorang dapat diulang setiap tahun, sehingga kita akan memiliki nilai tahunan untuk setiap siswa selama n tahun (dalam hal ini, faktor pengelompokan adalah siswa seperti yang kita miliki n jumlah pengamatan "bersarang" di dalam siswa). Dalam kasus Anda, Anda telah mengulangi langkah-langkah setiap saluran sensus decade. Dengan demikian, Anda bisa menggunakan TRTID10variabel Anda sebagai faktor pengelompokan untuk menangkap "antara varian satu dekade". Hal ini menyebabkan 3142 pengamatan bersarang di 635 traktat, dengan sekitar 4 dan 5 observasi per trus sensus.

Seperti disebutkan dalam komentar, menggunakan decadesebagai faktor pengelompokan tidak terlalu tepat, karena Anda hanya memiliki sekitar 5 dekade untuk setiap saluran sensus, dan efeknya dapat ditangkap dengan lebih baik dengan diperkenalkan decadesebagai kovariat.

Kedua, untuk menentukan apakah data Anda harus dimodelkan menggunakan poisson atau model binomial negatif (atau pendekatan inflasi nol). Pertimbangkan jumlah penayangan berlebih di data Anda. Karakteristik mendasar dari distribusi Poisson adalah equidispersion, yang berarti bahwa rata-rata sama dengan varian distribusi. Melihat data Anda, cukup jelas bahwa ada banyak penyebaran berlebihan. Variansi jauh lebih besar daripada sarana.

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 

Meskipun demikian, untuk menentukan apakah binomial negatif lebih tepat secara statistik, metode standar adalah melakukan uji rasio kemungkinan antara Poisson dan model binomial negatif, yang menunjukkan bahwa negbin lebih sesuai.

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

Setelah menetapkan ini, tes selanjutnya dapat mempertimbangkan apakah pendekatan multilevel (model campuran) dijamin menggunakan pendekatan yang sama, yang menunjukkan bahwa versi multilevel memberikan kesesuaian yang lebih baik. (Tes serupa dapat digunakan untuk membandingkan fit glmer dengan asumsi distribusi poisson ke objek glmer.nb, asalkan modelnya sama.)

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

Mengenai perkiraan model poisson dan nb mereka sebenarnya diharapkan sangat mirip satu sama lain, dengan perbedaan utama adalah kesalahan standar, yaitu jika ada penyebaran berlebihan, model poisson cenderung memberikan kesalahan standar yang bias. Mengambil data Anda sebagai contoh:

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

Perhatikan bagaimana estimasi koefisien semuanya sangat mirip, perbedaan utamanya hanya signifikansi dari salah satu kovariat Anda, serta perbedaan dalam varian efek acak, yang menunjukkan bahwa varian level-unit yang ditangkap oleh parameter overdispersion di nb model ( thetanilai dalam objek glmer.nb) menangkap beberapa varians antara saluran yang ditangkap oleh efek acak.

Mengenai koefisien eksponensial (dan interval kepercayaan terkait), Anda dapat menggunakan yang berikut ini:

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

Pikiran terakhir, tentang nol inflasi. Tidak ada implementasi multilevel (setidaknya yang saya ketahui) dari model poisson atau negbin yang meningkat nol yang memungkinkan Anda menentukan persamaan untuk komponen campuran nol yang meningkat. yang glmmADMBmodel yang memungkinkan Anda memperkirakan parameter inflasi konstan nol. Pendekatan alternatif akan menggunakan zeroinflfungsi dalam psclpaket, meskipun ini tidak mendukung model bertingkat. Dengan demikian, Anda dapat membandingkan kecocokan binomial negatif level tunggal, dengan binomial negatif tingkat tunggal nol. Kemungkinannya adalah jika inflasi nol tidak signifikan untuk model tingkat tunggal, ada kemungkinan bahwa itu tidak akan signifikan untuk spesifikasi multilevel.

Tambahan

Jika Anda khawatir tentang autokorelasi spasial, Anda bisa mengendalikan ini menggunakan beberapa bentuk regresi geografis tertimbang (meskipun saya percaya ini menggunakan data titik, bukan area). Atau, Anda dapat mengelompokkan saluran sensus Anda dengan faktor pengelompokan tambahan (negara bagian, kabupaten) dan memasukkan ini sebagai efek acak. Terakhir, dan saya tidak yakin apakah ini sepenuhnya layak, dimungkinkan untuk menggabungkan ketergantungan spasial dengan menggunakan, misalnya, jumlah rata-rata R_VACtetangga urutan pertama sebagai kovariat. Dalam kasus apa pun, sebelum pendekatan semacam itu, akan masuk akal untuk menentukan apakah autokorelasi spasial memang ada (menggunakan Global Moran's I, tes LISA, dan pendekatan serupa).