I would like to conduct a Bayesian meta-analysis on some plant population growth rates that I have already calculated. The structure of the data is nested because I have several plant species within multiple sites. For example, I have data on Species_1 and Species_2 in Site_1, Species_1 and Species_3 in Site_2, Species_3 and Species_4 in Site_3 and so on. I calculated the growth rates of each species population (i.e. Species_1 in Site_1, etc.) using a Bayesian hierarchical model using JAGS because I was interested in including the observation error during the sampling process. Thus, for each population I have a posterior distribution of the average growth rate.
Now, I would like to summarize the results into an overall growth rate and I was thinking about doing a Bayesian meta-analysis following the example of this paper:
Haase, P. et al. The recovery of European freshwater biodiversity has come to a halt. Nature 620, 582–588 (2023).
I was able to easily replicate their analysis using the R package brms, which uses Stan, but I would like to to it in JAGS to be consistent with the previous analysis. I have used the following model in brms:
set.seed(5345)
data_brms <- list(Site = paste("Site", sample(1:5, replace = T, size = 100), sep = '_'),
Species = paste("Species", sample(1:10, replace = T, size = 100), sep = '_'),
growth_rate_mean = rnorm(100, mean = 0, sd = 1),
growth_rate_sd = rgamma(100, shape = 2, scale = 1/2))
brms_prior = c(set_prior("normal(0, 3)", class = "Intercept"),
set_prior("gamma(0.001, 0.001)", class = "sd"))
mod_brms <- brm(growth_rate_mean|se(growth_rate_sd) ~ 1 + (1|Site/Species),
data = dummy_data,
prior = brms_prior,
cores = 4,
chains = 4,
iter = 10000,
warmup = 1000)
Which results in the following model:
functions {
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
vector<lower=0>[N] se; // known sampling error
// data for group-level effects of ID 1
int<lower=1> N_1; // number of grouping levels
int<lower=1> M_1; // number of coefficients per level
int<lower=1> J_1[N]; // grouping indicator per observation
// group-level predictor values
vector[N] Z_1_1;
// data for group-level effects of ID 2
int<lower=1> N_2; // number of grouping levels
int<lower=1> M_2; // number of coefficients per level
int<lower=1> J_2[N]; // grouping indicator per observation
// group-level predictor values
vector[N] Z_2_1;
int prior_only; // should the likelihood be ignored?
}
transformed data {
vector<lower=0>[N] se2 = square(se);
}
parameters {
real Intercept; // temporary intercept for centered predictors
vector<lower=0>[M_1] sd_1; // group-level standard deviations
vector[N_1] z_1[M_1]; // standardized group-level effects
vector<lower=0>[M_2] sd_2; // group-level standard deviations
vector[N_2] z_2[M_2]; // standardized group-level effects
}
transformed parameters {
real sigma = 0; // dispersion parameter
vector[N_1] r_1_1; // actual group-level effects
vector[N_2] r_2_1; // actual group-level effects
real lprior = 0; // prior contributions to the log posterior
r_1_1 = (sd_1[1] * (z_1[1]));
r_2_1 = (sd_2[1] * (z_2[1]));
lprior += normal_lpdf(Intercept | 0, 3);
lprior += cauchy_lpdf(sd_1 | 0, 1)
- 1 * cauchy_lccdf(0 | 0, 1);
lprior += cauchy_lpdf(sd_2 | 0, 1)
- 1 * cauchy_lccdf(0 | 0, 1);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = rep_vector(0.0, N);
mu += Intercept;
for (n in 1:N) {
// add more terms to the linear predictor
mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_2_1[J_2[n]] * Z_2_1[n];
}
target += normal_lpdf(Y | mu, se);
}
// priors including constants
target += lprior;
target += std_normal_lpdf(z_1[1]);
target += std_normal_lpdf(z_2[1]);
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept;
}
I would like to replicate this model in JAGS, but I cannot get the nested random effects right. I tried followin this answer: Nested Random effect in JAGS/ WinBUGS but it did not work so I followed the second answer from this post: Three-level nested mixed-effects model but the results differ a lot from the ones I get from brms. My JAGS model is:
set.seed(5345)
data_jags <- list(N = 100,
Site = as.factor(paste("Site", sample(1:5, replace = T, size = 100), sep = '_')),
Species = as.factor(paste("Species", sample(1:10, replace = T, size = 100), sep = '_')),
growth_rate_mean = rnorm(100, mean = 0, sd = 1),
growth_rate_sd = rgamma(100, shape = 2, scale = 1/2),
nSites = 5,
nSpecies = 10)
# Set the interaction term:
interaction <- cbind(data_jags$Species, data_jags$Site)
interaction <- interaction[!duplicated(interaction),]
data_jags[[8]] <- as.factor(interaction[,2])
names(data_jags)[8] <- 'interaction'
# Define model:
sink("JAGSmodel.txt")
cat("
model{
for(i in 1:N){
P[i] <- 1/(growth_rate_sd[i]^2)
growth_rate_mean[i] ~ dnorm(regression_fitted[i], P[i])
regression_fitted[i] <- intercept + Species_randomeffect[Species[i]] + Site_randomeffect[Site[i]]
}
# priors:
regression_precision ~ dgamma(0.001, 0.001)
intercept ~ dnorm(0, 1/9)
# random effects
for(i in 1:nSites){Site_randomeffect[i] ~ dnorm(0, Site_precision)}
for(i in 1:nSpecies){Species_randomeffect[i] ~ dnorm(Site_randomeffect[interaction[i]], Species_precision)}
Species_precision ~ dgamma(0.001, 0.001)
Site_precision ~ dgamma(0.001, 0.001)
}
}
",fill = TRUE)
sink()
parameters <- c('intercept', 'Site_sd', 'Species_sd')
mod_jags <- run.jags(model = 'JAGSmodel.txt',
data = dummy_data,
n.chains = 4,
burnin = 1000,
sample = 10000,
thin = 1,
parameters)
The results from each model are very different and I do not know if this is because brms is having trouble converging (it gives a warning with the dummy data) or because the JAGS model is incorrectly defined. If I use my original data brms works perfectly but the results from JAGS do not resemble the results from brms.
