Convex and non-convex regularization methods for spatial point processes intensity estimation

TL;DR

This paper introduces regularization techniques for spatial point process intensity estimation, ensuring consistency and sparsity, with practical implementation and application to forestry data.

Contribution

It develops new regularized estimation methods for spatial point processes, providing theoretical guarantees and demonstrating their effectiveness through simulations and real data.

Findings

Methods achieve consistency and sparsity.

Numerical implementation is feasible and effective.

Application to forestry data demonstrates practical utility.

Abstract

This paper deals with feature selection procedures for spatial point processes intensity estimation. We consider regularized versions of estimating equations based on Campbell theorem derived from two classical functions: Poisson likelihood and logistic regression likelihood. We provide general conditions on the spatial point processes and on penalty functions which ensure consistency, sparsity and asymptotic normality. We discuss the numerical implementation and assess finite sample properties in a simulation study. Finally, an application to tropical forestry datasets illustrates the use of the proposed methods.

Tables

Table 1. Table 1: The first and the second derivatives of several penalty functions.

Penalty	$p_{\lambda }^{\prime }\left(\theta \right)$	$p_{\lambda }^{\prime \prime }\left(\theta \right)$
${\ell }_2$	$\lambda \theta$	$\lambda$
${\ell }_1$	$\lambda$	$0$
Elastic net	$\lambda \left\{\left(1-\gamma \right)\theta +\gamma \right\}$	$\lambda \left(1-\gamma \right)$
SCAD
MC+

Table 2. Table 2: Details of some regularization methods.

Method	${\sum }_{j=1}^p{p}_{{\lambda }_j}\left(\left|{\beta }_j\right|\right)$
Ridge	${\sum }_{j=1}^p\frac{1}{2}\lambda {\beta }_j^2$
Lasso	${\sum }_{j=1}^p\lambda \left|{\beta }_j\right|$
Enet*	${\sum }_{j=1}^p\lambda \left\{\gamma \right|{\beta }_j|+\frac{1}{2}\left(1-\gamma \right){\beta }_j^2\}$
AL*	${\sum }_{j=1}^p{\lambda }_j\left|{\beta }_j\right|$
Aenet*	${\sum }_{j=1}^p{\lambda }_j\left\{\gamma \right|{\beta }_j|+\frac{1}{2}\left(1-\gamma \right){\beta }_j^2\}$
SCAD
MC+	${\sum }_{j=1}^p\left\{\left(\lambda \left|{\beta }_j\right|-\frac{{\beta }_j^2}{2\gamma }\right)𝕀\left(\left|{\beta }_j\right|\le \gamma \lambda \right)+\frac{1}{2}\gamma {\lambda }^2𝕀\left(\left|{\beta }_j\right|\ge \gamma \lambda \right)\right\}$

Table 3. Table 3: Details of the sequences a n subscript 𝑎 𝑛 a_{n} , b n subscript 𝑏 𝑛 b_{n} and c n subscript 𝑐 𝑛 c_{n} for a given regularization method.

Method	$a_n$	$b_n$	$c_n$
Ridge	${\lambda }_n\underset{j=1,\dots s}{\max }\left\{\left|{\beta }_{0j}\right|\right\}$	$0$	${\lambda }_n$
Lasso	${\lambda }_n$	${\lambda }_n$	0
Enet	${\lambda }_n\left[\left(1-\gamma \right)\underset{j=1,\dots s}{\max }\left\{\left|{\beta }_{0j}\right|\right\}+\gamma \right]$	$\gamma {\lambda }_n$	$\left(1-\gamma \right){\lambda }_n$
AL	$\underset{j=1,\dots s}{\max }\left\{{\lambda }_{n,j}\right\}$	$\underset{j=s+1,\dots p}{\min }\left\{{\lambda }_{n,j}\right\}$	0
Aenet	$\underset{j=1,\dots s}{\max }\left\{{\lambda }_{n,j}\left(\left(1-\gamma \right)\left|{\beta }_{0j}\right|+\gamma \right)\right\}$	$\gamma \underset{j=s+1,\dots p}{\min }\left\{{\lambda }_{n,j}\right\}$	$\left(1-\gamma \right)\underset{j=1,\dots ,s}{\max }\left\{{\lambda }_{n,j}\right\}$
SCAD	$0\text{<sup>*</sup>}$	${\lambda }_n\text{<sup>**</sup>}$	$0\text{<sup>*</sup>}$
MC+	$0\text{<sup>*</sup>}$	${\lambda }_n-\frac{K_1}{\gamma {\left|D_n\right|}^{1/2}}\text{<sup>**</sup>}$	$0\text{<sup>*</sup>}$

Table 4. Table 4: Empirical selection properties (TPR, FPR, and PPV in % percent \% ) based on 2000 replications of Thomas processes on the domain D ⊖ R symmetric-difference 𝐷 𝑅 D\ominus R ( μ = 400 𝜇 400 \mu=400 ) for different values of κ 𝜅 \kappa and for the three different scenarios. Different penalty functions are considered as well as two estimating equations, the regularized Poisson likelihood (PL) and the regularized weighted Poisson likelihood (WPL).

Method	$\kappa =5\times {10}^{-4}$						$\kappa =5\times {10}^{-5}$
	Regularized PL			Regularized WPL			Regularized PL			Regularized WPL
	TPR	FPR	PPV	TPR	FPR	PPV	TPR	FPR	PPV	TPR	FPR	PPV
	Scenario 1
Ridge	100	100	10	100	100	10	100	100	10	100	100	10
Lasso	100*	27	35	56	0*	98	89	35	34	33	0*	62
Enet	100*	59	18	39	4	36	91	60	21	31	0*	57
AL	100*	1	93	58	0*	100*	88	7	72	35	0*	67
Aenet	100*	6	72	59	0*	99	89	12	61	34	0*	64
SCAD	100*	18	41	66	0*	98	90	17	46	31	0*	56
MC+	100*	21	36	68	0*	96	90	21	42	30	0*	54
	Scenario 2
Ridge	100	100	10	100	100	10	100	100	10	100	100	10
Lasso	100*	25	35	52	1	88	90	38	29	31	0*	55
Enet	100*	52	19	49	4	62	90	60	20	24	1	38
AL	99	4	80	52	0*	100*	87	9	67	36	0*	67
Aenet	99	8	65	53	0*	99	88	14	54	35	0*	65
SCAD	100*	17	43	64	0*	92	88	17	45	28	0*	50
MC+	100*	18	41	59	1	87	88	21	41	27	0*	50
	Scenario 3
Ridge	100	100	13	100	100	13	100	100	13	100	100	13
Lasso	100*	56	24	52	2	87	98	89	15	13	2	20
Enet	100*	76	18	47	4	63	99	94	14	8	2	11
AL	100*	29	42	52	0*	100*	95	77	17	18	2	30
Aenet	100*	38	33	54	0*	99	96	82	16	15	1	25
SCAD	100*	34	33	58	0*	85	95	71	18	13	1	22
MC+	100*	35	32	56	0*	84	95	71	18	13	1	23

Table 5. Table 5: Empirical selection properties (TPR, FPR, and PPV in % percent \% ) based on 2000 replications of Thomas processes on the domain D 𝐷 D ( μ = 1600 𝜇 1600 \mu=1600 ) for different values of κ 𝜅 \kappa and for the three different scenarios. Different penalty functions are considered as well as two estimating equations, the regularized Poisson likelihood (PL) and the regularized weighted Poisson likelihood (WPL).

Method	$\kappa =5\times {10}^{-4}$						$\kappa =5\times {10}^{-5}$
	Regularized PL			Regularized WPL			Regularized PL			Regularized WPL
	TPR	FPR	PPV	TPR	FPR	PPV	TPR	FPR	PPV	TPR	FPR	PPV
	Scenario 1
Ridge	100	100	10	100	100	10	100	100	10	100	100	10
Lasso	100	26	35	52	0*	100*	98	48	22	56	0*	96
Enet	100	64	16	55	6	50	99	76	14	50	5	45
AL	100	0*	98	50	0	100	96	6	77	55	0*	98
Aenet	100	4	79	54	0*	100*	97	11	60	57	0*	96
SCAD	100	17	50	60	0*	100*	98	18	47	52	0*	90
MC+	100	22	47	60	0*	97	98	23	42	44	0*	79
	Scenario 2
Ridge	100	100	10	100	100	10	100	100	10	100	100	10
Lasso	100	26	33	51	0*	97	98	43	24	52	1	91
Enet	100	56	18	51	5	55	99	69	15	49	4	62
AL	100	1	92	51	0	100	96	10	67	53	0*	99
Aenet	100	4	78	51	0*	100*	97	15	52	53	0*	98
SCAD	100	21	37	53	1	85	96	16	50	45	1	77
MC+	100	24	35	47	2	76	97	19	47	42	2	72
	Scenario 3
Ridge	100	100	13	100	100	13	100	100	13	100	100	13
Lasso	100	69	19	52	1	96	100	95	14	48	4	75
Enet	100	85	16	52	5	71	100	97	14	43	5	62
AL	100	43	32	51	0*	100*	99	86	15	51	2	86
Aenet	100	49	27	52	0*	99	99	89	15	50	3	82
SCAD	100	47	27	43	2	72	99	78	17	40	2	63
MC+	100	48	26	44	2	75	99	79	17	37	2	61

Table 6. Table 6: Empirical prediction properties (Bias, SD, and RMSE) based on 2000 replications of Thomas processes on the domain D ⊖ R symmetric-difference 𝐷 𝑅 D\ominus R ( μ = 400 𝜇 400 \mu=400 ) for different values of κ 𝜅 \kappa and for the three different scenarios. Different penalty functions are considered as well as two estimating equations, the regularized Poisson likelihood (PL) and the regularized weighted Poisson likelihood (WPL).

Method	$\kappa =5\times {10}^{-4}$						$\kappa =5\times {10}^{-5}$
	Regularized PL			Regularized WPL			Regularized PL			Regularized WPL
	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	Scenario 1
Oracle	0.11	0.18	0.21	0.64	0.20	0.67	0.29	0.81	0.86	0.57	0.54	0.78
Ridge	0.11	0.38	0.40	0.72	0.69	1.00	0.28	1.26	1.29	0.98	1.03	1.42
Lasso	0.28	0.32	0.42	1.06	0.32	1.11	0.47	0.99	1.10	1.40	0.73	1.58
Enet	0.24	0.38	0.44	1.28	0.28	1.31	0.45	1.04	1.13	1.59	0.58	1.70
AL	0.10	0.29	0.31	0.87	0.32	0.92	0.38	0.96	1.03	1.18	0.93	1.50
Aenet	0.14	0.30	0.33	0.93	0.39	1.01	0.40	0.96	1.04	1.29	0.82	1.53
SCAD	0.26	0.27	0.38	1.06	0.37	1.12	0.46	0.79	0.91	1.49	0.67	1.64
MC+	0.28	0.28	0.39	1.04	0.38	1.11	0.47	0.78	0.92	1.48	0.70	1.64
	Scenario 2
Oracle	0.12	0.23	0.26	0.71	0.26	0.76	0.30	0.78	0.84	0.59	0.62	0.84
Ridge	0.14	0.46	0.48	0.69	0.93	1.16	0.32	1.23	1.27	0.92	1.15	1.47
Lasso	0.34	0.33	0.48	1.20	0.37	1.26	0.45	0.96	1.06	1.50	0.69	1.65
Enet	0.38	0.40	0.55	1.40	0.35	1.44	0.44	1.03	1.12	1.78	0.49	1.85
AL	0.20	0.33	0.39	0.85	0.32	0.91	0.37	0.93	1.00	1.17	0.86	1.45
Aenet	0.25	0.33	0.42	0.96	0.34	1.02	0.40	0.94	1.02	1.29	0.78	1.51
SCAD	0.38	0.30	0.48	0.95	0.48	1.06	0.44	0.80	0.91	1.53	0.70	1.68
MC+	0.39	0.30	0.49	1.01	0.49	1.13	0.44	0.80	0.92	1.52	0.71	1.68
	Scenario 3
Oracle	0.12	0.46	0.48	0.70	0.26	0.75	0.65	1.14	1.31	0.87	0.88	1.24
Ridge	0.13	1.03	1.04	0.71	1.45	1.62	0.52	3.10	3.14	0.90	2.86	3.00
Lasso	0.20	0.69	0.71	1.26	0.40	1.32	0.51	2.91	2.95	1.93	0.68	2.04
Enet	0.21	0.83	0.86	1.53	0.40	1.58	0.52	2.94	2.99	2.03	0.60	2.12
AL	0.18	0.57	0.60	0.91	0.33	0.97	0.52	2.80	2.85	1.77	0.84	1.96
Aenet	0.22	0.61	0.65	1.04	0.36	1.10	0.52	2.80	2.85	1.86	0.73	2.00
SCAD	0.27	0.61	0.67	1.18	0.59	1.32	0.48	2.49	2.54	1.91	0.64	2.02
MC+	0.27	0.62	0.68	1.20	0.58	1.33	0.48	2.49	2.54	1.89	0.67	2.00

Table 7. Table 7: Empirical prediction properties (Bias, SD, and RMSE) based on 2000 replications of Thomas processes on the domain D 𝐷 D ( μ = 1600 𝜇 1600 \mu=1600 ) for different values of κ 𝜅 \kappa and for the three different scenarios. Different penalty functions are considered as well as two estimating equations, the regularized Poisson likelihood (PL) and the regularized weighted Poisson likelihood (WPL).

Method	$\kappa =5\times {10}^{-4}$						$\kappa =5\times {10}^{-5}$
	Regularized PL			Regularized WPL			Regularized PL			Regularized WPL
	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
	Scenario 1
Oracle	0.05	0.11	0.12	0.33	0.15	0.37	0.16	0.45	0.48	0.41	0.22	0.46
Ridge	0.04	0.21	0.21	0.70	0.55	0.90	0.13	0.72	0.73	0.74	0.58	0.94
Lasso	0.14	0.19	0.24	1.03	0.20	1.05	0.23	0.60	0.64	0.99	0.43	1.08
Enet	0.11	0.22	0.24	1.14	0.29	1.17	0.20	0.62	0.65	1.12	0.43	1.20
AL	0.04	0.18	0.18	0.87	0.18	0.89	0.16	0.58	0.60	0.87	0.42	0.96
Aenet	0.05	0.18	0.18	0.96	0.22	0.99	0.17	0.58	0.60	0.90	0.48	1.02
SCAD	0.19	0.18	0.26	1.30	0.34	1.34	0.14	0.53	0.55	1.37	0.51	1.46
MC+	0.20	0.18	0.27	1.33	0.28	1.36	0.15	0.53	0.55	1.38	0.52	1.48
	Scenario 2
