Wie passe ich ein Mehrebenenmodell für übermäßig verteilte Poisson-Ergebnisse an?

Ich möchte ein mehrstufiges GLMM mit einer Poisson-Verteilung (mit Überdispersion) unter Verwendung von R ausrüsten. Im Moment verwende ich lme4, aber mir ist aufgefallen, dass die quasipoissonFamilie kürzlich entfernt wurde.

Ich habe an anderer Stelle gesehen, dass Sie additive Überdispersion für Binomialverteilungen modellieren können, indem Sie einen zufälligen Achsenabschnitt mit einer Ebene pro Beobachtung hinzufügen. Gilt dies auch für die Poissonverteilung?

Gibt es einen besseren Weg, es zu tun? Gibt es andere Pakete, die Sie empfehlen würden?

Sie können Mehrebenen-GLMM mit einer Poisson-Verteilung (mit Überdispersion) unter Verwendung von R auf mehrere Arten anpassen. Nur wenige RPakete sind: lme4, MCMCglmm, armetc. Eine gute Referenz , um zu sehen ist Gelman und Hill (2007)

Ich werde ein Beispiel dafür mit rjagspackage in geben R. Es ist eine Schnittstelle zwischen Rund JAGS(wie OpenBUGSoder 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 … ich

$$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$

$\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)
}

Hier ist der RCode zu implementieren Verwendung es (sagen , es heißt: 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])

Sie können mit den Nachbauten Ihrer Parameter herumspielen und weitere Parameter einführen, um die Modellierung präziser zu gestalten ( wir denken gerne darüber nach ). Grundsätzlich kommt man auf die Idee.

Weitere Informationen zur Verwendung von rjagsund JAGSfinden Sie auf der Seite von John Myles White

Das lme4-Paket muss nicht verlassen werden, um eine Überdispersion zu berücksichtigen. Fügen Sie einfach einen zufälligen Effekt für die Beobachtungsnummer ein. Die genannten BUGS / JAGS-Lösungen sind wahrscheinlich zu viel für Sie, und wenn nicht, sollten Sie die leicht zu vergleichenden lme4-Ergebnisse haben.

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)

Dies wird hier diskutiert: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 informell und akademisch von Elston et al. (2001) .

Ich denke, dass das glmmADMB-Paket genau das ist, was Sie suchen.

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

Aber aus bayesianischer Sicht können Sie das MCMCglmm- Paket oder die BUGS / JAGS- Software verwenden, sie sind sehr flexibel und Sie können diese Art von Modell anpassen . (und die Syntax ist der von R sehr ähnlich)

EDIT dank @randel

Wenn Sie die Pakete glmmADMBund installieren möchten, R2admbist es besser, Folgendes zu tun:

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

Bisher gute Vorschläge. Hier ist noch einer. Mit der rhierNegbinRwFunktion des bayesmPakets können Sie ein hierarchisches negatives binomiales Regressionsmodell anpassen.