Modeling dragonfly population data with a Bayesian bivariate geometric mixed-effects model

Abstract

We develop a generalized linear mixed model (GLMM) for bivariate count responses for statistically analyzing dragonfly population data from the Northern Netherlands. The populations of the threatened dragonfly species Aeshna viridis were counted in the years 2015-2018 at 17 different locations (ponds and ditches). Two different widely applied population size measures were used to quantify the population sizes, namely the number of found exoskeletons ('exuviae') and the number of spotted egg-laying females were counted. Since both measures (responses) led to many zero counts but also feature very large counts, our GLMM model builds on a zero-inflated bivariate geometric (ZIBGe) distribution, for which we show that it can be easily parameterized in terms of a correlation parameter and its two marginal medians. We model the medians with linear combinations of fixed (environmental covariates) and random (location-specific intercepts) effects. Modeling the medians yields a decreased sensitivity to overly large counts; in particular, in light of growing marginal zero inflation rates. Because of the relatively small sample size (n = 114) we follow a Bayesian modeling approach and use Metropolis-Hastings Markov Chain Monte Carlo (MCMC) simulations for generating posterior samples.

Keywords: Aeshna viridis; Bayesian modeling; bivariate geometric distribution; count data; generalized linear model (GLM); mixed effects.

Figure 1.

Graphical representation of the hierarchical Bayesian model specified by (10)–(17). Arrows indicated dependencies, with distributional relationships indicated by labeled squares (‘Norm’, ‘Wis’, and ‘Gam’ indicate the normal, Wishart, and Gamma distributions, respectively). Observed quantities are indicated in gray and hyperparameters are given as plain nodes. The top-right panel contains the hierarchical distribution that makes up the $\mathrm{M}\mathrm{G}\mathrm{H}-t(\lambda ,v,d=2)$ prior for $\mathbf{\Sigma }$.

Figure 2.

Scatter plots of simulated data (cf. Section 4). For three $(M,N)$ combinations, we overlaid the scatter plots for $\theta =0.1$ (circles) and $\theta =0.9$ (dots). It can be seen that the dots ( $\theta =0.9$) are more strongly correlated than the circles ( $\theta =0.1$). To visualize overlaid points as clusters, we added jitter and reduced the opacity of points closer to $(0,0)$. The theoretical and the empirical correlations are provided in Table 1.

Figure 3.

Aeshna viridis count distribution. Left: Scatter plot of the dragonfly count data in log-log scale; the points have been slightly jittered to reveal clusters of observations. The sample Pearson correlation coefficient is $\hat{\rho }=0.299$. Right: Overlaid histograms showing the marginal count distributions along with fitted univariate zero-inflated geometric distributions in semi-log scale. The fitted zero-inflation parameter is $\pi =0.233$ ( $\pi =0.137$) and the fitted continuitized median is M = 5.68 (M = 2.74) for the exuviae (egg-laying female) counts.

Figure 4.

Graphical comparison of the covariate effects on the two population measures. Left: Scatter plot of the fractions of positive posterior samples for all covariates (egg-laying females vs. exuviae). A covariate effect is consistent across both population measures if the point is close to the diagonal. The dashed lines separate the positive and negative effects. Right: A grouped bar chart of the positive fractions. Values close to 1 (0) indicate significant positive (negative) effects. Covariate effects are consistent when the two bars point in the same direction and have approximately the same height.

Figure 5.

Location effects on the dragonfly population sizes. For each of the m = 17 locations (x-axis), the figure shows 95% confidence intervals for the random intercept parameters for exuviae (in dark gray) and egg-laying females (in light gray). The triangles and dots mark the posterior medians.

Figure 6.

Diagnostic plots. Top: Q-Q plots of response-specific randomized quantile residual (RQR) quantiles w.r.t. $\mathrm{U}\mathrm{n}\mathrm{i}\mathrm{f}\left(\right[0,1\left]\right)$ quantiles. Bottom: Histograms of response-specific Pearson residuals, where the observations' expectations and variances have been approximated using the respective model's posterior sample. Left: Model diagnostics for the proposed ZIBGe-GLMM model. Right: Model diagnostics for the model from [22], which is very akin to our ZIBGe-GLMM model except that it uses a bivariate Poisson distribution to build the GLMM likelihood.