Jak dopasować model wielopoziomowy do nadmiernie rozproszonych efektów Poissona?

32

Chcę dopasować wielopoziomowy GLMM z rozkładem Poissona (z nadmierną dyspersją) za pomocą R. W tej chwili używam lme4, ale zauważyłem, że ostatnio quasipoissonrodzina została usunięta.

Widziałem gdzie indziej, że można modelować nadmierną dyspersję addytywną dla rozkładów dwumianowych, dodając losowy punkt przecięcia z jednym poziomem na obserwację. Czy dotyczy to również rozkładu Poissona?

Czy jest na to lepszy sposób? Czy są inne pakiety, które poleciłbyś?

George Michaelides
źródło

Odpowiedzi:

22

Możesz dopasować wielopoziomowy GLMM z rozkładem Poissona (z nadmierną dyspersją), używając R na wiele sposobów. Kilka Rpakiety są: lme4, MCMCglmm, arm, itd. Dobrym odniesienia, aby zobaczyć to Gelman i Hill (2007)

Podam przykład wykonania tego przy użyciu rjagspakietu w R. Jest to interfejs pomiędzy Ri JAGS(jak OpenBUGSlub WinBUGS).

log θ I J = β 0 + β 1 T R E t m e n t i + δ I J δ I J ~ N ( 0 , σ 2 ϵ ) i = 1 I ,

nijPoisson(θij)
logθij=β0+β1 Treatmenti+δij
δijN(0,σϵ2)
T r e a t m e n t i = 0  lub  1 , , J - 1,  jeżeli   obserwacja i t h należy do grupy leczenia  1 lub,  2 , , J
i=1I,j=1J
Treatmenti=0 or 1,,J1 if the ith observation belongs to treatment group 1, or, 2,,J

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

Oto Rkod do wdrożenia go używać (mówią, że nazywa się: 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])

Możesz bawić się z tylnymi parametrami i możesz wprowadzić więcej parametrów, aby uczynić modelowanie bardziej precyzyjnym ( lubimy to uważać ). Zasadniczo masz pomysł.

Aby uzyskać więcej informacji na temat używania rjagsi JAGS, zobacz stronę Johna Mylesa White'a

suncoolsu
źródło
Dzięki!! Dopiero niedawno zacząłem przyglądać się analizie bayesowskiej i wciąż trudno mi to pojąć. Myślę, że jest to okazja, aby dowiedzieć się nieco więcej na ten temat.
George Michaelides,
1
Dlaczego nie dyspersja gamma?
Patrick McCann,
2
@Patrick na pewno możesz to zrobić. Ale ponieważ biorę dziennik średniej, wolę normalny efekt disp. Rozkład normalny dziennika jest kolejnym sposobem modelowania rozkładów podobnych do rozkładu gamma. HTH.
suncoolsu,
20

Nie trzeba opuszczać pakietu lme4, aby uwzględnić nadmierną dyspersję; wystarczy dołączyć losowy efekt dla liczby obserwacji. Wymienione rozwiązania BŁĘDÓW / JAGSÓW są prawdopodobnie dla ciebie przesadą, a jeśli nie są, powinieneś mieć łatwe do dopasowania wyniki lme4 do porównania.

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)

Jest to omówione tutaj: http://article.gmane.org/gmane.comp.lang.r.lme4.devel/4727 nieformalnie i naukowo przez Elston i in. (2001) .

Patrick McCann
źródło
Co jeśli model składa się z dwóch zmiennych nominalnych, jednej zmiennej ciągłej (wszystkie jako efekty ustalone) i jednej zmiennej grupującej (efekt losowy) z interakcjami trzeciego rzędu, a ponadto liczba mierzonych obiektów jest równa liczbie obserwacji lub zapisów w zestaw danych? Jak mam to uwzględnić w modelu?
Ladislav Naďo
7

Myślę, że pakiet glmmADMB jest dokładnie tym, czego szukasz.

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

Ale z punktu widzenia bayesowskiego można użyć pakietu MCMCglmm lub oprogramowania BUGS / JAGS , są one bardzo elastyczne i można dopasować do tego rodzaju modelu. (a składnia jest zbliżona do R)

EDYCJA dzięki @randel

Jeśli chcesz zainstalować pakiety glmmADMBi R2admb, lepiej zrobić:

install.packages("glmmADMB", repos="http://glmmadmb.r-forge.r-project.org/repos"‌​)   
install.packages("R2admb")
dickoa
źródło
Uważam, że obecnie pakiet powinien zostać zainstalowany za pomocą install.packages("glmmADMB",repos="http://glmmadmb.r-forge.r-project.org/repos")plusa install.packages('R2admb').
Randel
5

Good suggestions so far. Here's one more. You can fit a hierarchical negative binomial regression model using the rhierNegbinRw function of the bayesm package.

Ari B. Friedman
źródło