Oracle	0.05	0.15	0.16	0.36	0.17	0.40	0.18	0.46	0.49	0.39	0.26	0.47
Ridge	0.05	0.27	0.27	0.69	0.62	0.94	0.17	0.74	0.80	0.78	0.64	1.01
Lasso	0.16	0.20	0.25	1.16	0.24	1.18	0.23	0.60	0.64	1.14	0.43	1.22
Enet	0.17	0.23	0.29	1.24	0.24	1.26	0.23	0.63	0.67	1.33	0.42	1.40
AL	0.07	0.18	0.20	0.85	0.18	0.87	0.18	0.58	0.61	0.83	0.41	0.93
Aenet	0.09	0.19	0.21	0.94	0.20	0.96	0.20	0.59	0.62	0.92	0.41	1.01
SCAD	0.26	0.20	0.33	1.26	0.51	1.36	0.19	0.51	0.55	1.31	0.60	1.44
MC+	0.26	0.20	0.33	1.31	0.55	1.42	0.19	0.51	0.55	1.32	0.61	1.46
	Scenario 3
Oracle	0.13	0.31	0.34	0.43	0.18	0.47	0.31	0.96	1.01	0.75	0.35	0.83
Ridge	0.11	0.84	0.86	0.70	0.96	1.19	0.23	2.50	2.51	1.02	1.43	1.76
Lasso	0.12	0.64	0.65	1.14	0.29	1.17	0.22	2.41	2.42	1.40	0.61	1.52
Enet	0.13	0.71	0.73	1.35	0.30	1.39	0.23	2.42	2.43	1.63	0.56	1.73
AL	0.14	0.55	0.57	0.89	0.18	0.91	0.22	2.37	2.38	1.12	0.67	1.31
Aenet	0.15	0.56	0.58	1.00	0.22	1.03	0.22	2.36	2.37	1.26	0.64	1.41
SCAD	0.24	0.58	0.62	1.41	0.40	1.47	0.24	2.09	2.10	1.50	0.68	1.65
MC+	0.24	0.58	0.63	1.44	0.42	1.50	0.24	2.09	2.10	1.49	0.71	1.65

Table 8. Table 8: Empirical selection properties (TPR, FPR, and PPV in % percent \% ) based on 2000 replications of Thomas processes on the domain D 𝐷 D ( μ = 1600 𝜇 1600 \mu=1600 ) for κ = 5 × 10 − 5 𝜅 5 superscript 10 5 \kappa=5\times 10^{-5} , for two different scenarios, and for three different numbers of dummy points. Different estimating equations are considered, the regularized (un)weighted Poisson and (un)weighted logistic regression likelihoods, employing adaptive lasso regularization method.

Method	nd	Scenario 2						Scenario 3
		Unweighted			Weighted			Unweighted			Weighted
		TPR	FPR	PPV	TPR	FPR	PPV	TPR	FPR	PPV	TPR	FPR	PPV
Poisson	20	96	35	32	53	0*	96	98	82	16	47	2	79
	40	95	6	77	52	0*	95	98	83	16	46	2	77
	80	95	4	83	50	0*	94	98	83	16	43	2	74
Logistic	20	94	11	60	49	0*	91	98	72	20	41	2	73
	40	94	8	67	50	0*	93	99	81	16	43	2	74
	80	94	5	77	50	0*	93	99	83	16	42	2	73

Table 9. Table 9: Empirical prediction properties (Bias, SD, and RMSE) based on 2000 replications of Thomas processes on the domain D 𝐷 D ( μ = 1600 𝜇 1600 \mu=1600 ) for κ = 5 × 10 − 5 𝜅 5 superscript 10 5 \kappa=5\times 10^{-5} , for two different scenarios, and for three different numbers of dummy points. Different estimating equations are considered, the regularized (un)weighted Poisson and (un)weighted logistic regression likelihoods, employing adaptive lasso regularization method.

Method	nd	Scenario 2						Scenario 3
		Unweighted			Weighted			Unweighted			Weighted
		Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE	Bias	SD	RMSE
		No regularization
Poisson	20	0.37	0.64	0.74	0.29	0.74	0.79	0.28	2.15	2.16	0.42	2.06	2.11
	40	0.14	0.63	0.65	0.16	0.73	0.75	0.33	2.47	2.50	0.42	2.32	2.35
	80	0.17	0.64	0.66	0.11	0.75	0.76	0.26	2.57	2.58	0.43	2.40	2.43
Logistic	20	0.03	0.69	0.69	0.32	1.34	1.37	0.20	2.31	2.32	0.36	2.95	2.97
	40	0.07	0.60	0.61	0.12	0.96	0.97	0.23	2.31	2.32	0.37	2.56	2.58
	80	0.10	0.60	0.61	0.14	0.81	0.82	0.25	2.36	2.38	0.42	2.38	2.42
		Adaptive lasso
Poisson	20	0.30	0.59	0.67	0.86	0.47	0.98	0.30	2.00	2.03	1.14	0.68	1.33
	40	0.20	0.58	0.61	0.86	0.49	0.99	0.33	2.33	2.35	1.18	0.70	1.37
	80	0.18	0.59	0.62	0.88	0.51	1.02	0.28	2.41	2.43	1.22	0.71	1.41
Logistic	20	0.19	0.50	0.53	0.95	0.55	1.09	0.23	2.06	2.07	1.26	0.73	1.45
	40	0.18	0.52	0.55	0.89	0.52	1.03	0.23	2.15	2.16	1.22	0.72	1.42
	80	0.18	0.55	0.58	0.89	0.52	1.03	0.25	2.21	2.22	1.24	0.71	1.43

Table 10. Table 10: Barro Colorado Island data analysis: Parameter estimates of the regression coefficients for Beilschmiedia pendula Lauraceae trees applying regularized (un)weighted Poisson and logistic regression likelihoods with adaptive lasso regularization.

	Unweighted method		Weighted method
	Poisson estimates	Logistic estimates	Poisson estimates	Logistic estimates
Elev	0.39	0.40	0.41	0.45
Slope	0.26	0.32	0.51	0.60
Al	0	0	0	0
B	0.30	0.30	0	0
Ca	0.10	0.15	0	0
Cu	0.10	0.12	0	0
Fe	0.05	0	0	0
K	0	0	0	0
Mg	-0.17	-0.18	0	0
Mn	0.12	0.13	0.23	0.24
P	-0.60	-0.60	-0.50	-0.52
Zn	-0.43	-0.46	-0.35	-0.37
N	0	0	0	0
N.min	-0.12	-0.10	0	0
pH	-0.14	-0.14	0	0
Nb of cov.	12	11	5	5

Equations

$$\rho (u;\mathbf{\beta })=\exp ({\mathbf{\beta }}^{\top }\mathbf{z}(u)),$$

$$-\log L(\mathbf{\beta })+\lambda \Vert \mathbf{\beta }{\Vert }_1$$

$$-\log L(\mathbf{\beta })+\lambda \Vert \mathbf{\beta }{\Vert }_0.$$

$$\mu (B)=\mathbb{E}N(B),\text{mbox}B\subseteq {\mathbb{R}}^d.$$

$$\mu (B)={\int }_B\rho (u)\mathrm{d}u,\text{mbox}B\subseteq {\mathbb{R}}^d,$$

$${\alpha }^{(2)}(C)=\mathbb{E}\sum _{u,v\in \mathbf{X}}^{\ne }\mathbb{I}[(u,v)\in C],\text{mbox}C\subseteq {\mathbb{R}}^d\times {\mathbb{R}}^d.$$

$${\alpha }^{(2)}(C)=\int \int \mathbb{I}[(u,v)\in C]{\rho }^{(2)}(u,v)\mathrm{d}u\mathrm{d}v,\text{mbox}C\subseteq {\mathbb{R}}^d\times {\mathbb{R}}^d,$$

$$\mathbb{E}\sum _{u\in \mathbf{X}}k(u)=\int k(u)\rho (u)\mathrm{d}u$$

$$\mathbb{E}\sum _{u,v\in \mathbf{X}}^{\ne }k(u,v)=\int \int k(u,v){\rho }^{(2)}(u,v)\mathrm{d}u\mathrm{d}v.$$

$$g(u,v)=\frac{{\rho }^{(2)}(u,v)}{\rho (u)\rho (v)}$$

$${\rho }_c(u;\mathbf{\beta })=\exp ({\mathbf{\beta }}^{\top }\mathbf{z}(u))k(u-c;\omega )/\kappa ,$$

$$\rho (u;\mathbf{\beta })=\exp ({\mathbf{\beta }}^{\top }\mathbf{z}(u)).$$

$$k(u)=(2\pi {\omega }^2{)}^{-d/2}\exp (-\Vert u{\Vert }^2/(2{\omega }^2))$$

$$\log \mathbf{\Lambda }(u)={\mathbf{\beta }}^{\top }\mathbf{z}(u)+\mathbf{\varphi }(u)-{\sigma }^2/2,$$

$$\rho (u;\mathbf{\beta })=\exp ({\mathbf{\beta }}^{\top }\mathbf{z}(u)).$$

$$R(v-u;\zeta )=\exp (-\Vert u-v\Vert /\zeta ),\text{mbox}for\zeta >0.$$

$$\displaystyle L(\boldsymbol{\beta})={\prod_{u\in\mathbf{X}\cap D}\rho(u;\boldsymbol{\beta})}\exp\left({\int_{D}\big{(}1-\rho(u;\boldsymbol{\beta})\big{)}\mathrm{d}u}\right),$$

$$\ell (\mathbf{\beta })=\sum _{u\in \mathbf{X}\cap D}\log \rho (u;\mathbf{\beta })-{\int }_D\rho (u;\mathbf{\beta })\mathrm{d}u,$$

$$\ell (\mathbf{\beta })=\sum _{u\in \mathbf{X}\cap D}{\mathbf{\beta }}^{\top }\mathbf{z}(u)-{\int }_D\exp ({\mathbf{\beta }}^{\top }\mathbf{z}(u))\mathrm{d}u.$$

$$\mathbb{E}{\ell }^{(1)}(\mathbf{\beta })$$

$$={\int }_D\mathbf{z}(u)\exp ({\mathbf{\beta }}_0^{\top }\mathbf{z}(u))\mathrm{d}u-{\int }_D\mathbf{z}(u)\exp ({\mathbf{\beta }}^{\top }\mathbf{z}(u))\mathrm{d}u$$

$$={\int }_D\mathbf{z}(u)(\exp ({\mathbf{\beta }}_0^{\top }\mathbf{z}(u))-\exp ({\mathbf{\beta }}^{\top }\mathbf{z}(u)))\mathrm{d}u=0$$

$$\ell (w;\mathbf{\beta })=\sum _{u\in \mathbf{X}\cap D}w(u)\log \rho (u;\mathbf{\beta })-{\int }_{D}w(u)\rho (u;\mathbf{\beta })\mathrm{d}u,$$

$$w(u)=\{1+\rho (u)f(u){\}}^{-1},$$

$$\hat{K}(r)=\sum _{u,v\in \mathbf{X}\cap D}^{\ne }\frac{\mathbb{I}[\Vert u-v\Vert \le r]}{\hat{\rho }(u)\hat{\rho }(v)|D\cap D_{u-v}|}.$$

$$\ell (w;\mathbf{\beta })=$$

$$-{\int }_{D}w(u)\delta (u)\log \left(\frac{\rho (u;\mathbf{\beta })+\delta (u)}{\delta (u)}\right)\mathrm{d}u,$$

$$w(u)=\frac{\rho (u)+\delta (u)}{\delta (u)\{1+\rho (u)f(u)\}},$$

$$Q(w;\mathbf{\beta })=\ell (w;\mathbf{\beta })-|D|\sum _{j=1}^p{p}_{{\lambda }_j}(|{\beta }_j|),$$

$${\int }_D\rho (u;\mathbf{\beta })\mathrm{d}u\approx \sum _i^M{v}_i\rho (u_i;\mathbf{\beta })$$

$$\ell (\mathbf{\beta })\approx \sum _i^M{v}_i\{y_i\log \rho (u_i;\mathbf{\beta })-\rho (u_i;\mathbf{\beta })\},\text{mbox}where{y}_i=v_i^{-1}{\Delta }_i.$$

$$\ell (w;\mathbf{\beta })\approx \sum _i^M{w}_i{v}_i\{y_i\log \rho (u_i;\mathbf{\beta })-\rho (u_i;\mathbf{\beta })\}.$$

Full text

1 Introduction

Spatial point pattern data arise in many contexts where interest lies in describing the distribution of an event in space. Some examples include the locations of trees in a forest, gold deposits mapped in a geological survey, stars in a cluster star, animal sightings, locations of some specific cells in retina, or road accidents (see e.g. Møller and Waagepetersen, 2004; Illian et al., 2008; Baddeley et al., 2015). Interest in methods for analyzing spatial point pattern data is rapidly expanding accross many fields of science, notably in ecology, epidemiology, biology, geosciences, astronomy, and econometrics.

One of the main interests when analyzing spatial point pattern data is to estimate the intensity which characterizes the probability that a point (or an event) occurs in an infinitesimal ball around a given location. In practice, the intensity is often assumed to be a parametric function of some measured covariates (e.g. Waagepetersen, 2007; Guan and Loh, 2007; Møller and Waagepetersen, 2007; Waagepetersen, 2008; Waagepetersen and Guan, 2009; Guan and Shen, 2010; Coeurjolly and Møller, 2014). In this paper, we assume that the intensity function $\rho$ is parameterized by a vector $\boldsymbol{\beta}$ and has a log-linear specification

[TABLE]

where $\mathbf{z}(u)=\{z_{1}(u),\ldots,z_{p}(u)\}^{\top}$ are the $p$ spatial covariates measured at location $u$ and $\boldsymbol{\beta}=\{\beta_{1},\ldots,\beta_{p}\}^{\top}$ is a real $p$-dimensional parameter. When the intensity is a function of many variables, covariates selection becomes inevitable.

Variable selection in regression has a number of purposes: provide regularization for good estimation, obtain good prediction, and identify clearly the important variables (e.g. Fan and Lv, 2010; Mazumder et al., 2011). Identifying a set of relevant features from a list of many features is in general combinatorially hard and computationally intensive. In this context, convex relaxation techniques such as lasso (Tibshirani, 1996) have been effectively used for variable selection and parameter estimation simultaneously. The lasso procedure aims at minimizing:

[TABLE]

where $L(\boldsymbol{\beta})$ is the likelihood function for some model of interest. The $\ell_{1}$ penalty shrinks coefficients towards zero, and can also set coefficients to be exactly zero. In the context of variable selection, the lasso is often thought of as a convex surrogate for the best-subset selection problem:

[TABLE]

