Spatial distribution of invasive species: an extent of occurrence approach

Ecological Risk Assessment faces the challenge of determining the impact of invasive species on biodiversity conservation. Although many statistical methods have emerged in recent years in order to model the evolution of the spatio-temporal distribution of invasive species, the notion of extent of occurrence, formally deﬁned by the International Union for the Conservation of Nature, has not been properly handled. In this work, a novel and ﬂexible reconstruction of the extent of occurrence from occurrence data will be established from nonparametric support estimation theory. Mathemati-cally, given a random sample of points from some unknown distribution, we establish a new data-driven method for estimating its probability support S in general dimension. Under the mild geometric assumption that S is r − convex, the smallest r − convex set which contains the sample points is the natural estimator. A stochastic algorithm is proposed for determining an optimal estimate of r from the data under regularity conditions on the density function. The performance of this estimator is studied by reconstructing the extent of occurrence of an assemblage of invasive plant species in the Azores archipelago.


Introduction
in the computational geometry literature for providing reasonable global reconstructions if the sample points are (approximately) uniformly distributed on the set S (see Edelsbrunner 2014). In fact, despite being r −convexity a more general condition than convexity, C r (X n ) can achieve the same convergence rates than H (X n ) as proved by Rodríguez-Casal (2007). However, this estimator presents an important disadvantage in practice: it depends on the commonly unknown parameter r . Although the influence of r can be considerable, it must be specified by the practitioner (see Joppa et al. 2016) or selected through practical procedures without theoretical guarantees (see Burgman and Fox 2003). For the example of invasive species in Azorean islands, Fig. 5 shows C r (X 740 ) for different values of r . Small values of r provide fragmented estimators (many isolated points and connected components) leading to an EOO reconstruction which resembles X n (Fig. 5: second row,right). For an intermediate value of r , a realistic reconstruction of the EOO is obtained since sea areas are not inside the estimator (Fig. 5: third row,left). However, if large values of r are considered, then C r (X n ) basically coincides with H (X n ) (Fig. 5: third row,right). Therefore, arbitrary choices of r may provide incongruous EOO estimations.
Most of the available results in the literature about support estimation make special emphasis on asymptotic properties, especially consistency and convergence rates, but they do not usually give any criterion for selecting the unknown parameter r in C r (X n ) from the sample. The aim of this paper is to overcome this drawback and present a method for selecting the parameter r for the r −convex hull estimator from the available data. This problem has scarcely been studied in the statistical literature with just a couple of references available on the topic. First, Mandal and Murthy (1997) proposed a selector for r based on the concept of minimum spanning tree, but only consistency of the method was provided without considering optimality issues. Later, Rodríguez-Casal and Saavedra-Nieves (2016) proposed an automatic selection criterion based on a very intuitive idea for the selection of r but under the important restriction that the sample distribution is uniform. The idea for selecting r is as follows. According to Fig. 5 (bottom,right), sea areas are contained in C r (X n ) if the selected r is too large. So, the estimator contains a large ball empty of sample points, see gray balls in Fig. 5 (top,left) and (bottom, right). Janson (1987) calibrated the size of this maximal ball (or spacing) when the sample distribution is uniform on S. Berrendero et al. (2012) used this result to test uniformity when the support is unknown. However, Rodríguez-Casal and Saavedra-Nieves (2016) followed an opposite approach. They assume that X n comes from a uniform distribution on S and if a big enough spacing is found in C r (X n ), then it is incompatible with the assumption that data are uniform. As a consequence, it is concluded that r is too large. Therefore, it seems natural to select the largest value of r compatible with the uniformity assumption on C r (X n ).
Recently, Aaron et al. (2017) extended the results by Janson (1987) to the case where the data are generated from a density f that is bounded from below and Lipschitz continuous restricted to its bounded support. Here, we will use this extension in order to derive a test to decide, given a fixed r > 0, whether the unknown support S is r −convex with no more information apart from X n . In this case, if a large enough spacing is found in C r (X n ), then the null hypothesis of r −convexity will be rejected. A new data-driven selector for the shape index r will be established from this test. Following the scheme in Rodríguez-Casal and Saavedra-Nieves (2016), it is proposed to choose the largest value of r compatible with the r −convexity assumption. Once the parameter r is estimated from X n , a new data-driven support reconstruction, based on the estimator of r , will be proposed. As a consequence, a flexible reconstruction for the EOO, based on available data, will be obtained. Furthermore, when the support is convex, our EOO estimator will be similar to H (X n ). Therefore, the EOO definition given by IUCN is generalized. This paper is organized as follows. Mathematical tools are introduced in Sect. 2. First, the geometric assumptions on S and the optimal value of the parameter r to be estimated are introduced. Then, the maximal spacing and its estimator are formally defined. Some regularity assumptions on f are also established. In Sect. 3, we propose a procedure for testing the null hypothesis that S is r −convex for a given r > 0. This test will play a key role in the definition of the consistent estimator of r . Then, a new data-driven estimator for the support S is proposed in Sect. 4 and it will be seen that it achieves the same convergence rates as the convex hull for estimating convex sets. The main numerical features involving the practical application of the algorithm are exposed in Sect. 5. Section 6 contains a simulation study in order to analyze the performances of the r −convexity test and the proposed estimator of the parameter r . In Sect. 7, the behavior of the new support reconstruction will be analyzed estimating the EOO of an assemblage of terrestrial plant species in two Azorean islands. Conclusions are exposed in Sect. 8. Finally, proofs of theoretical results are deferred to Sect. 9.

