Anda dapat menyesuaikan multilevel GLMM dengan distribusi Poisson (dengan dispersi berlebih) menggunakan R dalam berbagai cara. Beberapa R
paket 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 rjags
paket di R
. Ini adalah antarmuka antara R
dan JAGS
(suka OpenBUGS
atau WinBUGS
).
nsaya j∼ P o i s s o n ( θsaya j)
δ i j ~ N ( 0 , σ 2logθsaya j= β0+ β1 T r e a t me n tsaya+ δsaya j
i=1...Saya,δsayaj∼ N( 0 , σ2ϵ)
T r e a t m e n t i = 0 atau 1 , … , J - 1 jika i t hi = 1 ... I,j = 1 ... J
T r e a t me n tsaya= 0 atau 1 , ... , J- 1 jika it h observasi termasuk dalam kelompok perlakuan 1 , atau, 2,…,J
Bagian δsaya j 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 JAGS
kode 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 R
kode 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 rjags
dan JAGS
, silakan lihat halaman John Myles White
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.
Ini dibahas di sini: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 secara informal dan akademis oleh Elston et al. (2001) .
sumber
Saya pikir paket glmmADMB persis seperti yang Anda cari.
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
glmmADMB
danR2admb
paket, lebih baik dilakukan:sumber
install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")
plusinstall.packages('R2admb')
.Saran yang bagus sejauh ini. Ini satu lagi. Anda bisa memasukkan model regresi binomial negatif hierarkis menggunakan
rhierNegbinRw
fungsibayesm
paket.sumber