The $\ell_{0}$ penalty $\|\boldsymbol{\beta}\|_{0}={\sum_{i=1}^{p}\mathbb{I}(|\beta_{i}|>0)}$ penalizes the number of nonzero coefficients in the model.

Since lasso can be suboptimal in model selection for some cases (e.g. Fan and Li, 2001; Zou, 2006; Zhang and Huang, 2008), many regularization methods then have been developped, motivating to go beyond $\ell_{1}$ regime to more aggressive non-convex penalties which bridges the gap between $\ell_{1}$ and $\ell_{0}$ such as SCAD (Fan and Li, 2001) and MC+ (Zhang, 2010).

More recently, there were several works on implementing variable selection for spatial point processes in order to reduce variance inflation from overfitting and bias from underfitting. Thurman and Zhu (2014) focused on using adaptive lasso to select variables for inhomogeneous Poisson point processes. This study then later was extended to the clustered spatial point processes by Thurman et al. (2015) who established the asymptotic properties of the estimates in terms of consistency, sparsity, and normality distribution. They also compared their results employing adaptive lasso to SCAD and adaptive elastic net in the simulation study and application, using both regularized weighted and unweighted estimating equations derived from the Poisson likelihood. Yue and Loh (2015) considered modelling spatial point data with Poisson, pairwise interaction point processes, and Neyman-Scott cluster models, incorporated lasso, adaptive lasso, and elastic net regularization methods into generalized linear model framework for fitting these point models. Note that the study by Yue and Loh (2015) also used an estimating equation derived from the Poisson likelihood. However, Yue and Loh (2015) did not provide the theoretical study in detail. Although, in application, many penalty functions have been employed to regularization methods for spatial point processes intensity estimation, the theoretical study is still restricted to some specific penalty functions.

In this paper, we propose regularized versions of estimating equations based on Campbell formula derived from the Poisson and the logistic regression likelihoods to estimate the intensity of the spatial point processes. We consider both convex and non-convex penalty functions. We provide general conditions on the penalty function to ensure an oracle property and a central limit theorem. Thus, we extend the work by Thurman et al. (2015) and obtain the theoretical results for more general penalty functions and under less restrictive assumptions on the asymptotic covariance matrix (see Remark 3). The logistic regression method proposed by Baddeley et al. (2014) is as easy to implement as the Poisson likelihood method, but is less biased since it does not require deterministic numerical approximation. We prove that the estimates obtained by regularizing the logistic regression likelihood can also satisfy asymptotic properties (see Remark 2). Our procedure is straightforward to implement since we only need to combine the spatstat R package with the two R packages glmnet and ncvreg.

The remainder of the paper is organized as follows. Section 2 gives backgrounds on spatial point processes. Section 3 describes standard parameter estimation methods when there is no regularization, while regularization methods are developed in Section 4. Section 5 develops numerical details induced by the methods introduced in Sections 3-4. Asymptotic properties following the work by Fan and Li (2001) for generalized linear models are presented in Section 6. Section 7 investigates the finite-sample properties of the proposed method in a simulation study, followed by an application to tropical forestry datasets in Section 8, and finished by conclusion and discussion in Section 9. Proofs of the main results are postponed to Appendices A-C.

2 Spatial point processes

Let $\mathbf{X}$ be a spatial point process on $\mathbb{R}^{d}$. Let $D\subset\mathbb{R}^{d}$ be a compact set of Lebesgue measure $|D|$ which will play the role of the observation domain. We view $\mathbf{X}$ as a locally finite random subset of $\mathbb{R}^{d}$, i.e. the random number of points of $\mathbf{X}$ in $B$, $N(B)$, is almost surely finite whenever $B\subset\mathbb{R}^{d}$ is a bounded region. Suppose $\mathbf{x}=\{x_{1},x_{2},\ldots,x_{m}\}$ denotes a realization of $\mathbf{X}$ observed within a bounded region $D$, where $x_{i},i=1,\ldots,m$ represent the locations of the observed points, and $m$ is the number of points. Note that $m$ is random and $0\leq m<\infty$. If $m=0$ then $\mathbf{x}=\emptyset$ is the empty point pattern in D. For further background material on spatial point processes, see for example Møller and Waagepetersen (2004).

2.1 Moments

The first and second-order properties of a point process are described by intensity measure and second-order factorial moment measure. First-order properties of a point process indicate the spatial distribution events in domain of interest. The intensity measure $\mu$ on $\mathbb{R}^{d}$ is given by

[TABLE]

If the intensity measure $\mu$ can be written as

[TABLE]

where $\rho$ is a nonnegative function, then $\rho$ is called the intensity function. If $\rho$ is constant, then $\mathbf{X}$ is said to be homogeneous or first-order stationary with intensity $\rho$. Otherwise, it is said to be inhomogeneous. We may interpret $\rho(u)\mathrm{d}u$ as the probability of occurence of a point in an infinitesimally small ball with centre $u$ and volume $\mathrm{d}u$.

Second-order properties of a point process indicate the spatial coincidence of events in the domain of interest. The second-order factorial moment measure $\alpha^{(2)}$ on $\mathbb{R}^{d}\times\mathbb{R}^{d}$ is given by

[TABLE]

where the $\neq$ over the summation sign means that the sum runs over all pairwise different points $u,v$ in $\mathbf{X}$, and $\mathbb{I}[.]$ is the indicator function. If the second-order factorial moment measure $\alpha^{(2)}$ can be written as

[TABLE]

where $\rho^{(2)}$ is a nonnegative function, then $\rho^{(2)}$ is called the second-order product density. Intuitively, $\rho^{(2)}(u,v)\mathrm{d}u\mathrm{d}v$ is the probability for observing a pair of points from $\mathbf{X}$ occuring jointly in each of two infinitesimally small balls with centres $u,v$ and volume $\mathrm{d}u,\mathrm{d}v$. Fore more detail description of moment measures of any order, see appendix C in Møller and Waagepetersen (2004).

Suppose $\mathbf{X}$ has intensity function $\rho$ and second-order product density $\rho^{(2)}$. Campbell theorem (see e.g. Møller and Waagepetersen, 2004) states that, for any function $k:\mathbb{R}^{d}\to[0,\infty)$ or $k:\mathbb{R}^{d}\times\mathbb{R}^{d}\to[0,\infty)$

[TABLE]

In order to study whether a point process deviates from independence (i.e., Poisson point process), we often consider the pair correlation function given by

[TABLE]

when both $\rho$ and $\rho^{(2)}$ exist with the convention $0/0=0$. For a Poisson point process (Section 2.2.1), we have $\rho^{(2)}(u,v)=\rho(u)\rho(v)$ so that $g(u,v)=1$. 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 with the same intensity function as $\mathbf{X}$. In the same spirit, we can define $\rho^{(k)}$ the $k$-th order intensity function (see Møller and Waagepetersen, 2004, for more details). If for any $u,v$, $g(u,v)$ depends only on $u-v$, the point process $\mathbf{X}$ is said to be second-order reweighted stationary.

2.2 Modelling the intensity function

We discuss spatial point process models specified by deterministic or random intensity function. Particularly, we consider two important model classes, namely Poisson and Cox processes. Poisson point processes serve as a tractable model class for no interaction or complete spatial randomness. Cox processes form major classes for clustering or aggregation. For conciseness, we focus on the two later classes of models. We could also have presented determinantal point processes (e.g. Lavancier et al., 2015) which constitute an interesting class of repulsive point patterns with explicit moments. This has not been further investigated for sake of brevity. In this paper, we focus on log-linear models of the intensity function given by (1.1).

2.2.1 Poisson point process

A point process $\mathbf{X}$ on $D$ is a Poisson point process with intensity function $\rho$, assumed to be locally integrable, if the following conditions are satisfied:

for any $B\subseteq D$ with $0\leq\mu(B)<\infty$, $N(B)\sim Poisson(\mu(B))$, 2. 2.

conditionally on $N(B)$, the points in $\mathbf{X}\cap B$ are i.i.d. with joint density proportional to $\rho(u)$, $u\in B$.

A Poisson point process with a log-linear intensity function is also called a modulated Poisson point process (e.g. Møller and Waagepetersen, 2007; Waagepetersen, 2008). In particular, for Poisson point processes, $\rho^{(2)}(u,v)=\rho(u)\rho(v)$, and $g(u,v)=1,\forall u,v\in D$.

2.2.2 Cox processes

A Cox process is a natural extension of a Poisson point process, obtained by considering the intensity function of the Poisson point process as a realization of a random field. Suppose that $\mathbf{\Lambda}=\{\mathbf{\Lambda}(u):u\in D\}$ is a nonnegative random field. If the conditional distribution of $\mathbf{X}$ given $\mathbf{\Lambda}$ is a Poisson point process on $D$ with intensity function $\mathbf{\Lambda}$, then $\mathbf{X}$ is said to be a Cox process driven by $\mathbf{\Lambda}$ (see e.g. Møller and Waagepetersen, 2004). There are several types of Cox processes. Here, we consider two types of Cox processes: a Neyman-Scott point process and a log Gaussian Cox process.

Neyman-Scott point processes. Let $\mathbf{C}$ be a stationary Poisson process (mother process) with intensity $\kappa>0$. Given $\mathbf{C}$, let $\mathbf{X}_{c},c\in\mathbf{C}$, be independent Poisson processes (offspring processes) with intensity function

[TABLE]

where $k$ is a probability density function determining the distribution of offspring points around the mother points parameterized by $\omega$. Then $\mathbf{X}=\cup_{c\in\mathbf{C}}\mathbf{X}_{c}$ is a special case of an inhomogeneous Neyman-Scott point process with mothers $\mathbf{C}$ and offspring $\mathbf{X}_{c},c\in\mathbf{C}$. The point process $\mathbf{X}$ is a Cox process driven by $\mathbf{\Lambda}(u)=\exp(\boldsymbol{\beta}^{\top}\mathbf{z}(u)){\sum_{c\in\mathbf{C}}k(u-c,\omega})/\kappa$ (e.g. Waagepetersen, 2007; Coeurjolly and Møller, 2014) and we can verify that the intensity function of $\mathbf{X}$ is indeed

[TABLE]

One example of Neyman-Scott point process is the Thomas process where

[TABLE]

is the density for $N_{d}(0,\omega^{2}\mathbf{I}_{d})$. Conditionally on a parent event at location $c$, children events are normally distributed around $c$. Smaller values of $\omega$ correspond to tighter clusters, and smaller values of $\kappa$ correspond to fewer number of parents. The parameter vector $\boldsymbol{\psi}=(\kappa,\omega)^{\top}$ is referred to as the interaction parameter as it modulates the spatial interaction (or, dependence) among events.

Log Gaussian Cox process. Suppose that $\log\mathbf{\Lambda}$ is a Gaussian random field. Given $\mathbf{\Lambda}$, the point process $\mathbf{X}$ follows Poisson process. Then $\mathbf{X}$ is said to be a log Gaussian Cox process driven by $\mathbf{\Lambda}$ (Møller and Waagepetersen, 2004). If the random intensity function can be written as

[TABLE]

where $\boldsymbol{\phi}$ is a zero-mean stationary Gaussian random field with covariance function $c(u,v;\boldsymbol{\psi})=\sigma^{2}R(v-u;\zeta)$ which depends on parameter $\boldsymbol{\psi}=(\sigma^{2},\zeta)^{\top}$ (Møller and Waagepetersen, 2007; Coeurjolly and Møller, 2014). The intensity function of this log Gaussian Cox process is indeed given by

[TABLE]

One example of correlation function is the exponential form (e.g. Waagepetersen and Guan, 2009)

[TABLE]

Here, $\boldsymbol{\psi}=(\sigma^{2},\zeta)^{\top}$ constitutes the interaction parameter vector, where $\sigma^{2}$ is the variance and $\zeta$ is the correlation scale parameter.

3 Parametric intensity estimation

One of the standard ways to fit models to data is by maximizing the likelihood of the model for the data. While maximum likelihood method is feasible for parametric Poisson point process models (Section 3.1), computationally intensive Markov chain Monte Carlo (MCMC) methods are needed otherwise (Møller and Waagepetersen, 2004). As MCMC methods are not yet straightforward to implement, estimating equations based on Campbell theorem have been developed (see e.g. Waagepetersen, 2007; Møller and Waagepetersen, 2007; Waagepetersen, 2008; Guan and Shen, 2010; Baddeley et al., 2014). We review the estimating equations derived from the Poisson likelihood in Section 3.2-3.3 and from the logistic regression likelihood in Section 3.4.

3.1 Maximum likelihood estimation

For an inhomogeneous Poisson point process with intensity function $\rho$ parameterized by $\boldsymbol{\beta}$, the likelihood function is

[TABLE]

and the log-likelihood function of $\boldsymbol{\beta}$ is

[TABLE]

where we have omitted the constant term ${\int_{D}1\mathrm{d}u}=|D|$. As the intensity function has log-linear form (1.1), (3.1) reduces to

[TABLE]

Rathbun and Cressie (1994) showed that the maximum likelihood estimator is consistent, asymptotically normal and asymptotically efficient as the sample region goes to $\mathbb{R}^{d}$.

3.2 Poisson likelihood

Let $\boldsymbol{\beta}_{0}$ be the true parameter vector. By applying Campbell theorem (2.1) to the score function, i.e. the gradient vector of $\ell(\boldsymbol{\beta})$ denoted by $\ell^{(1)}(\boldsymbol{\beta})$, we have

[TABLE]

when $\boldsymbol{\beta}=\boldsymbol{\beta}_{0}$. So, the score function of the Poisson log-likelihood appears to be an unbiased estimating equation, even though $\mathbf{X}$ is not a Poisson point process. The estimator maximizing (3.1) is referred to as the Poisson estimator. The properties of the Poisson estimator have been carefully studied. Schoenberg (2005) showed that the Poisson estimator is still consistent for a class of spatio-temporal point process models. The asymptotic normality for a fixed observation domain was obtained by Waagepetersen (2007) while Guan and Loh (2007) established asymptotic normality under an increasing domain assumption and for suitable mixing point processes.

Regarding the parameter $\boldsymbol{\psi}$ (see Section 2.2.2), Waagepetersen and Guan (2009) studied a two-step procedure to estimate both $\boldsymbol{\beta}$ and $\boldsymbol{\psi}$, and they proved that, under certain mixing conditions, the parameter estimates $(\boldsymbol{\hat{\beta}},\boldsymbol{\hat{\psi}})$ enjoy the properties of consistency and asymptotic normality.

3.3 Weighted Poisson likelihood

Although the estimating equation approach derived from the Poisson likelihood is simpler and faster to implement than maximum likelihood estimation, it potentially produces a less efficient estimate than that of maximum likelihood (Waagepetersen, 2007; Guan and Shen, 2010) because information about interaction of events is ignored. To regain some lack of efficiency, Guan and Shen (2010) proposed a weighted Poisson log-likelihood function given by

