Devil in the details: how can we avoid potential pitfalls of CATS regression when our data do not follow a Poisson distribution?

Modeling meta-community effect via offset

Beyond local trait selection, larger-scale effects can also influence species’ local abundance. High propagule pressure can increase the local abundance of species that are abundant in the surroundings. On the other hand, locally well-adapted species may be missing from the local community due to propagule limitation. The unique property of CATS is that it can measure the relative importance of local and meta-community-scale (dispersal) processes (Shipley, 2014).

If species survival, growth, and reproduction were independent of their traits, local abundances would differ from the meta-community-level means due to demographic stochasticity only (including stochasticity of dispersal). In this case, local abundances could be predicted well from mean abundances at the meta-community level, while using traits as predictors would not improve the fit. At the other extreme, when species abundances are independent of dispersal processes (i.e., no mass effect or propagule limitation) and fully determined by local processes, knowledge on the meta-community-level abundance would not be able to improve our ability to predict local abundances (Shipley, 2014). Therefore, the heart of CATS is the fitting of models with and without information on abundances in the meta-community level. Shipley (2014) called this information the “neutral prior.” However, Warton, Shipley & Hastie (2015) called attention to the term “prior” as being associated with Bayesian statistics, where it has a specific meaning. To avoid confusion, we will refer to it as “abundances expected a priori,” where “a priori” means “before knowing local conditions.”

The abundances expected a priori can be included into CATS regression models via offset terms. Warton, Shipley & Hastie (2015) suggested using the logarithm of relative abundances at the meta-community level in a Poisson regression to reproduce the original CATS model. This study discusses whether this suggestion is generally valid irrespective of the link function.

The aim of CATS models is to predict relative abundances of species, not to explain differences in total abundances among sites, which may be caused by differences in sampling intensity. Therefore, they always contain an intercept. The relative abundance predicted by a model containing only an intercept and offset should be equal to the relative abundances expected a priori (π): (1)${\displaystyle \frac{{\hat{y}}_i}{\sum _{i=1}^S{\hat{y}}_i}={\pi }_i}$ where ${\hat{y}}_i$ is the predicted abundance of species i, S is the number of species, and π_i is relative abundances expected a priori of species i. For most of the fitted models (but not for zero-inflated and two-stage models), ${\hat{y}}_i={\mu }_i$, where μ_i is the location parameter of the fitted distribution. Therefore, μ_i will be used instead of ${\hat{y}}_i$ where appropriate. If a canonical link is applied and the model contains an intercept term, the sum of predicted values is equal to the sum of observed values (y_tot). Thus, requirement Eq. (1) could be written in the following form: (2)${\displaystyle {\mu }_i={\pi }_i{y}_{tot}.}$ The GLM with an intercept and offset but no predictors can be written in the following general form: (3)${\displaystyle h\left({\mu }_i\right)={\beta }_0+O_i}$ where h() is the link function, β₀ is the intercept, and O_i is the offset for species i. Substituting Eq. (2) into Eq. (3), we obtain the following system of linear equations (note that the left side of equations can be replaced by numbers calculated from data on hand): (4)${\displaystyle h\left({\pi }_i{y}_{tot}\right)={\beta }_0+O_i.}$ This system contains S +1 variables (O ₁, O ₂, …, O_S and β₀), but only S equations, so it has no unique solution. It can be solved by choosing an arbitrary value for O ₁. Then, the other offsets can be calculated with the following formula: (5)${\displaystyle O_i=h\left({\pi }_i{y}_{tot}\right)-h\left({\pi }_1{y}_{tot}\right)+O_1.}$ In Poisson and negative binomial regression, the canonical link is the natural logarithm (Dobson, 2002). Thus, (6)${\displaystyle O_i=ln\left({\pi }_i{y}_{tot}\right)-ln\left({\pi }_1{y}_{tot}\right)+O_1=ln\left({\pi }_i\right)-ln\left({\pi }_1\right)+O_1.}$

