Comment puis-je adapter un modèle à plusieurs niveaux pour les résultats de poissons trop dispersés?

32

Je veux adapter un GLMM à plusieurs niveaux avec une distribution de Poisson (avec sur-dispersion) en utilisant R. En ce moment j'utilise lme4 mais j'ai remarqué que récemment la quasipoissonfamille a été supprimée.

J'ai vu ailleurs que vous pouvez modéliser une surdispersion additive pour les distributions binomiales en ajoutant une interception aléatoire avec un niveau par observation. Est-ce que cela s'applique également à la distribution poisson?

Y a-t-il une meilleure façon de le faire? Y a-t-il d'autres forfaits que vous recommanderiez?

r mixed-model poisson-distribution lme4-nlme overdispersion George Michaelides
la source

Réponses:

22

Vous pouvez adapter GLMM multiniveau avec une distribution de Poisson (avec sur-dispersion) en utilisant R de plusieurs manières. Peu de Rpaquets sont: lme4, MCMCglmm, arm, etc. Une bonne référence pour voir est Gelman et Hill (2007)

Je vais vous donner un exemple de cela en utilisant rjagspackage in R. C'est une interface entre Ret JAGS(comme OpenBUGSou 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 ,

njej~Pojesson(θjej)
bûche⁡θjej=β0+β1 Treunetmentje+δjej
δjej~N(0,σε2)
je=1…je,j=1…J
Treunetmentje=0 ou 1,…,J-1 si la jeth l'observation appartient au groupe de traitement 1, ou, 2,…,J

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

Voici le Rcode pour mettre en œuvre l' utiliser (dire qu'il est nommé: 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])

Vous pouvez jouer avec les postérieurs de vos paramètres et vous pouvez introduire plus de paramètres pour vous rendre la modélisation plus précise ( nous aimons penser cela ). Fondamentalement, vous avez l'idée.

Pour plus de détails sur l'utilisation de rjagset JAGS, veuillez consulter la page de John Myles White

suncoolsu
la source
Merci!! J'ai récemment commencé à me pencher sur l'analyse bayésienne et je la trouve toujours un peu difficile à comprendre. Je suppose que c'est l'occasion d'en apprendre un peu plus à ce sujet.
George Michaelides
1
Pourquoi pas la dispersion gamma?
Patrick McCann
2
@Patrick, vous pouvez certainement le faire. Mais comme je prends un log de la moyenne, je préfère un effet de disp normal. Une distribution log-normale est une autre façon de modéliser des distributions similaires à la distribution gamma. HTH.
suncoolsu
20

Pas besoin de quitter le paquet lme4 pour tenir compte de la surdispersion; il suffit d'inclure un effet aléatoire pour le numéro d'observation. Les solutions BUGS / JAGS mentionnées sont probablement exagérées pour vous, et si elles ne le sont pas, vous devriez avoir les résultats lme4 faciles à adapter pour la comparaison.

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)

Ceci est discuté ici: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 de manière informelle et académique par Elston et al. (2001) .

Patrick McCann
la source
Que se passe-t-il si un modèle se compose de deux variables nominales, une variable continue (tous sous forme d'effets fixes) et une variable de regroupement (effet aléatoire) avec des interactions de troisième ordre et, de plus, le nombre de sujets mesurés est égal au nombre d'observations ou d'enregistrements dans le jeu de données? Comment dois-je couvrir cela dans le modèle?
Ladislav Naďo
7

Je pense que le package glmmADMB est exactement ce que vous recherchez.

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

Mais d'un point de vue bayésien, vous pouvez utiliser le package MCMCglmm ou le logiciel BUGS / JAGS , ils sont très flexibles et vous pouvez adapter ce type de modèle. (et la syntaxe est proche de celle de R)

EDIT grâce à @randel

Si vous souhaitez installer les packages glmmADMBet R2admbil est préférable de faire:

install.packages("glmmADMB", repos="http://glmmadmb.r-forge.r-project.org/repos"‌​)   
install.packages("R2admb")
Dickoa
la source
Je pense qu'actuellement le paquet devrait être installé via install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")plus install.packages('R2admb').
Randel
5

Bonnes suggestions jusqu'à présent. En voici un de plus. Vous pouvez adapter un modèle de régression binomiale négative hiérarchique à l'aide de la rhierNegbinRwfonction du bayesmpackage.

Ari B. Friedman
la source