[TABLE]

where $w(\cdot)$ is a weight surface. By regarding (3.2), we see that a larger weight $w(u)$ makes the observations in the infinitesimal region $\mathrm{d}u$ more influent. By Campbell theorem, $\ell^{(1)}(w;\boldsymbol{\beta})$ is still an unbiased estimating equation. In addition, Guan and Shen (2010) proved that, under some conditions, the parameter estimates are consistent and asymptotically normal.

Guan and Shen (2010) showed that a weight surface $w(\cdot)$ that minimizes the trace of the asymptotic variance-covariance matrix of the estimates maximizing (3.2) can result in more efficient estimates than Poisson estimator. In particular, the proposed weight surface is

[TABLE]

where $f(u)={\int_{D}\{g(\|v-u\|;\boldsymbol{\psi})-1\}\mathrm{d}u}$ and $g(\cdot)$ is the pair correlation function. For a Poisson point process, note that $f(u)=0$ and hence $w(u)=1$, which reduces to maximum likelihood estimation. For general point processes, the weight surface depends on both the intensity function and the pair correlation function, thus incorporates information on both inhomogeneity and dependence of the spatial point processes. When clustering is present so that $g(v-u)>1$, then $f(u)>0$ and hence the weight decreases with $\rho(u)$. The weight surface can be achieved by setting $\hat{w}(u)=\{1+\hat{\rho}(u)\hat{f}(u)\}^{-1}$. To get the estimate $\hat{\rho}(u)$, $\boldsymbol{\beta}$ is substituted by $\boldsymbol{\tilde{\beta}}$ given by Poisson estimates, that is, $\hat{\rho}(u)=\rho(u;\boldsymbol{\tilde{\beta}})$. Alternatively, $\hat{\rho}(u)$ can also be computed nonparametrically by kernel method. Furthermore, Guan and Shen (2010) suggessted to approximate $f(u)$ by $K(r)-\pi r^{2}$, where $K(\cdot)$ is the Ripley’s $K-$function estimated by

[TABLE]

Guan et al. (2015) extended the study by Guan and Shen (2010) and considered more complex estimating equations. Specifically, $w(u)\mathbf{z}(u)$ is replaced by a function $h(u;\boldsymbol{\beta})$ in the derivative of (3.2) with respect to $\boldsymbol{\beta}$. The procedure results in a slightly more efficient estimate than the one obtained from (3.2). However, the computational cost is more important and since we combine estimating equations and penalization methods (see Section 4.1), we have not considered this extension.

3.4 Logistic regression likelihood

Although the estimating equations discussed in Section 3.2 and 3.3 are unbiased, these methods do not, in general, produce unbiased estimator in practical implementations. Waagepetersen (2008) and Baddeley et al. (2014) proposed another estimating function which is indeed close to the score of the Poisson log-likelihood but is able to obtain less biased estimator than Poisson estimates. In addition, their proposed estimating equation is in fact the derivative of the logistic regression likelihood.

Following Baddeley et al. (2014), we define the weighted logistic regression log-likelihood function by

[TABLE]

where $\delta(u)$ is a nonnegative real-valued function. Its role as well as an explanation of the name ’logistic method’ will be explained further in Section 5.2. Note that the score of (3.4) is an unbiased estimating equation. Waagepetersen (2008) showed asymptotic normality for Poisson and certain clustered point processes for the estimator obtained from a similar procedure. Furthermore, the methodology and results were studied by Baddeley et al. (2014) considering spatial Gibbs point processes.

To determine the optimal weight surface $w(\cdot)$ for logistic method, we follow Guan and Shen (2010) who minimized the trace of the asymptotic covariance matrix of the estimates. We obtain the weight surface defined by

[TABLE]

where $\rho(u)$ and $f(u)$ can be estimated as in Section 3.3.

4 Regularization techniques

This section discusses convex and non-convex regularization methods for spatial point process intensity estimation.

4.1 Methodology

Regularization techniques were introduced as alternatives to stepwise selection for variable selection and parameter estimation. In general, a regularization method attempts to maximize the penalized log-likelihood function $\ell(\boldsymbol{\theta})-\eta{\sum_{j=1}^{p}p_{\lambda_{j}}(|\theta_{j}|)}$, where $\ell(\boldsymbol{\theta})$ is the log-likelihood function of $\boldsymbol{\theta}$, $\eta$ is the number of observations, and $p_{\lambda}(\theta)$ is a nonnegative penalty function parameterized by a real number $\lambda\geq 0.$

Let $\ell(w;\boldsymbol{\beta})$ be either the weighted Poisson log-likelihood function (3.2) or the weighted logistic regression log-likelihood function (3.4). In a similar way, we define the penalized weighted log-likelihood function given by

[TABLE]

where $|D|$ is the volume of the observation domain, which plays the same role as the number of observations $\eta$ in our setting, $\lambda_{j}$ is a nonnegative tuning parameter corresponding to $\beta_{j}$ for $j=1,\ldots,p$, and $p_{\lambda}$ is a penalty function described in details in the next section.

4.2 Penalty functions and regularization methods

For any $\lambda\geq 0$, we say that $p_{\lambda}(\cdot):\mathbb{R}^{+}\to\mathbb{R}$ is a penalty function if $p_{\lambda}$ is a nonnegative function with $p_{\lambda}(0)=0$. Examples of penalty function are the

• •

$\ell_{2}\mbox{ norm: }p_{\lambda}(\theta)=\frac{1}{2}\lambda\theta^{2}$,

• •

$\ell_{1}\mbox{ norm: }p_{\lambda}(\theta)=\lambda\theta$,

• •

$\mbox{Elastic net: }\mbox{for }0<\gamma<1,p_{\lambda}(\theta)=\lambda\{\gamma\theta+\frac{1}{2}(1-\gamma)\theta^{2}\}$,

• •

$\mbox{SCAD: }\mbox{for any }\gamma>2,p_{\lambda}(\theta)=\begin{cases}\lambda\theta&\mbox{if }\theta\leq\lambda\\ \frac{\gamma\lambda\theta-\frac{1}{2}(\theta^{2}+\lambda^{2})}{\gamma-1}&\mbox{if }\lambda\leq\theta\leq\gamma\lambda\\ \frac{\lambda^{2}(\gamma^{2}-1)}{2(\gamma-1)}&\text{if }\theta\geq\gamma\lambda,\end{cases}$

• •

$\mbox{MC+: }\mbox{for any }\gamma>1,p_{\lambda}(\theta)=\begin{cases}\lambda\theta-\frac{\theta^{2}}{2\gamma}&\mbox{if }\theta\leq\gamma\lambda\\ \frac{1}{2}\gamma\lambda^{2}&\text{if }\theta\geq\gamma\lambda.\\ \end{cases}$

The first and second derivatives of the above functions are given by Table 1. It is to be noticed that $p^{\prime}_{\lambda}$ is not differentiable at $\theta=\lambda,\gamma\lambda$ (resp. $\theta=\gamma\lambda$) for SCAD (resp. for MC+) penalty.

As a first penalization technique to improve ordinary least squares, ridge regression (e.g. Hoerl and Kennard, 1988) works by minimizing the residual sum of squares subject to a bound on the $\ell_{2}$ norm of the coefficients. As a continuous shrinkage method, ridge regression achieves its better prediction through a bias-variance trade-off. Ridge can also be extended to fit generalized linear models. However, the ridge cannot reduce model complexity since it always keeps all the predictors in the model. Then, it was introduced a method called lasso (Tibshirani, 1996), where it employs $\ell_{1}$ penalty to obtain variable selection and parameter estimation simultaneously. Despite lasso enjoys some attractive statistical properties, it has some limitations in some senses (Fan and Li, 2001; Zou and Hastie, 2005; Zou, 2006; Zhang and Huang, 2008; Zhang, 2010), making huge possibilities to develop other methods. In the scenario where there are high correlations among predictors, Zou and Hastie (2005) proposed an elastic net technique which is a convex combination between $\ell_{1}$ and $\ell_{2}$ penalties. This method is particularly useful when the number of predictors is much larger than the number of observations since it can select or eliminate the strongly correlated predictors together.

The lasso procedure suffers from nonnegligible bias and does not satisfy an oracle property asymptotically (Fan and Li, 2001). Fan and Li (2001) and Zhang (2010), among others, introduced non-convex penalties to get around these drawbacks. The idea is to bridge the gap between $\ell_{0}$ and $\ell_{1}$, by trying to keep unbiased the estimates of nonzero coefficients and by shrinking the less important variables to be exactly zero. The rationale behind the non-convex penalties such as SCAD and MC+ can also be understood by considering its first derivative (see Table 1). They start by applying the similar rate of penalization as the lasso, and then continuously relax that penalization until the rate of penalization drops to zero. However, employing non-convex penalties in regression analysis, the main challenge is often in the minimization of the possible non-convex objective function when the non-convexity of the penalty is no longer dominated by the convexity of the likelihood function. This issue has been carefully studied. Fan and Li (2001) proposed the local quadratic approximation (LQA). Zou and Li (2008) proposed a local linear approximation (LLA) which yields an objective function that can be optimized using least angle regression (LARS) algorithm (Efron et al., 2004). Finally, Breheny and Huang (2011) and Mazumder et al. (2011) investigated the application of coordinate descent algorithm to non-convex penalties.

In (4.1), it is worth emphasizing that we allow each direction to have a different regularization parameter. By doing this, the $\ell_{1}$ and elastic net penalty functions are extended to the adaptive lasso (e.g. Zou, 2006) and adaptive elastic net (e.g. Zou and Zhang, 2009). Table 2 details the regularization methods considered in this study.

5 Numerical methods

We present numerical aspects in this section. For nonregularized estimation, there are two approaches that we consider. Weighted Poisson regression is explained in Section 5.1, while logistic regression is reviewed in Section 5.2. Penalized estimation procedure is done by employing coordinate descent algorithm (Section 5.3). We separate the use of the convex and non-convex penalties in Section 5.3.1 and 5.3.2.

5.1 Weighted Poisson regression

Berman and Turner (1992) developed a numerical quadrature method to approximate maximum likelihood estimation for an inhomogeneous Poisson point process. They approximated the likelihood by a finite sum that had the same analytical form as the weighted likelihood of generalized linear model with Poisson response. This method was then extended to Gibbs point processes by Baddeley and Turner (2000). Suppose we approximate the integral term in (3.1) by Riemann sum approximation

[TABLE]

where $u_{i},i=1,\ldots,M$ are points in $D$ consisting of the $m$ data points and $M-m$ dummy points. The quadrature weights $v_{i}>0$ are such that ${\sum_{i}v_{i}}=|D|$. To implement this method, the domain is firstly partitioned into $M$ rectangular pixels of equal area, denoted by $a$. Then one dummy point is placed in the center of the pixel. Let $\Delta_{i}$ be an indicator whether the point is an event of point process ($\Delta_{i}=1$) or a dummy point ($\Delta_{i}=0$). Without loss of generality, let $u_{i},\ldots,u_{m}$ be the observed events and $u_{m+1},\ldots,u_{M}$ be the dummy points. Thus, the Poisson log-likelihood function (3.1) can be approximated and rewritten as

[TABLE]

Equation (5.1) corresponds to a quasi Poisson log-likelihood function. Maximizing (5.1) is equivalent to fitting a weighted Poisson generalized linear model, which can be performed using standard statistical software. Similarly, we can approximate the weighted Poisson log-likelihood function (3.2) using numerical quadrature method by

[TABLE]

where $w_{i}$ is the value of the weight surface at point $i$. The estimate $\hat{w}_{i}$ is obtained as suggested by Guan and Shen (2010). The similarity beetween $(\ref{eq:qua})$ and $(\ref{eq:wqua})$ allows us to compute the estimates using software for generalized linear model as well. This fact is in particular exploited in the ppm function in the $\mathtt{spatstat}$ R package (Baddeley and Turner, 2005; Baddeley et al., 2015) with option method="mpl". To make the presentation becomes more general, the number of dummy points is denoted by $\texttt{nd}^{2}$ for the next sections.

5.2 Logistic regression

To perform well, the Berman-Turner approximation often requires a quite large number of dummy points. Hence, fitting such generalized linear models can be computationally intensive, especially when dealing with a quite large number of points. When the unbiased estimating equations are approximated using deterministic numerical approximation as in Section 5.1, it does not always produce unbiased estimator. To achieve unbiased estimator, we estimate (3.4) by

[TABLE]

where $\mathcal{D}$ is dummy point process independent of $\mathbf{X}$ and with intensity function $\delta$. The form (5.3) is related to the estimating equation defined by Baddeley et al. (2014, eq. 7). Besides that, we consider this form since if we apply Campbell theorem to the last term of (5.3), we obtain

[TABLE]

which is exactly what we have in the last term of (3.4). In addition, conditional on $\mathbf{X}\cup\mathcal{D}$, (5.3) is the weighted likelihood function for Bernoulli trials, $y(u)=1\{u\in\mathbf{X}\}$ for $u\in\mathbf{X}\cup\mathcal{D}$, with

[TABLE]

Precisely, (5.3) is a weighted logistic regression with offset term $-\log\delta$. Thus, parameter estimates can be straightforwardly obtained using standard software for generalized linear models. This approach is in fact provided in the $\mathtt{spatstat}$ package in R by calling the ppm function with option method="logi" (Baddeley et al., 2014, 2015).

In spatstat, the dummy point process $\mathcal{D}$ generates $\texttt{nd}^{2}$ points in average in $D$ from a Poisson, binomial, or stratified binomial point process. Baddeley et al. (2014) suggested to choose $\delta(u)=4m/|D|$, where $m$ is the number of points (so, $\texttt{nd}^{2}=4m$). Furthermore, to determine $\delta$, this option can be considered as a starting point for a data-driven approach (see Baddeley et al., 2014, for further details).

5.3 Coordinate descent algorithm

LARS algorithm (Efron et al., 2004) is a remarkably efficient method for computing an entire path of lasso solutions. For linear models, the computational cost is of order $O(Mp^{2})$, which is the same order as a least squares fit. Coordinate descent algorithm (Friedman et al., 2007, 2010) appears to be a more competitive algorithm for computing the regularization paths by costs $O(Mp)$ operations. Therefore we adopt cyclical coordinate descent methods, which can work really fast on large datasets and can take advantage of sparsity. Coordinate descent algorithms optimize a target function with respect to a single parameter at a time, iteratively cycling through all parameters until convergence criterion is reached. We detail this for some convex and non-convex penalty functions in the next two sections. Here, we only present the coordinate descent algorithm for fitting generalized weighted Poisson regression. A similar approach is used to fit penalized weighted logistic regression.

5.3.1 Convex penalty functions