In this case, it is appropriate to choose $O_1=ln\left({\pi }_1\right)$, which leads to the offset recommended by Warton, Shipley & Hastie (2015). For other link functions, the simplest choice is O ₁ = 0. However, this approach has two limitations. First, π_iy_tot has to be within the domain of the link function. For example, for binomial and beta distributions, where the canonical link is logit(x), offsets can be calculated only if π_iy_tot < 1; otherwise, $ln\left({\pi }_i{y}_{tot}/\left(1-{\pi }_i{y}_{tot}\right)\right)$ cannot be calculated. The second limitation is that a canonical link is not always the most appropriate link function, and sometimes, another link function has to be chosen. For example, for a Tweedie distribution with power parameter 1<p<2, the canonical link would be (Ohlsson & Johansson, 2006): (7)${\displaystyle h\left(\mu \right)=\frac{-1}{p-1}{\mu }^{-\left(p-1\right)}.}$

When the exact value of power parameter p is unknown, it can be estimated from data during model fitting, but in this case, the log-link is applied in R packages mgcv (Wood, 2017) and glmmTMB (Brooks et al., 2017). If a canonical link is not used, the sum of expected values may differ from the sum of observed values. But Eq. (1) remains true irrespective of the link function, and it can be converted to: (8)${\displaystyle \frac{{\mu }_i}{{\mu }_1}=\frac{{\pi }_i}{{\pi }_1}.}$ Combining Eqs. (3) and (8) results in: (9)${\displaystyle \frac{h^{-1}\left(\beta +O_i\right)}{h^{-1}\left(\beta +O_1\right)}=\frac{{\pi }_i}{{\pi }_1}}$ where $h^{-1}\left(x\right)$ is the inverse of link function.

Setting O₁ to an arbitrary value leads to a nonlinear equation system with S equations and S variables. Solving such a system is often a hard task. The situation would be much simpler if the following were true: (10)${\displaystyle h^{-1}\left(\beta +O_i\right)=h^{-1}\left(\beta \right)h^{-1}\left(O_i\right).}$ In this case, choosing $O_i=h\left({\pi }_i\right)$ satisfies condition Eq. (9). Since $h\left(x\right)=ln\left(x\right)$ and $h^{-1}\left(x\right)=e^x$ satisfy condition Eq. (10), it is reasonable to use log-link even if it is not the canonical link (for example, instead of logit in a binomial model). In medical statistics, binomial GLM with log-link is called relative risk regression, which is often recommended due to the easier interpretation of proportions than odds ratios (Marschner, 2015). If log-link is used with a binomial distribution, iteratively reweighted least squares (the standard method for fitting GLMs) may fail to converge to the maximum likelihood estimate (Marschner & Gillett, 2012). Therefore, alternative estimation procedures were developed and are implemented in the logbin R package (Donoghoe & Marschner, 2018). Log-link for beta regression is also available in the betareg R package (Cribari-Neto & Zeileis, 2010).

Using the mentioned distributions, we suppose that all data come from the same distribution, and only their parameters depend on species. In this case, the relative abundances at the meta-community level can be estimated by: (11)${\displaystyle {\pi }_i=\frac{{\overline{m}}_i}{\sum {\overline{m}}_i}}$ where $\overline{m_i}$ is the mean abundance of species i in plots representing the meta-community. Note that since only the ratio of relative abundances expected a priori are used, the offset could be simply $O_i=ln\left(\overline{m_i}\right)$ instead of $O_i=ln\left({\pi }_i\right)$.

If there are excessive zeros, it could be supposed that some of the zeros do not come from the distribution of “normal” abundances. Zero-inflated and two-part (hurdle) models are based on this assumption. In these models, there are two equations for two location-type parameters.

A zero-inflated model supposes that positive counts and some of the zeros comes from a Poisson or negative binomial process, while some zeros are “structural zeros” (i.e., species cannot occur there). The two parameters in this case are the probability of structural zeros (p) and the expected value of the Poisson or negative binomial process (μ) (Zuur et al., 2009). For estimating offsets, we must know the probability of structural zeros at the meta-community. Therefore, it seems that zero-inflated models have low practical relevance when a meta-community effect has to be modeled.

