Spatial point processes intensity estimation with a diverging number of covariates
Achmad Choiruddin Jean-François Coeurjolly Frédérique Letué
stat.ME math.ST stat.TH







1]Achmad Choiruddin
2, 3]Jean-François Coeurjolly
3]Frédérique Letué
[1]Department of Mathematical Sciences, Aalborg University, Denmark
[2]Department of Mathematics, Université du Québec à Montréal (UQAM), Canada
[3]Laboratory Jean Kuntzmann, Department of Probability and Statistics, Université Grenoble Alpes, France
Spatial point processes intensity estimation 
 with a diverging number of covariates
    [
    Received <date>; Accepted <date>
====================================================================================






Abstract: Feature selection procedures for spatial point processes parametric intensity estimation have been recently developed since more and more applications involve a large number of covariates. In this paper, we investigate the setting where the number of covariates diverges as the domain of observation increases. In particular, we consider estimating equations based on Campbell theorems derived from Poisson and logistic regression likelihoods regularized by a general penalty function. We prove that, under some conditions, the consistency, the sparsity, and the asymptotic normality are valid for such a setting. We support the theoretical results by numerical ones obtained from simulation experiments and an application to forestry datasets.



Key words and phrases:
Campbell formula, estimating equation, high dimensional regression, regularization method, variable selection.


figure
table


1214pt plus.8pt minus .6pt



§ INTRODUCTION
 


 §.§ Background

Spatial point pattern data arise in many contexts, e.g. ecology <cit.>, epidemiology <cit.>, criminology <cit.>, biology <cit.> and astronomy <cit.>, where interest lies in describing the distribution of an event in space.  Stochastic models generating spatial point patterns are called spatial point processes <cit.>.
	
Usually, the first step to analyze spatial point pattern data is to investigate the intensity. The intensity serves as the first-order characteristics of a spatial point process and often becomes the main interest in many studies, especially when the intensity is suspected to depend on spatial covariates. Examples include the study of spatial variation of specific disease risk related to pollution sources <cit.>, crime rate analysis in a city related to some demographical information <cit.>, and modeling of the spatial distribution of trees species in a forest related to some environmental factors <cit.>.

We focus in this study on the log-linear model for the intensity function of an inhomogeneous spatial point process defined by

    ρ (u;β)=exp(𝐳(u)^⊤β), u ∈ D ⊂ℝ^d,

where 𝐳(u)={ z_1(u),…,z_p(u)}^⊤ are the p spatial covariates measured at location u, d represents the state space of the spatial point processes (usually d=2,3) and β={β_1,…,β_p}^⊤ is a real p-dimensional parameter. Hence, our main concern is to assess the magnitudes of the vector β. For parametric estimation, while maximum likelihood estimation <cit.> has been widely implemented for Poisson point processes models, estimating equation-based methods <cit.> are simpler to implement for more general spatial point processes models, overcoming the possible drawback of MCMC methods which are usually computational expensive <cit.>. However, when the number of covariates is relatively large, maximum likelihood estimation  and estimating equation-based methods become undesirable: all covariates are selected yielding an increasing standard error for parameter estimates.



 §.§ Feature selection techniques


To select significant covariates, one may consider a traditional procedure such as a stepwise method. This technique starts with an initial set of covariates, then considers adding or deleting a covariate from the current set at each iteration using a criterion such as an F-statistic or AIC. However, such procedure has a number of limitations: it can be numerically unstable and exhibits high variance due to its discrete procedure <cit.>. It is even computationally unfeasible  especially when the number of covariates is too large <cit.>. 
	
To overcome this drawback, regularization techniques have recently been developed for spatial point processes intensity estimation. Such methods are able to perform variable selection while keeping interesting properties in terms of prediction. For Poisson point process models, the idea is to penalize the Poisson likelihood by a penalty function such as l_1 penalty <cit.>. For more general point process models, instead of employing the likelihood of the processes which often requires computational intensive MCMC methods <cit.>, penalized versions of estimating equations based on Campbell theorem derived both from Poisson and logistic regression likelihoods have been developed <cit.>. Furthermore, <cit.> and <cit.> show that, under some conditions, the estimates obtained from such procedures are consistent, sparse, and asymptotically normal.



 §.§ Issues in high dimensional data


The motivation of our paper comes from a study of biodiversity in a 50-hectare region (D=1,000m×  500m) of the tropical moist forest of Barro Colorado Island (BCI) in central Panama, where censuses have been carried out such that all free-standing woody stems at least 10 mm diameter at breast height were identified, tagged, and mapped, resulting in maps of over 350,000 individual trees with around 300 species <cit.>. In the same region, many environmental covariates such as topographical attributes and soil properties have been also collected. In particular, we are interested to study the spatial distribution of 3,604 locations of Beilschmiedia pendula Lauraceae (BPL) trees and to model its intensity as a parametric function of 93 covariates consisting of 2 topological attributes, 13 soil properties and 78 interactions between two soil nutrients.
	
Although it seems that the number of covariates is not very large with respect to the number of data points, two hours are required to estimate the parameters and select covariates using a standard stepwise procedure. To do this, we use the  function in  which intrinsically assumes that 𝐗 is a Poisson point process since we use the AIC in the stepwise procedure. For a general point process, other criteria could be investigated such as the one based on the F statistic, but they require to estimate the asymptotic covariance matrix of the estimates at each step: even if we know the right covariates, it is known as a difficult task (see e.g. 
<cit.>) especially when the number of parameters is large. That would easily triple the time of this estimation/selection procedure. To evaluate the performance of such a selection/estimation procedure, a simulation would be required which is unrealistic (1000 replications of a single model would take 250 days). This motivates us to consider regularization methods.
	
<cit.> and <cit.> are the first two theoretical works. Both these works have the important limitation that the number of covariates p is finite. We extend this in the present paper. Asymptotic properties which consider a diverging number of parameters for M-estimators have a long story <cit.> but have recently been investigated for penalized regression estimators by <cit.>. In particular, as argued by <cit.>, even though the asymptotic properties (i.e., consistency, sparsity, and asymptotic normality) proposed by <cit.> for penalized generalized linear models under the assumption that p is finite, are encouraging, there are many naive and simple model selection procedures which possess those properties. Establishing the validity of these asymptotic properties in a diverging number of parameters setting is, therefore, a major importance. We study this type of asymptotic properties in the spatial point processes framework. Hence, our work can be regarded as an extension of the study conducted by <cit.>.
	
A standard way of measuring asymptotic for spatial point process is the increasing domain asymptotic. Therefore, we investigate the problem where p=p_n grows with |D_n| the volume of the observation domain. In our setting, |D_n| plays the same role as n, the number of observations, in standard problems such as in linear models or generalized linear models. We obtain  consistency, sparsity, and asymptotic normality for our estimator. One of our main assumptions is that p_n^3/|D_n| → 0 as n →∞, which is similar to the one required by <cit.> when |D_n| is simply replaced by n (the sample size in their context).
	
Our results are general: (1) a large choice of penalty functions (either convex or non-convex function) and methods (e.g. ridge, lasso, elastic net, SCAD, and MC+) are available; (2) we include a large class of mixing spatial point processes. The implementation is done by combining the   <cit.>  package with the two  packages implementing penalized methods for generalized linear models:  <cit.> and  <cit.>.



 §.§ Outline of the paper

In Section <ref>, we introduce brief background on spatial point processes as well as regularization methods for spatial point processes intensity estimation. Section <ref> presents our asymptotic results. We investigate in Section <ref> the finite sample performance of our estimates in a simulation study and in an application to tropical forestry datasets. Conclusion and discussion are presented in Section <ref>. Proofs of the main results are postponed to Appendices <ref>-<ref>.



§ REGULARIZATION METHODS FOR SPATIAL POINT PROCESSES


This section gives brief introduction on spatial point processes and reviews regularization methods for spatial point processes intensity estimation previously studied by <cit.>  when the number of parameters is finite.

Let 𝐗 be a spatial point process on ℝ^d. We view 𝐗 as a locally finite random subset of ℝ^d. Let D ⊂ℝ^d be a compact set of Lebesgue measure |D| which will play the role of the observation domain. A realization of 𝐗 in D is thus a set 𝐱={x_1, x_2, …, x_m}, where x ∈ D and m is the observed number of points in D.  Suppose 𝐗 has intensity function ρ and second-order product density ρ^(2). Campbell theorem <cit.> states that, for any function k: ℝ^d → [0,∞) or k: ℝ^d ×ℝ^d → [0,∞)

    𝔼( ∑_u ∈𝐗 k(u) ) =∫_ℝ^d k(u) ρ (u)du
       𝔼(∑_u,v ∈𝐗^≠ k(u,v))=∫_ℝ^d∫_ℝ^d k(u,v) ρ^(2) (u,v)du dv.

We may interpret ρ(u) du as the probability of occurence of a point in an infinitesimally small ball with centre u and volume du. Intuitively, ρ^(2)(u,v)dudv is the probability for observing a pair of distinct points from 𝐗 occuring jointly in each of two infinitesimally small balls with centres u,v and volume du, dv. For further background materials on spatial point processes, see for example <cit.>.

In our study, we assume that the intensity function depends on parameter β, ρ= ρ(·;β). The standard parametric methods for estimating β are by maximizing the weighted Poisson likelihood <cit.> or the weighted logistic regression likelihood <cit.> given respectively by

    ℓ_PL(w;β)  =    ∑_u ∈𝐗∩ D w(u)logρ(u;β) - ∫_D w(u)ρ(u; β)du, 
    ℓ_LRL(w; β) =    ∑_u ∈𝐗∩ D w(u) log ( ρ(u;β)/δ(u)+ρ(u;β) )
       - ∫_D w(u) δ(u) log (ρ(u;β) + δ(u)/δ(u) ) du,

where w(·) is a weight non-negative function depending on the first and the second-order characterictics of 𝐗 and δ(·) is a non-negative real-valued function. The solution of maximizing (<ref>) (resp. (<ref>)) is called Poisson estimator (resp. the logistic regression estimator). We refer readers to <cit.> for further details on the weight function w(·) and to <cit.> for the role of function δ(·).

These standard methods cannot perform variable selection. To do so, <cit.> and <cit.> suggest to maximize a penalized version of (<ref>)-(<ref>)

    Q(w; β)=ℓ(w; β) - |D|∑_j=1^p p_λ_j(|β_j|),

where ℓ(w; β) is either the Poisson likelihood (<ref>) or the logistic regression likelihood (<ref>). We refer the second term of (<ref>) to a penalization term. In this term, we have mainly two parts: (1) a penalty function p_λ parameterized by λ≥ 0 and (2) the volume of the observation domain |D| which plays the same role as the sample size in the spatial point process framework.

For any non-negative λ, we say that p_λ(·) is a penalty function if p_λ is a non-negative function with p_λ(0)=0. Some examples, described in Table <ref>, include  l_2 penalty <cit.>, l_1 penalty <cit.>, elastic net <cit.>, SCAD <cit.>, and MC+ <cit.>. Note that, as indicated by (<ref>), we allow each direction to have different tuning parameters λ_j, j=1,…,p. Such a method is called an adaptive method (e.g. adaptive lasso <cit.> and adaptive elastic net <cit.>). For further backgrounds about penalty function and regularization methods, see, for example, <cit.>.









§ ASYMPTOTIC PROPERTIES


In this section, we present asymptotic properties of the regularized Poisson estimator when both |D_n| →∞ and p_n →∞ as n →∞. In particular, we consider 𝐗 as a d-dimensional point process observed over a sequence of observation domain D=D_n, n=1,2, … which expands to ℝ^d as n →∞. We assume that 𝐗 has a log-linear form  given by (<ref>) for which the dimension of parameter β, denoted now by p_n, diverges to ∞ as n →∞.  In Section <ref>, we provide notation and conditions, and discuss the differences   from the setting where p is fixed. Our main results are presented in Section <ref>.  For sake of conciseness, we do not present the asymptotic results for the regularized logistic regression estimator. The results are very similar. The main difference is lying in the conditions (𝒞.<ref>) and (𝒞.<ref>) for which the matrices  𝐀_n, 𝐁_n, and 𝐂_n have a different expression (see Remark <ref>).



 §.§ Notation and conditions
 
Throughout this section and Appendices <ref>-<ref>, let 

    ℓ_n (w;β)=   ℓ_n,PL(w;β) 
    
     =   ∑_u ∈𝐗∩ D_n w(u)logρ(u;β) - ∫_D_n w(u)ρ(u; β)du,  
    
    Q_n(w;β)=    ℓ_n (w;β) -  |D_n| ∑_j=1^p_n p_λ_n,j(|β_j|),

be respectively the weighted Poisson likelihood and its penalized version.

Let β_0 = {β_01,…,β_0s,β_0(s+1),…,β_0p_n}^⊤={β^⊤_01,β^⊤_02}^⊤ = (β_01^⊤, 0^⊤)^⊤ denote the p_n-dimensional vector to estimate, where β_01 is the s-dimensional vector of non-zero coefficients and β_02 is the (p_n-s)-dimensional vector of zero coefficients. We assume that the number of non-zero coefficients, s, does not depend on n.  Let 𝐳_01 and 𝐳_02 denote the corresponding s-dimensional and (p_n - s)-dimensional vectors of spatial covariates. We denote the regularized Poisson estimator by β̂=(β̂_1^⊤, β̂_2^⊤)^⊤. 




We recall the classical definition of strong mixing coefficients adapted to spatial point processes <cit.>: for k,l ∈ℕ∪{∞} and q ≥ 1, define

    α_k,l(q)=sup{   |P(A ∩ B)-P(A)P(B)|: A∈ℱ (Λ_1), B ∈ℱ (Λ_2), 
       Λ_1 ∈ℬ (ℝ^d), Λ_2 ∈ℬ (ℝ^d), |Λ_1| ≤ k, |Λ_2| ≤ l, d(Λ_1, Λ_2) ≥ q },

where ℱ is the σ-algebra generated by 𝐗∩Λ_i, i=1,2, d(Λ_1, Λ_2) is the minimal distance between sets Λ_1 and Λ_2, and ℬ(ℝ^d) denotes the class of Borel sets in ℝ^d.

We  define the p_n × p_n matrices 𝐀_n(w;β_0), 𝐁_n(w;β_0) and 𝐂_n(w;β_0) by

    𝐀_n(w;β_0)   =∫_D_n w(u)𝐳(u)𝐳(u)^⊤ρ(u;β_0)du, 
    𝐁_n(w;β_0)   =∫_D_n w(u)^2 𝐳(u)𝐳(u)^⊤ρ(u;β_0)du, 
    𝐂_n(w;β_0)   =∫_D_n∫_D_n w(u) w(v)𝐳(u)𝐳(v)^⊤{g(u,v)-1}ρ(u;β_0) ρ(v;β_0) du dv,

where g(u,v) is the classical pair correlation function <cit.> given by

    g(u,v)=ρ^(2)(u,v)/ρ(u)ρ(v),

when both ρ and ρ^(2) exist with the convention 0/0=0. For a Poisson point process, we have g(u,v)=1 since ρ^(2)(u,v)=ρ(u)ρ(v). If, for example, g(u,v)>1 (resp. g(u,v)<1), this indicates that pair of points are more likely (resp. less likely) to occur at locations u,v than for a Poisson point process. 

We denote the s × s top-left corner of 𝐀_n(w;β_0)  (𝐁_n(w;β_0), 𝐂_n(w;β_0)) by 𝐀_n,11(w;β_0) (𝐁_n,11 (w;β_0), 𝐂_n,11(w;β_0)). It is worth noticing that 𝐀_n,11(w;β_0), 𝐁_n,11(w;β_0) and 𝐂_n,11(w;β_0) depend on n only through D_n and not through p_n. In what follows, for a squared symmetric matrix 𝐌_n, ν_min(𝐌_n) and ν_max(𝐌_n) denote respectively the smallest and largest eigenvalue of 𝐌_n.

Under the conditions (𝒞.<ref>)-(𝒞.<ref>), we define the sequences a_n, b_n and c_n by

    a_n    =max_j=1,…,s | p'_λ_n,j(|β_0j|)| , 
    
    b_n    =inf_j=s+1,…,p_ninf_|θ| ≤ϵ_n 
    θ≠ 0 p'_λ_n,j(θ) , ϵ_n=K_1 √(p_n/|D_n|),	
    
    c_n    =  max_j=1,…,s |p^''_λ_n,j(|β_0j|) |,

where K_1 is any positive constant. 

Consider the following conditions (𝒞.<ref>)-(𝒞.<ref>) which are required to derive our asymptotic results:



  * For every n ≥ 1, D_n=nE={ne: e ∈ E}, where E ⊂ℝ^d is convex, compact, and contains the origin of ℝ^d in its interior. 

  * The intensity function has the log-linear specification given by (<ref>) where β∈Θ and Θ is an open convex bounded set of ℝ^p_n. Furthermore, we assume that there exists a neighborhood Ξ(β_0) of β_0 such that 

    sup_n≥ 1 sup_β∈Ξ(β_0) sup_u ∈ℝ^dρ(u; β) < ∞.
  

  * The covariates 𝐳 and the weight function w satisfy

    sup_n≥ 1 sup_i=1,…,p_n sup_u ∈ℝ^d |z_i(u)| < ∞,      sup_u ∈ℝ^d w(u)<∞.
 

  * There exists an integer t ≥ 1 such that for k=2, …, 2+t, the product density ρ^(k) exists and satisfies ρ^(k)< ∞. 

  * For the strong mixing coefficients (<ref>), we assume that there exists some t̃ > d(2+t)/ t such that α_2,∞(q)=O(q^-t̃). 

  * lim inf_nν_min(|D_n|^-1{𝐁_n,11(w;β_0)+𝐂_n,11(w;β_0)})>0. 

  * lim inf_nν_min(|D_n|^-1𝐀_n(w;β_0))> 0.  

  * The penalty function p_λ(·) is non-negative on ℝ^+, continuously differentiable on ℝ^+ ∖{0} with derivative p_λ' assumed to be a Lipschitz function on ℝ^+∖{0}.
Furthermore, given (λ_n,j)_n ≥ 1,  j=1, …, s, we assume that there exists (r̃_n,j)_n ≥ 1, where r̃_n,j√(|D_n|/p_n)→∞ as n →∞, such that, for n sufficiently large, p_λ_n,j is thrice continuously differentiable in the ball centered at |β_0j| with radius r̃_n,j and we assume that the third derivative is uniformly bounded. 

  * p_n^3 / |D_n| → 0 as n →∞. 


Conditions (𝒞.<ref>)-(𝒞.<ref>) are quite similar to the ones required by <cit.> in the setting when the number of parameters to estimate is fixed. Condition (𝒞.<ref>) is slightly stronger since we have to ensure that ρ(u; β) is finite for β in the neighborhood of β_0. Note that sup_u ∈ℝ^dρ(u; β_0) < ∞ follows directly from condition (𝒞.<ref>). We derive asymptotic properties when both |D_n| and p_n tend to infinity with n. However, to obtain an estimator which is consistent and has two other properties: sparsity and asymptotic normality, we need that the number of covariates does not grow too fast with respect to the volume of the observation domain. This condition is  stated by condition (𝒞.<ref>) which is similar to the one required by <cit.> when |D_n| is simply replaced by n (the sample size in their context).




 §.§ Main results
 
We state our main results here. Proofs are relegated to Appendices <ref>-<ref>.

We first show in Theorem <ref> that the regularized Poisson estimator converges in probability and exhibits its rate of convergence.



Assume the conditions (𝒞.<ref>)-(𝒞.<ref>) and (𝒞.<ref>)-(𝒞.<ref>) hold. Let a_n and c_n be given respectively by (<ref>) and (<ref>). If a_n=O(|D_n|^-1/2) and c_n=o(1), then there exists a local maximizer β̂ of Q_n(w;β)  such that β̂ -β_0=O_P(√(p_n) (|D_n|^-1/2+a_n)).


This implies that, the regularized Poisson estimator is root-(|D_n|/p_n) consistent. Note that, as expected, the convergence rate is √(p_n) times the convergence rate of the estimator obtained when p is fixed <cit.>. In addition, when we compare our results with the ones obtained by <cit.>, who also considered a diverging number of parameters setting, our estimator has the same rate of convergence when we replace |D_n| by n to their context. This rate of convergence also appears in other contexts considering diverging number of parameters setting <cit.>.



Now, we demonstrate in Theorem <ref> that such a root-(|D_n|/p_n) consistent estimator ensures the sparsity of β̂; that is, the estimate will correctly set β_2 to zero with probability tending to 1 as n →∞, and β̂_1 is asymptotically normal.



Assume the conditions (𝒞.<ref>)-(𝒞.<ref>) are satisfied. If a_n √(|D_n|)→⁠ 0, b_n √(|D_n|/p_n^2)→∞ and c_n √(p_n)→ 0 as n→∞, the root-(|D_n|/p_n) consistent local maximizer β̂=(β̂_1^⊤, β̂_2^⊤)^⊤ in Theorem <ref> satisfies:


  * Sparsity: P(β̂_2=0) → 1 as n →∞,

  * Asymptotic Normality: |D_n|^1/2Σ_n(w;β_0)^-1/2(β̂_1- β_01)𝒩(0, 𝐈_s),

where

    Σ_n(w;β_0)=     |D_n|{𝐀_n,11(w;β_0)+|D_n| Π_n }^-1{𝐁_n,11(w;β_0)+𝐂_n,11(w;β_0)}
       {𝐀_n,11(w;β_0)+|D_n| Π_n }^-1, 
    Π_n =    diag{p”_λ_n,1(|β_01|),…,p”_λ_n,s(|β_0s|)}.

As a consequence, Σ_n(w;β_0) is the asymptotic covariance matrix of β̂_1. Here, Σ_n(w;β_0)^-1/2 is the inverse of Σ_n(w;β_0)^1/2, where Σ_n(w;β_0)^1/2 is any square matrix with Σ_n(w;β_0)^1/2(Σ_n(w;β_0)^1/2)^⊤=Σ_n(w;β_0).




For lasso and adaptive lasso, Π_n= 0. For other penalties, since c_n=o(1), then Π_n=o(1). Since 𝐀_n,11(w;β_0)=O(|D_n|) from conditions (𝒞.<ref>)-(𝒞.<ref>), |D_n|  Π_n is asymptotically negligible with respect to 𝐀_n,11(w;β_0).




Theorems <ref> and <ref> remain true for the regularized logistic regression estimator if we replace in the expression of the matrices 𝐀_n, 𝐁_n, and 𝐂_n, w(u) by w(u) δ(u)/ (ρ(u; β_0) + δ(u)), u ∈ D_n and extend the condition (𝒞.<ref>) by adding sup_u ∈ℝ^dδ(u)<∞.
	
The proofs of Theorems <ref> and <ref> for this estimator are slightly different mainly because unlike the Poisson likelihood for which we have ℓ_n^2(w;β)=-𝐀_n(w;β), for the regularized logistic regre ℓ_n^2(w;β) is now stochastic and we only have 𝔼(ℓ_n^2(w;β))=-𝐀_n(w;β). Despite the additional difficulty, we maintain that no additional assumption is required. 


We show in Theorem <ref> that the sparsity and asymptotic normality are still valid when the number of parameters diverges. By Remark <ref>, when n is large enough, Σ_n(w;β_0) in (<ref>) becomes approximately

    |D_n|{𝐀_n,11(w;β_0)}^-1{𝐁_n,11(w;β_0)+𝐂_n,11(w;β_0)}{𝐀_n,11(w;β_0) }^-1,

which is precisely the asymptotic covariance matrix of the estimator of β_01 obtained by maximizing the likelihood function or solving estimating equations based on the submodel knowing that β_02=0. This shows that when n is sufficiently large, our estimator is as efficient as the oracle one.

To satisfy Theorem <ref>, we require that a_n √(|D_n|)→ 0, b_n √(|D_n|/p_n^2)→∞ and c_n √(p_n)→ 0 as n →∞ simultaneously. In particular, conditions on a_n and c_n ensure the asymptotic normality of β̂_1 while condition on b_n is used to prove the sparsity. Conditions regarding a_n and c_n are similar to the ones imposed by <cit.> when |D_n| is replaced by n in their context. However, we require a slightly stronger condition on b_n than the one required by <cit.> which in the present setting could be written as b_n √(|D_n|/p_n)→∞. As compensation, we do not need to impose, as <cit.> did, for any 0 < K_2 < ∞,  ν_max(|D_n|^-1𝐀_n(w;β_0)) < K_2. Such a condition is not straightforwardly satisfied in our setting since the other conditions only imply that ν_max(|D_n|^-1𝐀_n(w;β_0))=O(p_n).

Further details regarding a_n, b_n and c_n for each method are presented in Table <ref>. For the ridge regularization method, b_n=0, preventing from applying Theorem <ref> for this penalty. For lasso and elastic net, a_n=K_3 b_n for some constant K_3>0 (K_3=1 for lasso). The two conditions a_n √(|D_n|)→ 0 and b_n √(|D_n|/p_n^2)→∞ as n →∞ cannot be satisfied simultaneously. This is different for the adaptive versions where a compromise can be found by adjusting the λ_n,j's, as well as the two non-convex penalties SCAD and MC+, for which λ_n can be adjusted. For the regularization methods we consider in this study, the condition c_n √(p_n)→ 0 is implied by the condition a_n √(|D_n|)→ 0 as n→∞ and condition (𝒞.<ref>).







§ NUMERICAL RESULTS


This section is devoted to present numerical results. More precisely, we conduct simulation experiments in Section <ref> to assess the finite sample peformance of our estimates and apply our method to an application in ecology in Section <ref>. We apply the regularized Poisson likelihood (PL) and the regularized weighted Poisson likelihood (WPL) to select covariates and estimate their coefficients. Similar approach can be straightforwardly used for the regularized versions using logistic regression likelihood.

To numerically evaluate the parameters estimates, we apply Berman-Turner method <cit.> combined with coordinate descent algorithm  <cit.> to perform variable selection and parameter estimation. Berman-Turner device allows to show that maximizing (<ref>) is equivalent to fitting a weighted Poisson generalized linear model, so the standard software for generalized linear models (GLMs) can be used. This  has been exploited by the    package <cit.>. As we make links between spatial point processes intensity estimation and GLMs, we only have to deal with feature selection procedures for GLMs. Hence, we clearly have many advantages: the various computational strategies are carefully studied, and, in particular, efficiently implemented in . In this study, to compute the regularization path solutions, we employ coordinate descent algorithm <cit.>. This is implemented in the  <cit.> for regularization methods for GLMs using some convex penalties (i.e., ridge, lasso, elastic net, adaptive lasso and adaptive elastic net) and in the  <cit.> for regularization methods for GLMs using some non-convex penalties (i.e., SCAD and MC+). More details for computational strategies are discussed in detail  by <cit.>.

Our methods rely on the tuning parameter λ. Some previous studies  <cit.> suggest to use a modified BIC criterion to select the tuning parameter. We follow the literature and choose λ by minimizing WQBIC(λ), a modified version of the BIC criterion, defined by

    WQBIC(λ)=-2 ℓ(w;β̂ (λ))+ s(λ)log|D|,

where s(λ)=∑_j=1^p 𝕀{β̂_j(λ) ≠ 0} is the number of selected covariates with non-zero regression coefficients and |D| is the volume of observation domain. To implement the adaptive methods (i.e., adaptive lasso and adaptive elastic net), we follow <cit.> and define λ_j=λ / | β̃_̃j̃(ridge)|, j =1, ⋯,p, where β̃(ridge) is the estimates obtained from ridge regression and λ is a tuning parameter chosen by WQBIC(λ) criterion as described above. Following <cit.>, we fix γ=0.5 for elastic net and its adaptive version, γ=3.7 for SCAD, and γ=3 for MC+. For further discussion regarding the selection of γ for SCAD and MC+, see e.g. <cit.> and <cit.>.



 §.§ Simulation study


In this section, we investigate the behavior of our estimators in a simulation experiment in different situations when a large number of covariates for fitting spatial point process intensity estimation is involved. We intend to extend the setting considered by <cit.>. We start with relatively complex situation where strong multicollinearity is present (Scenarios <ref> and <ref>) and we then consider a more complex setting using real datasets (Scenarios <ref> and <ref>). We have two different scenarios (Scenarios <ref> and <ref>) for which the number of true covariates as well as their coefficients are different.

The spatial domain we consider is D=[0,1000] ×[0,500]. The true intensity function has the form ρ(u; β_0)=exp(𝐳(u)^⊤β_0), where 𝐳(u)={1, z_1(u),…,z_50(u)}^⊤ and β_0={β_0, β_01, ⋯, β_050}. We set  β_0 such that the mean number of points over D is equal to 1600. We consider two different scenarios described as follows.


  * We define the true vector β_0={β_0, 2, 0.75, 0, ⋯,0}. To define the covariates, we center and scale the 201 × 101 pixel images of elevation (x_1) and gradient of elevation (x_2) contained in the  datasets of  library in  
 and use them as two true covariates. In addition, we create two settings to define extra covariates: 


  * First, we generate 48 201 × 101 pixel images of covariates as a standard Gaussian white noise and denote them by x_3, …, x_50. Second, we transform them, together with x_1 and x_2, to have multicollinearity. In particular, we define 𝐳̃(u)=𝐕 ^⊤𝐱(u), where 𝐱(u)={x_1(u),…,x_50(u)}^⊤. More precisely, 𝐕 is such that Ω=𝐕 ^⊤𝐕, and (Ω)_ij=(Ω)_ji=0.7^|i-j| for i,j=1,…,50, except (Ω)_12=(Ω)_21=0, to preserve the correlation between x_1 and x_2.    
In this setting, 𝐳(u)={1,𝐳̃(u)}.


  * We center and scale the 13 50 × 25 pixel images of soil nutrients obtained from the study in tropical forest of Barro Colorado Island (BCI) in central Panama <cit.> and convert them to be  201 × 101 pixel images as x_1 and x_2. In addition, we consider the interaction between two soil nutrients such that we have 50 covariates in total. We use 48 covariates (13 soil nutrients and 35 interactions between them) as the extra covariates. Together with x_1 and x_2, we keep the structure of the covariance matrix to preserve the complexity of the situation. In this setting, we have 𝐳(u)=𝐱(u)={1, x_1(u),…,x_50(u)}^⊤.    



  * In this setting, we consider five true covariates out of 50 covariates. In addition of elevation (x_1) and gradient of elevation (x_2), we convert 50 × 25 pixel images of concentration of Aluminium (x_3), Boron (x_4) and Calcium (x_5) in the soil to be 201 × 101 pixel images as x_1 and x_2 and set them to be other three true covariates. All five covariates are centered and scaled. We define the true coefficient vector β_0={β_0, 5, 4,  3, 2, 1, 0, ⋯,0}. As in Scenario <ref>, we make two settings to define 45 extra covariates: 


  * This setting is similar to that of Scenario <ref>. We generate 45 201 × 101 pixel images of covariates as standard Gaussian white noise, denote them by x_6, …, x_50, and define 𝐳̃(u)=𝐕 ^⊤𝐱(u), where 𝐕 is such that Ω=𝐕 ^⊤𝐕, and (Ω)_ij=(Ω)_ji=0.7^|i-j| for i,j=1,…,50, except (Ω)_kl=(Ω)_lk=0, for k,l=1,⋯,5, k ≠ l, to preserve the correlation among x_1 - x_5. We still define 𝐳(u)={1,𝐳̃(u)}. 


  * We use the real dataset as in Scenario <ref> and consider similar setting. In this setting, we define 5 true covariates which have different regression coefficients as in  Scenario <ref>. 







With these scenarios, we simulate 2000 spatial point patterns from a Thomas point process using the 𝚛𝚃𝚑𝚘𝚖𝚊𝚜 function in the 𝚜𝚙𝚊𝚝𝚜𝚝𝚊𝚝 package. We set the interaction parameter κ to be κ=5 × 10^-4, κ=5 × 10^-5 and let ω=20. Briefly, smaller values of ω correspond to tighter clusters, and smaller values of κ correspond to a fewer number of parents <cit.>. For each scenario with different κ, we fit the intensity to the simulated point pattern realizations.

We report the performances of our estimates in terms of two characteristics: selection and prediction properties. We present the selection properties in Table <ref> and the prediction properties in Table <ref>







To evaluate the selection properties of the estimates, we consider the true positive rate (TPR), the false positive rate (FPR), and the positive predictive value (PPV). We want to find the methods which have a TPR close to 100% meaning that it can select correctly all the true covariates, a FPR close to 0 showing that it can remove all the extra covariates from the model, and a PPV close to 100% indicating that, for Scenario <ref> (resp. Scenario <ref>), it can keep exactly the two (resp. five) true covariates and remove all the 48 (resp. 45) extra covariates. In general, for both regularized PL and regularized WPL, the best selection properties are obtained from larger κ (5 × 10^-4) which indicates weaker spatial dependence. To compare the regularization methods, we emphasize here that the main difference between regularization methods which satisfy (adaptive lasso, adaptive elastic net, SCAD, and MC+) and which cannot satisfy (lasso, elastic net) our theorems is that the methods which cannot satisfy our theorems tend to over-select covariates, leading to suffering from larger FPR and smaller PPV in general. Among all regularization methods considered in this study, adaptive lasso and adaptive elastic net seem to outperform the other methods in most cases. Although adaptive lasso and adaptive elastic net perform quite similarly, the adaptive lasso is slightly better.

In this simulation study, we are still able to show that even when the strong multicollinearity exists such as in Scenario <ref>, our proposed methods work well for the penalization methods satisfying our theorems. However, as probably expected, our methods are getting difficult to distinguish between the important and the noisy covariates as the setting becomes more and more complex. In the experiments we conduct, we find that the regularized PL and WPL (with adaptive lasso) perform quite similar for the easiest (Scenario <ref>) and the toughest (Scenario <ref>) setting. For Scenarios <ref> and <ref>, the regularized WPL with adaptive lasso seems to be more favorable. From Table <ref>, we would recommend in general to combine the regularized WPL with the adaptive lasso to perform variable selection.

Table <ref> gives the prediction properties of the estimates (except for β_0 which is excluded) in terms of biases, standard deviations (SD), and square root of mean squared errors (RMSE), some criteria we define by

    Bias   = [ ∑_j=1^50{𝔼̂(β̂_j)-β_0j}^2 ]^1/2,
    SD= [ ∑_j=1^50σ̂_j^2 ]^1/2,
    RMSE= [ ∑_j=1^50𝔼̂(β̂_j-β_0j)^2 ]^1/2,

where 𝔼̂(β̂_j) and σ̂_j^2 are respectively the empirical mean and variance of the estimates β̂_j, for j=1,…,50.

In general, the properties improve with larger κ due to weaker spatial dependence. Regarding the regularization methods considered in this study, adaptive lasso and adaptive elastic net perform best. Adaptive elastic net becomes more preferable than adaptive lasso for a clustered process (κ= 5 × 10^-5) and for a structured spatial data (Scenarios <ref> and <ref>). The  adaptive elastic net is more efficient than the adaptive lasso especially in the complex situation: large number of covariates, strong multicollinearity, clustered processes, and complex spatial structure due to the advantage of combining l_1 and l_2 penalties.

By employing regularized WPL, we have potentially more efficient estimates than that of the regularized PL, especially for the more clustered process. However, this does not mean that the regularized WPL is able to improve the RMSE since it usually introduces extra biases. Regularized WPL seems more appropriate for the case having covariates with complex spatial structure (Scenarios <ref> and <ref>). Otherwise, regularized PL is more favorable. From Table <ref>, when the focus is on prediction, we would recommend to apply adaptive elastic net as a general advice, and we would combine with regularized WPL if the covariates have complex spatial structure (e.g. Scenarios <ref> and <ref>) or combine with regularized PL if there is no evidence of complex spatial structure in the covariates (e.g. Scenarios <ref> and <ref>).

Note that, from Table <ref>, the adaptive lasso is more preferable if the focus is on variable selection while, from Table <ref>, the adaptive elastic net is more favorable if the focus is for prediction. To have a more general recommendation, we would recommend applying adaptive elastic net when we are faced with a complex situation: a large number of covariates, strong multicollinearity, clustered processes and complex spatial structure. By combining l_1 and l_2 penalties, the adaptive elastic net provides a nice balance between selection and  prediction properties. This is why in most complex cases (Scenario <ref> with κ=5 × 10^-5), adaptive elastic net decides to choose more covariates than adaptive lasso (which includes true and noisy covariates) to suffer from slightly less appropriate properties for the selection performance but to be able to improve significantly the prediction properties.




 §.§ Application to forestry datasets


We now consider the study of ecology in a tropical rainforest in Barro Corrolado Island (BCI), Panama, described previously in Section <ref>. In particular, we are interested  in studying the spatial distribution of 3,604 locations of Beilschmiedia pendula Lauraceae (BPL) trees by modeling its intensity as a log-linear function of 93 covariates consisting of 2 topological attributes, 13 soil properties, and 78 interactions between two soil nutrients.Regarding the relatively large number of covariates, we apply our proposed methods to select few covariates among them and estimate their coefficients. In particular, we use the regularized Poisson methods with the lasso, adaptive lasso, and adaptive elastic net. Note that we center and scale all the covariates to observe which covariates owing relatively large effect on the intensity.



 







We present in Table <ref> the number of selected and non-selected covariates by each method. Out of 93 covariates, more than 50% from the total number of covariates are selected by regularized PL while much fewer covariates are selected by regularized WPL. The regularized PL seems to overfit the model. 

Regarding lasso method, 77 covariates are selected by regularized PL method while 20 covariates are selected by regularized WPL. Compared to the two adaptive methods (i.e., adaptive lasso and adaptive elastic net), lasso tends to keep less important covariates. This may explain why lasso cannot satisfy our Theorem <ref>. In terms of selection properties, adaptive lasso and adaptive elastic net perform similarly when regularized WPL is applied.

Table <ref> gives the information regarding nine covariates commonly selected among the six methods. Although the magnitudes of the estimates can be slightly different, the signs all agree with each other.

These results suggest that BPL trees favor to live in the areas of higher elevation and slope with a high concentration of Copper and Manganese in the soil. Furthermore, BPL trees prefer to live in the areas with lower concentration levels of Phosphorus and Zinc in the soil. The interaction between Aluminum and Phosphorus gives a negative association with the appearance of BPL trees while the interaction between Magnesium and Phosphorus and the interaction between Nitrogen mineralization and pH show a positive association with the occurrence of BPL trees. The maps of 3,604 locations of BPL trees, as well as the nine commonly selected covariates, are depicted in Figure <ref>.


















§ CONCLUSION AND DISCUSSION


We consider feature selection techniques for spatial point processes intensity estimation by regularizing  estimating equations derived from Poisson and logistic regression likelihoods in a setting where the number of parameters diverges as the volume of observation domain increases. Under some conditions, we prove that the estimates obtained from such setting satisfy consistency, sparsity, and asymptotic normality. Our results are available for large classes of spatial point processes and for many penalty functions.

We conduct  simulation experiments to evaluate the finite sample properties of the regularized Poisson estimator and regularized weighted Poisson estimator. From the results, we would recommend in general the combination between regularized WPL and adaptive lasso if the concern is on variable selection. Furthermore, when the focus is on prediction, the regularized WPL combined with the adaptive elastic net is more preferable for the situation where there is a complex spatial structure in the covariates. For more general advice, we would recommend using the adaptive elastic net rather than the adaptive lasso since the adaptive elastic net is able to balance the selection and the prediction properties by combining the l_1 and the l_2 penalties.



To implement our methods, we combine the   package and the two  packages  and  dealing with penalized generalized linear models. This results in a computationally fast procedure even when the number of covariates is large. It is worth noticing that, as other regularization methods, our methods also rely on the selection of the tuning parameter. As the study in a classical regression analysis, the BIC-type methods are proposed to obtain selection consistent estimator <cit.>. We have numerical evidence from simulation studies that this criterion can satisfy the selection consistency.  Theoretical justification in this spatial point process framework is the purpose of a future research.

We apply our methods to the Barro Corrolado Island study to estimate the intensity of Beilschmiedia pendula Lauraceae (BPL) tree as a log-linear function of 93 environmental covariates. Regularized weighted Poisson likelihood combined with adaptive elastic net performs similarly to adaptive lasso. Among 93 covariates, we find nine spatial covariates which may have a high influence to the appearance of BPL trees, including two topological attributes: elevation and slope, four soil nutrients: Copper, Manganese, Phosphorus and Zinc, and three interaction between two soil properties: the interaction between  Aluminum and Phosphorus, between Magnesium and Phosphorus, and between Nitrogen mineralisation and pH.

A further work would consider to include the 296 other species of trees, which were surveyed in the same observation region, to study the existence of any competition between BPL and other species of trees in the forest. In such a situation, the methods used in this study may face some computational issues. The Dantzig selector <cit.> might be a good alternative since the implementation for linear models (and generalized linear models) results in a linear programming. Thus, more competitive algorithms are available. It would be interesting to bring this approach to spatial point process framework.

14pt
Acknowledgements
4pt

We thank A. L. Thurman who kindly shared the  code used for the simulation study in <cit.> and P. Breheny who kindly provided his code used in   package. We also thank R. Drouilhet for technical help. The research of A. Choiruddin is supported by The Danish Council for Independent Research – Natural Sciences, grant DFF – 7014-00074 "Statistics for point processes in space and beyond", and by the ”Centre for Stochastic Geometry and Advanced Bioimaging”, funded by grant 8721 from the Villum Foundation. The research of J.-F. Coeurjolly is supported by the Natural Sciences and Engineering Research Council of Canada. The research of F. Letué is supported by ANR-11-LABX-0025 Persyval-lab (project Persyvact2).

The BCI soils data sets were collected and analyzed by J. Dalling, R. John, K. Harms, R. Stallard and J. Yavitt with support from NSF DEB021104,021115, 0212284,0212818 and OISE 0314581, and STRI Soils Initiative and CTFS and assistance from P. Segre and J. Trani. Datasets are available at the CTFS website 
 .



plainnat










§ AUXILIARY LEMMA
  

The following lemma is used in the proof of Theorem <ref> and Lemma <ref> (which includes Lemma <ref> and Theorem <ref>). Throughout the proofs, the notation 𝐗_n = O_P (x_n) or 𝐗_n = o_P (x_n) for a random vector 𝐗_n and a sequence of real numbers x_n means that 𝐗_n=O_P(x_n) and 𝐗_n=o_P(x_n). In the same way for a vector 𝐕_n or a squared matrix 𝐌_n, the notation 𝐕_n=O(x_n) and  𝐌_n=O(x_n) mean that 𝐕_n=O(x_n) and  𝐌_n=O(x_n).

 
Under conditions (𝒞.<ref>)-(𝒞.<ref>), the following result holds as n→∞
	
    ℓ^(1)_n(w;β_0)  =O_P( √( p_n |D_n| ))  .




Using Campbell Theorems (<ref>)-(<ref>), the score vector ℓ^(1)_n(w;β_0) has variance

    Var[ ℓ^(1)_n(w; β_0)]= 𝐁_n(w; β_0)+𝐂_n(w; β_0).

Conditions (𝒞.<ref>)-(𝒞.<ref>) allow us to obtain that sup_u ∈ℝ^d∫_ℝ^d{g(u,v)-1}dv < ⁠∞. We then deduce using conditions (𝒞.<ref>)-(𝒞.<ref>) that

    𝐁_n(w; β_0)+𝐂_n(w; β_0)=O(p_n |D_n|).

The result is proved since for any centered real-valued stochastic process Y_n with finite variance Var[Y_n], Y_n=O_P(√(Var[Y_n])).









§ PROOF OF THEOREM <REF>
 

In the proof of this result and the following ones, the notation κ stands for a generic constant which may vary from line to line. In particular this constant is independent of n, β_0 and 𝐤.


Let d_n = √(p)_n(|D_n|^-1/2+a_n), and 𝐤={k_1, k_2, …, k_p_n}^⊤. We remind the reader that the estimate of β_0 is defined as the maximum of the function Q_n (given by (<ref>)) over Θ, an open convex bounded set of ℝ^p_n for any n ≥ 1. For any 𝐤 such that 𝐤≤ K<∞, β_0 + d_n 𝐤∈Θ for n  sufficiently large. Assume this is valid in the following.
To prove Theorem <ref>, we aim at proving that for any given ϵ>0, there exists sufficiently large K>0 such that for n sufficiently large 

    P(sup_𝐤 = KΔ_n(𝐤)>0)≤ϵ,
      Δ_n(𝐤) = Q_n(w;β_0+d_n𝐤)-Q_n(w;β_0).

Equation (<ref>) will imply that with probability at least 1-ϵ, there exists a local maximum in the ball {β_0+d_n𝐤:𝐤≤ K}, and therefore  a local maximizer β̂ is such that β̂-β_0=O_P(d_n).  We decompose Δ_n(𝐤) as Δ_n(𝐤)= T_1+T_2 where

    T_1     = ℓ_n(w;β_0+d_n𝐤)-ℓ_n(w; β_0) 
    
    	T_2     = |D_n|∑_j=1^p_n( p_λ_n,j(|β_0j|)- p_λ_n,j(|β_0j+d_nk_j|)).

Since ρ(u;·) is infinitely continuously differentiable and ℓ_n^(2)(w;β) =-𝐀_n(w;β), then using a second-order Taylor expansion there exists t∈ (0,1) such that

    T_1 =     d_n 𝐤^⊤ℓ_n^(1)(w;β_0) - 1/2d_n^2𝐤^⊤𝐀_n(w;β_0) 𝐤
       + 1/2d_n^2𝐤^⊤( 𝐀_n(w;β_0) -𝐀_n(w;β_0 + td_n 𝐤) ) 𝐤 .

By conditions (𝒞.<ref>)-(𝒞.<ref>), there exists a non-negative constant κ such that

    1/2𝐀_n(w;β_0) -𝐀_n(w;β_0 + td_n 𝐤) ≤κ d_n |D_n| p_n.

Now, denote ν̌:= lim inf_n→∞ν_min(|D_n|^-1𝐀_n(w;β_0)). By condition (𝒞.<ref>), we have  that for any 𝐤

    0<ν̌≤𝐤^⊤( |D_n|^-1𝐀_n(w;β_0)) 𝐤/𝐤^2.

Therefore, we have

    T_1 ≤ d_n ℓ_n^(1)(w;β_0) 𝐤  - ν̌/2 d_n^2 |D_n| 𝐤^2 +κ p_n d_n^3  |D_n|𝐤^2.

Now by the condition (𝒞.<ref>) and by assumption that a_n=O(|D_n|^-1/2), we obtain p_n d_n  = o(1), so κ p_n d_n^3  |D_n|𝐤^2 = o(1) d_n^2 |D_n|𝐤^2. Hence, for n sufficiently large

    T_1 ≤ d_n ℓ_n^(1)(w;β_0) 𝐤  - ν̌/4 d_n^2 |D_n| 𝐤^2.


Regarding the term T_2,

    T_2≤ T_2^' := |D_n|∑_j=1^s( p_λ_n,j(|β_0j|)- p_λ_n,j(|β_0j+d_nk_j|))

since for any j the penalty function p_λ_n,j is non-negative and p_λ_n,j(|β_0j|)=0 for j=s+1,…,p_n.

From (𝒞.<ref>), for n sufficiently large, p_λ_n,j is twice continuously differentiable for every β_j = β_0j+t d_n k_j with t∈ (0,1). Therefore using a third-order Taylor expansion, there exist t_j ∈ (0,1), j=1,…,s such that -T_2^'= T_2,1^' + T_2,2^' + T_2,3^', where

    T_2,1^'   =d_n|D_n|∑_j=1^s k_j p_λ_n,j^'(|β_0j|) (β_0,j)  ≤√(s) a_n d_n |D_n|  𝐤≤  d_n^2 |D_n|  𝐤, 
    
    T_2,2^'   = 1/2 d_n^2 |D_n|∑_j=1^s k_j^2 p^''_λ_n,j  (|β_0j|) ≤  c_n d_n^2|D_n|  𝐤^2, 
    
    T_2,3^'   = 1/6 d_n^3|D_n| ∑_j=1^s k_j^3 p^'''_λ_n,j  (|β_0j+t_j d_n k_j|) ≤κ d_n^3 |D_n|.

The three inequalities above are obtained using the definitions of a_n and c_n, condition (𝒞.<ref>) and Cauchy-Schwarz inequality. We deduce that for n sufficiently large

    T_2 ≤ | T_2^'|    ≤ 2 d_n^2 |D_n| 𝐤,

and then

    Δ_n(𝐤) ≤ d_n ℓ_n^(1)(w;β_0) 𝐤 -ν̌/4 d_n^2 |D_n| 𝐤^2 + 2 d_n^2 |D_n| 𝐤.

We now return to (<ref>): for n sufficiently large

    P(sup_𝐤= KΔ_n(𝐤)>0) ≤P(
    ℓ_n^(1)(w;β_0) > ν̌/4 d_n|D_n| K - 2 d_n |D_n|
     ).

Since d_n |D_n|=O(√(p_n |D_n|)), by choosing K large enough, there exists κ such that for n  sufficiently large

    P( sup_𝐤=KΔ_n(𝐤) >0) ≤P( ℓ_n^(1)(w;β_0) >κ√(p_n |D_n|)) ≤ϵ

for any given ϵ>0 from (<ref>) in Lemma <ref>.













§ PROOF OF THEOREM <REF>
 
Before proving Theorem <ref>, we present Lemmas <ref>-<ref>. Lemma <ref> is used to prove Theorem <ref>(i) while  Lemma <ref> is used to derive Theorem <ref>(ii).


Assume the conditions (𝒞.<ref>)-(𝒞.<ref>) hold. If a_n=O(|D_n|^-1/2) and b_n √(|D_n|/p_n^2)→∞ as n→∞, then with probability tending to 1, for any β_1 satisfying β_1 - β_01=O_P(√(p_n/|D_n|)), and for any constant K_1 > 0,

    Q_n(w;(β_1^⊤,0^⊤)^⊤)
    = max_β_2≤ K_1 √(p_n/|D_n|)
    Q_n(w;(β_1^⊤,β_2^⊤)^⊤).



Let ε_n= K_1 √(p_n/|D_n|). It is sufficient to show that with probability tending to 1 as n→∞, for any β_1 satisfying β_1 -β_01=O_P(√(p_n/|D_n|)), we have for any j=s+1, …, p_n


    ∂ Q_n(w;β)/∂β_j<0    0<β_j<ε_n,



    ∂ Q_n(w;β)/∂β_j>0    -ε_n<β_j<0.

From (<ref>),

    ∂ℓ_n(w;β)/∂β_j = ∂ℓ_n(w;β_0)/∂β_j + R_n,

where R_n = ∫_D_n w(u) z_j(u)(ρ(u;β)-ρ(u;β_0)) du. Using similar arguments used in the proof of Lemma <ref>, we can prove that

    ∂ℓ_n(w;β_0)/∂β_j = O_P(√(|D_n|)).

Let u ∈ℝ^d. By Taylor expansion, there exists t∈ (0,1), such that 

    ρ(u;β) = ρ(u;β_0) + (β-β_0)^⊤𝐳(u) ρ(u;β_0 + t(β-β_0 )).

For n sufficiently large, β_0 + t(β-β_0 ) ∈Ξ(β_0) defined in condition (𝒞.<ref>). Therefore, for n sufficiently large, we have by Cauchy-Schwarz inequality and conditions (𝒞.<ref>)-(𝒞.<ref>)

    |R_n| ≤κ∫_D_nβ-β_0 𝐳(u)du = O_P(√(|D_n| p_n^2)).

We therefore deduce that for any j=s+1, …, p_n

    ∂ℓ_n(w;β)/∂β_j = O_P(√(|D_n| p_n ^2) ) .


Now, we want to prove (<ref>). Let 0<β_j<ε_n and b_n be the sequence given by (<ref>). By condition (𝒞.<ref>), b_n is well-defined and since by the assumption b_n √(|D_n|/p_n^2)→∞, in particular,  b_n>0 for n sufficiently large. Therefore, for n sufficiently large,

    P ( ∂ Q_n(w;β)/∂β_j<0 )   =P ( ∂ℓ_n(w;β)/∂β_j - |D_n|p'_λ_n,j(|β_j|)(β_j)<0 )
       =P ( ∂ℓ_n(w;β)/∂β_j< |D_n|p'_λ_n,j(|β_j|) )
       ≥P ( ∂ℓ_n(w;β)/∂β_j< |D_n|b_n )
       = P ( ∂ℓ_n(w;β)/∂β_j< √(|D_n|p_n^2) √(|D_n|/p_n^2)b_n ).

The assertion (<ref>) is therefore deduced from (<ref>) and from the assumption that b_n √(|D_n|/p_n^2)→∞ as n  →∞. We proceed similarly to prove (<ref>).






 
Under the conditions (𝒞.<ref>)-(𝒞.<ref>) and the conditions required in Lemma <ref>, the following convergence holds in distribution as n→∞
	
    {𝐁_n,11(w; β_01)+𝐂_n,11(w; β_01)}^-1/2ℓ^(1)_n,1(w; β_01) 𝒩(0,𝐈_s) ,

where ℓ^(1)_n,1(w;β_0) is the first s components of ℓ^(1)_n(w;β_0) and 𝐁_n,11(w;β_0) (resp. 𝐂_n,11(w;β_0)) is the s × s top-left corner of 𝐁_n(w;β_0) (resp  𝐂_n(w;β_0)).



By Lemma <ref> and by using Campbell Theorems (<ref>)-(<ref>),

    Var[ ℓ^(1)_n,1(w; β_0)]= 𝐁_n,11(w; β_0)+𝐂_n,11(w; β_0).

The remainder of the proof follows <cit.>. Let C_i=i+(-1/2,1/2]^d be the unit box centered at i ∈ℤ^d and define ℐ_n={i ∈ℤ^d, C_i ∩ D_n ≠∅}. Set D_n=⋃_i ∈ℐ_n C_i,n, where C_i,n=C_i ∩ D_n. We have

    ℓ^(1)_n,1(w;β_0)=∑_i ∈ℐ_n Y_i,n

where

    Y_i,n=∑_u ∈𝐗∩ C_i,n w(u)𝐳_01(u) - ∫_ C_i,n w(u)𝐳_01(u)exp(β_01^⊤𝐳_01(u) )du.

For any n ≥ 1 and any i ∈ℐ_n, Y_i,n has zero mean, and by condition (𝒞.<ref>),

    sup_n ≥ 1sup_i ∈ℐ_n𝔼(Y_i,n^2+δ) < ∞.


If we combine (<ref>) with conditions  (𝒞.<ref>)-(𝒞.<ref>), we can apply <cit.>, a central limit theorem for triangular arrays of random fields.










 We now focus on the proof of Theorem <ref>. Since Theorem <ref>(i) is proved by Lemma <ref>, we only need to prove Theorem <ref>(ii), which is the asymptotic normality of β̂_1. As shown in Theorem <ref>, there is a root-(|D_n|/p_n) consistent local maximizer β̂ of Q_n(w;β), and it can be shown that there exists an estimator β̂_1 in Theorem <ref> that is a root-(|D_n|/p_n) consistent local maximizer of Q_n (w;(β_1^⊤,0^⊤)^⊤), which is regarded as a function of  β_1, and that satisfies

    ∂ Q_n(w;β̂)/∂β_j=0    j=1,…,s β̂=( β̂_1^⊤,0^ ⊤)^⊤.

There exists t∈ (0,1) and β̃= β̂ + t(β_0-β̂)  such that for j=1,⋯,s

    0
    =   ∂ℓ_n(w;β̂)/∂β_j-|D_n|p'_λ_n,j(|β̂_j|)(β̂_j) 
    
    =   ∂ℓ_n(w;β_0)/∂β_j+∑_l=1^s∂^2 ℓ_n(w; β̃)/∂β_j ∂β_l(β̂_l-β_0l)-|D_n|p'_λ_n,j(|β̂_j|)(β̂_j) 
    
    =   ∂ℓ_n(w;β_0)/∂β_j+∑_l=1^s∂^2 ℓ_n(w; β_0)/∂β_j ∂β_l(β̂_l-β_0l)+∑_l=1^sΨ_n,jl(β̂_l-β_0l)
       -|D_n|p'_λ_n,j(|β_0j|)(β_0j)-|D_n|ϕ_n,j,

where 

    Ψ_n,jl=∂^2 ℓ_n(w;β̃)/∂β_j ∂β_l-∂^2 ℓ_n(w;β_0)/∂β_j ∂β_l

and ϕ_n,j=p'_λ_n,j(|β̂_j|)(β̂_j)-p'_λ_n,j(|β_0j|)(β_0j). Since p'_λ is a Lipschitz function by condition (𝒞.<ref>), there exists κ≥ 0 such that by condition on a_n

    ϕ_n,j    =p'_λ_n,j(|β̂_j|)(β̂_j)-p'_λ_n,j(|β_0j|)(β_0j) 
        = (p'_λ_n,j(|β̂_j|)-p'_λ_n,j(|β_0j|)) (β̂_j) + p'_λ_n,j(|β_0j|) ((β̂_j)-(β_0j)) 
       ≤κ| |β̂_j| - |β_0j| | + 2a_n 
       ≤κ | β̂_j - β_0j | + 2a_n .

We now decompose ϕ_n,j as ϕ_n,j=T_1+T_2 where

    T_1=ϕ_n,j𝕀(|β̂_j- β_0j| ≤r̃_n,j) 
    	     
    	T_2=ϕ_n,j𝕀(|β̂_j-β_0j| > r̃_n,j)

and where r̃_n,j is the sequence defined in the condition (𝒞.<ref>). Under this condition, the following Taylor expansion can be derived for the term T_1: there exists t∈ (0,1) and β̌_j= β̂_j+t(β_0j-β̂_j) such that 

    T_1     =p”_λ_n,j(|β_0j|)(β̂_j- β_0j) 𝕀(|β̂_j- β_0j| ≤r̃_n,j) 
          + 1/2 (β̂_j - β_0j)^2 p”'_λ_n,j(|β̃_j|) sign(β̌_j) 𝕀(|β̂_j- β_0j| ≤r̃_n,j) 
       =p”_λ_n,j(|β_0j|)(β̂_j- β_0j) 𝕀(|β̂_j- β_0j| ≤r̃_n,j) + O_P(p_n/|D_n|)

where the latter equation ensues from Theorem <ref> and condition (𝒞.<ref>). Again,
from Theorem <ref>,  𝕀(|β̂_j- β_0j| ≤r̃_n,j)  1  which implies that 𝕀(|β̂_j- β_0j| ≤r̃_n,j)  1,  so T_1=p”_λ_n,j(|β_0j|)(β̂_j- β_0j) (1+o_P(1))+  O_P(p_n/|D_n|).

Regarding the term T_2, we have by (<ref>)

    T_2    ≤{κ | β̂_j - β_0j | + 2a_n } 𝕀(|β̂_j-β_0j| > r̃_n,j) 
       = κ | β̂_j - β_0j |  𝕀(|β̂_j-β_0j| > r̃_n,j) + o(|D_n|^-1/2).

We want to prove that T_2 = o_P(|D_n|^-1/2). Define S_n=|β̂_j-β_0j|  𝕀(|β̂_j-β_0j| > r̃_n,j) and T_n=𝕀(S_n> δ|D_n|^-1/2) for some δ>0. We claim that the result is proved if we prove that 𝔼T_n → 0 for any δ>0.  Condition (𝒞.<ref>) implies in particular that for n large enough, r̃_n,j > √(p_n/|D_n|) > √(1/|D_n|). Using this, it can be checked that the binary random variable T_n reduces to T_n=𝕀(|β̂_j-β_0j| > r̃_n,j)L^1→0 as n→∞.

Then, we  deduce that 

    ϕ_n,j = p”_λ_n,j(|β_0j|)(β̂_j- β_0j) (1+o_P(1))+ O_P(p_n/|D_n|) + o_P(|D_n|^-1/2).

	

Let ℓ^(1)_n,1(w;β_0) (resp. ℓ^(2)_n,1(w;β_0)) be the first s components (resp. s × s top-left corner) of ℓ^(1)_n(w;β_0) (resp. ℓ^(2)_n(w;β_0)). Let also Ψ_n be the s × s matrix containing Ψ_n,jl, j,l=1,…,s. Finally, let the vector 𝐩'_n, the vector ϕ_n and the s × s matrix 𝐌_n be defined by

    𝐩'_n   ={p'_λ_n,1(|β_01|)(β_01),…,p'_λ_n,s(|β_0s|)(β_0s)}^⊤, 
    ϕ_n   ={ϕ_n,1,…,ϕ_n,s}^⊤, 
    𝐌_n   ={𝐁_n,11(w; β_0)+𝐂_n,11(w; β_0)}^-1/2.

We rewrite both sides of (<ref>) as

    ℓ^(1)_n,1(w;β_0)+ℓ^(2)_n,1(w;β_0)(β̂_1-β_01)+ Ψ_n (β̂_1-β_01) -|D_n| 𝐩'_n-|D_n| ϕ_n   =0.

By definition of Π_n given by (<ref>) and from (<ref>), we obtain ϕ_n=Π_n (β̂_1-β_01)(1+o_P(1))+ O_P(p_n/|D_n|) +  o_P(|D_n|^-1/2). Using this, we deduce, by premultiplying both sides of (<ref>) by 𝐌_n, that

    𝐌_n ℓ^(1)_n,1(w;β_0)-   𝐌_n (𝐀_n,11(w;β_0)+ |D_n| Π_n)(β̂_1-β_01) 
        =O(|D_n|  𝐌_n 𝐩'_n ) + o_P(|D_n|  𝐌_nΠ_n (β̂_1-β_01) )  
         + O_P (𝐌_n  p_n) +  o_P(𝐌_n  |D_n|^1/2) 
         + O_P (𝐌_n Ψ_n (β̂_1-β_01) ) .



Now, 𝐌_n=O(1/ √(|D_n|))  by condition (𝒞.<ref>), Ψ_n=O_P(√(p_n |D_n|))  by conditions (𝒞.<ref>)-(𝒞.<ref>) and Theorem <ref>, and β̂_1-β_01=O_P(√(p_n/|D_n|)) by Theorem <ref> and Theorem <ref>(i). Finally, since by assumptions that a_n √(|D_n|)→ 0 and c_n √(p_n)→ 0 as n →∞, we deduce that

    |D_n|  𝐌_n 𝐩'_n   = O(a_n √(D_n)) =o(1),
    
    |D_n|  𝐌_nΠ_n (β̂_1-β_01)    = O_P(√(|D_n|) c_n √(p_n/|D_n|)) =o_P(1)  ,
    𝐌_n  √(|D_n|)   =  O(1),
    𝐌_n   p_n    =  O( √(p_n^2/|D_n|))=o(1),
    𝐌_n Ψ_n (β̂_1-β_01)   =O_P( √(p_n^2/|D_n|))=o_P(1).

The last two lines are obtained from (𝒞.<ref>). Therefore, we have that

    𝐌_n ℓ^(1)_n,1(w;β_0)-𝐌_n (𝐀_n,11(w;β_0)+ |D_n| Π_n)(β̂_1-β_01) =o_P(1).

By (<ref>) in Lemma <ref> and by Slutsky's Theorem, we deduce that

    {𝐁_n,11(w; β_0)+𝐂_n,11(w; β_0)}^-1/2×
    {𝐀_n,11(w;β_0)+|D_n| Π_n}(β̂_1-β_01)   𝒩(0,𝐈_s)

as n →∞, which can be rewritten, in particular under (𝒞.<ref>), as 

    |D_n|^1/2Σ_n(w;β_0)^-1/2(β̂_1-β_01)𝒩(0,𝐈_s)

where Σ_n(w,β_0) is given by (<ref>).















=1.7pc
=2pt
914pt plus.8pt minus .6pt

natexlab#1#1url<#>1urlprefixURL



.65cm
Department of Mathematical Sciences, Aalborg University, Denmark
1pt
E-mail: achmad@math.aau.dk
6pt

Department of Mathematics, Université du Québec à Montréal (UQAM), Canada and

Department of Probability and Statistics, Université Grenoble Alpes, France
1pt
E-mail: coeurjolly.jean-francois@uqam.ca
6pt

Department of Probability and Statistics, Université Grenoble Alpes, France
1pt
E-mail: frederique.letue@univ-grenoble-alpes.fr 






