Bagaimana saya menyesuaikan model multilevel untuk hasil poisson yang terlalu tersebar?

32

Saya ingin mencocokkan GLMM bertingkat dengan distribusi Poisson (dengan dispersi berlebih) menggunakan R. Saat ini saya menggunakan lme4 tetapi saya perhatikan bahwa baru-baru ini quasipoissonkeluarga telah dihapus.

Saya telah melihat di tempat lain bahwa Anda dapat memodelkan over-dispersi aditif untuk distribusi binomial dengan menambahkan intersepsi acak dengan satu level per observasi. Apakah ini berlaku untuk distribusi poisson juga?

Apakah ada cara yang lebih baik untuk melakukannya? Apakah ada paket lain yang akan Anda rekomendasikan?

George Michaelides
sumber

Jawaban:

22

Anda dapat menyesuaikan multilevel GLMM dengan distribusi Poisson (dengan dispersi berlebih) menggunakan R dalam berbagai cara. Beberapa Rpaket adalah: lme4, MCMCglmm, arm, dll Sebuah referensi yang baik untuk melihat adalah Gelman dan Hill (2007)

Saya akan memberikan contoh untuk melakukan hal ini menggunakan rjagspaket di R. Ini adalah antarmuka antara Rdan JAGS(suka OpenBUGSatau WinBUGS).

nijPoisson(θij)
δ i j ~ N ( 0 , σ 2
logθij=β0+β1 Treatmenti+δij
i=1...Saya,
δijN(0,σϵ2)
T r e a t m e n t i = 0  atau  1 , , J - 1  jika  i t h
i=1I,j=1J
Treatmenti=0 or 1,,J1 if the ith observation belongs to treatment group 1, or, 2,,J

Bagian δij dalam kode di atas model overdispersi. Tetapi tidak ada yang menghentikan Anda dari pemodelan korelasi antara individu (Anda tidak percaya bahwa individu benar-benar independen) dan di dalam individu (tindakan berulang). Juga, parameter laju dapat diskalakan oleh beberapa konstanta lain seperti pada rate models. Silakan lihat Gelman dan Hill (2007) untuk referensi lebih lanjut. Berikut adalah JAGSkode untuk model sederhana:

data{
        for (i in 1:I){         
            ncount[i,1] <- obsTrt1[i]
            ncount[i,2] <- obsTrt2[i]
                ## notice I have only 2 treatments and I individuals 
    }                               
}

model{
    for (i in 1:I){ 
        nCount[i, 1] ~ dpois( means[i, 1] )
        nCount[i, 2] ~ dpois( means[i, 2] )

        log( means[i, 1] ) <- mu + b * trt1[i] + disp[i, 1]
        log( means[i, 2] ) <- mu + b * trt2[i] + disp[i, 2]

        disp[i, 1] ~ dnorm( 0, tau)
        disp[i, 2] ~ dnorm( 0, tau)

    }

    mu  ~ dnorm( 0, 0.001)
    b   ~ dnorm(0, 0.001)
    tau ~ dgamma( 0.001, 0.001)
}

Berikut adalah Rkode untuk menerapkan penggunaan itu (mengatakan itu bernama: overdisp.bug)

dataFixedEffect <- list("I"       = 10,
                        "obsTrt1" = obsTrt1 , #vector of n_i1
                        "obsTrt2" = obsTrt2,  #vector of n_i2
                        "trt1"    = trt1,     #vector of 0
                        "trt2"    = trt2,     #vector of 1
                       )

initFixedEffect <- list(mu = 0.0 , b = 0.0, tau = 0.01)

simFixedEffect <- jags.model(file     = "overdisp.bug",
                             data     = dataFixedEffect,
                             inits    = initFixedEffect,
                             n.chains = 4,
                             n.adapt  = 1000)

sampleFixedEffect <- coda.samples(model          = simFixedEffect,
                                  variable.names = c("mu", "b", "means"),
                                  n.iter         = 1000)

meansTrt1 <- as.matrix(sampleFixedEffect[ , 2:11])
meansTrt2 <- as.matrix(sampleFixedEffect[ , 12:21])

Anda dapat bermain-main dengan eksterior parameter Anda dan Anda dapat memperkenalkan lebih banyak parameter untuk membuat pemodelan Anda lebih tepat ( kami suka memikirkannya ). Pada dasarnya, Anda mendapatkan idenya.

Untuk detail lebih lanjut tentang penggunaan rjagsdan JAGS, silakan lihat halaman John Myles White

suncoolsu
sumber
Terima kasih!! Saya baru saja mulai melihat ke dalam analisis bayesian dan saya masih merasa agak sulit untuk dipahami. Saya kira ini adalah kesempatan untuk belajar lebih banyak tentangnya.
George Michaelides
1
Mengapa tidak dispersi gamma?
Patrick McCann
2
@ Patrick Anda pasti bisa melakukannya. Tapi karena saya mengambil log berarti saya lebih suka efek disp biasa. Distribusi log normal adalah cara lain untuk memodelkan distribusi yang mirip dengan distribusi gamma. HTH.
suncoolsu
20

Tidak perlu meninggalkan paket lme4 untuk memperhitungkan penyebaran berlebihan; cukup sertakan efek acak untuk nomor observasi. Solusi BUGS / JAGS yang disebutkan mungkin berlebihan bagi Anda, dan jika tidak, Anda harus memiliki hasil lme4 yang mudah untuk dicocokkan sebagai perbandingan.

data$obs_effect<-1:nrow(data)
overdisp.fit<-lmer(y~1+obs_effect+x+(1|obs_effect)+(1+x|subject_id),data=data,family=poisson)

Ini dibahas di sini: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 secara informal dan akademis oleh Elston et al. (2001) .

Patrick McCann
sumber
Bagaimana jika suatu model terdiri dari dua variabel nominal, satu variabel kontinu (semua sebagai efek tetap) dan satu variabel pengelompokan (efek acak) dengan interaksi tingkat ketiga dan, terlebih lagi, jumlah subjek yang diukur sama dengan jumlah pengamatan atau catatan dalam dataset? Bagaimana saya harus membahas ini dalam model?
Ladislav Naďo
7

Saya pikir paket glmmADMB persis seperti yang Anda cari.

install.packages ("glmmADMB", repos = "http://r-forge.r-project.org")

Tetapi dalam sudut pandang bayesian Anda dapat menggunakan paket MCMCglmm atau perangkat lunak BUGS / JAGS , mereka sangat fleksibel dan Anda dapat menggunakan model seperti ini. (dan sintaksnya dekat dengan R)

EDIT terima kasih kepada @randel

Jika Anda ingin menginstal glmmADMBdan R2admbpaket, lebih baik dilakukan:

install.packages("glmmADMB", repos="http://glmmadmb.r-forge.r-project.org/repos"‌​)   
install.packages("R2admb")
dickoa
sumber
Saya percaya saat ini paket harus diinstal melalui install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")plus install.packages('R2admb').
Randel
5

Saran yang bagus sejauh ini. Ini satu lagi. Anda bisa memasukkan model regresi binomial negatif hierarkis menggunakan rhierNegbinRwfungsi bayesmpaket.

Ari B. Friedman
sumber