Saya sedang membangun model Bayesian hierarkis yang agak rumit untuk meta-analisis menggunakan R dan JAGS. Menyederhanakan sedikit, dua tingkat kunci dari model memiliki mana adalah th pengamatan titik akhir (dalam hal ini, hasil panen GM vs non-GM) dalam studi , adalah efek untuk studi , s adalah efek untuk berbagai variabel tingkat studi (status perkembangan ekonomi negara tempat studi dilakukan, spesies tanaman, metode studi, dll) diindeks oleh keluarga fungsi , dan
Saya terutama tertarik untuk memperkirakan nilai-nilai s. Ini berarti bahwa menjatuhkan variabel tingkat studi dari model bukanlah pilihan yang baik.
Ada korelasi tinggi di antara beberapa variabel tingkat studi, dan saya pikir ini menghasilkan autokorelasi besar dalam rantai MCMC saya. Plot diagnostik ini menggambarkan lintasan rantai (kiri) dan autokorelasi yang dihasilkan (kanan):
Sebagai konsekuensi dari autokorelasi, saya mendapatkan ukuran sampel efektif 60-120 dari 4 rantai masing-masing 10.000 sampel.
Saya punya dua pertanyaan, yang satu jelas obyektif dan yang lain lebih subyektif.
Selain menipis, menambahkan lebih banyak rantai, dan menjalankan sampler lebih lama, teknik apa yang dapat saya gunakan untuk mengelola masalah autokorelasi ini? Dengan "mengelola", maksud saya, "menghasilkan perkiraan yang cukup baik dalam jumlah waktu yang wajar." Dalam hal daya komputasi, saya menjalankan model ini pada MacBook Pro.
Seberapa serius tingkat autokorelasi ini? Diskusi di sini dan di blog John Kruschke menunjukkan bahwa, jika kita menjalankan modelnya cukup lama, "autokorelasi yang kasar mungkin semuanya telah dirata-ratakan" (Kruschke) dan jadi itu bukan masalah besar.
Berikut kode JAGS untuk model yang menghasilkan plot di atas, kalau-kalau ada orang yang cukup tertarik untuk menelusuri detail:
model {
for (i in 1:n) {
# Study finding = study effect + noise
# tau = precision (1/variance)
# nu = normality parameter (higher = more Gaussian)
y[i] ~ dt(alpha[study[i]], tau[study[i]], nu)
}
nu <- nu_minus_one + 1
nu_minus_one ~ dexp(1/lambda)
lambda <- 30
# Hyperparameters above study effect
for (j in 1:n_study) {
# Study effect = country-type effect + noise
alpha_hat[j] <- gamma_countr[countr[j]] +
gamma_studytype[studytype[j]] +
gamma_jour[jourtype[j]] +
gamma_industry[industrytype[j]]
alpha[j] ~ dnorm(alpha_hat[j], tau_alpha)
# Study-level variance
tau[j] <- 1/sigmasq[j]
sigmasq[j] ~ dunif(sigmasq_hat[j], sigmasq_hat[j] + pow(sigma_bound, 2))
sigmasq_hat[j] <- eta_countr[countr[j]] +
eta_studytype[studytype[j]] +
eta_jour[jourtype[j]] +
eta_industry[industrytype[j]]
sigma_hat[j] <- sqrt(sigmasq_hat[j])
}
tau_alpha <- 1/pow(sigma_alpha, 2)
sigma_alpha ~ dunif(0, sigma_alpha_bound)
# Priors for country-type effects
# Developing = 1, developed = 2
for (k in 1:2) {
gamma_countr[k] ~ dnorm(gamma_prior_exp, tau_countr[k])
tau_countr[k] <- 1/pow(sigma_countr[k], 2)
sigma_countr[k] ~ dunif(0, gamma_sigma_bound)
eta_countr[k] ~ dunif(0, eta_bound)
}
# Priors for study-type effects
# Farmer survey = 1, field trial = 2
for (k in 1:2) {
gamma_studytype[k] ~ dnorm(gamma_prior_exp, tau_studytype[k])
tau_studytype[k] <- 1/pow(sigma_studytype[k], 2)
sigma_studytype[k] ~ dunif(0, gamma_sigma_bound)
eta_studytype[k] ~ dunif(0, eta_bound)
}
# Priors for journal effects
# Note journal published = 1, journal published = 2
for (k in 1:2) {
gamma_jour[k] ~ dnorm(gamma_prior_exp, tau_jourtype[k])
tau_jourtype[k] <- 1/pow(sigma_jourtype[k], 2)
sigma_jourtype[k] ~ dunif(0, gamma_sigma_bound)
eta_jour[k] ~ dunif(0, eta_bound)
}
# Priors for industry funding effects
for (k in 1:2) {
gamma_industry[k] ~ dnorm(gamma_prior_exp, tau_industrytype[k])
tau_industrytype[k] <- 1/pow(sigma_industrytype[k], 2)
sigma_industrytype[k] ~ dunif(0, gamma_sigma_bound)
eta_industry[k] ~ dunif(0, eta_bound)
}
}
sumber
Jawaban:
Mengikuti saran dari user777, sepertinya jawaban untuk pertanyaan pertama saya adalah "gunakan Stan." Setelah menulis ulang model di Stan, berikut adalah lintasannya (4 rantai x 5.000 iterasi setelah burn-in):
Dan plot autokorelasi:
Jauh lebih baik! Untuk kelengkapan, inilah kode Stan:
sumber