Since $\ell(w;\boldsymbol{\beta})$ given by $(\ref{eq:wqua})$ is a concave function of the parameters, the Newton-Raphson algorithm used to maximize the penalized log-likelihood function can be done using the iteratively reweighted least squares (IRLS) method. If the current estimate of the parameters is $\boldsymbol{\tilde{\beta}}$, we construct a quadratic approximation of the weighted Poisson log-likelihood function using Taylor’s expansion:

[TABLE]

where $C(\boldsymbol{\tilde{\beta}})$ is a constant, $y_{i}^{*}$ are the working response values and $\nu_{i}$ are the weights,

[TABLE]

Regularized Poisson linear model works by firstly identifying a decreasing sequence of $\lambda\in[\lambda_{\min},\lambda_{\max}]$, for which starting with minimum value of $\lambda_{\max}$ such that the entire vector $\boldsymbol{\hat{\beta}}=0$. For each value of $\lambda$, an outer loop is created to compute $\ell_{Q}(w;\boldsymbol{\beta})$ at $\boldsymbol{\tilde{\beta}}$. Secondly, a coordinate descent method is applied to solve a penalized weighted least squares problem

[TABLE]

The coordinate descent method is explained as follows. Suppose we have the estimate $\tilde{\beta}_{l}$ for $l\neq j$, $l,j=1,\ldots,p$. The method consists in partially optimizing (5.5) with respect to $\beta_{j}$, that is

[TABLE]

Friedman et al. (2007) have provided the form of the coordinate-wise update for penalized regression using several penalties such as nonnegative garrote (Breiman, 1995), lasso, elastic net, fused lasso (Tibshirani et al., 2005), group lasso (Yuan and Lin, 2006), Berhu penalty (Owen, 2007), and LAD-lasso (Wang et al., 2007a). For instance, the coordinate-wise update for the elastic net, which embraces the ridge and lasso regularization by setting respectively $\gamma$ to 0 or 1, is

[TABLE]

where ${\tilde{y}_{i}}^{(j)}={\tilde{\beta}_{0}}+{\sum_{l\neq j}z_{il}{\tilde{\beta}_{l}}}$ is the fitted value excluding the contribution from covariate $z_{ij}$, and $S(z,\lambda)$ is the soft-thresholding operator with value

[TABLE]

The update (5.6) is repeated for $j=1,\ldots,p$ until convergence. Coordinate descent algorithm for several convex penalties is implemented in the R package glmnet (Friedman et al., 2010). For (5.6), we can set $\gamma=0$ to implement ridge and $\gamma=1$ to lasso, while we set $0<\gamma<1$ to apply elastic net regularization. For adaptive lasso, we follow Zou (2006), take $\gamma=1$ and replace $\lambda$ by $\lambda_{j}=\lambda/|{\tilde{\beta}_{j}}|^{\tau}$, where $\boldsymbol{\tilde{\beta}}$ is an initial estimate, say $\boldsymbol{\tilde{\beta}}(ols)$ or $\boldsymbol{\tilde{\beta}}(ridge)$, and $\tau$ is a positive tuning parameter. To avoid the computational evaluation for choosing $\tau$, we follow Zou (2006, Section 3.4) and Wasserman and Roeder (2009) who also considered $\tau=1$, so we choose $\lambda_{j}=\lambda/|{\tilde{\beta}_{j}}(ridge)|$, where $\boldsymbol{\tilde{\beta}}(ridge)$ is the estimates obtained from ridge regression. Implementing adaptive elastic net follows along similar lines.

5.3.2 Non-convex penalty functions

Breheny and Huang (2011) have investigated the application of coordinate descent algorithm to fit penalized generalized linear model using SCAD and MC+, for which the penalty is non-convex. Mazumder et al. (2011) also studied the coordinate-wise optimization algorithm in linear models considering more general non-convex penalties.

Mazumder et al. (2011) concluded that, for a known current estimate $\tilde{\theta}$, the univariate penalized least squares function $Q_{u}(\theta)=\frac{1}{2}(\theta-\tilde{\theta})^{2}+p_{\lambda}(|\theta|)$ should be convex to ensure that the coordinate-wise procedure converges to a stationary point. Mazumder et al. (2011) found that this turns out to be the case for SCAD and MC+ penalty, but it cannot be satisfied by bridge (or power) penalty and some cases of log-penalty.

Breheny and Huang (2011) derived the solution of coordinate descent algorithm for SCAD and MC+ in generalized linear models cases, and it is implemented in the $\mathtt{ncvreg}$ package of $\mathtt{R}$. Let $\tilde{\boldsymbol{\beta}_{l}}$ be a vector containing estimates $\tilde{\beta}_{l}$ for $l\neq j$, $l,j=1,\ldots,p$, and we wish to partially optimize (5.5) with respect to $\beta_{j}$. If we define $\tilde{g}_{j}={\sum_{i=1}^{M}\nu_{j}z_{ij}(y_{i}-{\tilde{y}_{i}}^{(j)})}$ and $\tilde{\eta}_{j}={\sum_{i=1}^{M}\nu_{j}z_{ij}^{2}}$, the coordinate-wise update for SCAD is

[TABLE]

for any $\gamma>\max_{j}(1+1/\tilde{\eta}_{j})$. Then, for $\gamma>\max_{j}(1/\tilde{\eta}_{j})$ and the same definition of $\tilde{g}_{j}$ and $\tilde{\eta}_{j}$, the coordinate-wise update for MC+ is

[TABLE]

where $S(z,\lambda)$ is the soft-thresholding operator given by (5.7).

5.4 Selection of regularization or tuning parameter

It is worth noticing that coordinate descent procedures (and other computation procedures computing the penalized likelihood estimates) rely on the tuning parameter $\lambda$ so that the choice of $\lambda$ is also becoming an important task. The estimation using a large value of $\lambda$ tends to have smaller variance but larger biases, whereas the estimation using a small value of $\lambda$ leads to have zero biases but larger variance. The trade-off between the biases and the variances yields an optimal choice of $\lambda$ (Fan and Lv, 2010).

To select $\lambda$, it is reasonable to identify a range of $\lambda$ values extending from a maximum value of $\lambda$ for which all penalized coefficients are zero to $\lambda=0$ (e.g. Friedman et al., 2010; Breheny and Huang, 2011). After that, we select a $\lambda$ value which optimizes some criterion. By fixing a path of $\lambda\geq 0$, we select the tuning parameter $\lambda$ which minimizes $\mathrm{WQBIC}(\lambda)$, a weighted version of the BIC criterion, defined by

[TABLE]

where $s(\lambda)={\sum_{j=1}^{p}\mathbb{I}\{{\hat{\beta}_{j}}(\lambda)\neq 0\}}$ is the number of selected covariates with nonzero regression coefficients and $|D|$ is the observation volume which represents the sample size. For linear regression models, $\mathbf{Y}=\mathbf{X}^{\top}\boldsymbol{\hat{\beta}}+\boldsymbol{\epsilon}$, Wang et al. (2007b) proposed BIC-type criterion for choosing $\lambda$ by

[TABLE]

where $\eta$ is the number of observations and $\mathrm{DF}(\lambda)$ is the degree of freedom. This criterion is consistent, meaning that, it selects the correct model with probability approaching 1 in large samples when a set of candidate models contains the true model. Their findings is in line with the study of Zhang et al. (2010) for which the criterion was presented in more general way, called generalized information criterion (GIC). The criterion WQBIC is the specific form of GIC proposed by Zhang et al. (2010).

The selection of $\gamma$ for SCAD and MC+ is another task, but we fix $\gamma=3.7$ for SCAD and $\gamma=3$ for MC+, following Fan and Li (2001) and Breheny and Huang (2011) respectively, to avoid more complexities.

6 Asymptotic theory

In this section, we present the asymptotic results for the regularized weighted Poisson likelihood estimator when considering $\mathbf{X}$ as a $d$-dimensional point process observed over a sequence of observation domain $D=D_{n},n=1,2,\ldots$ which expands to $\mathbb{R}^{d}$ as $n\to\infty$. The regularization parameters $\lambda_{j}=\lambda_{n,j}$ for $j=1,\dots,p$ are now indexed by $n$. These asymptotic results also hold for the regularized unweighted Poisson likelihood estimator. For sake of conciseness, we do not present the asymptotic results for the regularized logistic regression estimate. The results are very similar. The main difference is lying in the conditions ($\mathcal{C}$.6) and ($\mathcal{C}$.7) for which the matrices $\mathbf{A}_{n},\mathbf{B}_{n},$ and $\mathbf{C}_{n}$ have a different expression (see Remark 2).

6.1 Notation and conditions

We recall the classical definition of strong mixing coefficients adapted to spatial point processes (e.g. Politis et al., 1998): for $k,l\in\mathbb{N}\cup\{\infty\}$ and $q\geq 1$, define

[TABLE]

where $\mathscr{F}$ is the $\sigma$-algebra generated by $\mathbf{X}\cap\Lambda_{i},i=1,2,d(\Lambda_{1},\Lambda_{2})$ is the minimal distance between sets $\Lambda_{1}$ and $\Lambda_{2}$, and $\mathscr{B}(\mathbb{R}^{d})$ denotes the class of Borel sets in $\mathbb{R}^{d}$.

Let $\boldsymbol{\beta}_{0}=\{\beta_{01},\ldots,\beta_{0s},\beta_{0(s+1)},\ldots,\beta_{0p}\}^{\top}=\{\boldsymbol{\beta}^{\top}_{01},\boldsymbol{\beta}^{\top}_{02}\}^{\top}=(\boldsymbol{\beta}_{01}^{\top},\mathbf{0}^{\top})^{\top}$ denote the $p$-dimensional vector of true coefficient values, where ${\boldsymbol{\beta}_{01}}$ is the $s$-dimensional vector of nonzero coefficients and $\boldsymbol{\beta}_{02}$ is the (p-s)-dimensional vector of zero coefficients.

We define the $p\times p$ matrices $\mathbf{A}_{n}(w;\boldsymbol{\beta}_{0}),\mathbf{B}_{n}(w;\boldsymbol{\beta}_{0}),$ and $\mathbf{C}_{n}(w;\boldsymbol{\beta}_{0})$ by

[TABLE]

Consider the following conditions ($\mathcal{C}$.1)-($\mathcal{C}$.8) which are required to derive our asymptotic results, where $o$ denotes the origin of $\mathbb{R}^{d}$:

1. ($\mathcal{C}$.1)

For every $n\geq 1,D_{n}=nE=\{ne:e\in E\}$, where $E\subset\mathbb{R}^{d}$ is convex, compact, and contains $o$ in its interior. 2. ($\mathcal{C}$.2)

We assume that the intensity function has the log-linear specification given by (1.1) where $\beta\in\Theta$ and $\Theta$ is an open convex bounded set of $\mathbb{R}^{p}$. 3. ($\mathcal{C}$.3)

The covariates $\mathbf{z}$ and the weight function $w$ satisfy

[TABLE] 4. ($\mathcal{C}$.4)

There exists an integer $t\geq 1$ such that for $k=2,\ldots,2+t$, the product density $\rho^{(k)}$ exists and satisfies $\rho^{(k)}<\infty$. 5. ($\mathcal{C}$.5)

For the strong mixing coefficients (6.1), we assume that there exists some $\tilde{t}>d(2+t)/t$ such that $\alpha_{2,\infty}(q)=O(q^{-\tilde{t}})$. 6. ($\mathcal{C}$.6)

There exists a $p\times p$ positive definite matrix $\mathbf{I}_{0}$ such that for all sufficiently large $n$, $|D_{n}|^{-1}\{\mathbf{B}_{n}(w;\boldsymbol{\beta}_{0})+\mathbf{C}_{n}(w;\boldsymbol{\beta}_{0})\}\geq\mathbf{I}_{0}$. 7. ($\mathcal{C}$.7)

There exists a $p\times p$ positive definite matrix $\mathbf{I}_{0}^{\prime}$ such that for all sufficiently large $n$, we have $|D_{n}|^{-1}\mathbf{A}_{n}(w;\boldsymbol{\beta}_{0})\geq\mathbf{I}_{0}^{\prime}$. 8. ($\mathcal{C}$.8)

The penalty function $p_{\lambda}(\cdot)$ is nonnegative on $\mathbb{R}+$, continuously differentiable on $\mathbb{R}^{+}\setminus\{0\}$ with derivative ${p}_{\lambda}^{\prime}$ assumed to be a Lipschitz function on $\mathbb{R}^{+}\setminus\{0\}$. Furthermore, given $(\lambda_{n,j})_{n\geq 1},\mbox{ for }j=1,\ldots,s,$ we assume that there exists $(\tilde{r}_{n,j})_{n\geq 1}$, where $|D_{n}|^{1/2}\tilde{r}_{n,j}\to\infty$ as $n\to\infty$, such that, for $n$ sufficiently large, $p_{\lambda_{n,j}}$ is thrice continuously differentiable in the ball centered at $|\beta_{0j}|$ with radius $\tilde{r}_{n,j}$ and we assume that the third derivative is uniformly bounded.

Under the condition ($\mathcal{C}$.8), we define the sequences $a_{n}$, $b_{n}$ and $c_{n}$ by

[TABLE]

These sequences $a_{n}$, $b_{n}$ and $c_{n}$, detailed in Table 3 for the different methods considered in this paper, play a central role in our results. Even if this will be discussed later in Section 6.3, we specify right now that we require that $a_{n}|D_{n}|^{1/2}\to 0$, $b_{n}|D_{n}|^{1/2}\to\infty$ and $c_{n}\to 0$.

6.2 Main results

We state our main results here. Proofs are relegated to Appendices A-C.

We first show in Theorem 1 that the penalized weighted Poisson likelihood estimator converges in probability and exhibits its rate of convergence.

Theorem 1**.**

Assume the conditions ($\mathcal{C}$.1)-($\mathcal{C}$.8) hold and let $a_{n}$ and $c_{n}$ be given by (6.2) and (6.4). If $a_{n}=O(|D_{n}|^{-1/2})$ and $c_{n}=o(1)$, then there exists a local maximizer ${\boldsymbol{\hat{\beta}}}$ of $Q(w;\boldsymbol{\beta})$ such that ${\bf\|\boldsymbol{\hat{\beta}}-\boldsymbol{\beta}_{0}\|}=O_{\mathrm{P}}(|D_{n}|^{-1/2}+a_{n})$.

This implies that, if $a_{n}=O(|D_{n}|^{-1/2})$ and $c_{n}=o(1)$, the penalized weighted Poisson likelihood estimator is root-$|D_{n}|$ consistent. Furthermore, we demonstrate in Theorem 2 that such a root-$|D_{n}|$ consistent estimator ensures the sparsity of $\boldsymbol{\hat{\beta}}$; that is, the estimate will correctly set $\boldsymbol{\beta}_{2}$ to zero with probability tending to 1 as $n\to\infty$, and $\boldsymbol{\hat{\beta}}_{1}$ is asymptotically normal.

