I am analyzing ecological data in R, where I aim to understand the impact of urbanization on species trends. My response variable is the coefficient of species trends (estimate), and my main predictor is the coefficient of urbanization trends (urb_trend).

enter image description here

I'm using the glmmTMB package to fit a mixed effects model, accounting for random effects by SPECIES and nested within Country.Name / SITE_ID. My model looks like this:

mod_h1 <- glmmTMB(estimate ~ urb_trend + (1|SPECIES) + (1|Country.Name/SITE_ID),
                  data = final_df,
                  weights = inverse_variance_weights,
                  family = gaussian)

I've encountered issues with heteroscedasticity, as evident from diagnostic plots generated by the DHARMa package:

sim_res <- simulateResiduals(fittedModel = mod_h1)

enter image description here

Despite trying several approaches to mitigate this—like transforming the response variable, adding various predictors to the model (including polynomial terms), adjusting weights, and experimenting with the dispformula—the issue persists.

(A) What strategies within glmmTMB or other packages might better address heteroscedasticity in the context of ecological data analysis?

(B) How can I more effectively determine which variables might be influencing the variability of errors in my model?

(C) Have I overlooked an alternative approach, transformation, or model specification that could help resolve this issue?

  • $\begingroup$ Your question is quite open-ended and it might be challenging to get any specific advice without access to the data. In any case, the question seems better suited to Cross Validated. $\endgroup$
    – dipetkov
    Apr 15 at 8:13
  • 1
    $\begingroup$ I've shared the dataset (final_df.csv) and the R script (Hypotheses.R) that contains the code and detailed information about the data and the analysis approach on GitHub. You can access them here: github.com/Paucom9/URBAN_TRENDS. The data is public, and I welcome any advice or insights. However, please do not publish, share, or use this data without permission. $\endgroup$ Apr 15 at 9:33
  • $\begingroup$ If it's really difficult to substantially improve problems with heteroscedasticity, you could also try to calculate HC (robust) standard errors, to account for the increased uncertainty due to heteroscedasticity. E.g., parameters::model_parameters(model, vcov = "HC3") $\endgroup$
    – Daniel
    Apr 15 at 11:49
  • 1
    $\begingroup$ It may be that your response variable is bounded between -1 and 1, and so the assumption that error variance is constant probably doesn't hold (Gaussian errors can range between -Inf and Inf). Maybe find an error family that is more appropriate for your data. Some error families (e.g., betabinomials) allow you to define a dispersion family to capture extra variation - in the case of betabinomials, this is extra-binomial variation. $\endgroup$
    – André.B
    Apr 15 at 16:19
  • $\begingroup$ Thank you for yout suggestion @André.B , but the response variable consists of coefficients from linear models representing species trends, and these are not inherently bounded between -1 and 1. Instead, they can theoretically assume any value, making them unbounded and continuous. $\endgroup$ Apr 15 at 17:03

1 Answer 1


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.

enter image description here

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

enter image description here

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.

enter image description here enter image description here

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:

  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(
      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(
      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(

M3 <- mixture$sample(
  data = data,
  parallel_chains = 4


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

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, ])
  • 1
    $\begingroup$ PS: I'm not sure weights = log(1/(std.error)) is the right way to incorporate the estimated std. errors with {glmmTMB}. I think the weights should be the inverse of the variance, so weights = 1/(std.error^2). It may be worth it to double check that the weights are used as intended. $\endgroup$
    – dipetkov
    May 6 at 10:20
  • $\begingroup$ Thank you for you answer @dipetkov $\endgroup$ May 13 at 8:01
  • $\begingroup$ We calculate estimates (i.e., butterfly population trends) and standard errors using the following GLS model: gls_model <- gls(log_sindex ~ YEAR, data = data_subset, correlation = corAR1(form = ~ YEAR)). In this model, log_sindex represents the logarithm of an annual abundance index. The code is available in the Butterfly population trends.R script. Additionally, I've uploaded a matrix of normalized distances (0-1 values) between sites (distance_sites_normalized_matrix.csv) and phylogenetic data for European butterfly species (EUROPEAN_BUTTERFLIES_FULLMCC_DROPTIPED.nwk). $\endgroup$ May 13 at 8:42
  • $\begingroup$ Phylogenetic data can be written as Tree <- read.tree(file_path, "EUROPEAN_BUTTERFLIES_FULLMCC_DROPTIPED.nwk")) Tree$tip.label <- gsub("_", " ", Tree$tip.label) # Replace underscores with spaces in the tip labels cor_phylo <- corBrownian(phy = Tree, form = ~ SPECIES). How can I account in the model for both spatial and phylogenetic autocorrelation? All species are butterflies but it is true that not all are found accross all sites. However this can be accounted by including SPECIES as random factor, right? @dipetkov $\endgroup$ May 13 at 8:48
  • $\begingroup$ Climate variability can be accounted using final_df$genzname, which classifies site climate regions following Metzger et al.2013. genzname levels include: Cold and wet (0.12% of sites); Cold and mesic (28.34%); Cool temperate and dry (10.1%); Cool temperate and xeric (0.52%); Cool temperate and moist (22.3%); Warm temperate and mesic (34.7%); Warm temperate and xeric (3.9%). We plan to analyze species variability by using species traits, detailed in lines 31-40 of Hypotheses.R. @dipetkov $\endgroup$ May 13 at 8:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.