A two-part (hurdle) model fits two separate models: a binomial model for presence/absence data and a truncated Poisson or negative binomial model for positive abundances. In this case, the two parameters are the probability of presence (p) and mean of the Poisson or negative binomial distribution (μ), from which the fitted zero-truncated distribution is deduced (not the mean of the truncated distribution itself) (Zuur et al., 2009). For simplicity, let us imagine that we really fit a two-part model as two separate GLMs. The first GLM is a binomial model for binary data. Therefore, we estimate offsets using the standard procedure (but mean abundances at the meta-community level have to be calculated from binary data). In the second GLM, we fit a truncated Poisson or truncated negative binomial distribution for the non-zero abundances. Similar to Poisson regression, it is assumed that $ln\left({\mu }_i\right)$ is a linear combination of trait values and the offset, but the expected abundance is the following for a Poisson distribution: (12)${\displaystyle {\hat{y}}_i=\frac{{\mu }_i}{1-exp\left(-{\mu }_i\right)}.}$ For a negative binomial distribution, the expected abundance is: (13)${\displaystyle {\hat{y}}_i=\frac{{\mu }_i}{1-{\left(\frac{{\mu }_i+\theta }{\theta }\right)}^{-\theta }}}$ ${\hat{y}}_i$ is the expected or mean abundance of species i when it present, while μ_i is the expected or mean abundance when only structural zeros are excluded. For setting the offset, we need a priori expectation for the latter.

The meta-community level mean of species’ abundance when present (${\overline{m}}_i^{+}$) can be easily estimated. Assuming a Poisson process, the mean abundance excluding structural zeros (${\tilde{m}}_i$) can be estimated by solving the following nonlinear equation: (14)${\displaystyle {\overline{m}}_i^{+}=\frac{{\tilde{m}}_i}{1-exp\left(-{\tilde{m}}_i\right)}}$ If a negative binomial distribution is assumed, a similar approach can be applied if θ is known. Then, $ln\left({\tilde{m}}_i\right)$ could be used as an offset.

Relative importance of environmental selection and dispersal processes

The relative importance of local and meta-community-level processes can be calculated from variation explained by models containing only traits (as independent variables), only offset (calculated from meta-community-level abundances; see above), or both traits and offset (Shipley, 2014). Using the classic R² as a measure of explained variance is suitable only in OLS regression. Different generalizations of R² are suggested for GLMs (Cameron & Windmeijer, 1996; Cameron & Windmeijer, 1997; Menard, 2000; Nakagawa & Schielzeth, 2013). Shipley (2014) proposed using a generalization based on Kullback–Leibler divergence (Cameron & Windmeijer, 1997).

I will show below that formula Eq. (4) from Shipley (2014) is valid for only a Poisson-model, and different formulas have to be used for other distributions. A definition of Kullback–Leibler R² is not available for models with an offset; therefore, solution for this case by Shipley (2014) and its alternatives will be discussed. Increasing the number of predictors (traits) always improves the fit of the model (i.e., increases the Kullback–Leibler R²). Therefore, R² values of models with different numbers of predictors (traits) cannot be compared. Shipley (2014) proposed an “adjustment” procedure based on randomization of traits. Although this procedure is correct, it is time consuming for large datasets . Thus, an alternative deterministic adjustment is proposed.

R² for models without offset

Kullback–Leibler R² is a generalization of the classic R² used in ordinary least squares regression: (15)${\displaystyle R^2=1-\frac{\sum {\left(y_i-{\hat{y}}_i\right)}^2}{\sum {\left(y_i-\overline{y}\right)}^2}=\frac{\sum {\left(y_i-\overline{y}\right)}^2-\sum {\left(y_i-{\hat{y}}_i\right)}^2}{\sum {\left(y_i-\overline{y}\right)}^2}}$ where $\sum {\left(y_i-{\hat{y}}_i\right)}^2$ and $\sum {\left(y_i-\overline{y}\right)}^2$ are the squared Euclidean distances between observed values and predictions of models with and without predictors, respectively. Therefore, R² is a proportional decrease of distance between model prediction and observed values due to the inclusion of predictors in the model. For other distributions, the squared Euclidean distance can be replaced with Kullback–Leibler divergence with the same interpretation: (16)${\displaystyle R_{KL}^2=1-\frac{K\left(\mathbf{y};\mu \right)}{K\left(\mathbf{y};{\mu }^0\right)}=\frac{K\left(\mathbf{y};{\mu }^0\right)-K\left(\mathbf{y};\mu \right)}{K\left(\mathbf{y};{\mu }^0\right)}}$ where y is the vector of observed values, and μ and μ⁰ are vectors of values predicted by the evaluated and intercept-only models, respectively.

