Je construis un modèle bayésien hiérarchique assez complexe pour une méta-analyse utilisant R et JAGS. Pour simplifier un peu, les deux niveaux clés du modèle ont où est la ème observation de la (dans ce cas, les rendements des cultures GM vs non GM) dans l'étude , est l'effet pour l'étude , les s sont des effets pour diverses variables au niveau de l'étude (le statut de développement économique du pays où le étude a été réalisée, espèces cultivées, méthode d'étude, etc.) indexées par une famille de fonctions , et le
Je suis principalement intéressé par l'estimation des valeurs des . Cela signifie que supprimer des variables de niveau étude du modèle n'est pas une bonne option.
Il existe une forte corrélation entre plusieurs des variables au niveau de l'étude, et je pense que cela produit de grandes autocorrélations dans mes chaînes MCMC. Ce tracé de diagnostic illustre les trajectoires de la chaîne (à gauche) et l'autocorrélation qui en résulte (à droite):
À la suite de l'autocorrélation, j'obtiens des tailles d'échantillons efficaces de 60 à 120 à partir de 4 chaînes de 10 000 échantillons chacune.
J'ai deux questions, l'une clairement objective et l'autre plus subjective.
Outre l'amincissement, l'ajout de chaînes et le fonctionnement de l'échantillonneur plus longtemps, quelles techniques puis-je utiliser pour gérer ce problème d'autocorrélation? Par «gérer», je veux dire «produire des estimations raisonnablement bonnes dans un délai raisonnable». En termes de puissance de calcul, j'exécute ces modèles sur un MacBook Pro.
Quelle est la gravité de ce degré d'autocorrélation? Les discussions ici et sur le blog de John Kruschke suggèrent que, si nous exécutons le modèle assez longtemps, "l'autocorrélation grumeleuse a probablement toutes été moyennée" (Kruschke) et donc ce n'est pas vraiment un gros problème.
Voici le code JAGS pour le modèle qui a produit l'intrigue ci-dessus, juste au cas où quelqu'un serait suffisamment intéressé pour parcourir les détails:
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)
}
}
la source
Réponses:
Suite à la suggestion de user777, il semble que la réponse à ma première question soit "utilisez Stan". Après avoir réécrit le modèle dans Stan, voici les trajectoires (4 chaînes x 5000 itérations après rodage):
Et les tracés d'autocorrélation:
Bien mieux! Pour être complet, voici le code Stan:
la source