Extending null scenarios with Faddy distributions in a probabilistic randomization protocol for presence-absence data

Abstract

Navarro and Manly (Popul Ecol 51:505–512, 2009) (NM) have proposed a randomization protocol for null model analysis of species occurrences at discrete locations based on probability distributions and generalized linear models. In the NM method, presences-absences are governed by independent Bernoulli random variables. In addition, a non-observable non-negative random variable (“quasi-abundance”) from either Poisson, Binomial or Negative Binomial distributions are log-linearly related to the qualitative effects of species and location. By connecting the probability of occurrence of each species on each location and the quasi-abundance distributions, one generalized linear model for the observed presences-absences is selected by profile deviance, and the resulting fitted probabilities of the null model with minimum deviance is used to generate random matrices via parametric bootstrap. This work contributes with a unified theoretical formulation of the NM method, based on Faddy distributions, to allow general distributions of over-dispersed and under-dispersed discrete random variables. For a subset of the Faddy models, the log concave property of the inverse link function guarantees convergence to a global minimum deviance thus providing unique estimates for the linear parameters of the models. The method is illustrated using presence-absence data of island lizard communities. Interpretations of this combined GLM-parametric bootstrap protocol are discussed, highlighting the way fitted probabilities under the chosen null model are related to the row and column totals of the observed table. Additional properties of the probabilistic NM protocol, with possible avenues of future research, are also discussed.

Appendices

Appendix

The basic null model is given in Eq. (1) in terms of the expected value of the unobserved quasi-abundance $\mu_{kl} = E(N_{kl} )$, k = 1,.., t; l = 1,…,s:

$$ \mu_{kl} = \exp (\gamma + \alpha_{k} + \beta_{l} ),k = {1},..,t;l = {1}, \ldots ,s $$

where ${\alpha }_{1}={\beta }_{1}=0$. These unobserved means are connected to probabilities of observed Bernoulli-distributed outcomes ${Y}_{kl}$, $E({Y}_{kl})={\pi }_{kl}$, using different generalized linear models. One collection of models defines link functions in terms of the diffusion approximation to the mean quasi-abundance for Faddy distributions with parameters b and c, as described in Sects. 3.2 (Eq. 20) and Sect. 3.5.1. Another generalized linear model based on Faddy distributions is given by Eq. 21. This “Appendix” focuses mostly on models given by Eq. 20 (indicated here as Eq. 30):

