Summary: I fit three Bayesian models to understand the problem of over-/under-dispersion when modeling this dataset about species trends across Europe. A mixture of two gaussians seems to handle the heterogeneity better than a model with a single normal component (at least according to the DHARMa residual QQ plot) but the mixture model has issues of its own.
Ultimately, I suspect that it would be more productive to re-examine the available data and modify/expand the model in the light of what you know about the scientific problem. For example: How are the std. errors of the species trends calculated; can you improve the estimation procedure? Can you specify one model that simultaneously models species trends, urbanization trends and correlations due to geography, so that all uncertainties are reflected appropriately? Can you incorporate additional information about the country/site climate and other features relevant to the species biology. (The species are all butterflies? I notice that not all species are found across all the sites. This may be important as well.)
Now to the analysis.
In this case the outcome variable, estimate, is derived from a preliminary analysis. A comment explains that estimates are "coefficients of linear models of abundance versus year". As a result, each estimate comes with an estimate of its standard error.
Let's start with a model with known measurement error. Here is the specification in {brms}; the formula is similar to the glmmTMB formula except the measurement error is specified via se(...), not via a weights argument.
M1 <- brm(
estimate | se(std.error, sigma = FALSE) ~ urb_trend + (1 | SPECIES) + (1 | SITE_ID),
data = final_df
)
Note: I simplify the random effects due to geographic location; this has little effect on how the model fits the data. (The R code to reproduce the analysis is attached at the end.)
And here are two posterior predictive (PP) checks that compare the empirical distribution of the data y (ie. the species trends) to the distribution of data yrep simulated from the fitted model. On the left is a LOO-PIT QQ plot (LOO = leave-one-out cross validation, PIT = probability integral transform); it's very similar to a DHARMa residual QQ plot as DHARMa also uses PIT by default. On the right the kernel density estimate of y is overlaid with the densities of 100 replicates yrep sampled from the fitted model.
Based on these two PP checks, the overall fit is not bad. I see three issues indicating lack of fit: (a) the model predicts more ys close to 0 that are observed in the data (that's why the yrep density curves peak above the y density); (b) there is an asymmetric bump at around y = 0.3 which the model doesn't predict (this can only be resolved by adding a predictor that explains this "feature" of the species trends), and (c) the model predicts large (in absolute value) ys outside of the observed [-1,1] range.
Next, let's try a model that ignores the "known" std.error and instead estimates it. This model assumes all data points have the same residual std. deviation, σ.
M2 <- brm(
estimate ~ urb_trend + (1 | SPECIES) + (1 | SITE_ID),
data = final_df
)
The PP checks indicate that M2 fits the data worse than M1. This is expected as M2 ignores information (the estimated std. error of the species trends). As a result the model is overdispersed: it predicts too few small-size outcomes (close to 0) and too many medium-size outcomes (close to -0.3 or +0.3).
So far we have two unsatisfactory models but one is under-dispersed and the other over-dispersed. Next I combine them in a mixture of two gaussians: one component corresponds to M1 and has known std. errors, the other component corresponds to M2 and has unknown residual std. deviation, σ. The mixing proportion, λ, is the weight of M1 in the mixture (and 1-λ is the weight of M2); λ another parameter to be estimated.
I don't know how to specify this model with {brms}, so I implement it directly in Stan (see code below). Here are the (final) PP checks. I also followed the DHARMa for Bayesians instructions to make a DHARMa residual plot.
The diagnostics indicate that M3 = mixture(M1,M2) handles heterogeneity in the uncertainties about species trends better than either M1 or M2. However, some lack of fit at 0 remains. All three models make predictions outside of the observed range [-1,1] for the species trend estimate. As I wrote in the first paragraph, I believe it will be easier to address these issues by reviewing the procedure to obtain the estimates and their std.errors as well as by thinking about what other information may be available and pertinent to ecological problem being studied.
R code:
library("brms")
library("DHARMa")
library("tidyverse")
library("bayesplot")
library("posterior")
library("cmdstanr")
options(
brms.backend = "cmdstanr",
brms.threads = 4,
mc.cores = 4
)
final_df <- read_csv(here::here("URBAN_TRENDS-main", "final_df.csv"))
# Y ~ normal(μ, SE) with SE known
M1 <- brm(
estimate | se(std.error, sigma = FALSE) ~ urb_trend + (1 | SPECIES) + (1 | SITE_ID),
data = final_df
)
# Y ~ normal(μ, σ)
M2 <- brm(
estimate ~ urb_trend + (1 | SPECIES) + (1 | SITE_ID),
data = final_df
)
pp_check(M1, type = "loo_pit_qq")
pp_check(M1, type = "dens_overlay", ndraws = 100)
# Y ~ mixture(normal(μ, SE), normal(μ, σ))
# This model is fitted in Stan using the cmdstanr interface.
code <- "
data {
int N;
int N_A;
int N_B;
vector[N] X;
vector[N] Y;
vector<lower=0>[N] S;
array[N] int<lower=1, upper=N_A> A;
array[N] int<lower=1, upper=N_B> B;
}
parameters {
real alpha;
real beta;
real<lower=0> sigma;
real<lower=0> tau_a;
real<lower=0> tau_b;
real<lower=0, upper=1> lambda;
vector[N_A] z_a;
vector[N_B] z_b;
}
transformed parameters {
vector[N] mu = alpha + beta * X + tau_a * z_a[A] + tau_b * z_b[B];
}
model {
alpha ~ std_normal();
beta ~ std_normal();
tau_a ~ normal(0, 1);
tau_b ~ normal(0, 1);
sigma ~ std_normal();
z_a ~ std_normal();
z_b ~ std_normal();
lambda ~ beta(9, 1);
// There is (currently) no way to vectorize mixture models at the observation level
for (n in 1:N) {
target += log_mix(
lambda,
normal_lpdf(Y[n] | mu[n], S[n]),
normal_lpdf(Y[n] | mu[n], sigma)
);
}
}
generated quantities {
vector[N] log_lik;
array[N] real Yrep;
for (n in 1:N) {
log_lik[n] = log_mix(
lambda,
normal_lpdf(Y[n] | mu[n], S[n]),
normal_lpdf(Y[n] | mu[n], sigma)
);
int u = bernoulli_rng(lambda);
real sd = (u == 1) ? S[n] : sigma;
Yrep[n] = normal_rng(mu[n], sd);
}
}
"
data <- list(
Y = final_df$estimate,
X = final_df$urb_trend,
S = final_df$std.error,
A = as.numeric(as.factor(final_df$SPECIES)),
B = as.numeric(as.factor(final_df$SITE_ID))
)
data$N <- length(data$Y)
data$N_A <- n_distinct(data$A)
data$N_B <- n_distinct(data$B)
mixture <- cmdstan_model(
write_stan_file(code=code)
)
M3 <- mixture$sample(
data = data,
parallel_chains = 4
)
M3$diagnostic_summary()
M3$summary()
draws <- as_draws_matrix(M3)
Yrep <- subset_draws(draws, "Yrep")
log_lik <- subset_draws(draws, "log_lik")
sim <- createDHARMa(
simulatedResponse = t(Yrep),
observedResponse = data$Y,
integerResponse = FALSE
)
plot(sim)
loo_object <- loo(log_lik, save_psis = TRUE)
ppc_loo_pit_qq(y = data$Y, yrep = Yrep, psis_object = loo_object$psis_object)
ppc_dens_overlay(y = data$Y, yrep = Yrep[1:100, ])
parameters::model_parameters(model, vcov = "HC3")$\endgroup$