Как мне соответствовать многоуровневой модели для перераспределенных результатов по пуассону?

32

Я хочу установить многоуровневый GLMM с распределением Пуассона (с избыточной дисперсией), используя R. В настоящее время я использую lme4, но я заметил, что недавно quasipoissonсемейство было удалено.

В другом месте я видел, что вы можете моделировать аддитивную избыточную дисперсию для биномиальных распределений, добавляя случайный перехват с одним уровнем на наблюдение. Относится ли это и к распределению Пуассона?

Есть ли лучший способ сделать это? Есть ли другие пакеты, которые вы бы порекомендовали?

Джордж Михаэлидис
источник

Ответы:

22

Вы можете установить многоуровневый 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яJ~пояssоN(θяJ)
журналθяJзнак равноβ0+β1 TреaTмеNTя+δяJ
δяJ~N(0,σε2)
язнак равно1...я,Jзнак равно1...J
TреaTмеNTязнак равно0 или 1,...,J-1 если яTчас наблюдение относится к группе лечения 1, или, 2,...,J

Часть в приведенном выше коде моделирует избыточную дисперсию. Но никто не мешает вам моделировать корреляцию между индивидами (вы не верите, что индивиды действительно независимы) и внутри индивидов (повторные измерения). Кроме того, параметр скорости может быть масштабирован некоторой другой константой, как в . Пожалуйста, смотрите Gelman and Hill (2007) для получения дополнительной ссылки. Вот код для простой модели:δяJrate modelsJAGS

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, пожалуйста, см . Страницу Джона Майлса Уайта

suncoolsu
источник
Благодарность!! Я только недавно начал изучать байесовский анализ, и мне все еще трудно его понять. Я предполагаю, что это возможность узнать немного больше об этом.
Джордж Михаэлидес
1
Почему не гамма-дисперсия?
Патрик Макканн
2
@ Патрик, ты определенно можешь это сделать. Но так как я беру среднее значение, я предпочитаю нормальный эффект дисп. Логарифмическое нормальное распределение является еще одним способом моделирования распределений, подобных гамма-распределению. НТН.
Suncoolsu
20

Нет необходимости оставлять пакет lme4 для учета чрезмерной дисперсии; просто включите случайный эффект для числа наблюдения. Упомянутые решения BUGS / JAGS для вас, вероятно, излишни, а если нет, то для сравнения вам нужно легко подобрать результаты lme4.

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)

Это обсуждается здесь: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 неофициально и академически Elston et al. (2001) .

Патрик Макканн
источник
Что, если модель состоит из двух номинальных переменных, одной непрерывной переменной (все как фиксированные эффекты) и одной группирующей переменной (случайный эффект) с взаимодействиями 3-го порядка, и, кроме того, количество измеряемых объектов равно количеству наблюдений или записей в набор данных? Как я должен покрыть это в модели?
Ладислав Наго
7

Я думаю, что пакет glmmADMB - это именно то, что вы ищете.

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

Но с байесовской точки зрения вы можете использовать пакет MCMCglmm или программное обеспечение BUGS / JAGS , они очень гибкие, и вы можете соответствовать этой модели. (и синтаксис близок к R)

РЕДАКТИРОВАТЬ благодаря @randel

Если вы хотите установить glmmADMBи R2admbпакеты, лучше сделать:

install.packages("glmmADMB", repos="http://glmmadmb.r-forge.r-project.org/repos"‌​)   
install.packages("R2admb")
dickoa
источник
Я считаю, что в настоящее время пакет должен быть установлен через install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")плюс install.packages('R2admb').
Рандел
5

Хорошие предложения пока. Вот еще один. Вы можете установить иерархическую модель отрицательной биномиальной регрессии, используя rhierNegbinRwфункцию bayesmпакета.

Ари Б. Фридман
источник