Theorem 2**.**

Assume the conditions ($\mathcal{C}$.1)-($\mathcal{C}$.8) hold. If $a_{n}|D_{n}|^{1/2}\to 0$, $b_{n}|D_{n}|^{1/2}\to\infty$ and $c_{n}\to 0$ as $n\to\infty$, the root-$|D_{n}|$ consistent local maximizers ${\boldsymbol{\hat{\beta}}}=(\boldsymbol{\hat{\beta}}_{1}^{\top},\boldsymbol{\hat{\beta}}_{2}^{\top})^{\top}$ in Theorem 1 satisfy:

1. (i)

Sparsity: $\mathrm{P}(\boldsymbol{\hat{\beta}}_{2}=0)\to 1$ as $n\to\infty$, 2. (ii)

Asymptotic Normality: $|D_{n}|^{1/2}\boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0})^{-1/2}(\boldsymbol{\hat{\beta}}_{1}-\boldsymbol{\beta}_{01})\xrightarrow{d}\mathcal{N}(0,\mathbf{I}_{s})$,

where

[TABLE]

and where $\mathbf{A}_{n,11}(w;\boldsymbol{\beta}_{0})$ $(\mbox{resp. }\mathbf{B}_{n,11}(w;\boldsymbol{\beta}_{0}),\mathbf{C}_{n,11}(w;\boldsymbol{\beta}_{0}))$ is the $s\times s$ top-left corner of $\mathbf{A}_{n}(w;\boldsymbol{\beta}_{0})$ $(\mbox{resp. }\mathbf{B}_{n}(w;\boldsymbol{\beta}_{0}),\mathbf{C}_{n}(w;\boldsymbol{\beta}_{0}))$.

As a consequence, $\boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0})$ is the asymptotic covariance matrix of $\boldsymbol{\hat{\beta}}_{1}$. Note that $\boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0})^{-1/2}$ is the inverse of $\boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0})^{1/2}$, where $\boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0})^{1/2}$ is any square matrix with \boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0})^{1/2}\big{(}\boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0})^{1/2}\big{)}^{\top}=\boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0}).

Remark 1**.**

For lasso and adaptive lasso, $\boldsymbol{\Pi}_{n}=\mathbf{0}$. For other penalties, since $c_{n}=o(1)$, then $\boldsymbol{\|}\Pi_{n}\|=o(1)$. Since $\|\mathbf{A}_{n,11}(w;\boldsymbol{\beta}_{0})\|=O(|D_{n}|)$ from conditions ($\mathcal{C}$.2) and ($\mathcal{C}$.3), $|D_{n}|\,\|\boldsymbol{\Pi}_{n}\|$ is asymptotically negligible with respect to $\|\mathbf{A}_{n,11}(w;\boldsymbol{\beta}_{0})\|$.

Remark 2**.**

Theorems 1 and 2 remain true for the regularized weighted logistic regression likelihood estimates if we extend the condition ($\mathcal{C}$.3) by replacing in the expression of the matrices $\mathbf{A}_{n},\mathbf{B}_{n},$ and $\mathbf{C}_{n}$, $w(u)$ by ${w(u)\delta(u)}/({\rho(u;\boldsymbol{\beta})+\delta(u)}),u\in D_{n}$ and by adding $\sup_{u\in\mathbb{R}^{d}}\delta(u)<\infty$.

Remark 3**.**

We want to highlight here the main theoretical differences with the work by Thurman et al. (2015). First, the methodology and results are available for the logistic regression likelihood. Second, we consider very general penalty function while Thurman et al. (2015) only considered the adaptive lasso method. Third, we do not assume, as in Thurman et al. (2015), that $|D_{n}|^{-1}\mathbf{M}_{n}\to\mathbf{M}$ as $n\to\infty$ (where $\mathbf{M}_{n}$ is $\mathbf{A}_{n},\mathbf{B}_{n},$ or $\mathbf{C}_{n}$), when $\mathbf{M}$ is a positive definite matrix. Instead we assume sharper condition assuming $\lim_{n\to\infty}\nu_{\min}(|D_{n}|^{-1}\mathbf{M}_{n})>0$, where $\mathbf{M}_{n}$ is either $\mathbf{A}_{n}$ or $\mathbf{B}_{n}+\mathbf{C}_{n}$ and $\nu_{\min}(\mathbf{M^{\prime}})$ is the smallest eigenvalue of a positive definite matrix $\mathbf{M^{\prime}}$. This makes the proofs a little bit more technical.

6.3 Discussion of the conditions

We adopt the conditions ($\mathcal{C}$.1)-($\mathcal{C}$.6) based on the paper from Coeurjolly and Møller (2014). In condition ($\mathcal{C}$.1), the assumption that $E$ contains $o$ in its interior can be made without loss of generality. If instead $u$ is an interior point of $E$, then condition ($\mathcal{C}$.1) could be modified to that any ball with centre $u$ and radius $r>0$ is contained in $D_{n}=nE$ for all sufficiently large $n$. Condition ($\mathcal{C}$.3) is quite standard. From conditions ($\mathcal{C}$.2)-($\mathcal{C}$.5), the matrices $\mathbf{A}_{n}(w;\boldsymbol{\beta}_{0})$, $\mathbf{B}_{n}(w;\boldsymbol{\beta}_{0})$ and $\mathbf{C}_{n}(w;\boldsymbol{\beta}_{0})$ are bounded by $|D_{n}|$ (see e.g. Coeurjolly and Møller, 2014).

Combination of conditions ($\mathcal{C}$.1)-($\mathcal{C}$.6) are used to establish a central limit theorem for $|D_{n}|^{-1/2}\ell_{n}^{(1)}(w;\boldsymbol{\beta}_{0})$ using a general central limit theorem for triangular arrays of nonstationary random fields obtained by Karácsony (2006), which is an extension from Bolthausen (1982), then later extended to nonstationary random fields by Guyon (1995). As pointed out by Coeurjolly and Møller (2014), condition ($\mathcal{C}$.6) is a spatial average assumption like when establishing asymptotic normality of ordinary least square estimators for linear models. This condition is also useful to make sure that the matrix $|D_{n}|^{-1}\{\mathbf{B}_{n}(w;\boldsymbol{\beta}_{0})+\mathbf{C}_{n}(w;\boldsymbol{\beta}_{0})\}$ is invertible. Conditions ($\mathcal{C}$.6)-($\mathcal{C}$.7) ensure that the matrix $\boldsymbol{\Sigma}_{n}(w;\boldsymbol{\beta}_{0})$ is invertible for sufficiently large $n$. Conditions ($\mathcal{C}$.1)-($\mathcal{C}$.6) are discussed in details for several models by Coeurjolly and Møller (2014). They are satisfied for a large class of intensity functions and a large class of models including Poisson and Cox processes discussed in Section 2.2.

Condition ($\mathcal{C}$.8) controls the higher order terms in Taylor expansion of the penalty function. Roughly speaking, we ask the penalty function to be at least Lipschitz and thrice differentiable in a neighborhood of the true parameter vector. As it is, the condition looks technical, however, it is obviously satisfied for ridge, lasso, elastic net (and the adaptive versions). According to the choice of $\lambda_{n}$, it is satisfied for SCAD and MC+ when $|\beta_{0j}|$, for $j=1,\ldots,s$, is not equal to $\gamma\lambda_{n}$ and/or $\lambda_{n}$.

Theorem 2 requires the conditions $a_{n}|D_{n}|^{1/2}\to 0$, $b_{n}|D_{n}|^{1/2}\to\infty$ and $c_{n}\to 0$ as $n\to\infty$ simultaneously. By requiring these assumptions, the corresponding penalized weighted Poisson likelihood estimators possess the oracle property and perform as well as weighted Poisson likelihood estimator which estimates $\boldsymbol{\beta}_{1}$ knowing the fact that $\boldsymbol{\beta}_{2}=\mathbf{0}$.

For the ridge regularization method, $b_{n}=0$, preventing from applying Theorem 2 for this penalty. For lasso and elastic net, $a_{n}=K_{2}b_{n}$ for some constant $K_{2}>0$ ($K_{2}$=1 for lasso). The two conditions $a_{n}|D_{n}|^{1/2}\to 0$ and $b_{n}|D_{n}|^{1/2}\to\infty$ as $n\to\infty$ cannot be satisfied simultaneously. This is different for the adaptive versions where a compromise can be found by adjusting the $\lambda_{n,j}$’s, as well as the two non-convex penalties SCAD and MC+, for which $\lambda_{n}$ can be adjusted. For the regularization methods considered in this paper, the condition $c_{n}\to 0$ is implied by the condition $a_{n}|D_{n}|^{1/2}\to 0$ as $n\to\infty$.

7 Simulation study

We conduct a simulation study with three different scenarios, described in Section 7.1, to compare the estimates of the regularized Poisson likelihood (PL) and that of the regularized weighted Poisson likelihood (WPL). We also want to explore the behaviour of the estimates using different regularization methods. Empirical findings are presented in Section 7.2. Furthermore, we compare, in Section 7.3, the estimates of the regularized (un)weighted logistic likelihood and the ones of the regularized (un)weighted Poisson likelihood.

7.1 Simulation set-up

The setting is quite similar to that of Waagepetersen (2007) and Thurman et al. (2015). The spatial domain is $D=[0,1000]\times[0,500]$. We center and scale the $201\times 101$ pixel images of elevation ($x_{1}$) and gradient of elevation ($x_{2}$) contained in the bei datasets of spatstat library in R (R Core Team, 2016), and use them as two true covariates. In addition, we create three different scenarios to define extra covariates:

1. Scenario 1.

We generate eighteen $201\times 101$ pixel images of covariates as standard Gaussian white noise and denote them by $x_{3},\ldots,x_{20}$. We define $\mathbf{z}(u)=\mathbf{x}(u)=\{x_{1}(u),\ldots,x_{20}(u)\}^{\top}$ as the covariates vector. The regression coefficients for $z_{3},\ldots,z_{20}$ are set to zero. 2. Scenario 2.

First, we generate eighteen $201\times 101$ pixel images of covariates as in the scenario 1. Second, we transform them, together with $x_{1}$ and $x_{2}$, to have multicollinearity. Third, we define $\mathbf{z}(u)=\mathbf{V}^{\top}\mathbf{x}(u)$, where $\mathbf{x}(u)=\{x_{1}(u),\ldots,x_{20}(u)\}^{\top}$. More precisely, $\mathbf{V}$ is such that $\boldsymbol{\Omega}=\mathbf{V}^{\top}\mathbf{V},$ and $(\Omega)_{ij}=(\Omega)_{ji}=0.7^{|i-j|}$ for $i,j=1,\ldots,20$, except $(\Omega)_{12}=(\Omega)_{21}=0$, to preserve the correlation between $x_{1}$ and $x_{2}$. The regression coefficients for $z_{3},\ldots,z_{20}$ are set to zero. 3. Scenario 3.

We consider a more complex situation. We center and scale the 13 soil nutrients covariates obtained from the study in tropical forest of Barro Colorado Island (BCI) in central Panama (see Condit, 1998; Hubbell et al., 1999, 2005), and use 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 $\mathbf{z}(u)=\mathbf{x}(u)=\{x_{1}(u),\ldots,x_{15}(u)\}^{\top}$. The regression coefficients for $z_{3},\ldots,z_{15}$ are set to zero.

The different maps of the covariates obtained from scenarios 2 and 3 are depicted in Appendix D. Except for $z_{3}$ which has high correlation with $z_{2}$, the extra covariates obtained from scenario 2 tend to have a constant value (Figure 3). This is completely different from the ones obtained from scenario 3 (Figure 4).

The mean number of points over the domain $D$, $\mu$, is chosen to be 1600. We set the true intensity function to be $\rho(u;\boldsymbol{\beta}_{0})=\{\beta_{0}+\beta_{1}z_{1}(u)+\beta_{2}z_{2}(u)\}$, where $\beta_{1}=2$ represents a relatively large effect of elevation, $\beta_{2}=0.75$ reflects a relatively small effect of gradient, and $\beta_{0}$ is selected such that each realization has 1600 points in average. Furthermore, we erode regularly the domain $D$ such that, with the same intensity function, the mean number of points over the new domain $D\ominus R$ becomes 400. The erosion is used to observe the convergence of the procedure as the observation domain expands. We consider the default number of dummy points for the Poisson likelihood, denoted by $\texttt{nd}^{2}$, as suggested in the spatstat R package, i.e. $\texttt{nd}^{2}\approx 4m$, where $m$ is the number of points. With these scenarios, we simulate 2000 spatial point patterns from a Thomas point process using the $\mathtt{rThomas}$ function in the $\mathtt{spatstat}$ package. We also consider two different $\kappa$ parameters $(\kappa=5\times 10^{-4},\kappa=5\times 10^{-5})$ as different levels of spatial interaction and let $\omega=20$. For each of the four combinations of $\kappa$ and $\mu$, we fit the intensity to the simulated point pattern realizations. We also fit the oracle model which only uses the two true covariates.

All models are fitted using modified internal function in $\mathtt{spatstat}$ (Baddeley et al., 2015), $\mathtt{glmnet}$ (Friedman et al., 2010), and $\mathtt{ncvreg}$ (Breheny and Huang, 2011). A modification of the $\mathtt{ncvreg}$ $\mathtt{R}$ package is required to include the penalized weighted Poisson and logistic likelihood methods.

7.2 Simulation results

To better understand the behaviour of Thomas processes designed in this study, Figure 1 shows the plot of the four realizations using different $\kappa$ and $\mu$. The smaller value of $\kappa$, the tighter the clusters since there are fewer parents. When $\mu=400$, i.e. by considering the realizations observed on $D\ominus R$, the mean number of points over the 2000 replications and standard deviation are 396 and 47 (resp. 400 and 137) when $\kappa=5\times 10^{-4}$ (resp. $\kappa=5\times 10^{-5}$). When $\mu=1600$, the mean number of points and standard deviation are 1604 and 174 (resp. 1589 and 529) when $\kappa=5\times 10^{-4}$ (resp. $\kappa=5\times 10^{-5}$).