Table 1 shows that distributions that could be used in CATS regression belong to the exponential family, so their density functions can be written in the following general form (McCullagh & Nelder, 1999): (17)${\displaystyle f\left(y;\psi ,\varphi \right)=exp\left\{\frac{y\psi -b\left(\psi \right)}{a\left(\varphi \right)}-c\left(y,\varphi \right)\right\}}$ where ψ is the natural or canonical parameter, ϕ is the dispersal parameter, and a, b, and c are specific functions. The mean and variance of y are: (18)${\displaystyle \begin{array}{c}E\left(\mathbf{y}\right)=\mu =b^{\prime }\left(\psi \right)\hfill \\ Var\left(\mathbf{y}\right)=a\left(\varphi \right)b^{\prime \prime }\left(\psi \right)\hfill \end{array}}$ where b′ and b″ are the first and second derivatives of function b. Note that negative binomial and Tweedie distributions belong to this family only if parameters θ and p are known constants. Function $h\left(\mu \right)=\psi$ is called a canonical link function.

For members of the exponential family, Kullback–Leibler divergence can be calculated as the difference between the likelihood of a full model (i.e., a model where predicted and observed values are equal) and a fitted model (Cameron & Windmeijer, 1997): (19)${\displaystyle K\left(\mathbf{y};\mu \right)=2\left[l\left({\mu }^{full};\mathbf{y}\right)-l\left(\mu ;\mathbf{y}\right)\right].}$ Thus, for members of the exponential family, $R_{KL}^2$ could be deduced as a corrected version of likelihood ratio R²or McFadden R²: (20)${\displaystyle R_L^2=1-\frac{l\left(\mu ;\mathbf{y}\right)}{l\left({\mu }^0;\mathbf{y}\right)}.}$ A drawback of $R_L^2$ is that its maximum is not 1, but $1-\left[l\left({\mu }^{full};\mathbf{y}\right)/l\left({\mu }^0;\mathbf{y}\right)\right]$. Since its minimum is zero, $R_L^2$ can be rescaled to the interval of 0–1 by dividing it by its maximum, which results in $R_{KL}^2$: (21)${\displaystyle \frac{R_L^2}{\left(1-\frac{l\left({\mu }^{full};\mathbf{y}\right)}{l\left({\mu }^0;\mathbf{y}\right)}\right)}=\frac{l\left({\mu }^0;\mathbf{y}\right)-l\left(\mu ;\mathbf{y}\right)}{l\left({\mu }^0;\mathbf{y}\right)}/\frac{l\left({\mu }^0;\mathbf{y}\right)-l\left({\mu }^{full};\mathbf{y}\right)}{l\left({\mu }^0;\mathbf{y}\right)}=\frac{l\left({\mu }^0;\mathbf{y}\right)-l\left(\mu ;\mathbf{y}\right)}{l\left({\mu }^0;\mathbf{y}\right)-l\left({\mu }^{full};\mathbf{y}\right)}=R_{KL}^2.}$ Substituting Eqs. (19) into (16) results in: (22)${\displaystyle R_{KL}^2=1-\frac{l\left({\mu }^{full};\mathbf{y}\right)-l\left(\mu ;\mathbf{y}\right)}{l\left({\mu }^{full};\mathbf{y}\right)-l\left({\mu }^0;\mathbf{y}\right)}=\frac{l\left(\mu ;\mathbf{y}\right)-l\left({\mu }^0;\mathbf{y}\right)}{l\left({\mu }^{full};\mathbf{y}\right)-l\left({\mu }^0;\mathbf{y}\right)}.}$

