Kesulitan menemukan model yang bagus cocok untuk data hitungan dengan efek campuran - ZINB atau yang lainnya?

Saya memiliki satu set data yang sangat kecil tentang kelimpahan lebah soliter yang saya kesulitan analisis. Ini menghitung data, dan hampir semua penghitungan berada dalam satu perawatan dengan sebagian besar nol dalam perawatan lainnya. Ada juga beberapa nilai yang sangat tinggi (masing-masing di dua dari enam situs), sehingga distribusi jumlah memiliki ekor yang sangat panjang. Saya bekerja di R. Saya telah menggunakan dua paket berbeda: lme4 dan glmmADMB.

Model campuran Poisson tidak cocok: model sangat overdispersi ketika efek acak tidak dipasang (model glm), dan underdispersi ketika efek acak dipasang (model glmer). Saya tidak mengerti mengapa ini terjadi. Desain eksperimental membutuhkan efek acak bersarang jadi saya harus memasukkannya. Distribusi kesalahan lognormal Poisson tidak meningkatkan kecocokan. Saya mencoba distribusi kesalahan binomial negatif menggunakan glmer.nb dan tidak bisa menyesuaikannya - batas iterasi tercapai, bahkan ketika mengubah toleransi menggunakan glmerControl (tolPwrss = 1e-3).

