Mengelola autokorelasi tinggi di MCMC

10

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

ysayaj=αj+ϵsaya
αj=hγh(j)+ϵj
ysayajsayajαjjγhϵs adalah istilah kesalahan. Perhatikan bahwa s bukan koefisien pada variabel dummy. Sebagai gantinya, ada variabel berbeda untuk nilai tingkat studi yang berbeda. Misalnya, ada untuk negara-negara berkembang dan untuk negara-negara maju. γγγdevelHaihalsayangγdevelHaihaled

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):
autokorelasi tinggi dalam output MCMC

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.

  1. 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.

  2. 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)
}
}
Dan Hicks
sumber
1
Untuk apa nilainya, model bertingkat kompleks cukup banyak alasan bahwa Stan ada, untuk alasan persis yang Anda identifikasi.
Sycorax berkata Reinstate Monica
Saya awalnya mencoba membangun ini di Stan, beberapa bulan yang lalu. Studi melibatkan sejumlah temuan yang berbeda, yang (setidaknya pada saat itu; saya belum memeriksa untuk melihat apakah ada perubahan) memerlukan penambahan lapisan kompleksitas pada kode dan berarti Stan tidak dapat memanfaatkan perhitungan matriks yang membuatnya begitu cepat.
Dan Hicks
1
Saya tidak terlalu memikirkan kecepatan, seperti efisiensi HMC yang mengeksplorasi posterior. Pemahaman saya adalah bahwa karena HMC dapat mencakup lebih banyak dasar, setiap iterasi memiliki autokorelasi yang lebih rendah.
Sycorax berkata Reinstate Monica
Oh, ya, itu poin yang menarik. Saya akan memasukkannya ke dalam daftar hal yang harus saya coba.
Dan Hicks

Jawaban:

6

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):
masukkan deskripsi gambar di sini Dan plot autokorelasi:
masukkan deskripsi gambar di sini

Jauh lebih baik! Untuk kelengkapan, inilah kode Stan:

data {                          // Data: Exogenously given information
// Data on totals
int n;                      // Number of distinct finding i
int n_study;                // Number of distinct studies j

// Finding-level data
vector[n] y;                // Endpoint for finding i
int study_n[n_study];       // # findings for study j

// Study-level data
int countr[n_study];        // Country type for study j
int studytype[n_study];     // Study type for study j
int jourtype[n_study];      // Was study j published in a journal?
int industrytype[n_study];  // Was study j funded by industry?

// Top-level constants set in R call
real sigma_alpha_bound;     // Upper bound for noise in alphas
real gamma_prior_exp;       // Prior expected value of gamma
real gamma_sigma_bound;     // Upper bound for noise in gammas
real eta_bound;             // Upper bound for etas
}

transformed data {
// Constants set here
int countr_levels;          // # levels for countr
int study_levels;           // # levels for studytype
int jour_levels;            // # levels for jourtype
int industry_levels;        // # levels for industrytype
countr_levels <- 2;
study_levels <- 2;
jour_levels <- 2;
industry_levels <- 2;
}

parameters {                    // Parameters:  Unobserved variables to be estimated
vector[n_study] alpha;      // Study-level mean
real<lower = 0, upper = sigma_alpha_bound> sigma_alpha;     // Noise in alphas

vector<lower = 0, upper = 100>[n_study] sigma;          // Study-level standard deviation

// Gammas:  contextual effects on study-level means
// Country-type effect and noise in its estimate
vector[countr_levels] gamma_countr;     
vector<lower = 0, upper = gamma_sigma_bound>[countr_levels] sigma_countr;
// Study-type effect and noise in its estimate
vector[study_levels] gamma_study;
vector<lower = 0, upper = gamma_sigma_bound>[study_levels] sigma_study;
vector[jour_levels] gamma_jour;
vector<lower = 0, upper = gamma_sigma_bound>[jour_levels] sigma_jour;
vector[industry_levels] gamma_industry;
vector<lower = 0, upper = gamma_sigma_bound>[industry_levels] sigma_industry;


// Etas:  contextual effects on study-level standard deviation
vector<lower = 0, upper = eta_bound>[countr_levels] eta_countr;
vector<lower = 0, upper = eta_bound>[study_levels] eta_study;
vector<lower = 0, upper = eta_bound>[jour_levels] eta_jour;
vector<lower = 0, upper = eta_bound>[industry_levels] eta_industry;
}

transformed parameters {
vector[n_study] alpha_hat;                  // Fitted alpha, based only on gammas
vector<lower = 0>[n_study] sigma_hat;       // Fitted sd, based only on sigmasq_hat

for (j in 1:n_study) {
    alpha_hat[j] <- gamma_countr[countr[j]] + gamma_study[studytype[j]] + 
                    gamma_jour[jourtype[j]] + gamma_industry[industrytype[j]];
    sigma_hat[j] <- sqrt(eta_countr[countr[j]]^2 + eta_study[studytype[j]]^2 +
                        eta_jour[jourtype[j]] + eta_industry[industrytype[j]]);
}
}

model {
// Technique for working w/ ragged data from Stan manual, page 135
int pos;
pos <- 1;
for (j in 1:n_study) {
    segment(y, pos, study_n[j]) ~ normal(alpha[j], sigma[j]);
    pos <- pos + study_n[j];
}

// Study-level mean = fitted alpha + Gaussian noise
alpha ~ normal(alpha_hat, sigma_alpha);

// Study-level variance = gamma distribution w/ mean sigma_hat
sigma ~ gamma(.1 * sigma_hat, .1);

// Priors for gammas
gamma_countr ~ normal(gamma_prior_exp, sigma_countr);
gamma_study ~ normal(gamma_prior_exp, sigma_study);
gamma_jour ~ normal(gamma_prior_exp, sigma_study);
gamma_industry ~ normal(gamma_prior_exp, sigma_study);
}
Dan Hicks
sumber