McCullagh & Nelder (1999) called $2\left[l\left({\mu }^{full};\mathbf{y}\right)-l\left(\mu ;\mathbf{y}\right)\right]$ the scaled deviance (D^∗), so $R_{KL}^2$ can also be calculated from scaled deviances (D^∗) or deviances (D) of fitted and intercept-only models: (23)${\displaystyle R_{KL}^2=1-\frac{D^{\ast }\left(\mu ;\mathbf{y}\right)}{D^{\ast }\left({\mu }^0;\mathbf{y}\right)}=1-\frac{D\left(\mu ;\mathbf{y}\right)/a\left(\varphi \right)}{D\left({\mu }^0;\mathbf{y}\right)/a\left(\varphi \right)}=1-\frac{D\left(\mu ;\mathbf{y}\right)}{D\left({\mu }^0;\mathbf{y}\right)}.}$

If there is no offset, in a generalized linear model fitted by ML estimation with a canonical link, the expectations in an intercept-only model is equal to the mean of observed values: ${\mu }^0=\overline{y}$. Formulas for this case are listed in Appendix S1. Appendix S2 shows that formula Eq. (4) from Shipley (2014) is equivalent to the formula given for a Poisson regression in Appendix S1.

These formulas assume that the likelihood is a function of μ only, and if there are other parameters, their values are constants known a priori (i.e., not estimated during regression). If these parameters are estimated in regression, we can obtain different estimates for the evaluated and the intercept-only models. For a negative binomial distribution with unknown dispersion (θ), Cameron & Windmeijer (1996) suggested using a parameter estimated for an evaluated model when the likelihood of full and intercept-only models is calculated. Applying this approach, $R_{KL}^2$ may decrease when a regressor is added to the model due to changes in estimated θ. The same approach can be applied for the power parameter of a Tweedie distribution. GLM fitting programs usually give the log-likelihood of the fitted model and the deviance of fitted and intercept-only models.

R² for models with offset

In the previous section, μ was the prediction of model containing an intercept and predictors (traits), while μ⁰ was the prediction of an intercept-only model. A possible solution for how we should include the offset is to define μ⁰ as the prediction of a model without predictors (i.e., a model containing only an intercept and offset). At first glance, it seems to be a natural generalization of R² shown in the previous section. However, there is a drawback in this approach: the effect of meta-community-level processes cannot be directly measured since R² for models containing only an intercept and offset (but no traits) would always be zero. Shipley (2014) suggested an indirect measure of a pure meta-community effect: (24)${\displaystyle \frac{R^2\left(\text{traits};\text{offset}\right)-R^2\left(\text{traits}\right)}{1-\overline{R^2\left(\text{random traits}\right)}}.}$

The nominator is an adjustment for removing bias (see next section on adjustment), so now, we should focus on the denominator. This subtraction is based on the assumption that $R^2\left(\text{traits}\right)$ is the variation explained by traits, while $R^2\left(\text{traits};\text{offset}\right)$ is the variation explained by traits and offset (i.e., meta-community effect) together. However, this assumption is not satisfied when in calculation of R-squared μ⁰ is the prediction of a model containing only an intercept and offset. To understand why, we should recall the geometric interpretation of Kullback–Leibler R² shortly mentioned above: R² is the proportional decrease of distance between observed and predicted values (or proportional improvement of fit) due to the inclusion of predictors in the model. Therefore, in Eq. (24), both $R^2\left(\text{traits}\right)$ and $R^2\left(\text{traits};\text{offset}\right)$ are proportional improvements of fit due to the inclusion of traits, but they are proportional to different original distances of observed and predicted values. Therefore, their difference has no simple interpretation and does not measure the pure meta-community effect.

Instead of defining μ⁰ as a prediction of a model without a predictor, it could be defined as a prediction of an intercept-only model, even if an offset is applied. This definition allows us to calculate a meaningful R² for models with offset but no predictors as a direct measure of the meta-community effect. This definition may result in negative R² when including an offset increases the distance between observed and predicted values instead of decreasing it. A negative value is nonsense if R² is interpreted as explained variation, but it is meaningful if R² is interpreted as a proportional change in the distance between observed and predicted values. This geometric interpretation seems more useful in CATS regression, where it has a meaning that includes meta-community-level relative abundances decreases the goodness-of-fit.

Note that in R environment, to avoid negative R² values, the following formula is applied instead of Eq. (15): (25)${\displaystyle R^2=1-\frac{\sum {\left(y_i-{\hat{y}}_i\right)}^2}{\sum {\left(y_i-{\hat{y}}_i\right)}^2+\sum {\left({\hat{y}}_i-\overline{y}\right)}^2}.}$