Karena banyak dari nol akan disebabkan oleh fakta bahwa saya tidak melihat lebah (mereka sering merupakan benda hitam kecil), saya selanjutnya mencoba model nol-meningkat. ZIP tidak pas. ZINB adalah model yang paling cocok sejauh ini, tapi saya masih tidak terlalu senang dengan model yang cocok. Saya bingung harus mencoba apa selanjutnya. Saya memang mencoba model rintangan tetapi tidak dapat mencocokkan distribusi yang terpotong dengan hasil yang tidak nol– saya pikir karena begitu banyak dari nol yang ada dalam perlakuan kontrol (pesan kesalahannya adalah “Kesalahan dalam model.frame.default (rumus = s.bee ~ tmt + lu +: panjang variabel berbeda (ditemukan untuk 'pengobatan') ”).

Selain itu, saya berpikir bahwa interaksi yang saya sertakan adalah melakukan sesuatu yang aneh pada data saya karena koefisiennya sangat kecil - walaupun model yang berisi interaksi itu paling baik ketika saya membandingkan model yang menggunakan AICctab dalam paket bbmle.

Saya menyertakan beberapa skrip R yang akan mereproduksi set data saya. Variabelnya adalah sebagai berikut:

d = Tanggal Julian, df = Tanggal Julian (sebagai faktor), d.sq = df kuadrat (jumlah lebah meningkat kemudian turun sepanjang musim panas), st = situs, s.bee = jumlah lebah, tmt = pengobatan, lu = jenis penggunaan lahan, hab = persentase habitat semi alami di lanskap sekitarnya, ba = batas bidang bundar.

Setiap saran bagaimana saya bisa mendapatkan model yang baik (distribusi kesalahan alternatif, berbagai jenis model dll) akan sangat diterima!

Terima kasih.

d <- c(80,  80,  121, 121, 180, 180, 86,  86,  116, 116, 144, 144, 74,  74, 143, 143, 163, 163, 71, 71,106, 106, 135, 135, 162, 162, 185, 185, 83,  83,  111, 111, 133, 133, 175, 175, 85,  85,  112, 112,137, 137, 168, 168, 186, 186, 64,  64,  95,  95,  127, 127, 156, 156, 175, 175, 91,  91, 119, 119,120, 120, 148, 148, 56, 56)
df <- as.factor(d)
d.sq <- d^2
st <- factor(rep(c("A", "B", "C", "D", "E", "F"), c(6,12,18,10,14,6)))
s.bee <- c(1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,4,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,3,0,0,0,0,5,0,0,2,0,50,0,10,0,4,0,47,3)
tmt <- factor(c("AF","C","C","AF","AF","C","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C",
"C","AF","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C","AF","C","AF","C",
"C","AF","AF","C","AF","C","AF","C","AF","C","C","AF","C","AF","C","AF","AF","C","AF","C",
"AF","C","AF","C","AF","C"))
lu <- factor(rep(c("p","a","p","a","p"), c(6,12,28,14,6)))
hab <- rep(c(13,14,13,14,3,4,3,4,3,4,3,4,3,4,15,35,37,35,37,35,37,35,37,0,2,1,2,1,2,1), 
        c(1,2,2,1,1,1,1,2,2,1,1,1,1,1,18,1,1,1,2,2,1,1,1,14,1,1,1,1,1,1))
ba <-  c(480,6520,6520,480,480,6520,855,1603,855,1603,1603,855,855,12526,855,5100,855,5100,2670,7679,7679,2670,
2670,7679,2670,7679,7679,2670,2670,7679,2670,7679,2670,7679,2670,7679,1595,3000,1595,3000,3000,1595,1595,3000,1595
,3000,4860,5460,4860,5460,5460,4860,5460,4860,5460,4860,4840,5460,4840,5460,3000,1410,3000,1410,3000,1410)
data <- data.frame(st,df,d.sq,tmt,lu,hab,ba,s.bee)
with(data, table(s.bee, tmt) )

# below is a much abbreviated summary of attempted models:

library(MASS)
library(lme4)
library(glmmADMB)
library(coefplot2)

###
### POISSON MIXED MODEL

    m1 <- glmer(s.bee ~ tmt + lu + hab + (1|st/df), family=poisson)
    summary(m1)

    resdev<-sum(resid(m1)^2)
    mdf<-length(fixef(m1)) 
    rdf<-nrow(data)-mdf
    resdev/rdf
# 0.2439303
# underdispersed. ???



###
### NEGATIVE BINOMIAL MIXED MODEL

    m2 <- glmer.nb(s.bee ~ tmt + lu + hab + d.sq + (1|st/df))
# iteration limit reached. Can't make a model work.



###
### ZERO-INFLATED POISSON MIXED MODEL

    fit_zipoiss <- glmmadmb(s.bee~tmt + lu + hab + ba + d.sq +
                    tmt:lu +
                    (1|st/df), data=data,
                    zeroInflation=TRUE,
                    family="poisson")
# has to have lots of variables to fit
# anyway Poisson is not a good fit



###
### ZERO-INFLATED NEGATIVE BINOMIAL MIXED MODELS

## BEST FITTING MODEL SO FAR:

    fit_zinb <- glmmadmb(s.bee~tmt + lu + hab +
                    tmt:lu +
                    (1|st/df),data=data,
                    zeroInflation=TRUE,
                    family="nbinom")
    summary(fit_zinb)
# coefficients are tiny, something odd going on with the interaction term
# but this was best model in AICctab comparison

# model check plots
    qqnorm(resid(fit_zinb))
    qqline(resid(fit_zinb))

    coefplot2(fit_zinb)

    resid_zinb <- resid(fit_zinb , type = "pearson")
    hist(resid_zinb)

    fitted_zinb <- fitted (fit_zinb)
    plot(resid_zinb ~ fitted_zinb)


## MODEL WITHOUT INTERACTION TERM - the coefficients are more realistic:

    fit_zinb2 <- glmmadmb(s.bee~tmt + lu + hab +
                    (1|st/df),data=data,
                    zeroInflation=TRUE,
                    family="nbinom")

# model check plots
    qqnorm(resid(fit_zinb2))
    qqline(resid(fit_zinb2))

    coefplot2(fit_zinb2)

    resid_zinb2 <- resid(fit_zinb2 , type = "pearson")
    hist(resid_zinb2)

    fitted_zinb2 <- fitted (fit_zinb2)
    plot(resid_zinb2 ~ fitted_zinb2)



# ZINB models are best so far
# but I'm not happy with the model check plots

Posting ini sudah empat tahun, tetapi saya ingin mengikuti apa yang dikatakan fsociety dalam komentar. Diagnosis residu dalam GLMM tidak mudah, karena plot residual standar dapat menunjukkan non-normalitas, heteroskedastisitas, dll., Bahkan jika modelnya ditentukan dengan benar. Ada paket R DHARMa, yang secara khusus cocok untuk mendiagnosis jenis model ini.

Paket ini didasarkan pada pendekatan simulasi untuk menghasilkan residu skala dari model campuran linier umum dan menghasilkan plot diagnostik yang mudah ditafsirkan. Berikut adalah contoh kecil dengan data dari pos asli dan model pas pertama (m1):

library(DHARMa)
sim_residuals <- simulateResiduals(m1, 1000)
plotSimulatedResiduals(sim_residuals)

Plot di sebelah kiri menunjukkan plot QQ dari residu yang diskalakan untuk mendeteksi penyimpangan dari distribusi yang diharapkan, dan plot di sebelah kanan mewakili residu vs nilai yang diprediksi saat melakukan regresi kuantil untuk mendeteksi penyimpangan dari keseragaman (garis merah harus horizontal dan pada 0,25 , 0,50 dan 0,75).

Selain itu, paket ini juga memiliki fungsi spesifik untuk pengujian dispersi berlebih / kurang dan inflasi nol, antara lain:

testOverdispersionParametric(m1)

Chisq test for overdispersion in GLMMs

data:  poisson
dispersion = 0.18926, pearSS = 11.35600, rdf = 60.00000, p-value = 1
alternative hypothesis: true dispersion greater 1

testZeroInflation(sim_residuals)

DHARMa zero-inflation test via comparison to expected zeros with 
simulation under H0 = fitted model


data:  sim_residuals
ratioObsExp = 0.98894, p-value = 0.502
alternative hypothesis: more