Como encaixo um modelo multinível para resultados de poisson super dispersos?

32

Quero ajustar um GLMM multinível com uma distribuição Poisson (com excesso de dispersão) usando R. No momento, estou usando o lme4, mas notei que recentemente a quasipoissonfamília foi removida.

Já vi em outro lugar que você pode modelar a super dispersão aditiva para distribuições binomiais adicionando uma interceptação aleatória com um nível por observação. Isso também se aplica à distribuição de poisson?

Tem algum jeito melhor de fazer isso? Existem outros pacotes que você recomendaria?

George Michaelides
fonte

Respostas:

22

Você pode ajustar o GLMM multinível com uma distribuição Poisson (com excesso de dispersão) usando R de várias maneiras. Alguns Rpacotes são: lme4, MCMCglmm, arm, etc. Uma referência boa para ver é Gelman e Hill (2007)

Vou dar um exemplo de fazer isso usando rjagspacote no R. É uma interface entre Re JAGS(como OpenBUGSou WinBUGS).

log θ i j = β 0 + β 1 t r e um t m e n t i + δ i j δ i j ~ N ( 0 , σ 2 ϵ ) i = 1 ... I ,

nEujPoEusson(θEuj)
registroθEuj=β0 0+β1 TreumatmentEu+δEuj
δEujN(0 0,σϵ2)
Eu=1...Eu,j=1...J
TreumatmentEu=0 0 ou 1,...,J-1 se o Euth observação pertence ao grupo de tratamento 1ou 2,...,J

δEujrate 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)
}

Aqui está o Rcódigo para implementar usá-lo (diga que é nomeado 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])

Você pode brincar com as partes posteriores dos seus parâmetros e introduzir mais parâmetros para tornar a modelagem mais precisa ( gostamos de pensar isso ). Basicamente, você entendeu a idéia.

Para mais detalhes sobre o uso rjagse JAGS, consulte a página de John Myles White

suncoolsu
fonte
Obrigado!! Só recentemente comecei a analisar a análise bayesiana e ainda acho um pouco difícil de entender. Eu acho que essa é uma oportunidade para aprender um pouco mais sobre isso.
George Michaelides
1
Por que não dispersão gama?
Patrick McCann
2
@ Patrick você definitivamente pode fazer isso. Mas como estou registrando a média, prefiro o efeito disp normal. Uma distribuição normal de log é outra maneira de modelar distribuições semelhantes à distribuição gama. HTH.
suncoolsu
20

Não é necessário deixar o pacote lme4 para dar conta de sobredispersão; inclua apenas um efeito aleatório para o número de observação. As soluções BUGS / JAGS mencionadas provavelmente são um exagero para você e, se não forem, você deverá ter os resultados do lme4 fáceis de ajustar para comparação.

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)

Isso é discutido aqui: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 informal e academicamente por Elston et al. (2001) .

Patrick McCann
fonte
E se um modelo consistir em duas variáveis ​​nominais, uma variável contínua (todas como efeitos fixos) e uma variável de agrupamento (efeito aleatório) com interações de terceira ordem e, além disso, o número de sujeitos medidos é igual ao número de observações ou registros no conjunto de dados? Como devo cobrir isso no modelo?
Ladislav Naďo 28/01
7

Eu acho que o pacote glmmADMB é exatamente o que você está procurando.

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

Mas, do ponto de vista bayesiano, você pode usar o pacote MCMCglmm ou o software BUGS / JAGS , eles são muito flexíveis e você pode ajustar esse tipo de modelo. (e a sintaxe é próxima da R)

EDIT graças a @randel

Se você deseja instalar os pacotes glmmADMBe R2admb, é melhor fazer:

install.packages("glmmADMB", repos="http://glmmadmb.r-forge.r-project.org/repos"‌​)   
install.packages("R2admb")
dickoa
fonte
Acredito que atualmente o pacote deve ser instalado via install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")plus install.packages('R2admb').
Randel
5

Boas sugestões até agora. Aqui está mais um. Você pode ajustar um modelo de regressão binomial negativa hierárquica usando a rhierNegbinRwfunção do bayesmpacote.

Ari B. Friedman
fonte