Tables 4 and 5 present the selection properties of the estimates using the penalized PL and the penalized WPL methods. Similarly to Bühlmann and Van De Geer (2011), the indices we consider are the true positive rate (TPR), the false positive rate (FPR), and the positive predictive value (PPV). TPR corresponds to the ratio of the selected true covariates over the number of true covariates, while FPR corresponds to the ratio of the selected noisy covariates over the number of noisy covariates. TPR explains how the model can correctly select both $z_{1}$ and $z_{2}$. Finally, FPR investigates how the model uncorrectly select among $z_{3}$ to $z_{p}$ ($p=20$ for scenarios 1 and 2 and $p=15$ for scenario 3). PPV corresponds to the ratio of the selected true covariates over the total number of selected covariates in the model. PPV describes how the model can approximate the oracle model in terms of selection. Therefore, we want to find the methods which have a TPR and a PPV close to 100$\%$, and a FPR close to 0.

Generally, for both the penalized PL and the penalized WPL methods, the best selection properties are obtained for a larger value of $\kappa$ which shows weaker spatial dependence. For a more clustered one, indicated by a smaller value of $\kappa$, it seems more difficult to select the true covariates. As $\mu$ increases from 400 (Table 4) to 1600 (Table 5), the TPR tends to improve, so the model can select both $z_{1}$ and $z_{2}$ more frequently.

Ridge, lasso, and elastic net are the regularization methods that cannot satisfy our theorems. It is firstly emphasized that all covariates are always selected by the ridge so that the rates are never changed whatever method used. For the penalized PL with lasso and elastic net regularization, it is shown that they tend to have quite large value of FPR, meaning that they wrongly keep the noisy covariates more frequently. When the penalized WPL is applied, we gain smaller FPR, but we suffer from smaller TPR at the same time. This smaller TPR actually comes from the unselection of $z_{2}$ which has smaller coefficient than that of $z_{1}$.

When we apply adaptive lasso, adaptive elastic net, SCAD, and MC+, we achieve better performance, especially for FPR which is closer to zero which automatically improves the PPV. Adaptive elastic net (resp. elastic net) has slightly larger FPR than adaptive lasso (resp. lasso). Among all regularization methods considered in this paper, adaptive lasso seems to outperform the other ones.

Considering scenarios 1 and 2, we observe best selection properties for the penalized PL combined with adaptive lasso. As the design is getting more complex for scenario 3, applying the penalized PL suffers from much larger FPR, indicating that this method may not be able to overcome the complicated situation. However, when we use the penalized WPL, the properties seem to be more stable for the different designs of simulation study. One more advantage when considering the penalized WPL is that we can remove almost all extra covariates. It is worth noticing that we may suffer from smaller TPR when we apply the penalized WPL, but we lose the only less informative covariates. From Tables 4 and 5, when we are faced with complex situation, we would recommend the use of the penalized WPL method with adaptive lasso penalty if the focus is on selection properties. Otherwise, the use of the penalized PL combined with adaptive lasso penalty is more preferable.

Tables 6 and 7 give the prediction properties of the estimates in terms of biases, standard deviations (SD), and square root of mean squared errors (RMSE), some criterions we define by

[TABLE]

where $\hat{\mathbb{E}}(\hat{\beta}_{j})$ and $\hat{\sigma}_{j}^{2}$ are respectively the empirical mean and variance of the estimates $\hat{\beta}_{j}$, for $j=1,\ldots,p$, where $p=20$ for scenarios 1 and 2, and $p=15$ for scenario 3.

In general, the properties improve with larger value of $\kappa$ and $\mu$ due to weaker spatial dependence and larger sample size. For the oracle model where the model contains only $z_{1}$ and $z_{2}$, the WPL estimates are more efficient than the PL estimates, particularly in the more clustered case, agreeing with the findings by Guan and Shen (2010).

When the regularization methods are applied, the bias increases in general, especially when we consider the penalized WPL method. The regularized WPL has a larger bias since this method does not select $z_{2}$ much more frequently. Furthermore, weighted method seems to introduce extra bias, even though the regularization is not considered as in the oracle model. For a low clustered process, the SD using the penalized WPL is similar to that of the penalized PL which may be because of the weaker dependence represented by larger $\kappa$, making weight surface $w(\cdot)$ closer to 1. However, a larger RMSE is obtained from the penalized WPL. When we observe the more clustered process, we obtain smaller SD using the penalized WPL which explains why in some cases (mainly scenario 3) the RMSE gets smaller.

For the ridge method, the bias is closest to that of the oracle model, but it has the largest SD. Among the regularization methods, the adaptive lasso method has the best performance in terms of prediction.

Considering scenarios 1 and 2, we obtain best properties when we apply the penalized PL with adaptive lasso penalty. As the design is getting much more complex for scenario 3, when we use the penalized PL with adaptive lasso, the SD is doubled and even quadrupled due to the overselection of many unimportant covariates. In particular, for the more clustered process, the better properties are even obtained by applying the regularized WPL combined with adaptive lasso. From Tables 6 and 7, when the focus is on prediction properties, we would recommend to apply the penalized WPL combined with adaptive lasso penalty when the observed point pattern is very clustered and when covariates have a complex stucture of covariance matrix. Otherwise, the use of the penalized PL combined with adaptive lasso penalty is more favorable. Our recommendations in terms of prediction support as what we recommend in terms of selection.

7.3 Logistic regression

Our concern here is to compare the estimates of the penalized (un)weighted logistic likelihood to that of the penalized (un)weighted Poisson likelihood with different number of dummy points. We remind that the number of dummy points comes up when we discretize the integral terms in (3.2) and in (3.4). In the following, to ease the presentation, we use the term Poisson estimates (resp. logistic estimates) for parameter estimates obtained using the regularized Poisson likelihood (resp. the regularized logistic regression likelihood).

We consider three different numbers of dummy points denoted by $\texttt{nd}^{2}$. By these different numbers of dummy points, we want to observe the properties with three different situations: (a) $\texttt{nd}^{2}<m$, (b) $\texttt{nd}^{2}\approx m$, and (c) $\texttt{nd}^{2}>m$, where $m$ is the number of points. In the following, $m\approx 1600$ and $\texttt{nd}^{2}=$ 400, 1600, and 6400. Note that the choice by default from the Poisson likelihood in spatstat corresponds to case (c). Baddeley et al. (2014) showed that for datasets with very large number of points and for very structured point processes, the logistic likelihood method is clearly preferable as it requires a smaller number of dummy points to perform quickly and efficiently. We want to investigate a similar comparison when these methods are regularized.

We only repeat the results for $\kappa=5\times 10^{-5}$ and $\mu=1600$, and for scenarios 2 and 3. We use the same selection and prediction indices examined in Section 7.2 and consider only the adaptive lasso method.

Table 8 presents selection properties for the Poisson and logistic likelihoods with adaptive lasso regularization. For unweighted versions of the procedure, the regularized logistic method outperforms the regularized Poisson method when $\texttt{nd}=20$, i.e. when the number of dummy points is much smaller than the number of points. When $\texttt{nd}^{2}\approx m$ or $\texttt{nd}^{2}>m$, the methods tend to have similar performances. When we consider weighted versions of the regularized logistic and Poisson likelihoods, the results do not change that much with nd and the regularized Poisson likelihood method slightly outperforms the regularized logistic likelihood method. In addition, for scenario 3 which considers a more complex situation, the methods tend to select the noisy covariates much more frequently.

Empirical biases, standard deviation and square root of mean squared errors are presented in Table 9. We include all empirical results for the standard Poisson and logistic estimates (i.e. no regularization is considered). Let us first consider the unweighted methods with no regularization. The logistic method clearly has smaller bias, especially when $\texttt{nd}=20$, which explains why in most situations the RMSE is smaller. However, for the weighted methods, although the logistic method has smaller bias in general, it produces much larger SD, leading to larger RMSE for all cases. When we compare the weighted and the unweighted methods for logistic estimates, in general, not only do we fail to reduce the SD, but we also have larger bias. When the adaptive lasso regularization is considered, combined with the unweighted methods, we can preserve the bias in general and simultaneously improve the SD, and hence improve the RMSE. The logistic likelihood method slightly outperforms the Poisson likelihood method. When the weighted methods are considered, we obtain smaller SD, but we have larger bias. For weighted versions of the Poisson and logistic likelihoods, the results do not change that much with nd and the weighted Poisson method slightly outperforms the weighted logistic method. From Tables 8 and 9, when the number of dummy points can be chosen as $\texttt{nd}^{2}\approx m$ or $\texttt{nd}^{2}>m$, we would recommend to apply the Poisson likelihood method. When the number of dummy points should be chosen as $\texttt{nd}^{2}<m$, the logistic likelihood method is more favorable. Our recommendations regarding whether weighted or unweighted methods follow the ones as in Section 7.2.

8 Application to forestry datasets

In a 50-hectare region ($D=1,000\mathrm{m}\times 500\mathrm{m}$) of the tropical moist forest of Barro Colorado Island (BCI) in central Panama, censuses have been carried out where 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 more than 300 species (see Condit, 1998; Hubbell et al., 1999, 2005). It is of interest to know how the very high number of different tree species continues to coexist, profiting from different habitats determined by e.g. topography or soil properties (see e.g. Waagepetersen, 2007; Waagepetersen and Guan, 2009). In particular, the selection of covariates among topological attributes and soil minerals as well as the estimation of their coefficients are becoming our most concern.

We are particularly interested in analyzing the locations of 3,604 Beilschmiedia pendula Lauraceae (BPL) tree stems. We model the intensity of BPL trees as a log-linear function of two topological attributes and 13 soil properties as the covariates. Figure 2 contains maps of the locations of BPL trees, elevation, slope, and concentration of Phosporus. BPL trees seem to appear in greater abundance in the areas of high elevation, steep slope, and low concentration of Phosporus. The covariates maps are depicted in Figure 4.

We apply the regularized (un)weighted Poisson and the logistic likelihoods, combined with adaptive lasso regularization to select and estimate parameters. Since we do not deal with datasets which have very large number of points, we can set the default number of dummy points for Poisson likelihood as in the spatstat package, i.e. the number of dummy points can be chosen to be larger than the number of points, to perform quickly and efficiently. It is worth emphasizing that we center and scale the 15 covariates to observe which one has the largest effect on the intensity. The results are presented in Table 10: 12 covariates for the Poisson likelihood and 11 for the logistic method are selected out of the 15 covariates using the unweighted methods while only 5 covariates (both for the Poisson and logistic methods) are selected using the weighted versions. The unweighted methods tend to overfit the model by overselecting unimportant covariates.

The weighted methods tend to keep out the uninformative covariates. Both Poisson and logistic estimates own similar selection and estimation results. First, we find some differences on estimation between the unweighted and the weighted methods, especially for slope and Manganese (Mn), for which the weighted methods have approximately two times larger estimators. Second, we may loose some nonzero covariates when we apply the weighted methods, even though it is only for the covariates which have relatively small coefficient. Boron (B) has high correlation with many of the other covariates, particularly with those which are not selected. This is possibly why Boron, which is selected and may have nonnegligible coefficient in the unweighted methods, is not chosen in the model. This may explain why the weighted methods introduce extra biases. However, since the situation appears to be quite close to the scenario 3 from the simulation study, the weighted methods are more favorable in terms of both selection and prediction.

In this application, we do not face any computational problem. Nevertheless, if we have to model a species of trees with much more points, the default value for nd will lead to numerical problems. In such a case, the logistic likelihood would be a good alternative.

These results suggest that BPL trees favor to live in areas of higher elevation and slope. This result is different from the findings by Waagepetersen (2007) and Guan and Loh (2007) which concluded based on standard error estimation that BPL trees do not really prefer either high or low altitudes. However, we have the same conclusion with the analysis by Guan and Shen (2010) and Thurman et al. (2015) that BPL trees prefer to live on higher altitudes. Further, higher levels of Manganese (Mn) and lower levels of both Phosporus (P) and Zinc (Zn) concentrations in soil are associated with higher appearance of BPL trees.

9 Conclusion and discussion

We develop regularized versions of estimating equations based on Campbell theorem derived from the Poisson and the logistic likelihoods. Our procedure is able to estimate intensity function of spatial point processes, when the intensity is a function of many covariates and has a log-linear form. Furthermore, our procedure is also generally easy to implement in R since we need to combine spatstat package with glmnet and ncvreg packages. We study the asymptotic properties of both regularized weighted Poisson and logistic estimates in terms of consistency, sparsity, and normality distribution. We find that, among the regularization methods considered in this paper, adaptive lasso, adaptive elastic net, SCAD, and MC+ are the methods that can satisfy our theorems.

We carry out some scenarios in the simulation study to observe selection and prediction properties of the estimates. We compare the penalized Poisson likelihood (PL) and the penalized weighted Poisson likelihood (WPL) with different penalty functions. From the results, when we deal with covariates having a complex covariance matrix and when the point pattern looks quite clustered, we recommend to apply the penalized WPL combined with adaptive lasso regularization. Otherwise, the regularized PL with adaptive lasso is more preferable. The further and more careful investigation to choose the tuning parameters may be needed to improve the selection properties. We note the bias increases quite significantly when the regularized WPL is applied. When the penalized WPL is considered, a two-step procedure may be needed to improve the prediction properties: (1) use the penalized WPL combined with adaptive lasso to chose the covariates, then (2) use the selected covariates to obtain the estimates. This post-selection inference procedure has not been investigated in this paper.

We also compare the estimates obtained from the Poisson and the logistic likelihoods. When the number of dummy points can be chosen to be either similar to or larger than the number of points, we recommend the use of the Poisson likelihood method. Nevertheless, when the number of dummy points should be chosen to be smaller than the number of points, the logistic method is more favorable.

A further work would consist in studying the situation when the number of the covariates is much larger than the sample size. In such a situation, the coordinate descent algorithm used in this paper may cause some numerical troubles. The Dantzig selector procedure introduced by Candes and Tao (2007) might be a good alternative as the implementaion for linear models (and for generalized linear models) results in a linear programming. It would be interesting to bring this approach to spatial point process setting.

Acknowledgements

We thank A. L. Thurman who kindly shared the R code used for the simulation study in Thurman et al. (2015) and P. Breheny who kindly provided his code used in ncvreg R package. We also thank R. Drouilhet for technical help. The research of J.-F. Coeurjolly is funded by ANR-11-LABX-0025 Persyval-lab (2011, project Oculo-Nimbus and Persyvact). The research of F. Letué is funded 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 Center for Tropical Forest Science website http://ctfs.si.edu/webatlas/datasets/bci/soilmaps/BCIsoil.html.

Appendix A Auxiliary Lemma

The following result is used in the proof of Theorems 1-2. Throughout the proofs, the notation $\mathbf{X}_{x}=O_{\mathrm{P}}(x_{n})$ or $\mathbf{X}_{b}=o_{\mathrm{P}}(x_{n})$ for a random vector $\mathbf{X}_{n}$ and a sequence of real numbers $x_{n}$ means that $\|\mathbf{X}_{n}\|=O_{\mathrm{P}}(x_{n})$ and $\|\mathbf{X}_{n}\|=o_{\mathrm{P}}(x_{n})$. In the same way for a vector $\mathbf{V}_{n}$ or a squared matrix $\mathbf{M}_{n}$, the notation $\mathbf{V}_{n}=O(x_{n})$ and $\mathbf{M}_{n}=O(x_{n})$ mean that $\|\mathbf{V}_{n}\|=O(x_{n})$ and $\|\mathbf{M}_{n}\|=O(x_{n})$.