$$ \pi_{kl} = \left\{ {\begin{array}{ll} {1 - \exp \left\{ {\frac{b}{1 - c}\left[ {1 - \left( {1 + \frac{{\mu_{kl} }}{b}} \right)^{1 - c} } \right]} \right\},} & {c < 1} \\ {1 - \left( {\frac{b}{{b + \mu_{kl} }}} \right)^{b} ,} & {c = 1} \\ \end{array} } \right. $$

(30)

where b > 0. The estimation of each generalized linear model is performed by maximum likelihood, and a model-selection process using profile deviance is conducted for varying values of b and c, the deviance function given by Eq. 25.

The following lemma applies to any model under the NM protocol:

Lemma

If $P(Y_{kl} = 0) = P(N_{kl} = 0);\quad P(Y_{kl} = 1) = P(N_{kl} \ge 1)$, then, for all probability mass functions on the non-negative integers $N_{kl}$:

$${\pi}_{\text{kl}}\left({\mu}_{\text{kl}}\right){\le}{\mu}_{\text{kl}}.$$

(31)

Proof

$$ \mu _{{kl}} = E\left( {N_{{kl}} } \right) = \sum\limits_{{n_{{kl}} = 0}}^{\infty } {n_{{kl}} } P\left( {N_{{kl}} {\text{ }} = {\text{ }}n_{{kl}} } \right) \ge \sum\limits_{{n_{{kl}} = 1}}^{\infty } P \left( {N_{{kl}} {\text{ }} = {\text{ }}n_{{kl}} } \right) = P\left( {N_{{kl}} \ge 1} \right) = P\left( {Y_{{kl}} {\text{ }} = {\text{ }}1} \right){\text{ }} = {\text{ }}\pi _{{kl}} . $$

The equality in ${\pi}_{\text{kl}}\left({\mu}_{\text{kl}}\right){\le}{\mu}_{\text{kl}}$ is achieved when ${\text{P}}\left({\text{N}}_{\text{kl}}\text{ = 0}\right)\text{ = 1}-{ \pi}_{\text{kl}}$, ${\text{P}}\left({\text{N}}_{\text{kl}}\text{ = 1}\right)\text{=}{ \pi}_{\text{kl}}$, and ${\text{P}}\left({\text{N}}_{\text{kl}}\text{ > 1}\right)\text{ = 0}$.

The inequality (31) is essential for the appropriate application of the “Faddy inspired” models (30) and will constrain the set of competing models in the profile deviance.

Theorem

Under the NM protocol, the generalized linear models connecting the probabilities of the observed Bernoulli-distributed outcomes $Y_{kl}$ with parameter $E(Y_{kl} ) = \pi_{kl}$, and the quasi-abundance whose link functions are defined in terms of the diffusion approximation to the mean for Faddy distributions with parameters b and c, are restricted to the set

R = {b > 0, c ≤ 1, b ≥ ‒c}.

Proof

For (31) to be true, the second derivative of ${\pi}_{\text{kl}}\left({\mu}_{\text{kl}}\right)$ evaluated at ${\mu}_{\text{kl}}$= 0 must be less than or equal to 0. And it turns out for Eq. (30) this is also a sufficient condition. Evaluating the derivative in Eq. (30) we have:

$${\left.\frac{d\pi ({\mu }_{kl})}{d{\mu }_{kl}}\right|}_{{\mu }_{kl}=0}={\left.\mathrm{exp}\left(\left(b\left(1-{\left({1+\mu }_{kl}/b\right)}^{1-c}\right)\right)/(1-c)\right){\left({1+\mu }_{kl}/b\right)}^{-c}\right|}_{{\mu }_{kl}=0}=1$$

In addition, ${\pi}_{\text{kl}}\left({\mu}_{\text{kl}}\right){\le}{\mu}_{\text{kl}}$ implies that $\frac{{d}^{2}\pi ({\mu }_{kl})}{{d}^{2}{\mu }_{kl}} \le 0$ and

$${\left.\frac{{d}^{2}\pi ({\mu }_{kl})}{d{{\mu }_{kl}}^{2}}\right|}_{{\mu }_{kl}=0}={\left.-\frac{\mathrm{exp}\left(\left(b\left(1-{\left({1+\mu }_{kl}/b\right)}^{1-c}\right)\right)/(1-c)\right){\left(\frac{b+{\mu }_{kl}}{b}\right)}^{-2c}\left(c{\left(\frac{b+{\mu }_{kl}}{b}\right)}^{c}+b+{\mu }_{kl}\right)}{b+{\mu }_{kl}}\right|}_{{\mu }_{kl}=0}\le 0$$

Note that the last term in the numerator is the only term that can be negative. Therefore, the last inequality implies that

$$-(b+c)\le 0$$

For $b>0$, and $0\le c\le 1$, this last expression is always valid. However, the validity of this expression forces the inequality:

$$b\ge - c.$$

In particular, for c < 0, we must have: $b>-c$.

The case b = ‒c in Eq. (30)

On the lower boundary of the valid region R, b = ‒c, c < 0, the Faddy model can be written as:

$$\pi \left({\mu }_{kl}\right)=1-\mathrm{exp}\left(1-\left(b/\left(1+b\right)\right){\left(1+{\mu }_{kl}/b\right)}^{1+b}\right)$$

In this case, the probability mass function of ${\text{N}}_{\text{kl}}$ is under-dispersed. The profile deviance of the lizard data suggested that the minimum deviance is possibly attained at this boundary as b goes to infinity. Thus, we were also interested in the “Faddy model” given by:

$$\pi \left({\mu }_{kl}\right)=\underset{b\to \infty }{\mathrm{lim}}\left(1-\mathrm{exp}\left(1-\left(b/\left(1+b\right)\right){\left(1+{\mu }_{kl}/b\right)}^{1+b}\right)\right)=1-\mathrm{exp}\left(1-\mathit{exp}\left({\mu }_{kl}\right)\right)$$

(32)

Putting model (32) to compete in the profile deviance for the lizard data, turns out to be the optimal model among all models of the form (30).

Log-concavity for models given by Eq. (30)

The proposed Faddy models (30) fitted by maximum likelihood have good properties in the permissible region R. For every data set and for every starting value in this region, the minimization procedures will converge to a global minimum. This statement can be proved by noticing that the “inverse links functions” defined by expression (30) for b and c in R, all have the property of log-concavity. Therefore, as far as maximum likelihood estimation is concerned, a subset of the link functions induced by Faddy models have the same good properties as the usual link functions for binary data (except for the exponential family properties). For parameters outside this region, convergence to the global minimum may or not may be assured; it will depend on the observed binary data and on the initial starting values. The model fitting process was also tried for the lizard data using Faddy models with b < ‒c and values of the deviances were computed but they were highlighted so that their interpretation should be made with caution. To prove the claim that the deviance function

$$D\left({\varvec{\Pi}}\text{,}{\text{y}}\right)\text{=}-2\left[\sum_{{y}_{kl}=0}\mathrm{log}{(1-\pi }_{kl})+\sum_{{y}_{kl}=1}\mathrm{log}{(\pi }_{kl})\right]$$

is concave when the inverse links functions ${\pi }_{kl}$ are defined by expression (30), we make use of a result given in Tang and Ye (2020) that the sum of concave functions is concave. Thus, a sufficient condition for the existence of a unique MLE estimator of the log-linear parameters is that $\mathrm{log}({\pi }_{kl}\left({\mu }_{kl}\right))$ and $\mathrm{log}(1-{\pi }_{kl}\left({\mu }_{kl}\right))$ be concave. This is equivalent to show that

$\frac{{\partial }^{2}\mathrm{log}({\pi }_{kl}\left({\mu }_{kl}\right))}{\partial {{\mu }_{kl}}^{2}}\le 0$ and $\frac{{\partial }^{2}\mathrm{log}({1-\pi }_{kl}\left({\mu }_{kl}\right))}{\partial {{\mu }_{kl}}^{2}}\le 0$

It can be shown that, for c < 1 and b > 0:

$$\frac{{\partial }^{2}\mathrm{log}({1-\pi }_{kl}\left({\mu }_{kl}\right))}{\partial {{\mu }_{kl}}^{2}}=\frac{{e}^{{\mu }_{kl}}{\left(\frac{b+{e}^{{\mu }_{kl}}}{b}\right)}^{-c-1}\left(\left(c-1\right){e}^{{\mu }_{kl}}-b\right)}{b}\le 0$$

The case c = 1, b > 0 is analogous:

$$\frac{{\partial }^{2}\mathrm{log}({1-\pi }_{kl}\left({\mu }_{kl}\right))}{\partial {{\mu }_{kl}}^{2}}=-\frac{b(b+1){\left(\frac{b}{b+{\mu }_{kl}}\right)}^{b}}{{(b+{\mu }_{kl})}^{2}}\le 0$$

Thus, $\mathrm{log}(1-{\pi }_{kl}\left({\mu }_{kl}\right))$ is log-concave. Showing log-concavity of $\mathrm{log}({\pi }_{kl}\left({\mu }_{kl}\right))$ by looking at the sign of $\frac{{\partial }^{2}\mathrm{log}({\pi }_{kl}\left({\mu }_{kl}\right))}{\partial {{\mu }_{kl}}^{2}}$ is more complex but the resulting expression can be inspected in a similar way as the previous case.

$$ \frac{{\partial ^{2} {\text{log}}(\pi _{{kl}} \left( {\mu _{{kl}} } \right))}}{{\partial \mu _{{kl}} ^{2} }} = \frac{{e^{{\frac{{b\left( { - 1 + \left( {\frac{{b + e\mu _{{kl}} }}{b}} \right)^{{1 - c}} } \right)}}{{ - 1 + c}}}} + \mu _{{kl}} \left( {\frac{{b + e^{{\mu _{{kl}} }} }}{b}} \right)^{{ - 1 - 2c}} \left( { - e^{{\mu _{{kl}} }} \left( {b + e^{{\mu _{{kl}} }} } \right) + \left( { - 1 + e^{{\frac{{b\left( { - 1 + \left( {\frac{{b + e\mu _{{kl}} }}{b}} \right)^{{1 - c}} } \right)}}{{ - 1 + c}}}} } \right)\left( {\frac{{b + e^{{\mu _{{kl}} }} }}{b}} \right)^{c} ( - b + ( - 1 + c)e^{{\mu _{{{\text{kl}}}} }} } \right)}}{{b\left( { - 1 + e^{{\frac{{b\left( { - 1 + \left( {\frac{{b + e\mu _{{kl}} }}{b}} \right)^{{1 - c}} } \right)}}{{ - 1 + c}}}} } \right)}} $$

With the exception of the last term in the numerator, all terms are positive for c < 1 and b > 0. Thus, we need to demonstrate that:

$$ f = \left( { - e^{{\mu _{{kl}} }} \left( {b + e^{{\mu _{{kl}} }} } \right) + \left( { - 1 + e^{{\frac{{b\left( { - 1 + \left( {\frac{{b + e\mu _{{kl}} }}{b}} \right)^{{1 - c}} } \right)}}{{ - 1 + c}}}} } \right)\left( {\frac{{b + e^{{\mu _{{kl}} }} }}{b}} \right)^{c} ( - b + ( - 1 + c)e^{{\mu _{{kl}} }} } \right) \le 0$$

One can show that the inequality does not hold for R₁ = {c < 0, 0 < b < ‒c} (Region 1). At first glance, it appears that the inequality holds for R₂ = {c < 0, b ≥ ‒c} (Region 2), and also for R₃ = {0 ≤ c < 1, b > 0} (Region 3). Therefore, numerically maximizing the function f over the Region

$${\text{R}}_{2}{ \cup } \, {\text{R}}_{3}=\mathrm{R}$$

always yields a non-positive value. To show that the expression fails over R₁, we can demonstrate that it fails in the neighbourhood of ${\mu }_{kl}=0$. It can be shown that the evaluation of f at 0 in the limit and its first two derivatives yields:

$$\underset{{\mu }_{\mathit{kl}}\to 0}{\mathrm{lim}}f=0$$

$$\underset{{\mu }_{\mathit{kl}}\to 0}{\mathrm{lim}}\frac{\partial f}{\partial {\mu }_{kl}}=0$$

$$\underset{{\mu }_{\mathit{kl}}\to 0}{\mathrm{lim}}\frac{{\partial }^{2}f}{\partial {{\mu }_{kl}}^{2}}=-c-b>0$$

By continuity of f and its first two derivatives, there exists $\varepsilon >0$, such that for $\varepsilon >{\mu }_{kl}>0$ and $-c>b$ we have ${f}^{{^{\prime}}{^{\prime}}}\left({\mu }_{kl}\right)>0$, then ${f}^{^{\prime}}\left({\mu }_{kl}\right)>0$ and, consequently, $f\left({\mu }_{kl}\right)>0.$ However, a necessary condition for log-concavity of $f$ is that $f\left({\mu }_{kl}\right)\le 0$ for all ${\mu }_{kl}>0.$ Thus, f fails to be log-concave in R₁ = {c < 0, 0 < b < ‒c}. By the same logic we can see that, for Regions 2 and 3, $f\left({\mu }_{kl}\right)\le 0$ in the neighbourhood of $0$. Although this is not an analytic proof that the function f(${\mu }_{kl},b,c)$ over Regions 2 and 3 is always negative, this “proof” by numerical maximization seems sufficient.

Remark

The most extreme function that is not in the “Faddy family” (30) but satisfies (31) is:

$${\pi }^{*}\left({\mu }_{kl}\right)=\mathrm{min}\left({\mu }_{kl},1\right).$$

We explored this model and we found that it yields a higher likelihood (smaller deviance) than any other model in the permissible region. However, in a “real-world” setting one would not be able to justify this kind of model.