Я строю довольно сложную иерархическую байесовскую модель для мета-анализа с использованием R и JAGS. Упрощенно, два ключевых уровня модели имеют α j = ∑ h γ h ( j ) + ϵ j, где y i j - i- е наблюдение за конечной точкой (в данном случае , GM против урожайности не-GM культур) в исследовании j , α j - эффект для исследования j
Меня в первую очередь интересует оценка значений s. Это означает, что удаление переменных уровня обучения из модели не является хорошим вариантом.
Существует высокая корреляция между несколькими переменными уровня исследования, и я думаю, что это приводит к большим автокорреляциям в моих цепочках MCMC. Этот диагностический график иллюстрирует траектории цепочки (слева) и результирующую автокорреляцию (справа):
В результате автокорреляции я получаю эффективные размеры выборки 60-120 из 4 цепочек по 10 000 выборок в каждой.
У меня два вопроса, один явно объективный, а другой более субъективный.
Какие методы я могу использовать для решения этой проблемы автокорреляции, кроме прореживания, добавления новых цепочек и более длительного запуска сэмплера? Под «управлять» я подразумеваю «производить достаточно хорошие оценки в разумные сроки». С точки зрения вычислительной мощности, я использую эти модели на MacBook Pro.
Насколько серьезна эта степень автокорреляции? Обсуждения как здесь, так и в блоге Джона Крушке показывают, что, если мы просто запустим модель достаточно долго, «комковатая автокорреляция, вероятно, все усреднена» (Крушке), и, следовательно, это не имеет большого значения.
Вот код JAGS для модели, которая создала сюжет выше, на тот случай, если кому-то будет интересно узнать подробности:
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)
}
}
источник
Ответы:
Следуя предложению user777, похоже, что ответ на мой первый вопрос «использовать Stan». После переписывания модели в Stan, вот траектории (4 цепочки x 5000 итераций после выгорания):
и графики автокорреляции:
Намного лучше! Для полноты вот код Стэна:
источник