Możesz dopasować wielopoziomowy GLMM z rozkładem Poissona (z nadmierną dyspersją), używając R na wiele sposobów. Kilka R
pakiety są: lme4
, MCMCglmm
, arm
, itd. Dobrym odniesienia, aby zobaczyć to Gelman i Hill (2007)
Podam przykład wykonania tego przy użyciu rjags
pakietu w R
. Jest to interfejs pomiędzy R
i JAGS
(jak OpenBUGS
lub WinBUGS
).
log θ I J = β 0 + β 1 T R E t m e n t i + δ I J δ I J ~ N ( 0 , σ 2 ϵ ) i = 1 … I ,
nij∼Poisson(θij)
logθij=β0+β1 Treatmenti+δij
δij∼N(0,σ2ϵ)
T r e a t m e n t i = 0 lub 1 , … , J - 1, jeżeli obserwacja i t h należy do grupy leczenia 1 lub, 2 , … , Ji=1…I,j=1…J
Treatmenti=0 or 1,…,J−1 if the ith observation belongs to treatment group 1, or, 2,…,J
δijrate models
JAGS
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)
}
Oto R
kod do wdrożenia go używać (mówią, że nazywa się: 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])
Możesz bawić się z tylnymi parametrami i możesz wprowadzić więcej parametrów, aby uczynić modelowanie bardziej precyzyjnym ( lubimy to uważać ). Zasadniczo masz pomysł.
Aby uzyskać więcej informacji na temat używania rjags
i JAGS
, zobacz stronę Johna Mylesa White'a
Nie trzeba opuszczać pakietu lme4, aby uwzględnić nadmierną dyspersję; wystarczy dołączyć losowy efekt dla liczby obserwacji. Wymienione rozwiązania BŁĘDÓW / JAGSÓW są prawdopodobnie dla ciebie przesadą, a jeśli nie są, powinieneś mieć łatwe do dopasowania wyniki lme4 do porównania.
Jest to omówione tutaj: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 nieformalnie i naukowo przez Elston i in. (2001) .
źródło
Myślę, że pakiet glmmADMB jest dokładnie tym, czego szukasz.
Ale z punktu widzenia bayesowskiego można użyć pakietu MCMCglmm lub oprogramowania BUGS / JAGS , są one bardzo elastyczne i można dopasować do tego rodzaju modelu. (a składnia jest zbliżona do R)
EDYCJA dzięki @randel
Jeśli chcesz zainstalować pakiety
glmmADMB
iR2admb
, lepiej zrobić:źródło
install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")
plusainstall.packages('R2admb')
.Good suggestions so far. Here's one more. You can fit a hierarchical negative binomial regression model using the
rhierNegbinRw
function of thebayesm
package.źródło