Jak analizować dane dotyczące liczby podłużnej: rozliczanie czasowej autokorelacji w GLMM?

16

Witaj guru statystyczni i kreatorzy programowania R,

Interesuje mnie modelowanie chwytów zwierząt jako funkcji warunków środowiskowych i dnia w roku. W ramach innego badania mam liczbę przechwyceń przez ~ 160 dni w ciągu trzech lat. Na każdy z tych dni mam temperaturę, opady, prędkość wiatru, wilgotność względną itp. Ponieważ dane były zbierane wielokrotnie z tych samych 5 wykresów, używam wykresu jako efektu losowego.

Rozumiem, że nlme może łatwo uwzględnić czasową autokorelację w resztkach, ale nie obsługuje funkcji łącza niegaussowskiego, takich jak lme4 (które nie mogą obsłużyć autokorelacji?). Obecnie myślę, że może zadziałać użycie pakietu nlme w R na logu (liczba). Więc moim rozwiązaniem w tej chwili byłoby uruchomienie czegoś takiego:

m1 <- lme(lcount ~ AirT + I(AirT^2) + RainAmt24 + I(RainAmt24^2) + RHpct + windspeed + 
      sin(2*pi/360*DOY) + cos(2*pi/360*DOY), random = ~1|plot, correlation =
      corARMA(p = 1, q = 1, form = ~DOY|plot), data = Data)

gdzie DOY = dzień roku. W ostatecznym modelu może być więcej interakcji, ale to jest mój ogólny pomysł. Mógłbym również potencjalnie spróbować modelować strukturę wariancji za pomocą czegoś podobnego

weights = v1Pow

Nie jestem pewien, czy jest lepszy sposób na regresję modelu mieszanego Poissona czy coś takiego? Właśnie znalazłem dyskusję matematyczną w rozdziale 4 „Modelów regresji do analizy szeregów czasowych” Kedema i Fokianosa. W tej chwili było to trochę poza mną, szczególnie w aplikacji (kodowanie w R). Widziałem także rozwiązanie MCMC w Zuur i in. Książka modeli efektów mieszanych (rozdz. 23) w języku BŁĘDÓW (używając winBUGS lub JAG). Czy to moja najlepsza opcja? Czy istnieje prosty pakiet MCMC w R, który by to obsługiwał? Naprawdę nie znam technik GAMM ani GEE, ale chętnie zbadam te możliwości, jeśli ludzie będą sądzić, że zapewnią lepszy wgląd.Moim głównym celem jest stworzenie modelu do przewidywania chwytania zwierząt w danych warunkach środowiskowych. Po drugie, chciałbym wyjaśnić, na co reagują zwierzęta pod względem ich aktywności.

Będziemy wdzięczni za wszelkie przemyślenia na temat najlepszego sposobu postępowania (filozoficznego), sposobu kodowania tego w R lub BŁĘDACH. Jestem dość nowy w R i BŁĘDACH (winBUGS), ale się uczę. To także pierwszy raz, kiedy próbowałem rozwiązać czasową autokorelację.

Dzięki, Dan

djhocking
źródło
1
Ogólnie jestem wielkim fanem GEE, ale unikałbym go tutaj, ponieważ masz tylko 5 klastrów (wykresów). Aby działać dobrze asymptotycznie, GEE zwykle wymaga większej (około 40) liczby klastrów.
StatsStudent
Jako właściciel Maca łatwiej mi było używać STAN niż WINBUGS.
eric_kernfeld

Odpowiedzi:

3

Log przekształcający twoją odpowiedź jest opcją, choć nie idealną. Framework GLM jest ogólnie preferowany. Jeśli nie znasz GLM, zacznij od ich przejrzenia, zanim przejrzysz rozszerzenia modeli mieszanych. W przypadku danych zliczeń prawdopodobne są odpowiednie założenia Poissona lub ujemnego dwumianowego rozkładu. Ujemna dwumianowa jest wskazywana, jeśli wariancja jest wyższa niż średnia wskazująca na nadmierną dyspersję ( https://en.wikipedia.org/wiki/Overdispersion ). Interpretacja oszacowań parametrów jest równoważna dla obu.

W R istnieje kilka opcji, a moje doświadczenie najczęściej cytowało lme4.

#glmer
library(lme4) 
glmer(count ~ AirT + I(AirT^2) + RainAmt24 + I(RainAmt24^2) + RHpct + windspeed + sin(2*pi/360*DOY) + cos(2*pi/360*DOY) + (1|plot), family=poisson, data = Data) 
# use glmer.nb with identical syntax but no family for negative binomial.

# glmmADMB with negative binomial
install.packages("glmmADMB", repos=c("http://glmmadmb.r-forge.r-project.org/repos", getOption("repos")),type="source") 
require(glmmADMB)
glmmadmb(count ~ AirT + I(AirT^2) + RainAmt24 + I(RainAmt24^2) + RHpct + windspeed + sin(2*pi/360*DOY) + cos(2*pi/360*DOY) + (1|plot), 
           family="nbinom", zeroInflation=FALSE, data=Data)

# glmmPQL, requires an estimate for theta which can be obtained from a 
# glm model in which the correlation structure is ignored.
library(MASS)
glmmPQL(count ~ AirT + I(AirT^2) + RainAmt24 + I(RainAmt24^2) + RHpct + windspeed + sin(2*pi/360*DOY) + cos(2*pi/360*DOY) , random = list(~1 | plot), data = Data,family = negative.binomial(theta = 4.22, link = log))

Te linki mogą być również pomocne:

https://udrive.oit.umass.edu/xythoswfs/webui/_xy-11096203_1-t_yOxYgf1s http://www.cell.com/trends/ecology-evolution/pdf/S0169-5347(09)00019-6.pdf

Oba są autorstwa Ben Bolkera, autora lme4.

Nie testowałem przykładów, ale powinny dać ci pojęcie, od czego zacząć. Podaj dane, jeśli chcesz zweryfikować ich wdrożenie.

t-student
źródło