Uncertainty analysis of species distribution models

The maximum entropy model, a commonly used species distribution model (SDM) normally combines observations of the species occurrence with environmental information to predict the geographic distributions of animal or plant species. However, it only produces point estimates for the probability of species existence. To understand the uncertainty of the point estimates, we analytically derived the variance of the outputs of the maximum entropy model from the variance of the input. We applied the analytic method to obtain the standard deviation of dengue importation probability and Aedes aegypti suitability. Dengue occurrence data and Aedes aegypti mosquito abundance data, combined with demographic and environmental data, were applied to obtain point estimates and the corresponding variance. To address the issue of not having the true distributions for comparison, we compared and contrasted the performance of the analytical expression with the bootstrap method and Poisson point process model which proved of equivalence of maximum entropy model with the assumption of independent point locations. Both Dengue importation probability and Aedes aegypti mosquito suitability examples show that the methods generate comparatively the same results and the analytic method we introduced is dramatically faster than the bootstrap method and directly apply to maximum entropy model.

The counts O = {o 1 , o 2 , ...o n } follow a multinomial distribution, whose true parameter P = {p 1 , p 2 , ..., p n } is unknown. Let m = n i=1 o i be the total number of observations. The maximum likelihood estimator of P = {p 1 , p 2 , ..., p n }, is and it follows a normal distribution [1], From the analysis in Kapur et al. [2], one can show that the maximum likelihood estimates of P from maximum entropy model (1) are achieved when For brevity, let a j =Ê(f j (X)). Based on equation (2), the vector of a j can be expressed as Let g(A) denote the maximum entropy optimization, model (1), as a function from R k to R n . In other words, the function takes as input the vector A, specifying right hand sides of the equality constraintsÊ(f j (X)), and outputs a probability estimate across the geographic region P . We would like to understand the uncertainty in the output g(A) as a function of the uncertainty of the input A. This can be done following steps similar to those in the delta method [4, p.75].
To understand the uncertainty in the output g(A), we begin by writing a first order Taylor expansion of g around E(A) where ∇g(·) is an n × k matrix of partial derivatives, with entry (i, j) specified by ∂pi ∂aj . If we can compute an expression for these partial derivatives, then everything on the right hand side above is constant, except [A − E(A)] whose distribution we know because we know the distribution of A. g(A) is an affine transformation of [A − E(A)], and can be approximated as To complete the analysis of the output uncertainty, we continue by deriving an expression for ∂pi ∂aj . We introduce some additional notation, following Kapur et al. [2]. Let λ be the Lagrange multiplier for constraint (1c), and µ j be the multiplier for constraint (1b) for j = 1, . . . , k. It can be shown [2] that the optimal p i have the expression PLOS 2/7 Using the constraint (1c), we have We can now substitute the expression for e λ+1 back into (7) to derive We now have an expression of the p i in terms of the dual multipliers µ j . But, we would like to compute ∂pi ∂aj , which we can do by first computing partial derivatives with respect to the µ j and using the expression What remains to be computed is ∂pi ∂µr and ∂µr ∂aj . From (9), we have where we derived the first equality from the chain rule, the second equality by substituting using expression (9), and the third equality by using the fact that a r = n t=1 f r (x t )p t because p i 's are feasible and optimal in the maximum entropy optimization.
Then, we want to get the value of ∂µr ∂aj . It is hard to get the expression of µ r in terms of a j , however, we can derive ∂aj ∂µr and use the Inverse Function Theorem [5] to calculate ∂µr ∂aj . To get the expression of a j in terms of µ r , we substitute the expression of the optimal p i (7) into constraint (1b), and then plug-in the expression of λ in terms of µ r (8). Then, we express a j as From 10, we have where we derived the first equality from the chain rule, the second equality by substituting using expression (9). The final equality comes from the definition of covariance, where we take the covariance of feature f r and f j with respect to the maximum entropy model results p i . If the determinant of the covariance is non-zero, following the Inverse Function Theorem [5], the inverse is differentiable. We denote the covariance matrix of features with respect to the maximum entropy model results as Ψ, where Ψ rj = cov P (f r , f j ). We denote the inverse covariance matrix as Ψ −1 and refer to its (r, j)th entry as (Ψ −1 ) rj . By the inverse function theorem, ∂µr ∂aj is equal to (Ψ −1 ) rj . Finally, we can express the ∂pi ∂aj as

Increasing programming speed
In the calculation of the relative probability for Aedes aegypti for a 1km 2 square grid, we have 933,680 grid cells in total. The problem of computing the covariance of the output mainly comes from Σ where The matrix is of size 933680 × 933680, which can cause out of memory errors. To figure out a way of speeding the calculation, we first split Σ into two parts, where The covariance of the output P can be estimated as First, to calculate ∇g · F·Σ1·F T m · (∇g) T , we define The i th diagonal element of ∇g · F · Σ 1 · F T · (∇g) T , denoted as d 1 i , can be calculated as To calculate ∇g · F·Σ2·F T m · (∇g) T , we define The i th diagonal element of ∇g · F · Σ 2 · F T · (∇g) T , denoted as d 1 i , is calculated by After we have calculated both d 1 i and d 2 i based on 12 and 13, we can calculate the variance of the p i as

Comparison between Analytic method and Poisson PPM
Poisson PPM was proved to be equivalent to maximum entropy model with hidden assumptions of independence data. We showed the plots of variance vs. point estimations and standard deviation comparison between Poisson PPM and analytic method for both Dengue importation probability and Aedes Aegypti suitability probability. Fig 1a shows the relationship between Dengue importation probability point estimates and estimated variance using analytic method, which indicating the possible improper use of Poisson PPM approach. Fig 1b shows a standard deviation comparison between the analytic and Poisson PPM method for Dengue importation probability. The regression line between two results is s a = 0.054s p with R 2 = 0.805, where s a and s p stand for the standard deviation estimates from the analytic and Poisson PPM methods, respectively. Poisson PPM gives a much larger standard deviation comparing to analytic and bootstrap method indicating the possible violate of the case location independence assumption. Fig 1c shows the relationship between Aedes Aegypti suitability point estimates and estimated variance using analytic method. There is more linear relationship comparing to Dengue importation cases. Fig 1d shows the standard deviation comparison between analytic method and Poisson PPM approach for Aedes Aegypti suitability. Each red dot represent the standard deviation estimates for each grid using Poisson PPM and analytic method respectively. The blue dot shows the diagonal line when two methods aligned well. We have the relationship s a = 0.0917s p with R 2 = 0.812, where s a and s p stand for the standard deviation estimates from the analytic and Poisson PPM methods, respectively. Similar as Dengue case, Poisson PPM doesn't seem to functional well with a much larger standard deviation estimation comparing to analytic and bootstrap method.