If there is no offset, Eqs. (15) and (25) result in the same value, but they differ if offset is applied. A generalization of Eq. (25) could be: (26)${\displaystyle R^2=1-\frac{K\left(\mathbf{y};\mu \right)}{K\left(\mathbf{y};\mu \right)+K\left(\mu ;{\mu }^0\right)}.}$

Adjusted R² and partitioning of explained variation

As a goodness-of-fit measure, a drawback of R² (both in the classic form and its generalization) is that including an additional predictor in the model always increases R², even if the predictor is independent of the dependent variable. Therefore, it has a positive expected value instead of zero when a dependent variable is not related to the predictors. To remove this bias, Fisher (1925) suggested using the following adjustment: (27)${\displaystyle R_{adj}^2=1-\left(1-R^2\right)\frac{n-1}{n-k-1}}$ where n is the number of data points, and k is the number of predictors.

Unfortunately, this adjustment is valid for only for “classic” R² of ordinary least squares regression. Since Shipley (2014) has not found a similar solution for Kullback–Leibler R², he proposed a procedure based on reshuffling trait values to remove the bias. Although the suggested procedure is correct, it has no unique result, in contrast to the correction using a closed form. Ricci (2010) has shown that there is a simple general adjustment for $R_{KL}^2$ when the dependent variable follows a distribution belonging to the exponential family. Let us write the formula for $R_{KL}^2$ using scaled deviances: (28)${\displaystyle R_{KL}^2=\frac{D\left({\mu }^0;\mathbf{y}\right)-D\left(\mu ;\mathbf{y}\right)}{D\left({\mu }^0;\mathbf{y}\right)}.}$

Recall that $R_{KL}^2$ is the proportional improvement of model fit. Let us focus on the denominator, which is the absolute improvement. In a bias-corrected version, the absolute improvement should be zero when predictors have no effect on the dependent variable. Under this condition, for members of the exponential family, $D\left({\mu }^0;\mathbf{y}\right)-D\left(\mu ;\mathbf{y}\right)$ approximately follows a Chi-square distribution with degree of freedom equals to the number of predictors (k). Since the expected value of the Chi-square distribution is its degree of freedom, the following is an approximately bias-free goodness-of-fit measure: (29)${\displaystyle R_{KL,adj}^2=\frac{D\left({\mu }^0;\mathbf{y}\right)-D\left(\mu ;\mathbf{y}\right)-k}{D\left({\mu }^0;\mathbf{y}\right)}.}$

Appendix S3 shows that Eq. (27) is a special case of Eq. (29) for a Gaussian distribution with dispersion estimated from the data. Note that Ricci (2010) applied an alternative derivation of adjustment by generalization of the shrinkage factor and used the deviance instead of the scaled deviance in the formulas.

For partitioning variation, we should fit models containing both traits and offset and models with only traits and only offset. Let us denote the corresponding adjusted R² values by R²(trait; offset), R²(trait), and R²(offset), respectively. R²(trait; offset) measures the whole variation explained by studied traits and relative abundances at the meta-community level. The pure trait effect (i.e., variation explained only by traits) is R²(trait; offset)- R²(offset), while the pure meta-community effect is R²(trait; offset)- R²(trait). Variation that can be explained by both traits and the meta-community effect is R²(trait)+R²(offset)- R²(trait; offset). These formulas are analogous to partitioning of the variation of community composition into environmental and spatial components (Borcard, Legendre & Drapeau, 1992; Peres-Neto et al., 2006).

Example 1: fitting Poisson model to over-dispersed counts

The first example illustrates the problems that arise when an inappropriate model is fitted. The type of abundance often clearly determines the type of model to be fitted. However, when abundance is measured by a number of individuals, a Poisson or negative binomial model should be fitted depending on whether there is a random or aggregated spatial pattern of individuals. If the spatial pattern is random, the number of individuals will follow a Poisson distribution, where the variance is equal to the mean. An aggregated spatial pattern leads to over-dispersed counts (i.e., the variance is higher than mean), which can be modeled by a negative binomial distribution.

