Zarządzanie wysoką autokorelacją w MCMC

10

Buduję raczej złożony hierarchiczny model bayesowski do metaanalizy przy użyciu R i JAGS. Upraszczając nieco, dwa kluczowe poziomy modelu mają gdzie jest obserwacją punkt końcowy (w tym przypadku plony GMO i GMO) w badaniu , jest efektem dla badania , są efektami dla różnych zmiennych na poziomie badania (status rozwoju gospodarczego kraju, w którym badanie zostało przeprowadzone, gatunki upraw, metoda badania itp.) indeksowane przez rodzinę funkcji i

yjajot=αjot+ϵja
αjot=hγh(jot)+ϵjot
yjajotjajotαjotjotγhϵs są terminami błędów. Zauważ, że nie są współczynnikami dla zmiennych fikcyjnych. Zamiast tego istnieją różne zmienne dla różnych wartości na poziomie badania. Na przykład istnieje dla krajów rozwijających się i dla krajów rozwiniętych. γγγremivmilopjansolγremivmilopmire

Interesuje mnie przede wszystkim szacowanie wartości . Oznacza to, że usunięcie zmiennych na poziomie badania z modelu nie jest dobrym rozwiązaniem. γ

Istnieje kilka korelacji między kilkoma zmiennymi na poziomie badania i myślę, że powoduje to dużą autokorelację w moich łańcuchach MCMC. Ten wykres diagnostyczny ilustruje trajektorie łańcucha (po lewej) i wynikową autokorelację (po prawej):
wysoka autokorelacja na wyjściu MCMC

W wyniku autokorelacji uzyskuję skuteczne wielkości próbek 60–120 z 4 łańcuchów po 10 000 próbek.

Mam dwa pytania, jedno wyraźnie obiektywne, a drugie bardziej subiektywne.

  1. Jakie techniki można wykorzystać do zarządzania tym problemem autokorelacji, oprócz przerzedzania, dodawania kolejnych łańcuchów i dłuższego uruchamiania samplera? Przez „zarządzaj” mam na myśli „przedstawianie dość dobrych oszacowań w rozsądnym czasie”. Pod względem mocy obliczeniowej używam tych modeli na MacBooku Pro.

  2. Jak poważny jest ten stopień autokorelacji? Dyskusje zarówno tutaj, jak i na blogu Johna Kruschke sugerują, że jeśli poprowadzimy ten model wystarczająco długo, „grudkowata autokorelacja prawdopodobnie wszystkie zostały uśrednione” (Kruschke), a więc nie jest to naprawdę wielka sprawa.

Oto kod JAGS dla modelu, który wyprodukował powyższą fabułę, na wypadek, gdyby ktoś był wystarczająco zainteresowany, aby przejrzeć szczegóły:

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
źródło
1
Ze względu na swoją wartość złożone modele wielopoziomowe są właściwie powodem, dla którego Stan istnieje, z dokładnie tych powodów, które można zidentyfikować.
Sycorax mówi Przywróć Monikę
Początkowo próbowałem zbudować to w Stanie, kilka miesięcy temu. Badania obejmowały różną liczbę ustaleń, które (przynajmniej wtedy; nie sprawdziłem, czy wszystko się zmieniło) wymagały dodania kolejnej warstwy złożoności do kodu i oznaczały, że Stan nie mógł skorzystać z obliczeń macierzowych dzięki czemu jest tak szybki.
Dan Hicks,
1
Nie myślałem o prędkości tak bardzo, jak o wydajności, z jaką HMC bada tył. Rozumiem, że ponieważ konsola HMC może objąć znacznie więcej gruntów, każda iteracja ma niższą autokorelację.
Sycorax mówi Przywróć Monikę
Och, tak, to interesujący punkt. Umieszczę to na mojej liście rzeczy do wypróbowania.
Dan Hicks,

Odpowiedzi:

6

Zgodnie z sugestią użytkownika777 wygląda na to, że odpowiedź na moje pierwsze pytanie brzmi „użyj Stana”. Po przepisaniu modelu w Stan, oto trajektorie (4 łańcuchy x 5000 iteracji po wypaleniu):
wprowadź opis zdjęcia tutaj I wykresy autokorelacji:
wprowadź opis zdjęcia tutaj

Dużo lepiej! Dla kompletności, oto kod 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
źródło