Lemma 1**.**

Under the conditions ($\mathcal{C}$.1)-($\mathcal{C}$.6), the following convergence holds in distribution as $n\to\infty$

[TABLE]

Moreover as $n\to\infty$,

[TABLE]

Proof.

Let us first note that using Campbell Theorems (2.1)-(2.2)

[TABLE]

The proof of (A.1) follows Coeurjolly and Møller (2014). Let $C_{i}=i+(-1/2,1/2]^{d}$ be the unit box centered at $i\in\mathbb{Z}^{d}$ and define $\mathscr{I}_{n}=\{i\in\mathbb{Z}^{d},C_{i}\cap D_{n}\neq\emptyset\}$. Set $D_{n}={\displaystyle\bigcup_{i\in\mathscr{I}_{n}}C_{i,n}}$, where $C_{i,n}=C_{i}\cap D_{n}$. We have

[TABLE]

where

[TABLE]

For any $n\geq 1$ and any $i\in\mathscr{I}_{n}$, $Y_{i,n}$ has zero mean, and by condition ($\mathcal{C}$.4),

[TABLE]

If we combine (A.3) with conditions ($\mathcal{C}$.1)-($\mathcal{C}$.6), we can apply Karácsony (2006, Theorem 4), a central limit theorem for triangular arrays of random fields, to obtain (A.1) which also implies that

[TABLE]

as $n\to\infty$. The second result (A.2) is deduced from condition ($\mathcal{C}$.6) which in particular implies that $|D_{n}|^{1/2}\{\mathbf{B}_{n}(w;\boldsymbol{\beta}_{0})+\mathbf{C}_{n}(w;\boldsymbol{\beta}_{0})\}^{-1/2}=O(1)$.

∎

Appendix B Proof of Theorem 1

In the proof of this result and the following ones, the notation $\kappa$ stands for a generic constant which may vary from line to line. In particular this constant is independent of $n$, $\boldsymbol{\beta}_{0}$ and $\mathbf{k}$.

Proof.

Let $d_{n}=|D_{n}|^{-1/2}+a_{n}$, and $\mathbf{k}=\{k_{1},k_{2},\ldots,k_{p}\}^{\top}\in\mathbb{R}^{p}$. We remind the reader that the estimate of $\boldsymbol{\beta}_{0}$ is defined as the maximum of the function $Q$ (given by (4.1)) over $\Theta$, an open convex bounded set of $\mathbb{R}^{p}$. For any $\mathbf{k}$ such that $\|\mathbf{k}\|\leq K<\infty$, $\boldsymbol{\beta}_{0}+d_{n}\mathbf{k}\in\Theta$ for $n$ sufficiently large. Assume this is valid in the following. To prove Theorem 1, we follow the main argument by Fan and Li (2001) and aim at proving that for any given $\epsilon>0$, there exists $K>0$ such that for $n$ sufficiently large

[TABLE]

Equation (B.1) will imply that with probability at least $1-\epsilon$, there exists a local maximum in the ball $\{\boldsymbol{\beta}_{0}+d_{n}\mathbf{k}:\|\mathbf{k}\|\leq K\}$, and therefore a local maximizer $\boldsymbol{\hat{\beta}}$ such that $\|{\boldsymbol{\hat{\beta}}-\boldsymbol{\beta}_{0}}\|=O_{\mathrm{P}}(d_{n})$. We decompose $\Delta_{n}(\mathbf{k})$ as $\Delta_{n}(\mathbf{k})=T_{1}+T_{2}$ where

[TABLE]

Since $\rho(u;\cdot)$ is infinitely continuously differentiable and $\ell_{n}^{(2)}(w;\boldsymbol{\beta})=-\mathbf{A}_{n}(w;\boldsymbol{\beta})$, then using a second-order Taylor expansion there exists $t\in(0,1)$ such that

[TABLE]

Since $\Theta$ is convex and bounded and since $w(\cdot)$ and $\mathbf{z}(\cdot)$ are uniformly bounded by conditions ($\mathcal{C}$.2)-($\mathcal{C}$.3), there exists a nonnegative constant $\kappa$ such that

[TABLE]

Let $\nu_{\min}(\mathbf{M})$ be the smallest eigenvalue of a squared matrix $\mathbf{M}$. By condition ($\mathcal{C}$.7),

[TABLE]

Hence

[TABLE]

Regarding the term $T_{2}$,

[TABLE]

since for any $j$ the penalty function $p_{\lambda_{n,j}}$ is nonnegative and $p_{\lambda_{n,j}}(|\beta_{0j}|)=0$ for $j=s+1,\dots,p$.

Since $d_{n}|D_{n}|^{1/2}=O(1)$, then by ($\mathcal{C}$.8), for $n$ sufficiently large, $p_{\lambda_{n,j}}$ is twice continuously differentiable for every $\beta_{j}=\beta_{0j}+td_{n}k_{j}$ with $t\in(0,1)$. Therefore using a third-order Taylor expansion, there exist $t_{j}\in(0,1)$, $j=1,\dots,s$ such that

[TABLE]

Now by definition of $a_{n}$ and $c_{n}$ and from condition ($\mathcal{C}$.8), we deduce that there exists $\kappa$ such that

[TABLE]

from Cauchy-Schwarz inequality. Since $c_{n}=o(1)$, $d_{n}=o(1)$ and $a_{n}d_{n}|D_{n}|=O(d_{n}^{2}|D_{n}|)$, then for $n$ sufficiently large

[TABLE]

We now return to (B.1): for $n$ sufficiently large

[TABLE]

Since $d_{n}|D_{n}|=O(|D_{n}|^{1/2})$, by choosing $K$ large enough, there exists $\kappa$ such that for $n$ sufficiently large

[TABLE]

for any given $\epsilon>0$ from (A.2).

∎

Appendix C Proof of Theorem 2

To prove Theorem 2(i), we provide Lemma 2 as follows.

Lemma 2**.**

Assume the conditions ($\mathcal{C}$.1)-($\mathcal{C}$.6) and condition ($\mathcal{C}$.8) hold. If $a_{n}=O(|D_{n}|^{-1/2})$ and $b_{n}|D_{n}|^{1/2}\to\infty$ as $n\to\infty$, then with probability tending to $1$, for any $\boldsymbol{\beta}_{1}$ satisfying $\|{\boldsymbol{\beta}_{1}-\boldsymbol{\beta}_{01}}\|=O_{\mathrm{P}}(|D_{n}|^{-1/2})$, and for any constant $K_{1}>0$,

[TABLE]

Proof.

It is sufficient to show that with probability tending to $1$ as ${n\to\infty}$, for any ${\boldsymbol{\beta}_{1}}$ satisfying $\|{\boldsymbol{\beta}_{1}-\boldsymbol{\beta}_{01}}\|=O_{\mathrm{P}}(|D_{n}|^{-1/2})$, for some small $\varepsilon_{n}=K_{1}|D_{n}|^{-1/2}$, and for $j=s+1,\ldots,p$,

[TABLE]

[TABLE]

First note that by (A.2), we obtain $\|\ell^{(1)}_{n}(w;\boldsymbol{\beta}_{0})\|=O_{\mathrm{P}}(|D_{n}|^{1/2})$. Second, by conditions ($\mathcal{C}$.2)-($\mathcal{C}$.3), there exists $t\in(0,1)$ such that

[TABLE]

Third, let $0<\beta_{j}<\varepsilon_{n}$ and $b_{n}$ the sequence given by (6.3). By condition ($\mathcal{C}$.8), $b_{n}$ is well-defined and since by assumption $b_{n}|D_{n}|^{1/2}\to\infty$, in particular, $b_{n}>0$ for $n$ sufficiently large. Therefore, for $n$ sufficiently large,

[TABLE]

$\mathrm{P}\left({\partial Q(w;\boldsymbol{\beta})}/{\partial\beta_{j}}<0\right)\xrightarrow{}1\mbox{ as }n\to\infty$ since ${\partial\ell_{n}(w;\boldsymbol{\beta})}/{\partial\beta_{j}}=O_{\mathrm{P}}(|D_{n}|^{1/2})$ and $b_{n}|D_{n}|^{1/2}\xrightarrow{}\infty$. This proves (C.1). We proceed similarly to prove (C.2). ∎

Proof.

We now focus on the proof of Theorem 2. Since Theorem 2(i) is proved by Lemma 2, we only need to prove Theorem 2(ii), which is the asymptotic normality of $\boldsymbol{\hat{\beta}}_{1}$. As shown in Theorem 1, there is a root-$|D_{n}|$ consistent local maximizer $\boldsymbol{\hat{\beta}}$ of $Q(w;\boldsymbol{\beta})$, and it can be shown that there exists an estimator $\boldsymbol{\hat{\beta}}_{1}$ in Theorem 1 that is a root-$(|D_{n}|)$ consistent local maximizer of Q\Big{(}w;({\boldsymbol{\beta}_{1}}^{\top},\mathbf{0}^{\top})^{\top}\Big{)}, which is regarded as a function of $\boldsymbol{\beta}_{1}$, and that satisfies

[TABLE]

There exists $t\in(0,1)$ and $\boldsymbol{\breve{\beta}}=\boldsymbol{\hat{\beta}}+t(\boldsymbol{\beta}_{0}-\boldsymbol{\hat{\beta}})$ such that

[TABLE]

where

[TABLE]

and $\phi_{n,j}=p^{\prime}_{\lambda_{n,j}}(|\hat{\beta}_{j}|)\operatorname{sign}(\hat{\beta}_{j})-p^{\prime}_{\lambda_{n,j}}(|\beta_{0j}|)\operatorname{sign}(\beta_{0j})$. We decompose $\phi_{n,j}$ as $\phi_{n,j}=T_{1}+T_{2}$ where

[TABLE]

and where $\tilde{r}_{n,j}$ is the sequence defined in the condition ($\mathcal{C}$.8). Under this condition, the following Taylor expansion can be derived for the term $T_{1}$: there exists $t\in(0,1)$ and $\check{\beta}_{j}=\hat{\beta}_{j}+t(\beta_{0j}-\hat{\beta}_{j})$ such that

[TABLE]

where the latter equation ensues from Theorem 1 and condition ($\mathcal{C}$.8). Again, from Theorem 1, $\mathbb{I}(|\hat{\beta}_{j}-\beta_{0j}|\leq\tilde{r}_{n,j})\xrightarrow{L^{1}}1$ which implies that $\mathbb{I}(|\hat{\beta}_{j}-\beta_{0j}|\leq\tilde{r}_{n,j})\xrightarrow{\mathrm{P}}1$, so T_{1}=p^{\prime\prime}_{\lambda_{n,j}}(|{\beta}_{0j}|)(\hat{\beta}_{j}-\beta_{0j})\big{(}1+o_{\mathrm{P}}(1)\big{)}+O_{\mathrm{P}}(|D_{n}|^{-1}).

Regarding the term $T_{2}$, since $p^{\prime}_{\lambda}$ is a Lipschitz function, there exists $\kappa\geq 0$ such that

[TABLE]

By Theorem 1, $|\hat{\beta}_{j}-\beta_{0j}|=O_{\mathrm{P}}(|D_{n}|^{-1/2})$ and $\mathbb{I}(|\hat{\beta}_{j}-\beta_{0j}|>\tilde{r}_{n,j})=o_{\mathrm{P}}(1)$, so $T_{2}=o_{\mathrm{P}}(|D_{n}|^{-1/2})$ and we deduce that

[TABLE]

Let $\ell^{(1)}_{n,1}(w;\boldsymbol{\beta}_{0})$ (resp. $\ell^{(2)}_{n,1}(w;\boldsymbol{\beta}_{0})$) be the first $s$ components (resp. $s\times s$ top-left corner) of $\ell^{(1)}_{n}(w;\boldsymbol{\beta}_{0})$ (resp. $\ell^{(2)}_{n}(w;\boldsymbol{\beta}_{0})$). Let also $\boldsymbol{\Psi}_{n}$ be the $s\times s$ matrix containing $\Psi_{n,jl},j,l=1,\ldots,s$. Finally, let the vector $\mathbf{p}^{\prime}_{n}$, the vector $\boldsymbol{\phi}_{n}$ and the $s\times s$ matrix $\mathbf{M}_{n}$ be

[TABLE]

We rewrite both sides of (C.3) as

[TABLE]

By definition of $\boldsymbol{\Pi}_{n}$ given by (6.6) and from (C.4), we obtain \boldsymbol{\phi}_{n}=\boldsymbol{\Pi}_{n}(\boldsymbol{\hat{\beta}}_{1}-\boldsymbol{\beta}_{01})\big{(}1+o_{\mathrm{P}}(1)\big{)}+o_{\mathrm{P}}(|D_{n}|^{-1/2}). Using this, we deduce, by premultiplying both sides of (C.5) by $\mathbf{M}_{n}$, that

[TABLE]

The condition ($\mathcal{C}$.6) implies that there exists an $s\times s$ positive definite matrix $\mathbf{I}_{0}^{\prime\prime}$ such that for all sufficiently large $n$, we have $|D_{n}|^{-1}(\mathbf{B}_{n,11}(w;\boldsymbol{\beta}_{0})+\mathbf{C}_{n,11}(w;\boldsymbol{\beta}_{0}))\geq\mathbf{I}_{0}^{\prime\prime}$, hence $\|\mathbf{M}_{n}\|=O(|D_{n}|^{-1/2})$.

Now, $\|\boldsymbol{\Psi}_{n}\|=O_{\mathrm{P}}(|D_{n}|^{1/2})$ by conditions ($\mathcal{C}$.2)-($\mathcal{C}$.3) and by Theorem 1, and $\|\boldsymbol{\hat{\beta}}_{1}-\boldsymbol{\beta}_{01}\|=O_{\mathrm{P}}(|D_{n}|^{-1/2})$ by Theorem 1 and by Theorem 2(i). Finally, since by assumption $a_{n}=o(|D_{n}|^{-1/2})$, we deduce that

[TABLE]

Therefore, we have that

[TABLE]

From (A.1), Theorem 2(i) and by Slutsky’s Theorem, we deduce that

[TABLE]

as $n\to\infty$, which can be rewritten, in particular under ($\mathcal{C}$.7), as

[TABLE]

where $\mathbf{\Sigma}_{n}(w,\boldsymbol{\beta}_{0})$ is given by (6.5). ∎

Appendix D Map of covariates