The presence of over-dispersion can be checked by comparing the Akaike Information Criteria of Poisson and negative binomial models or using diagnostic plots of residuals (Fig. 1). Dunn-Smyth (or randomized quantile) residuals (Dunn & Smyth, 1996) are especially useful for this purpose because if the model’s assumptions (specified distribution and log-linear relationship) are satisfied, they follow a standard normal distribution (Warton, Shipley & Hastie, 2015; Feng, Li & Sadeghpour, 2020).

This example uses simulated data. Abundances (y) of 20 species in a plot were simulated. Abundances follow a negative binomial distribution with a mean that has a log-linear relation to values of a trait. The dispersion parameter is constant, and trait values follow a normal distribution: (30)${\displaystyle \begin{array}{c}{y}_i\sim \text{NegBin}\left({\mu }_i=exp\left(0.5\ast x_i\right),\theta =1\right)\hfill \\ x_i\sim N\left(m=10,\sigma =3\right).\hfill \end{array}}$

The simulation was repeated 50 times, and Poisson and negative binomial models were fitted to each simulated plots separately. Figure 1 shows a typical diagnostic plot, which has a strong nonlinearity of the QQ plot indicating that the distributional assumption of a Poisson regression is not satisfied. The diagnostic plots of a negative binomial model do not indicate any problem.

The estimated slopes are dispersed around the real value (0.5) in both Poisson and negative binomial models (Fig. 2). The variation among estimates was higher in the Poisson distribution, while the confidence intervals were narrower due to the neglect of over-dispersion. These two facts together may result in over-interpretation of differences in the strength of selection among plots.

Example 2: choosing appropriate offset

The second example illustrates the importance of choosing an appropriate offset, without which fitted relative abundances may differ considerably from the relative abundances in the meta-community, even if no traits are included in the model. The example uses the dataset of Raevel, Violle & Munoz (2012), which contains the abundance (number of individuals) of 97 species at 52 sites. Data were transformed into a presence/absence scale, and then meta-community-level abundances were measured as the number of occurrences at the 52 sites.

For modeling of the presence/absence data, a binomial distribution has to be applied. The canonical link for this distribution is the logit link. The proposed way of calculating the offset for the logit link in this paper cannot be applied because the product of the number of species in the plot and relative abundance at the meta-community scale was larger than one in 39 species-site combinations (and it excludes 22 of 52 sites). Therefore, the logarithm of relative abundances in the meta-community was used.

Two link functions were tested: the canonical (logit) link and log link, as suggested in this paper. The latter was fitted using the logbin package (Donoghoe & Marschner, 2018).

A model containing only intercept and offset terms was fitted, so the predicted relative abundances in plots should equal to the relative abundance in the meta-community. This requirement was satisfied in the model using log link (not shown). However, when the canonical link was used, there is a non-linear relationship between two vectors of relative abundances (Fig. 3).

Example 3: comparing formulas for estimation of variation components

The aim of the next example is to compare variation components estimated by the method of Shipley (2014) and the new method proposed in this paper. Simulated data were used, where trends of variation components were predictable. Data were generated using the following model: (31)${\displaystyle \begin{array}{c}{y}_i\sim \text{Poisson}\left({\lambda }_i\right)\hfill \\ log{\lambda }_i=a+log{\pi }_i+s\ast t_i\hfill \end{array}}$ where π_i is the meta-community-level relative abundances, $\left|s\right|$ is the strength of selection, and t_i is the trait value. To remove changes due to total community size, the intercept (a) was set to: (32)${\displaystyle a=logA-log\sum _i{\pi }_i{e}^{s{t}_i}.}$

Thus, ∑λ_i = A for any value of s.

The species pool consists of 50 species, and their traits follow a standard normal distribution. The expected community size (A) was set to 2500, and the strength of selection (s) changes from 0 to 3. Pure selection and pure meta-community effects were calculated for each simulated community separately using formulas from Shipley (2014) and the method proposed in this paper.

There is a good agreement between variation components calculated by the two ways (Fig. 4). As expected, the pure meta-community effect decreases, while the pure selection effect increases with increasing strength of selection, and the later was about zero at s = 0 (Fig. 5). The joint effect and unexplained variation do not change considerably with changing strength of selection.