¿Cómo ajusto un modelo multinivel para resultados de Poisson sobredispersos?

Quiero ajustar un GLMM multinivel con una distribución de Poisson (con sobredispersión) usando R. En este momento estoy usando lme4 pero noté que recientemente la quasipoissonfamilia fue eliminada.

He visto en otras partes que puede modelar la sobredispersión aditiva para distribuciones binomiales agregando una intercepción aleatoria con un nivel por observación. ¿Esto se aplica también a la distribución de Poisson?

Hay una mejor manera de hacerlo? ¿Hay otros paquetes que recomendarías?

Puede ajustar GLMM multinivel con una distribución de Poisson (con sobredispersión) usando R de múltiples maneras. Unos Rpaquetes son: lme4, MCMCglmm, arm, etc. Una referencia buena para ver que es Gelman y Hill (2007)

Daré un ejemplo de esto usando el rjagspaquete in R. Es una interfaz entre Ry JAGS(me gusta OpenBUGSo WinBUGS).

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

$$n_{ij} \sim \mathrm{Poisson}(\theta_{ij})$$ $$\log \theta_{ij} = \beta_0 + \beta_1 \mbox{ } \mathtt{Treatment}_{i} + \delta_{ij}$$ $$\delta_{ij} \sim N(0, \sigma^2_{\epsilon})$$ $$i=1 \ldots I, \quad j = 1\ldots J$$ $\mathtt{Treatment}_i = 0 \mbox{ or } 1, \ldots, J-1 \mbox{ if the } i^{th} \mbox{ observation belongs to treatment group } 1 \mbox{, or, } 2, \ldots, J$

La parte en el código anterior modela la sobredispersión. Pero no hay nadie que le impida modelar la correlación entre individuos (no cree que los individuos sean realmente independientes) y dentro de los individuos (medidas repetidas). Además, el parámetro de velocidad puede ser escalado por alguna otra constante como en . Consulte Gelman y Hill (2007) para obtener más referencias. Aquí está el código para el modelo simple:$\delta_{ij}$rate 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)
}

Aquí está el Rcódigo para implementar el uso que (dicen que se denomina: 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])

Puede jugar con las partes posteriores de sus parámetros y puede introducir más parámetros para que su modelado sea más preciso ( nos gusta pensar esto ). Básicamente, te haces una idea.

Para obtener más detalles sobre el uso rjagsy JAGS, consulte la página de John Myles White

No es necesario dejar el paquete lme4 para tener en cuenta la sobredispersión; solo incluya un efecto aleatorio para el número de observación. Las soluciones BUGS / JAGS mencionadas son probablemente excesivas para usted, y si no lo son, debería tener los resultados lme4 fáciles de ajustar para comparar.

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)

Esto se discute aquí: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 de manera informal y académica por Elston et al. (2001) .

Creo que el paquete glmmADMB es exactamente lo que estás buscando.

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

Pero desde el punto de vista bayesiano, puede usar el paquete MCMCglmm o el software BUGS / JAGS , son muy flexibles y puede adaptarse a este tipo de modelo. (y la sintaxis está cerca de la R)

EDITAR gracias a @randel

Si desea instalar los paquetes glmmADMBy R2admb, es mejor hacerlo:

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

Buenas sugerencias hasta ahora. Aquí hay uno más. Puede ajustar un modelo jerárquico de regresión binomial negativa utilizando la rhierNegbinRwfunción del bayesmpaquete.