Mathemathical tools
Regularity conditions, namely shape assumptions on S, will be introduced next. In addition, we will discuss which is the optimal value of the shape index r to be estimated. Next, required conditions on the density function f will be also presented. Finally, basic notions on maximal spacings are established.

About geometric assumptions on S and the optimal r
In this work, S is assumed to be r −convex for some r > 0. Definition 1 establishes the formal definition of this geometric property. and B r (x), the open ball with center x and radius r , whereas D c denotes the complementary of D .
In practice, C r (X n ) can be computed as the intersection of the complements of all open balls of radius larger than or equal to r that do not intersect X n . To illustrate the importance of a good choice of r , Fig. 2 shows the computation of C r (X 740 ) for (X 740 ) (red color) and B r * (x) for r * ≥ r (gray color) such that B r * (x) ∩ X 740 = ∅ taking r = 0.3 (left) and r = 5 (right) r = 0.3 (left) and r = 5 (right) for the example in Azorean islands. Computations have been done using the alphahull package in R, see Pateiro-López and Rodríguez-Casal (2010). Note that C 0.3 (X 740 ) is an acceptable EOO reconstruction equal to the intersection of the complements of all gray open balls represented. However, if we select r = 5, marine areas are clearly inside the C 5 (X 740 ).
It is well-known that the concept of r −convex hull is closely related to the closing of A by B r (0) from the mathematical morphology, see Serra (1983). It can be shown that The problem of reconstructing a r −convex support S using a data-driven procedure could be easily solved if the parameter r is selected from the data set. The first step is to determine precisely the optimal value of r to be estimated, which is established in Definition 2: we propose to estimate the largest value of r which verifies that S is r −convex.
Definition 2 Let S ⊂ R d a compact, nonconvex and r −convex set for some r > 0. It is defined Remark 1 For simplicity in the exposition, it is assumed that S is not convex; otherwise, r 0 would be infinity, and the convex hull of the sample provides a good reconstruction.
Remark 2 If the supreme in (1) is a maximum, then S is r convex for r ≤ r 0 . In this case, if r < r 0 , C r (X n ) is a non-admissible estimator since it is always outperformed by C r 0 (X n ). This happens because, with probability one, C r (X n ) ⊂ C r 0 (X n ) ⊂ S.
is a more general condition (center). Circular ring with inner circle of radius r 0 (right) Proposition 2.4 in Rodríguez-Casal and Saavedra-Nieves (2016) ensures that under the shape restriction detailed below, the supreme in (1) is a maximum. The mild regularity condition we need is the following: (R) S an S c satisfy the rolling property with rolling positive constants r and λ, respectively.
Following Cuevas et al. (2012), it is said A satisfies the (outside) r −rolling condition if each boundary point a ∈ ∂ A is contained in a closed ball with radius r whose interior does not meet A. There exist interesting relationships between this property and r −convexity. In particular, Cuevas et al. (2012) proved that if A is compact and r −convex, then A fulfills the r −rolling condition. According to Fig. 3 (left), the reciprocal is not always true. Proposition 2.2 in Rodríguez-Casal and Saavedra-Nieves (2016) shows that (R) is a (mild) sufficient condition to ensure the r−rolling condition implies r−convexity. Condition (R) was essentially analyzed by Walther (1997Walther ( , 1999 but just the case r = λ was taken into account. In this work, the radius λ can be different from r, see Fig. 3 (center).

About maximal spacings
The optimal value of the shape index r to be estimated is just established in Definition 2. Some concepts on maximal spacings theory must be handled to propose a consistent estimate of r 0 .
The notion of maximal-spacing in several dimensions was introduced and studied by Deheuvels (1983) for uniformly distributed data on the unit cube. Later on, Janson (1987) extended these results to uniformly distributed data on any bounded set and derived the asymptotic distribution of different maximal-spacings notions without conditions on the shape of the support S. Aaron et al. (2017) generalized the results by Janson (1987) to the non-uniform case.
The shape of the considered spacings will be defined by a given set A ⊂ R d . For the validity of the theoretical results, it is sufficient to assume that A is a compact and convex set. For practical purposes, the usual choices are A = [0, 1] d or A = B 1 [0], being B r [x] the closed ball with center x and radius r . For a general dimension d, the first definition of maximal spacing is that used by Janson (1987) under the restriction of data are uniformly distributed: If the Lebesgue measure of the set A is one, Δ * n (X n ) d represents the Lebesgue measure of the largest set {x} ⊕ γ A ⊂ S\X n . The concept of maximal spacing can be related easily to the maximal inner radius when A = B 1 [0]. If Int(S) = ∅, the maximal inner radius of S is defined as Note that the value of the maximal spacing depends on S and also on X n . However, the definition of the maximal inner radius relies only on S. Aaron et al. (2017) extended the definition of maximal-spacing assuming that X n is drawn according to a density f with bounded support S, the Lebesgue measure of the set A is one and its barycenter is the origin of R d . In this more general setting, the maximal spacing is defined as The previous definition of maximal spacing relies also on density f . In this way, it distinguishes between low and high density regions. Throughout this paper, Janson (1987) calibrated the volume of the maximal spacing under uniformity assumptions without conditions on the shape of the support S. The corresponding extension established in Theorem 2 in Aaron et al. (2017) is shown in Theorem 1 modifying slightly the original hypotheses on f and on the shape of S. The result remains true if it is assumed that S is under (R) and the density function f satisfies ( f L 0,1 ): All through this paper, we assume that the random sample of points, X n , is generated from a density f that satisfies the regularity condition ( f L 0,1 ). Note that it includes the uniform distribution and also, more realistic scenarios with non-uniform sampling allowing to deal with observer biased data.
Theorem 1 (Aaron et al. (2017)) Let X n be a random and i.i.d sample drawn according to a density f that satisfies ( f L 0,1 ) with compact and nonempty support S under (R). Let U be a random variable with distribution and let β be a constant specified in Janson (1987). Then, we have that

Remark 3
The value of constant β does not depend on S. It is explicitly given in Janson (1987). Specifically, In particular, for the bidimensional case, β = 1.

About nonparametric estimation of maximal spacings
A plug-in estimator of the maximal spacing Δ n (X n ) will be introduced next. Note that the definition of Δ n (X n ) relies on the support S and also on the density function f (both are usually unknown). Under the assumption of r −convexity, S will be estimated as C r (X n ). As for the density function f , following the ideas in Aaron et al. (2017), a non-conventional density estimator will be introduced in Definition 3.

Definition 3 Let r > 0 and let
This nonparametric estimator only takes n different values: the evaluation of the usual kernel estimator on the sample points. In fact, for each point x ∈ S,f n (x) is equal to f n (X i ) where X i is the closest sample point to x. It can be checked that, with probability one and for n large enough, there exists a point in X n as close to x as desired. Therefore, this density estimator is only a simplification of the usual one with clear computational advantages for estimating Δ n and V n . This estimator is just a slight modification of the one proposed in Aaron et al. (2017), avoiding zero values.
Finally, some technical hypotheses on the kernel function must be established. Observe that this condition is satisfied, for instance, by the Gaussian kernel.
where p is a polynomial and φ is a bounded real function of bounded variation, verifying that c K = u K (u)du < ∞, K ≥ 0 and there exists r K and Then, we define the following plug-in estimator of Δ n (X n ) Given the definition off n and the assumption ( f L 0,1 ), it is expected thatf n does not go to zero on C r (X n ), see Lemma 1. This is important in formulae (2). For instance, if S is r −convex,δ(C r (X n )\X n ) should converge to zero as the sample size increases. However, if S C r (S), the plug-in estimator of Δ n (X n ) is expected to converge to a positive constant.

A new test for r−convexity
We will introduce a consistent hypothesis test based on X n drawn according to an unknown density f on the unknown support S, to assess r −convexity for a certain r > 0. This test is crucial for defining an estimator of r 0 that would allow the datadriven estimation of the support S.
Given r > 0, the null hypothesis that S is r −convex will be tested takingV n,r as statistic. The idea that supports this procedure is simple: Under ( f L 0,1 ) and (R), Theorem 1 allows us to detect which values of V n (X n ) are large enough to be incompatible with these two assumptions. A similar reasoning can be also applied if we consider V n,r , the test is based on the opposite approach: Under ( f L 0,1 ) and (R), if the test statistic takes large enough values, it will mean that the selected r is not appropriate and a smaller value of r should be considered.
Theorem 2 Let r > 0 and let X n be a random and i.i.d sample drawn according to a density f that satisfies ( f L 0,1 ) with compact and nonempty support S under (R). Let f n be the modified density estimator introduced in Definition 3 whose kernel function is supposed to satisfy condition (K p φ ) and the sequence h n of smoothing parameters fulfills h n = O(n −ζ ) for some 0 < ζ < 1/d. Letδ(C r (X n )\X n ) be the maximal spacing estimator established in equation (2). Given the statistical testing problem, has an asymptotic level not larger than α.
Remark 4 Note that the optimal kernel sequence size for estimating f , h n = h 0 n −1/(d+4) , satisfies the hypotheses under which Theorem 2 holds. Therefore, any reasonable bandwidth selector should be suitable for testing r −convexity.
Remark 5 Under (R) r 0 is the maximum of the set {γ > 0 : C γ (S) = S}. Hence, the hypotheses of the test introduced in Theorem 2 can be rewritten as follows: The performance of this test can be illustrated using the real database of invasive plants in Azorean islands. Given the sample X 740 , the practitioner could be interested in testing the null hypothesis that the EOO is r −convex, for instance, for r = 5. According to Fig. 5 (third row,right), it is clear that large Atlantic Ocean areas are inside C 5 (X 740 ) and the EOO is overestimated. Moreover,V 740,5 will be too large. In fact, although larger samples sizes were considered, its volume would take a constant value (see gray ball inside the EOO reconstruction). Therefore, the null hypothesis of 5−convexity should be rejected. Note that the situation is the opposite if testing r −convexity for r = 0.3 is the goal. In this case,V 740,0.3 should be clearly smaller. Furthermore, when the sample size increases, this volume tends to zero.

Selection and consistency results of the optimal smoothing parameter
The optimal estimation of the smoothing parameter r 0 from X n is based on the test previously proposed. Specifically, according to Definition 2, r 0 will be estimated bŷ That is, it is proposed to select the largest value of r compatible with the r −convexity assumption. Note that this choice depends on the significance level of the test, but its dependence is not explicitly given in the notation for the sake of clarity. The theoretical properties for the estimator of r 0 are considered next. First, the existence of the supreme defined in (3) must be guaranteed, a result which is proved in Theorem 3. In addition, it is also proved thatr 0 consistently estimates r 0 .
Theorem 3 Let f be a density function that satisfies ( f L 0,1 ) with compact, nonconvex and nonempty support S under (R). Letf n be the density estimator introduced in Definition 3 whose kernel function is supposed to satisfy condition (K p φ ) and the sequence h n of smoothing parameters fulfills h n = O(n −ζ ) for some 0 < ζ < 1/d. Let r 0 be the parameter defined in (1) andr 0 defined in (3). Let {α n } ⊂ (0, 1) be a sequence of significance levels converging to zero such that log(α n )/n → 0. Then,r 0 converges to r 0 in probability.

Remark 6
For the sake of clarity, S is assumed non-convex throughout the paper. However, if S is convex, it can be shown thatr 0 goes to infinity (which is the value of r 0 in this case) because, with high probability, the test is not rejected for all values of r .
We use again the example of invasive plants in Azorean islands in order to illustrate the behavior of this estimator. Under ( f L 0,1 ) and (R), ifV n,r is large enough, then the null hypothesis of r −convexity will be rejected. Therefore, a smaller value of r should be selected. This case corresponds to Fig. 5 (third row, right) taking r = 5. Observe that the null hypothesis of r −convexity would be also rejected for all r ≥ r because C r (X n ) ⊂ C r (X n ) and, consequently,V n,r ≥V n,r . However, the situation is completely opposite in Fig. 5 (second row, right) when r = 0.03. Here, the size of the maximal spacing found in C 0.03 (X 740 )\X 740 does not allow to reject that the support is 0.03−convex. As a consequence, a bigger r than 0.03 should be considered.

Consistency and convergence rates of resulting support estimator
The behavior of the random set Cr 0 (X n ) as an estimator of S can be studied once the consistency ofr 0 has been proved. Two metrics between sets are usually considered in order to assess the performance of a support estimator. Specifically, let A and C be two closed, bounded, nonempty subsets of R d . The Hausdorff distance between A and C is defined by Besides, if A and C are two bounded and Borel sets, then the distance in measure between A and C is defined by d μ (A, C) = μ(A C), where μ denotes the Lebesgue measure and , the symmetric difference, that is, A C = (A\C) ∪ (C\ A). Hausdorff distance quantifies the physical proximity between two sets, whereas the distance in measure is useful to quantify their similarity in content. However, neither of these distances are completely useful for measuring the similarity between the shape of two sets. The Hausdorff distance between boundaries, d H (∂ A, ∂C), can be also used to evaluate the performance of the estimators (see Baíllo and Cuevas 2001;Cuevas and Rodríguez-Casal 2004;Rodríguez-Casal 2007or Genovese et al. 2012. In particular, if lim r →r + 0 d H (S, C r (S)) = 0, then the consistency of Cr 0 (X n ) can be proved easily from Theorem 3. However, the consistency cannot be guaranteed if d H (S, C r (S)) does not go to zero as r goes to r 0 from above (asr 0 does, see Proposition 1 below). This problem can be solved by considering the estimator C r n (X n ) where r n = νr 0 with ν ∈ (0, 1) fixed. This ensures that, for n large enough, with high probability, C r n (X n ) ⊂ S. From the practical point of view, the selection of ν is not a major issue becauser 0 is numerically approximated and the computed estimator always satisfies this property without multiplying by ν.
Theorem 4 Let X n be a random and i.i.d sample drawn according to a density f that satisfies ( f L 0,1 ) with compact, nonconvex and nonempty support S under (R). Let r 0 be the parameter defined in (1) andr 0 defined in (3). Let {α n } ⊂ (0, 1) be a sequence converging to zero such that log(α n )/n → 0. Let be ν ∈ (0, 1) and r n = νr 0 . Then, eventually almost sure, for some positive constant D.
The same convergence rate holds for d H (∂ S, ∂C r n (X n )) and d μ (S C r n (X n )).

Numerical implementation
The main numerical aspects of the estimation algorithm of r 0 in (1) are detailed in what follows. Although the method proposed in this work is fully data-driven, its practical implementation depends on the specification of the significance level of the test α. Choosing this value is clearly a much simpler problem than the specification of the shape index r 0 . From Theorem 3, with probability one, for a large enough n, the existence of the estimatorr 0 defined in (3) is guaranteed. However, in practice, this estimator might not exist for a specific sample X n and a given value of the significance level α. Therefore, the influence of α must be taken into account from the practical point of view. The null hypothesis of r −convexity will be (incorrectly) rejected for 0 < r ≤ r 0 with probability at most α. This is not important from the theoretical point of view. Since we are assuming that α = α n goes to zero as the sample size increases this has not theoretical relevance. But, what should be done, for a given sample, if H 0 is rejected for all r (or at least all reasonable values of r )? In order to fix a minimum acceptable value of r , it is assumed that S (and, hence, its estimator) will have no connected components with probability content lesser than a value p ∈ (0, 1). From an empirical approach, the connected components of the estimator will contain at least a proportion p of sample points. If p is sufficiently close to zero, too fragmented estimators (for instance, with isolated points or insignificant clusters) will not be considered even in the case that we reject H 0 for all r . In this latter case, the minimum value that ensures that all connected components of the estimator contain at least a proportion p of sample points will be taken. Therefore, this parameter p can be interpreted as a geometric stopping criterion that does not appear in theoretical results because the sequence α n is assumed to tend to zero. Note that this shape assumption is very flexible. It does not limit too much the number of connected components. In fact, the estimator could present, as maximum number of clusters, the largest integer less than or equal to 1/p. In particular, if p is close enough to zero, the number of connected components can reach very high values. In fact, it is expected that, when the sample size increases, p = p n tending to zero. An alternative procedure, very similar and computationally simpler, is to establish a maximum number of connected components instead of p.
Although the definition ofr 0 depends on a test that must be applied several times in practice, remark that multiple testing does not play any role. It is possible to writer 0 = sup{r > 0 :V n,r ≤ c n,α }. SinceV n,r ≤V n,r , for r ≤ r , dichotomy algorithms can be used to computer 0 . The practitioner must select a maximum number of iterations I and two initial points r m and r M with r m < r M such that the null hypothesis of r M −convexity is rejected and the null hypothesis of r m −convexity is accepted. According to the previous comments, it is assumed that the proportion of sample points in each connected component of C r m (X n ) must be at least p. Choosing a value close enough to zero is usually sufficient to select r m . According to Fig. 5 (second row, right), the maximal spacing in C 0.03 (X n ) will be small enough to accept 0.03−convexity. Therefore, taking r m ≤ 0.03 will be a good choice. However, if selecting this r m is not possible because, for very low values of r , the hypothesis of r −convexity is still rejected, then r 0 is estimated as the positive closest value to zero r such that all connected components of C r (X n ) contain at least a proportion p of sample points. On the other hand, if a large enough spacing for having a statistically significant test cannot be found in H (X n ), then we propose H (X n ) as the estimator for the support.
To sum up, the following inputs should be given: the significance level α ∈ (0, 1), a maximum number of iterations I , a proportion p and two initial values r m and r M . Given these parametersr 0 will be computed as follows: 1. In each iteration and while the number of them is smaller than I : Some technical aspects related to the computation of the maximal spacings must be also mentioned. In the proposed procedure, the null hypothesis needs to be tested I times. Since it involves the calculation of the maximal spacing, one may be aware of computational cost of the method. Nevertheless, as noted by Rodríguez-Casal and Saavedra-Nieves (2016), this maximal spacing does not need to be specifically determined and it is enough to check if there exists a point x such that In this case,V n,r ≥ c n,α and, therefore, the null hypothesis of r −convexity will be rejected. Furthermore, note that if this disc exists, then x / ∈ B c x,w n,α (X k ) where c x,w n,α = c x) and X k denotes the sample point such that x ∈ V or(X k ). Therefore,f n (x) = f n (X k ).
Then, the centers of the possible maximal balls that belong to the Voronoi tile with nucleus X i (1, . . . , n) necessarily lie in B c X i ,w n,α (X i ) c ∩ V or(X i ). We will follow the next steps: 1. Determine the set of candidates for ball centers where E(m) ⊂ X n denotes the extremes of the m−shape of X n when m = min c X j ,w n,α : X j ∈ X n , see Edelsbrunner (2014). If x ∈ D(r ), then we can guar- 3. If M(r ) ≤ĉ n,α , then the null hypothesis of r −convexity is not rejected.
It should be noted that if, for all Therefore, these points can be discarded in order to determine D(r ). Furthermore, E(m), ∂C r (X n ) and ∂ Bĉ * n,α,r (X n ) can be easily computed (at least for the bidimensional case). See Pateiro-López and Rodríguez-Casal (2010) for further details.

Simulation results
The behavior of the test of r −convexity is established in Sect. 3, and the performance of the estimator of r 0 described in Sect. 5 will be analyzed through a simulation study.
Three  Fig. 4, r 0 is equal to 0.25 in the first two models and it takes the value 0.307 in the third model. Random samples have been generated by acceptance-rejection method on the three described supports. A standard multivariate normal distribution was considered for replicates generation in Models 1 and 2 and a mixture of two normal densities, for Model 3. The vector of weights of each normal component in the mixture is (1/2, 1/2). The vectors of means are (−0.5, 0.5) and (0.5, −0.5), respectively, and the covariance matrices are Σ = (σ i j ) 2×2 , verifying that σ i j = 1/9 if i = j and σ i j = 1/15 if i = j.
Regarding the performance of the r −convexity test, a total of 250 replicates of sample sizes 100, 500, 1000, and 2000 were generated for each of the simulation models. Tables 1, 2, and 3 contain the mean number of rejections when α = 0.1 for Models 1, 2, and 3, respectively. It should be noted that different values of the parameter r have been considered. The exact value of r 0 was represented using bold  numbers for Models 1 and 2. For Model 3, the closest value to r 0 was also using bold numbers. Simulation results indicate that the test of r -convexity is well calibrated, as it is established in Theorem 2. Under the null hypothesis, the mean number of rejections tends, as the sample size increases, to the nominal level α. However, the r −convexity test exhibits better consistency behavior for Models 2 and 3 than for Model 1 although the shape of the supports in scenarios 1 and 2 is quite similar. The result in Aaron et al. (2017) on which test rests deals with convergence in distribution of extrema. Following Hall (1991), this convergence in distribution can be extremely slow. Therefore, bigger sample sizes should be considered for Model 1. Bootstrap calibration could be also investigated.
Regarding the estimation of r 0 , the behavior of the algorithm proposed in Sect. 5 was analyzed when α = 0.01, for the 250 replicates of the three models previously introduced. A maximum number of four connected components were allowed. Table 4 contains the empirical means (M) and the standard deviations (SD) of these estimations for Models 1, 2, and 3, respectively. A total of 214 estimations of r 0 were equal to ∞ for Model 1 when n = 100. This situation is repeated for Model 2 when n = 100 (for 223 samples) and n = 500 (11 replicates). These estimations have not been taken into account for computing the averages shown in Table 4. Finally, it is worth to mention that the asymptotic calibration of the test does not preclude observing the consistency of the estimator for the parameter r 0 whose estimation is the main goal of this work.

Extent of occurrence estimation
As an illustrative example, the new support estimator introduced in this work will be used for reconstructing the EOO of an assemblage of terrestrial invasive plants in two islands of the Azores Archipelago, Terceira and São Miguel. For this real dataset, we have shown that convexity assumption is very restrictive. According to Fig. 5 (first and second rows, left), sea areas are inside the classical estimator of the EOO. Obviously, it is overestimated given that terrestrial invasive plants do not occupy the Atlantic Ocean. The goal here is to reconstruct the EOO in this application overcoming the described limitation. First, it is necessary to estimate the optimal value r 0 from the sample of 740 geographical locations. If we select the significance level α equal to 0.01 and p = 0.05, the resulting estimator isr 0 = 0.057. In Fig. 6, Cr 0 (X 740 ) is shown. According to the results obtained, the EOO reconstruction has two different connected components corresponding to the two Azorean islands. In this case,r 0 corresponds to the minimum value that guarantees that all connected components contain at least 37 sample points to avoid insignificant clusters of invasive plants. For smaller values of p (for instance, 0.02 that corresponds to a minimum number of geographical locations equal to 15 in each cluster), the number of connected components does not change. Unlike classical EOO estimator, sea areas are not inside the reconstruction. Therefore, a more realistic estimator of the EOO can be determined. It should be noted that the dataset may not In the third row (red color), C 0.3 (X 740 ) (left) and C 5 (X 740 ) (right) be independent. Several sample points are collected at the same day, and they are very close to each other. This does not happen all days, but some data are clearly clustered. This may cause thatr 0 underestimate r 0 . The condition on the size of the connected components prevents the estimator to take too small values in these situations, at it may happen in this illustrative example.
The new method, although designed for handling more complex situations, provides similar reconstructions to those corresponding to the convex hull in those cases where the classical reconstruction works appropriately. For showing this, we will focus on the geographical locations from São Miguel island. Separately, the EOO will be estimated from data corresponding to years 2015 and 2016. A total of 33 and 48 geographical locations are available in 2015 and 2016, respectively. Figure 7 contains the EOO estimator in 2015 (left) and 2016 (center). In 2015, the resulting reconstruction of the EOO is equal to H (X 33 ). In 2016,r 0 = 1.404; however, the estimation of the EOO obtained, C 1.404 (X 48 ), is not so different from the convex hull. This last illustration suggests that, if more amount of data are available by year, this kind of analysis could be useful for studying the temporal changes in the spatial pattern of organisms, including invasive plants, on an area of interest.

Conclusions and open problems
The main goal of this work is to propose a new data-driven method for reconstructing a r −convex support in a consistent way. The route designed to reach this goal can be summarized as follows: (1) Defining the optimal value of r , r 0 , to be estimated, (2) establishing a nonparametric test to assess the null hypothesis that S is r −convex for a given r > 0, (3) defining the estimator of r 0 that strongly relies on the previous test (4) checking that the estimator of r 0 and the resulting support reconstruction are consistent and (5) studying the performance of the r −convexity test and the estimation algorithm of r 0 through simulations.
The definition of the estimatorr 0 depends on the r −convexity test established that, of course, could be used in an independent way. In many practical situations where the support is completely unknown and only a sample of points is available, it can be interesting to test if the corresponding support distribution is r −convex.
Furthermore, the behavior of the proposed support estimator was illustrated through the estimation of the EOO of an assemblage of terrestrial invasive plants in two Azorean islands, providing a novel tool for ERA. In this particular case, where convexity assumption on the EOO is too restrictive, our support estimator provides a more realistic and sophisticated reconstruction. Besides, we have also shown that when the classical convex reconstruction works appropriately, our estimator offers similar reconstructions. Furthermore, we have shown that estimating the EOO from annual (or any other time period) occurrences could be useful for detecting temporal changes in the spatial pattern of organisms.
Note that the resulting support estimator is spatially flexible. In other words, it is able to distinguish the different disconnected components of the support. Therefore, it could be used for estimating the number of support connected components. Another relevant application deals with the integrated nested Laplace approximation (INLA) methodology introduced in Rue et al. (2009) and extended by Lindgren et al. (2011), establishing the stochastic partial differential equation (SPDE) approach. As noted by these authors, the application of a SPDE model requires the determination of a physical domain for the process where a discretization mesh is constructed. Different domains yield to different results, and our proposal can be viewed as a data-driven alternative to obtain a reasonable domain of interest.
Finally, another interesting problem and intimately related to the EOO reconstruction is to estimate the area of occupancy (AOO). The IUCN defined the AOO as the area within its extent of occurrence. Under r −convexity, we could estimate the AOO as the area of the r −convex hull of the sample points. However, this estimator suffers from the drawback of not being rate-optimal. Arias-Castro et al. (2019) proposed, under uniform distribution assumptions, an optimal volume estimator based on the sample r −convex hull using a sample splitting strategy that attains the minimax lower bound. Therefore, the problem of estimating the AOO could be studied from a different perspective in future.

Proofs
In this section, the proofs of the stated theorems are presented. Aaron et al. (2017) considers that f is Hölder continuous with respect to Lebesgue measure. Under ( f L 0,1 ), this condition is also satisfied. Furthermore, Aaron et al. (2017) also assumed that there exists k < d and C ∂ S > 0 such that N (∂ S, ) ≤ C ∂ S −k where N (∂ S, ) denotes the inner covering number of ∂ S. Under (R), Theorem 1 in Walther (1997) guaranteed that ∂ S is a C 1 (d − 1)−dimensional submanifold. In this case, it is well known that the previous assumption is fulfilled for k = d − 1. More details can be found in Aaron et al. (2017).

Proof of Theorem 2
First, we will prove (a) and then, (b). (a) Under H 0 (C r (S) = S), with probability one, C r (X n ) ⊂ S. Then, for some + n , see Lemma 1. If G n holds, we can ensure that Therefore, under G n , it is verified that Δ n (X n ) ≥ (1 − + n )δ(C r (X n )\X n ). Consequently, V n (X n ) ≥ (1 − + n ) dV n,r . Hence, ifV n,r > c n,α , it is satisfied that As a consequence, (b) An auxiliary result must be taken into account for completing the proof of (b).
Lemma 2 Let X n be a random and i.i.d sample drawn according to a density f that satisfies ( f L 0,1 ) with compact, nonconvex and nonempty support S under (R). Let r 0 be the parameter defined in (1). Then, for all r > r 0 , there exists an open ball B ρ (x) such that B ρ (x) ∩ S = ∅ and P B ρ (x) ⊂ C r (X n ), eventually = 1.
Proof Proof of Lemma 8.4 in Rodríguez-Casal and Saavedra-Nieves (2016) remains true if sample distribution is not uniform. Therefore, in this more general setting, it allows to guarantee that there exists a closed ball of radius ρ > 0 that, with probability one and for n large enough, is inside C r (X n )\S.
where λ 1 is a positive constant. Under H 1 (S is not r −convex, S C r (S)), we will prove that, R(C r (X n )\X n ) ≥ ρ > 0, e.a.s.
Lemma 2 allows to guarantee that there exists a closed ball of radius ρ > 0 that, with probability one and for n large enough, is inside C r (X n )\X n . Consequently, Then,δ The proof is finished taking into account that c n,α tends to zero.

Proof of Theorem 3
Some auxiliary results are necessary. First we will prove that, with probability tending to one,r 0 is at least as big as r 0 .
Proposition 1 Let f be a density function that fulfills condition ( f L 0,1 ) with compact, nonconvex and nonempty support S under (R). Letf n be the corresponding density estimator introduced in Definition 3 whose kernel function is supposed to satisfy condition (K p φ ) and the sequence h n of smoothing parameters fulfills h n = O(n −ζ ) for some 0 < ζ < 1/d. Let r 0 be the parameter defined in (1) andr 0 defined in (3). Let {α n } ⊂ (0, 1) be a sequence converging to zero. Then, lim n→∞ P(r 0 ≥ r 0 ) = 1.
Lemma 1 (ii) guarantees that P(G c n , i. o.) = 0. Therefore, the second term of the inequality is negligible and its is verified that lim sup P(V n,r 0 > c n,α n ) ≤ P(V n (X n ) > (1 − + n ) d c n,α n ).
According to Theorem 1, U (X n ) d → U when n → ∞. Furthermore, notice that U has a continuous distribution, so convergence in distribution implies that sup u |P (U (X n ) ≤ u) − P(U ≤ u)| → 0.
It remains to prove thatr 0 cannot be arbitrarily larger that r 0 . An auxiliary result must be proved.
Proposition 2 Let X n be a random and i.i.d sample drawn according to a density f that satisfies ( f L 0,1 ) with compact, nonconvex and nonempty support S under (R). Let r 0 be the parameter defined in (1) and {α n } ⊂ (0, 1) a sequence converging to zero such that log(α n )/n → 0. Then, for any > 0, P r 0 ≤ r 0 + , eventually = 1.
Since, with probability one, X n ⊂ S we have B ρ (x 0 ) ∩ X n = ∅. Then, {x 0 } ⊕ ρ B 1 [0] ⊂ C r (X n )\X n . According to Lemma 1 (i), with probability one and for n large enough, there exists a constant λ 1 > 0 such that for all x ∈ C r (X n ), it is verified thatf 1/d n (x) ≥ λ 1 . Then, let γ be the positive constant λ 1 ρw 1/d d . Then, it is trivial to check that, with probability one and for n large enough, Therefore,δ(C r (X n )\X n ) ≥ γ > 0 and, consequently,V n,r = c γ > 0. Furthermore, since C r (X n ) ⊂ C r (X n ) for all r ≥ r , it is satisfied thatV n,r ≥V n,r = c γ > 0. On the other hand, since −u α n / log(α n ) = log(− log(1−α n ))/ log(α n ) → 1, we have c n,α n → 0. Then, with probability one and for n large enough, we have c n,α n < c γ . Therefore, according the definition established in (3),r 0 ≤ r .

Proof of Theorem 4 Theorem 3 of Rodríguez-Casal
, and D is some constant. Under the hypothesis of Theorem 4, this holds for any r ≤ min{r, λ}. Fix one r ≤ min{r, λ} such that r < νr 0 and define R n = { r ≤ r n ≤ r 0 }. Since, by Theorem 3, r n = νr 0 converges in probability to νr 0 and r < νr 0 < r 0 , we have that P(R n ) → 1. If the events E n and R n hold (notice that P(E n ∩ R n ) → 1), we have C r (X n ) ⊂ C r n (X n ) ⊂ S and, therefore, d H (S, C r n (X n )) ≤ d H (S, C r (X n )) ≤ D log n n 2/(d+1) .
This completes the proof of the first statement of Theorem 4. Similarly, it is possible to prove the result for the other error criteria considered in Theorem 4.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.