Вы можете установить многоуровневый GLMM с распределением Пуассона (с избыточной дисперсией), используя R несколькими способами. Немногие R
пакеты: lme4
, MCMCglmm
, arm
и т.д. Хорошая ссылка , чтобы увидеть это Гельман и Hill (2007)
Я приведу пример этого с использованием rjags
пакета в R
. Это интерфейс между R
и JAGS
(как OpenBUGS
или WinBUGS
).
войти θ я J = β 0 + β 1 т т е т т е н т я + δ я J δ я J ~ N ( 0 , σ 2 ε ) я = 1 ... I ,
Nя ж~ Р о я ы ы ö п ( θя ж)
журналθя ж= β0+ β1 Т т е с т те п тя+ δя ж
δя ж~ N( 0 , σ2ε)
я = 1 ... я,j = 1 … J
Т г е с тте п тя= 0 или 1 , … , J- 1, если ят ч наблюдение относится к группе лечения 1 , или, 2,…,J
Часть в приведенном выше коде моделирует избыточную дисперсию. Но никто не мешает вам моделировать корреляцию между индивидами (вы не верите, что индивиды действительно независимы) и внутри индивидов (повторные измерения). Кроме того, параметр скорости может быть масштабирован некоторой другой константой, как в . Пожалуйста, смотрите Gelman and Hill (2007) для получения дополнительной ссылки. Вот код для простой модели:δя жrate 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)
}
Вот R
код для реализации использовать его (скажем , это называется: 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])
Вы можете поиграть с постерами ваших параметров, и вы можете ввести больше параметров, чтобы сделать моделирование более точным ( нам нравится так думать ). По сути, вы поняли идею.
Для получения более подробной информации об использовании rjags
и JAGS
, пожалуйста, см . Страницу Джона Майлса Уайта
Нет необходимости оставлять пакет lme4 для учета чрезмерной дисперсии; просто включите случайный эффект для числа наблюдения. Упомянутые решения BUGS / JAGS для вас, вероятно, излишни, а если нет, то для сравнения вам нужно легко подобрать результаты lme4.
Это обсуждается здесь: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 неофициально и академически Elston et al. (2001) .
источник
Я думаю, что пакет glmmADMB - это именно то, что вы ищете.
Но с байесовской точки зрения вы можете использовать пакет MCMCglmm или программное обеспечение BUGS / JAGS , они очень гибкие, и вы можете соответствовать этой модели. (и синтаксис близок к R)
РЕДАКТИРОВАТЬ благодаря @randel
Если вы хотите установить
glmmADMB
иR2admb
пакеты, лучше сделать:источник
install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")
плюсinstall.packages('R2admb')
.Хорошие предложения пока. Вот еще один. Вы можете установить иерархическую модель отрицательной биномиальной регрессии, используя
rhierNegbinRw
функциюbayesm
пакета.источник