Bagaimana menganalisis data jumlah longitudinal: akuntansi untuk autokorelasi temporal di GLMM?

16

Halo guru statistik dan ahli pemrograman R,

Saya tertarik dalam memodelkan tangkapan hewan sebagai fungsi dari kondisi lingkungan dan hari dalam setahun. Sebagai bagian dari penelitian lain, saya memiliki jumlah tangkapan pada ~ 160 hari selama tiga tahun. Pada setiap hari ini saya memiliki suhu, curah hujan, kecepatan angin, kelembaban relatif, dll. Karena data dikumpulkan berulang kali dari 5 plot yang sama, saya menggunakan plot sebagai efek acak.

Pemahaman saya adalah bahwa nlme dapat dengan mudah menjelaskan autokorelasi temporal dalam residual tetapi tidak menangani fungsi tautan non-Gaussian seperti lme4 (yang tidak dapat menangani autokorelasi?). Saat ini, saya pikir mungkin bekerja untuk menggunakan paket nlme di R on log (hitung). Jadi solusi saya sekarang adalah menjalankan sesuatu seperti:

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)

di mana DOY = Day of the Year. Mungkin ada lebih banyak interaksi dalam model akhir, tetapi ini adalah ide umum saya. Saya juga berpotensi mencoba memodelkan struktur varians lebih jauh dengan sesuatu seperti

weights = v1Pow

Saya tidak yakin apakah ada cara yang lebih baik untuk dilakukan dengan regresi model campuran Poisson atau apa pun? Saya baru saja menemukan diskusi matematis di Bab 4 "Model Regresi untuk Analisis Rangkaian Waktu" oleh Kedem dan Fokianos. Itu sedikit di luar saya saat ini, terutama dalam aplikasi (pengkodean dalam R). Saya juga melihat solusi MCMC di Zuur et al. Buku Mixed Effects Models (Bab 23) dalam bahasa BUGS (menggunakan winBUGS atau JAG). Apakah itu pilihan terbaik saya? Apakah ada paket MCMC mudah di R yang akan menangani ini? Saya tidak terlalu terbiasa dengan teknik GAMM atau GEE tetapi akan bersedia untuk mengeksplorasi kemungkinan ini jika orang berpikir mereka akan memberikan wawasan yang lebih baik.Tujuan utama saya adalah membuat model untuk memprediksi penangkapan hewan mengingat kondisi lingkungan. Kedua, saya ingin menjelaskan apa yang hewan-hewan menanggapi dalam hal aktivitas mereka.

Setiap pemikiran tentang cara terbaik untuk melanjutkan (secara filosofis), bagaimana kode ini dalam R, atau dalam BUGS akan dihargai. Saya cukup baru untuk R dan BUGS (winBUGS) tetapi saya sedang belajar. Ini juga pertama kalinya saya mencoba membahas autokorelasi temporal.

Terima kasih, Dan

djhocking
sumber
1
Saya penggemar berat GEE pada umumnya, tetapi saya akan menghindari menggunakannya di sini karena Anda hanya memiliki 5 cluster (plot). Untuk berkinerja baik tanpa gejala, GEE biasanya membutuhkan lebih banyak (sekitar 40 atau lebih) jumlah cluster.
StatsStudent
Sebagai pemilik Mac, saya lebih mudah dengan STAN daripada WINBUGS.
eric_kernfeld

Jawaban:

3

Log mengubah respons Anda adalah sebuah opsi walaupun tidak ideal. Kerangka kerja GLM umumnya lebih disukai. Jika Anda tidak terbiasa dengan GLM maka mulailah dengan meninjau sebelum melihat ekstensi model campuran. Untuk data hitungan Poisson atau asumsi distribusi Binomial Negatif mungkin akan cocok. Binomial negatif diindikasikan jika variansnya lebih tinggi dari rata-rata yang mengindikasikan dispersi berlebihan ( https://en.wikipedia.org/wiki/Overdispersion ). Interpretasi estimasi parameter setara untuk keduanya.

Beberapa opsi ada di R dengan lme4 yang paling sering dikutip dalam pengalaman saya.

#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))

Tautan ini juga dapat membantu:

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

Keduanya oleh Ben Bolker, penulis lme4.

Saya belum menguji contoh tetapi mereka harus memberi Anda ide dari mana harus memulai. Harap berikan data jika Anda ingin memverifikasi implementasinya.

siswa-t
sumber