# New Development of Bayesian Variable Selection Criteria for Spatial Point
Process with Applications

Guanyu Hu 111(to whom correspondence should be addressed) Email:
<EMAIL_ADDRESS>Department of Statistics, University of Connecticut., Fred
Huffer, and Ming-Hui Chen

###### Abstract

Selecting important spatial-dependent variables under the nonhomogeneous
spatial Poisson process model is an important topic of great current interest.
In this paper, we use the Deviance Information Criterion (DIC) and Logarithm
of the Pseudo Marginal Likelihood (LPML) for Bayesian variable selection under
the nonhomogeneous spatial Poisson process model. We further derive the new
Monte Carlo estimation formula for LPML in the spatial Poisson process
setting. Extensive simulation studies are carried out to evaluate the
empirical performance of the proposed criteria. The proposed methodology is
further applied to the analysis of two large data sets, the Earthquake Hazards
Program of United States Geological Survey (USGS) earthquake data and the
Forest of Barro Colorado Island (BCI) data.

Keywords: BCI Data, DIC, LPML, Monte Carlo Estimation, USGS Earthquake Data

##  1 Introduction

There are three different types of spatial data, including point-reference
data, areal data, and spatial point pattern data. Spatial point pattern data
are the random locations of events in space Cressie (1993); Banerjee . (2014).
Spatial point patterns assume the randomness is associated with the locations
of the points. Spatial point process models have been developed for analyzing
spatial point pattern data Moller Waagepetersen (2003); Diggle (2013). Spatial
point pattern data are routinely encountered in many different fields such as
seismology, ecology, environmental science, and epidemiology. A common goal of
these fields is to connect spatially dependent covariates to the occurrence of
events of interest in space. Spatial point process regression models are well-
suited for this goal. Spatial Poisson point processes are widely used for the
purpose of regression analysis due to the easy implementation and the
attractive theoretical properties.

For regression problem, variable selection methods are widely studied in the
statistical literature. The information criterion based variable selection
methods such as Akaike information criterion (AIC) Akaike (1974) or Bayesian
information criterion (BIC) Schwarz . (1978) are easy to implement. Other than
the information criterion based methods, the regularized regression methods
such as ridge regression Hoerl Kennard (1970), LASSO Tibshirani (1996), and
elastic net Zou Hastie (2005) are also important tools for variable selection.
Within the Bayesian framework, the deviance information criterion (DIC)
Spiegelhalter . (2002) and the logarithm of the pseudo marginal likelihood
(LPML) Gelfand Dey (1994) are two popular model selection assessment criteria.
Priors based approaches (spike and slab prior Ishwaran . (2005) and shrinkage
prior Polson Scott (2010)) are also developed for Bayesian variable selection.

In the existing literature, the variable selection problem has been
extensively studied for point-reference data or areal data (Wang Zhu (2009);
Reich . (2010); Beale . (2010)). However, for the spatial point process model,
the literature on variable selection is relatively sparse. Thurman Zhu (2014)
proposed a regularized method that allows the selection of covariates at
multiple pixel resolutions. Yue Loh (2015) incorporated least absolute
shrinkage and selection operator (LASSO), adaptive LASSO and elastic net
regularization methods into the generalized linear model framework to select
important covariates for spatial point process models. Thurman . (2015)
proposed a regularized estimating equation approach for the spatial point
process model. Leininger . (2017) proposed Bayesian model assessment methods
using the predictive mean square error, empirical coverage, and ranked
probability scores. But these methods do not account for the tradeoff between
the goodness of fit and the model complexity.

The major contribution of this paper is to develop a novel implementation of
DIC and LPML under spatial point process models. Compared with point-reference
data, the traditional Monte Carlo method for computing Conditional Predictive
Ordinate CPO Chen . (2012) is not applicable for spatial point pattern data.
Therefore, we develop a new Monte Carlo method of CPO for spatial point
pattern data. Compared with the regularized methods, it is easy to carry out
Bayesian inference using posterior samples based on our model selection
approach. Our simulation studies show the promising empirical performance of
the proposed Bayesian estimate criteria. In addition, our proposed Bayesian
approach reveals some interesting features of the earthquake data sets.
Furthermore, our proposed criteria overcome the limitations of regularized
methods in selecting the important covariates and pixel resolution of the
covariates simultaneously for BCI data.

The rest of this paper is organized as follows. In Section 2, we review
nonhomogeneous spatial Poisson process and provide the Bayesian formulation of
the spatial Poisson process model. The detailed development of DIC and LPML
for the spatial Poisson process model is given in Section 3. Simulation
studies are presented in Section 4. In Section 5, we carry out an in-depth
analysis of the two data sets, the Earthquake Hazards Program of United States
Geological Survey (USGS) earthquake data and Forest of Barro Colorado Island
(BCI) data. We conclude the paper with a brief discussion in Section 6.

##  2 Bayesian Spatial Point Process

A spatial point pattern is a data set $\bm{y}=(s_{1},s_{2},...,s_{\ell})$
consists the locations $(s_{1},s_{2},...,s_{\ell})$ of points that are
observed in a bounded region $\mathcal{B}\subseteq\mathcal{R}^{2}$, which is a
realization of spatial point process $\bm{Y}$.
$N_{\bm{Y}}(A)=\sum_{i=1}^{\ell}1(s_{i}\in A)$ is a counting process
associated with the spatial point process $\bm{Y}$, which counts the number of
points of $\bm{Y}$ for area $A\subseteq\mathcal{B}$. For the process $\bm{Y}$,
there are many parametric distributions for a finite set of count variables
like Poisson processes, Gibbs processes, and Cox processes. In this paper, our
focus is on Poisson processes. For the Poisson process $\bm{Y}$ over
$\mathcal{B}$, which has the intensity function $\lambda(\bm{s})$,
$N_{\bm{Y}}(A)\sim Po(\lambda(A))$, where $\lambda(A)=\int_{A}\lambda(s)ds$.
In addition, if $A_{1}$ and $A_{2}$ are disjoint, then $N_{\bm{Y}}(A_{1})$ and
$N_{\bm{Y}}(A_{2})$ are independent, where $A_{1}\subseteq\mathcal{B}$ and
$A_{2}\subseteq\mathcal{B}$. Based on the properties of the Poisson process,
it is easy to obtain $E(N_{\bm{Y}}(A))=\text{Var}(N_{\bm{Y}}(A))=\lambda(A)$.
When $\lambda(s)=\lambda$, we have the constant intensity over the space
$\mathcal{B}$ and in this special case, $\bm{Y}$ reduces to a homogeneous
Poisson process (HPP). For a more general case, $\lambda(s)$ can be spatially
varying, which leads to a nonhomogeneous Poisson process (NHPP). For the NHPP,
the log-likelihood on $\mathcal{B}$ is given by

$\displaystyle\ell=\sum_{i=1}^{k}\text{log}\lambda(s_{i})-\int_{\mathcal{B}}\lambda(s)ds,$
(1)

where $\lambda(s_{i})$ is the intensity function for location $s_{i}$. To
build up a spatially varying intensity, $\lambda(s_{i})$ is often assumed to
take the following regression form:

$\displaystyle\lambda(s_{i})=\lambda_{0}\exp(\bm{\beta}\bm{Z}(s_{i})),$ (2)

where $\lambda_{0}$ is the baseline intensity function, $\bm{\beta}$ is the
vector of regression coefficient, and $\bm{Z}(s_{i})$ is the vector of
spatially varying covariates on location $s_{i}$.

For the model defined in (1) and (2), a gamma distribution for $\lambda_{0}$
is a conjugate prior. For the regression coefficients $\bm{\beta}$, there are
no such conjugate priors for the log-likelihood in (1). We specify a non-
informative normal prior for $\bm{\beta}$. Plugging (2) in (1) and using the
above priors for $\lambda_{0}$ and $\bm{\beta}$, we have

$\displaystyle\begin{split}\ell&=\sum_{i=1}^{k}(\text{log}\lambda_{0}+\bm{\beta}\bm{Z}(s_{i}))-\int_{\mathcal{B}}\lambda_{0}\exp(\bm{\beta}\bm{Z}(s))ds,\\\
\lambda_{0}&\sim\text{G}(a_{1},b_{1}),\\\
\bm{\beta}&\sim\text{MVN}(0,\sigma_{0}^{2}I),\end{split}$ (3)

where $l$ is the log-likelihood, and “G”,“MVN”, and “IG” are the shorthand of
gamma distribution, multivariate normal distribution, and the inverse gamma
distribution, respectively. Thus, the posterior distribution of
$\bm{\theta}=(\lambda_{0},\beta_{1},\cdots,\beta_{p})$ is given as follow

$\displaystyle\prod_{i=1}^{n}\lambda(s_{i})\times\exp\left(-\int_{\mathcal{B}}\lambda(s)ds\right)\times\pi(\bm{\theta}),$
(4)

where $\pi(\bm{\theta})=\pi(\bm{\beta})\pi(\lambda_{0})$ is the joint prior.
The analytical evaluation of the posterior distribution of $\bm{\theta}$ is
not possible. However, a Metropolis-Hasting algorithm within the Gibbs sampler
can be developed to sample from the posterior distribution in (4). The
algorithm requires sampling the following parameters in turn from their
respective full conditional distributions:

$\displaystyle\begin{split}f(\bm{\beta}|\cdot)&\propto
L\times\prod_{k=1}^{p}\frac{1}{\sqrt{2\pi\sigma^{2}_{0}}}\exp(-\frac{\beta^{2}_{k}}{2\sigma^{2}_{0}}),\\\
f(\lambda_{0}|\cdot)&\propto
L\times\frac{1}{\Gamma(a_{1})b_{1}^{a_{1}}}{\lambda_{0}}^{a_{1}-1}\exp(-\frac{\lambda_{0}}{b_{1}}).\end{split}$
(5)

We choose $\sigma_{0}^{2}=100$ and $a_{1}=b_{1}=0.01$, while yields a non-
informative joint prior for $\bm{\theta}$. To sample from the posterior
distribution of $\bm{\theta}$ in (4), an Metropolis–Hasting within Gibbs
algorithm is facilitated by R package nimble de Valpine . (2017). The
loglikelihood function of the spatial Poisson point process model used in the
MCMC iteration is directly defined using the RW_llFunction() sampler.

##  3 Bayesian Criteria for Variable Selection

We first introduce two Bayesian model selection criteria, Deviance Information
Criterion (DIC) and logarithm of the Pseudo-marginal likelihood (LPML). The
deviance information criterion is defined as

$\displaystyle\text{DIC}=\text{Dev}(\bar{\bm{\theta}})+2p_{D},$ (6)

where $\text{Dev}(\bar{\bm{\theta}})$ is the deviance function,
$p_{D}=\overline{\text{Dev}}({\bm{\theta}})-\text{Dev}(\bar{\bm{\theta}})$ is
the effective number of model parameters, $\bar{\bm{\theta}}$ is the posterior
mean of parameters $\bm{\theta}$, and $\overline{\text{Dev}}(\bm{\theta})$ is
the posterior mean of $\text{Dev}(\bm{\theta})$.

The LPML is defined as

$\displaystyle\text{LPML}=\sum_{i=1}^{n}\text{log}(\text{CPO}_{i}),$ (7)

where $\text{CPO}_{i}$ is the conditional predictive ordinate (CPO) for the
$i$-th subject. CPO is based on the leave-one-out-cross-validation. CPO
estimates the probability of observing $y_{i}$ in the future after having
already observed $y_{1},\cdots,y_{i-1},y_{i+1},\cdots,y_{n}$. The CPO for the
$i$-th subject is defined as

$\displaystyle\text{CPO}_{i}=f(y_{i}|\bm{y}_{-i})\equiv\int
f(y_{i}|\bm{\theta})\pi(\bm{\theta}|\bm{y}_{(-i)})d\bm{\theta},$ (8)

where $\bm{y}_{-i}$ is $y_{1},\cdots,y_{i-1},y_{i+1}\cdots,y_{n}$,

$\displaystyle\pi(\bm{\theta}|\bm{y}_{-i})=\frac{\prod_{j\neq
i}f(y_{j}|\bm{\theta})\pi(\bm{\theta})}{c(\bm{y}_{-i})},$ (9)

and $c(\bm{y}_{-i})$ is the normalizing constant. The $\text{CPO}_{i}$ in (8)
can be expressed as

$\displaystyle\text{CPO}^{-1}_{i}=\frac{1}{\int\frac{1}{f(y_{i}|\bm{\theta})}\pi(\bm{\theta}|\bm{y})d\bm{\theta}}.$
(10)

Using (10), a Monte Carlo estimate of $\text{CPO}_{i}$ in (10) is given by

$\displaystyle\widehat{\text{CPO}}^{-1}_{i}=\frac{1}{B}\sum_{b=1}^{B}\frac{1}{f(y_{i}|\bm{\theta}_{b})},$
(11)

where $\bm{\theta}_{b}$ is the $b$-th MCMC sample of $\bm{\theta}$ from
$\pi(\bm{\theta}|\bm{y})$.

In the context of variable selection, we select a variable subset model which
has the smallest DIC and the largest LPML.

###  3.1 DIC for Spatial Point Process

Since our main objective is to assess the fit of the spatial point pattern, we
specify the following deviance function

$\displaystyle\text{Dev}\left(D,\bm{Z},\bm{\beta},\lambda_{0}\right)=-2\times\left(\sum_{i=1}^{k}\text{log}\lambda(s_{i})-\int_{\mathcal{B}}\lambda(s)ds\right),$
(12)

where $D=(\bm{s}_{1},\bm{s}_{2},\cdots,\bm{s}_{k})$ is the observed points,
$\lambda(s)=\lambda_{0}\exp(\beta_{1}Z_{1}(s)+\cdots+\beta_{p}Z_{p}(s))$ is
the intensity function on location $\bm{s}$, and
$Z_{1}(s),Z_{2}(s),\cdots,Z_{p}(s)$ are spatial covariates. Therefore, the DIC
for the spatial point pattern is given as follows

$\displaystyle\text{DIC}=2E[\text{Dev}(D,\bm{Z},\bm{\beta},\lambda_{0})]-\text{Dev}(D,\bm{Z},\hat{\bm{\beta}},\hat{\lambda}_{0}),$
(13)

where $\hat{\bm{\beta}}$ and $\hat{\lambda}_{0}$ are the posterior means of
$\bm{\beta}$ and $\lambda_{0}$ and the expectation is taken with respect to
the posterior distribution of $\bm{\beta}$ and $\lambda_{0}$.

###  3.2 LPML for Spatial Point Processes

In this section we propose a definition of LPML which is suitable for general
Poisson processes. This definition is derived from the definition of LPML in
(7) by a limiting argument. We begin by slightly altering the definition of
$\text{CPO}_{i}$ given in (8). Let $f_{0}$ denote a fixed reference model for
the data $\bm{y}=(y_{1},y_{2},\ldots,y_{n})$, and define

$\text{CPO}_{i}=\frac{f(y_{i}|\bm{y}_{-i})}{f_{0}(y_{i}|\bm{y}_{-i})}\,.$ (14)

This definition is essentially equivalent to that in (8) since the reference
model cancels out when comparisons are made between models using
$\text{LPML}=\sum_{i=1}^{n}\log(\text{CPO}_{i})$. We shall assume that the
reference model $f_{0}$ involves no parameters (or that any parameters are
kept fixed at specified default values) and that the observations
$y_{1},y_{2},\ldots,y_{n}$ are independent. In this case
$f_{0}(y_{i}|\bm{y}_{-i})=f_{0}(y_{i})$. Repeating the argument which leads to
(10) then yields

$\text{CPO}_{i}^{-1}=\int\frac{f_{0}(y_{i})}{f(y_{i}|\bm{\theta})}\pi(\bm{\theta}|\bm{y})\,d\bm{\theta}.$
(15)

Suppose we have a full Bayesian model (such as that described in (2) and (3))
for a Poisson process $\bm{Y}$ with intensity $\lambda_{\bm{\theta}}(\cdot)$
on a region $\mathcal{B}\subset\mathcal{R}^{2}$. We observe the realization
$Y=\\{s_{1},s_{2},\ldots,s_{k}\\}$.

Let $A_{1},A_{2},\cdots,A_{n}$ be a partition of $\mathcal{B}$, i.e., disjoint
subsets such that $\bigcup_{i=1}^{n}A_{i}=\mathcal{B}$. Suppose we were to
make inference on $\bm{\theta}$ based only on the counts
$(N_{\bm{Y}}(A_{1}),N_{\bm{Y}}(A_{2}),\ldots,N_{\bm{Y}}(A_{n}))$ where
$N_{\bm{Y}}(A_{i})=\sum_{j=1}^{k}1(s_{j}\in A_{i})$. Conditional on
$\bm{\theta}$, these counts are independent Poisson random variables which we
take to be the values $y_{1},y_{2},\ldots,y_{n}$ making up the data $\bm{y}$
in the generic notation of equations (14) and (15). For our reference model we
use the Poisson process with a constant intensity of one. With these choices,
we make the following substitutions in (15):

$\displaystyle f(y_{i}|\bm{\theta})$
$\displaystyle=\frac{\lambda(A_{i})^{N_{\bm{Y}}(A_{i})}}{N_{\bm{Y}}(A_{i})!}\exp\\{-\lambda(A_{i})\\},$
$\displaystyle f_{0}(y_{i})$
$\displaystyle=\frac{|A_{i}|^{N_{\bm{Y}}(A_{i})}}{N_{\bm{Y}}(A_{i})!}\exp(-|A_{i}|),$

where $\lambda(A_{i})=\int_{A_{i}}\lambda(u)\,du$ and
$|A_{i}|=\int_{A_{i}}du$. (Here, to lighten the notation, we omit the
dependence of $\lambda$ on $\bm{\theta}$.) Therefore

$\displaystyle\frac{f_{0}(y_{i})}{f(y_{i}|\bm{\theta})}$
$\displaystyle=\left(\frac{\lambda(A_{i})}{|A_{i}|}\right)^{-N_{\bm{Y}}(A_{i})}\exp(\lambda(A_{i})-|A_{i}|),$

and we may write (15) as

$\text{CPO}_{i}^{-1}=E_{n}\left[\left(\frac{\lambda(A_{i})}{|A_{i}|}\right)^{-N_{\bm{Y}}(A_{i})}\exp\\{\lambda(A_{i})-|A_{i}|\\}\right],$

where $E_{n}$ represents the expectation over the posterior distribution of
$\bm{\theta}$ given $(N_{\bm{Y}}(A_{i}),\ldots,N_{\bm{Y}}(A_{n}))$.

For the fixed realization of the Poisson process
$\bm{Y}=\\{s_{1},s_{2},\ldots,s_{k}\\}$, we consider (informally) taking a
limit as the partition $A_{1},A_{2},\ldots,A_{n}$ becomes progressively finer.
When the partition is sufficiently fine, no set $A_{i}$ will contain more than
one of the points $s_{j}$; exactly $k$ of the sets $A_{i}$ will contain one
point, and the others will contain zero points. When the set $A_{i}$ is
sufficiently small, we expect $e^{\lambda(A_{i})-|A_{i}|}\approx
1+(\lambda(A_{i})-|A_{i}|)$. Therefore, when all the sets $|A_{i}|$ are
sufficiently small, for those sets $A_{i}$ which contain zero points of
$\bm{Y}$, we expect that

$\displaystyle\text{CPO}_{i}^{-1}$
$\displaystyle=E_{n}\left[\exp\\{\lambda(A_{i})-|A_{i}|\\}\right]\approx
1+E_{n}(\lambda(A_{i})-|A_{i}|),$

and based on the first order of Taylor expansion of $\log(1+x)$, we have

$\log\text{CPO}_{i}\approx-E_{n}(\lambda(A_{i})-|A_{i}|).$ (16)

If a set $A_{i}$ is sufficiently small and contains a single point $s_{j}$ of
$\bm{Y}$, then (assuming the intensity $\lambda(u)$ is a continuous function
of $u$ for all $\bm{\theta})$ and $A_{i}$ is smaller enough such that have
negligible area) we expect $\lambda(A_{i})/|A_{i}|\approx\lambda(s_{j})$ and
therefore

$\displaystyle\text{CPO}_{i}^{-1}$
$\displaystyle=E_{n}\left[\left(\frac{\lambda(A_{i})}{|A_{i}|}\right)^{-1}e^{\lambda(A_{i})-|A_{i}|}\right]\approx
E_{n}\lambda(s_{j})^{-1}(1+(\lambda(A_{i})-|A_{i}|))\approx
E_{n}\left[\lambda(s_{j})^{-1}\right],$

and

$\displaystyle\log\text{CPO}_{i}$
$\displaystyle\approx\log\left(E_{n}\left[\lambda(s_{j})^{-1}\right]\right)^{-1}.$
(17)

Combining (16) and (17), we obtain

$\displaystyle\text{LPML}_{n}$
$\displaystyle=\sum_{i=1}^{n}\log\text{CPO}_{i}=\sum_{i:N_{\bm{Y}}(A_{i})=1}\log\text{CPO}_{i}+\sum_{i:N_{\bm{Y}}(A_{i})=0}\log\text{CPO}_{i}$
$\displaystyle\approx\sum_{j=1}^{k}\log\left(E_{n}\left[\lambda(s_{j})^{-1}\right]\right)^{-1}-\sum_{i:N_{\bm{Y}}(A_{i})=0}E_{n}(\lambda(A_{i})-|A_{i}|)$
$\displaystyle\approx\sum_{j=1}^{k}\log\left(E_{n}\left[\lambda(s_{j})^{-1}\right]\right)^{-1}-\int_{\mathcal{B}}E_{n}[\lambda(u)]\,du+|\mathcal{B}|.$

As the partition gets progressively finer, the posterior conditional on the
counts $(N_{\bm{Y}}(A_{1}),$ $N_{\bm{Y}}(A_{2}),\ldots,N_{\bm{Y}}(A_{n}))$
converges to the posterior conditional on the complete data $\bm{Y}$
Waagepetersen (2004), and the various approximations we had made above become
exact in the limit. We obtain

$\lim_{n\rightarrow\infty}\text{LPML}_{n}=\sum_{j=1}^{k}\log\left(E\left[\lambda(s_{j})^{-1}\right]\right)^{-1}-\int_{\mathcal{B}}E[\lambda(u)]\,du+|\mathcal{B}|,$

where $E$ denotes the posterior expectation. Dropping the constant term
$|\mathcal{B}|$, we take the remainder of this expression as our definition of
the LPML for Poisson processes:

$\text{LPML}=\sum_{j=1}^{k}\log\left(E\left[\lambda(s_{j})^{-1}\right]\right)^{-1}-\int_{\mathcal{B}}E[\lambda(u)]\,du,$
(18)

where $\bm{Y}=\\{s_{1},s_{2},\ldots,s_{k}\\}$ and $E$ denotes the posterior
expectation. There is a natural Monte Carlo estimate of the LPML given by

$\widehat{\text{LPML}}=\sum_{j=1}^{k}\log\widetilde{\lambda}(s_{j})-\int_{\mathcal{B}}\overline{\lambda}(u)\,du,$
(19)

where
$\widetilde{\lambda}(s_{j})=(\frac{1}{B}\sum_{b=1}^{B}\lambda(s_{j}|\bm{\theta_{b}})^{-1})^{-1}$,
$\overline{\lambda}(u)=\frac{1}{B}\sum_{b=1}^{B}\lambda(u|\bm{\theta_{b}})$,
and $\\{\theta_{1},\theta_{2},\cdots,\theta_{B}\\}$ is a sample from the
posterior.

##  4 Simulation Study

In this section, we have four scenarios to generate covariates. The intensity
function of data generation model (DGM) in scenario 1 is
$\lambda(s)=\lambda_{0}\exp(\beta_{1}\times Z_{1}(s)+\beta_{2}\times
Z_{2}(s)+\beta_{3}\times Z_{1}(s)Z_{2}(s))$, where $Z_{1}(s)$ and $Z_{2}(s)$
are $x$-coordinate and $y$-coordinate of location $s$. In this simulation we
choose $\beta_{1}=2$, $\beta_{2}=0$, $\beta_{3}=1$, and $\lambda_{0}=30$ and
generate data on $[0,1]\times[0,1]$. We generate 100 data sets under this
setting, and the average number of points in this scenario is around 130. In
this simulation study, we compare seven different models. We give
$N(0,10^{2})$ prior on $\beta$s and $G(1,1)$ prior on $\lambda_{0}$. For each
replicates, 20,000 MCMC samples are drawn and last 10,000 samples are kept for
Bayesian inference. The full regression model includes x-coordinate,
y-coordinate and their cross effects. The average DIC and LPML values and the
selection percentage and model informations are shown in Table 1. This table
indicates that (i) the true model has the smallest average DIC value and the
largest average LPML value; (ii) the true model also has the largest selection
percentages under both DIC and LPML; and (iii) Model 4 has the second smallest
DIC value, the second largest LPML value, and the second largest selection
percentages. A partial explanation for (iii) is that both $Z_{1}(s)$ and
$Z_{2}(s)$ are correlated with the interaction term $Z_{1}(s)Z_{2}(s)$.

Table 1: Selection Results of Scenario 1

Model | Average DIC | Average LPML | DIC Selection % | LPML Selection %
---|---|---|---|---
Model 1($\beta_{1}$) | -1147.300 | 573.639 | 12 | 12
Model 2 ($\beta_{2}$) | -1093.641 | 546.810 | 0 | 0
Model 3 ($\beta_{3}$) | -1133.599 | 566.790 | 0 | 0
Model 4 ($\beta_{1},\beta_{2}$) | -1151.719 | 575.840 | 24 | 25
DGM ($\beta_{1},\beta_{3}$) | -1152.609 | 576.284 | 53 | 52
Model 5 ($\beta_{2},\beta_{3}$) | -1145.189 | 572.573 | 5 | 5
Model 6 ($\beta_{1},\beta_{2},\beta_{3}$) | -1151.656 | 575.788 | 6 | 6

11footnotetext: True Model

Furthermore, we also show the boxplots of the DIC difference and the LPML
difference between the true model and each of the candidate models in Figure
1. The medians and the inter quartile ranges (IQRs) of the DIC differences and
the LPML differences are show in Table 2

Table 2: Median and inter quartile ranges (IQRs) of the DIC and the LPML differences for Scenario 1 | DIC | LPML
---|---|---
Model 1 | 4.033 (0.517,7.339) | -2.009 (-0.253, -3.654)
Model 2 | 57.98 (48.03,69.35) | -28.97 (-24.00, -34.66)
Model 3 | 17.07 (11.60,26.40) | -8.520 (-5.789, -13.18)
Model 4 | 0.658 ( -0.378,1.798) | -0.339 (0.197, -0.887)
Model 5 | 6.541 (3.48,10.95) | -3.281 (-1.735, -5.481)
Model 6 | 1.244 (0.576,1.771) | -0.635 (-0.298, -0.908)

Except for Models 4 the first quartiles of DIC differences are above 0 while
the third quartiles of the LPML difference are below 0, which indicate that
the true model outperform models 1, 2, 3, 5, and 6 according to the DIC and
LPML values. However, the medians of the DIC and LPML differences are close to
0 for both Models 4 since both $Z_{1}(s)$ and $Z_{2}(s)$ are correlated with
the interaction term $Z_{1}(s)Z_{2}(s)$.

Figure 1: DIC and LPML Differences of Scenario 1

Furthermore, the intensity of DGM is
$\lambda(s)=\lambda_{0}\exp(\beta_{1}\times Z^{2}_{1}(s))$, where $Z_{1}(s)$
is the x-coordinate. We choose $\lambda_{0}=50$ and $\beta_{1}=4$. We generate
data on $[0,1]\times[0,1]$ and the average number of points over 100
replicates is 412. We compare the data generation model with three different
models: Model 1: $\lambda(s)=\lambda_{0}\exp(\beta_{1}\times Z_{1}(s))$; Model
2: $\lambda(s)=\lambda_{0}\exp(\beta_{2}\times Z_{2}(s))$, where $Z_{2}(s)$ is
the y-coordinate; Model 3: $\lambda(s)=\lambda_{0}\exp(\beta_{1}\times
Z_{1}(s)+\beta_{2}\times Z_{2}(s))$. We have the same prior distributions and
same number of MCMC samples with scenario 1

The average DIC and LPML values and the selection percentage and model
informations are shown in Table 3

Table 3: Selection Results of Scenario 2

Model | Average DIC | Average LPML | DIC Selection % | LPML Selection %
---|---|---|---|---
DGM | -4702.209 | 2351.101 | 94 | 94
Model 1 | -4690.384 | 2345.186 | 6 | 6
Model 2 | -4145.1794 | 2072.586 | 0 | 0
Model 3 | -4689.554 | 2344.769 | 0 | 0

11footnotetext: True Model

Furthermore, the medians and the inter quartile ranges (IQRs) of the DIC
differences and the LPML differences are show in Table 4

Table 4: Median and inter quartile ranges (IQRs) of the DIC and the LPML differences for Scenario 2 | DIC | LPML
---|---|---
Model 1 | 10.729 (6.676,17.11) | -5.365 (-8.562, -3.335)
Model 2 | 551.8 (531.1,584.3) | -275.9 (-292.2, -1265.5)
Model 3 | 12.16 (7.26,17.83) | -6.079 (-8,941,-3.628)

From the results in Table 3 and Table 4, our proposed criteria can effectively
select the true model and the difference of DIC and LPML between true model
and candidate models is significant.

The intensity function of DGM in scenario 3 is
$\lambda(s)=\exp(\beta_{1}\times Z_{1}(s)+\beta_{2}\times Z_{2}(s))$.
$\bm{Z}_{i},i=1,\cdots,4$ are generated from Gaussian random fields based on R
package RandomFields. Our choice of Gaussian random fields with mean 1 and the
covariance with exponential kernel model with variance 1, scale parameter 1
and nugget variance 0.2. In this simulation we choose $\beta_{1}=2$,
$\beta_{2}=1$, and $\lambda_{0}=1$ and generate data on $[0,1]\times[0,1]$. We
have $100\times 100$ grids on $[0,1]\times[0,1]$. Based on intensity surface
on $[0,1]\times[0,1]$, we independently generate locations in each grids. We
generate 100 data sets under this setting, and the average number of points in
each data set is 544. We have the same prior distributions and the same number
of MCMC samples with previous simulation. The average DIC and LPML values and
the selection percentage are shown in Table 6. From this table, we see that
the true model has the smallest DIC value, the largest LPML value, and the
highest selection percentages according to both DIC and LPML.

Table 5: Selection Results of Scenario 3 Model | Average DIC | Average LPML | DIC Selection % | LPML Selection %
---|---|---|---|---
Model 1 ($\beta_{1}$) | -8877.358 | 4438.677 | 0 | 0
Model 2 ($\beta_{2}$) | -6402.960 | 3201.478 | 0 | 0
Model 3 ($\beta_{3}$) | -5590.344 | 2795.170 | 0 | 0
Model 4 ($\beta_{4}$) | -5521.848 | 2760.921 | 0 | 0
DGM ($\beta_{1},\beta_{2}$) | -9303.699 | 4651.845 | 74 | 74
Model 5 ($\beta_{1},\beta_{3}$) | -8890.624 | 4445.308 | 0 | 0
Model 6 ($\beta_{1},\beta_{4}$) | -8879.830 | 4439.910 | 0 | 0
Model 7 ($\beta_{2},\beta_{3}$) | -6439.253 | 3219.623 | 0 | 0
Model 8 ($\beta_{2},\beta_{4}$) | -6409.395 | 3204.694 | 0 | 0
Model 9 ($\beta_{3},\beta_{4}$) | -5577.677 | 2788.834 | 0 | 0
Model 10 ($\beta_{1},\beta_{2},\beta_{3}$) | -9302.598 | 4651.292 | 15 | 15
Model 11 ($\beta_{1},\beta_{2},\beta_{4}$) | -9302.500 | 4651.243 | 9 | 9
Model 12 ($\beta_{1},\beta_{3},\beta_{4}$) | -8891.348 | 4445.667 | 0 | 0
Model 13 ($\beta_{2},\beta_{3},\beta_{4}$) | -6443.625 | 3221.808 | 0 | 0
Model 14 ($\beta_{1},\beta_{2},\beta_{3},\beta_{4}$) | -9301.493 | 4650.738 | 2 | 2

Furthermore, we also show the boxplots of the DIC difference and the LPML
difference between the true model and each of the candidate models in Figure
2. The medians and the inter quartile ranges (IQRs) of the DIC differences in
the top panel of Figure 2 are 1.608 (0.651, 2.012), 1.698 (0.826, 1.968), and
2.944 (1.492, 3.501) for Models 10, 11 and 14, respectively; and the medians
and the inter quartile ranges (IQRs) of the LPML differences in the bottom
panel of Figure 2 are -0.805 (-0.412, -0.985), -0.849 (-0.412, -0.985), and
-1.475 (-0.754, -1.758) for Models 10, 11 and 14, respectively. Unlike
scenario 3, the whole boxes of the DIC differences are above 0 while the boxes
of the LPML differences are below 0. Thus, there is an overwhelming evidence
in favor of the true model according both DIC and LPML.

Figure 2: DIC and LPML Differences of Scenario 3

Finally, the intensity of DGM of scenario 4 is
$\lambda(s)=\exp(\beta_{1}\times Z_{1}(s)+\beta_{2}\times Z_{2}(s)+W(s))$,
$\bm{Z}_{i},i=1,\cdots,3$ are generated from uniform distribution $U(0,1)$.
$W(s)$ is generated from Gaussian random with mean 0 and the covariance with
exponential kernel model with variance 1, scale parameter 1 and nugget
variance 0.2, which is unobserved in this simulation. In this simulation we
choose $\beta_{1}=4$, $\beta_{2}=4$, and $\lambda_{0}=1$ and generate data on
$[0,1]\times[0,1]$. We generate 100 data sets under this setting, and the
average number of points in each data set is 182. The full model is this
simulation is $\lambda(s)=\exp(\beta_{1}\times Z_{1}(s)+\beta_{2}\times
Z_{2}(s)+\beta_{3}\times Z_{3}(s))$. We have $100\times 100$ grids on
$[0,1]\times[0,1]$. Based on intensity surface on $[0,1]\times[0,1]$, we
independently generate locations in each grids. We have the same prior
distributions and same number of MCMC samples with previous scenarios. The
average DIC and LPML values and the selection percentage are shown in Table 6.
From this table, we see that the true model has the smallest DIC value, the
largest LPML value, and the highest selection percentages according to both
DIC and LPML. From this table, we see that the true model has the smallest DIC
value, the largest LPML value, and the highest selection percentages according
to both DIC and LPML.

Table 6: Selection Results of Scenario 4 Model | Average DIC | Average LPML | DIC Selection % | LPML Selection %
---|---|---|---|---
Model 1 ($\beta_{1}$) | -1694.704 | 847.344 | 0 | 0
Model 2 ($\beta_{2}$) | -1696.554 | 848.269 | 0 | 0
Model 3 ($\beta_{3}$) | -1459.131 | 729.560 | 0 | 0
DGM ($\beta_{1},\beta_{2}$) | -1873.874 | 936.921 | 88 | 88
Model 4 ($\beta_{1},\beta_{3}$) | -1693.393 | 846.683 | 0 | 0
Model 5 ($\beta_{2},\beta_{3}$) | -1695.379 | 847.676 | 0 | 0
Model 6 ($\beta_{1},\beta_{2},\beta_{3}$) | -1872.952 | 936.454 | 12 | 12

##  5 Real Data Analysis

###  5.1 USGS Earthquake Data

The earthquake data from USGS, the Earthquake Hazards Program of United States
Geological Survey (USGS), can be accessed via
https://earthquake.usgs.gov/earthquakes/. The data set we consider is composed
of conterminous U.S. earthquakes which have magnitude over 2.5 from 10-30-2018
to 11-27-2018. From USGS, conterminous U.S. refers to a rectangular region
including the lower 48 states and surrounding areas which are outside the
Conterminous U.S.. The total number of earthquakes is 155. The map of the
locations of earthquakes is shown in Figure 3.

Figure 3: Map of Locations of the Earthquake

In order to properly examine the relationship of these locations, we transform
the latitude and the longitude of the earthquakes to UTM (Universal Transverse
Mercator) coordinate system and then scale UTM coordinate system to a
$[0,1]\times[0,1]$ square. The locations of earthquakes in a unit square are
shown in Figure 4.

Figure 4: Map of Locations of the Earthquake in Unit Square

The spatial covariates used in this analysis include $x$-coordinates which are
transferred by longitudes, $y$-coordinates which are transferred by latitudes,
and the distance to New Madrid Seismic Zone Tuttle . (2005) which has
coordinates (36.58 N, 89.58 W). Similarly, we transform the latitude and the
longitude of the New Madrid Seismic Zone to UTM system and scale it to the
same unit square. New Madrid Seismic Zone is a major seismic zone and a
prolific source of intraplate earthquakes in the southern and midwestern
United States. The full spatial intensity model for the spatial point process
is given by

$\displaystyle\lambda(s)=\lambda_{0}\exp(\beta_{1}Z_{1}(s)+\beta_{2}Z_{2}(s)+\beta_{3}Z_{3}(s)),$
(20)

where $Z_{1}(s)$ is $x$-coordinate, $Z_{2}(s)$ is $y$-coordinate, and
$Z_{3}(s)$ is the distance to the New Madrid Seismic Zone. Including the
homogeneous spatial Poisson process model, we have 8 candidate models. With
the thinning interval to be 10, 2,000 samples are kept for calculation after a
burn-in of 30,000 samples using NIMBLE in R for each models. The prior for
$\beta$s is $N(0,10^{2})$ and the prior for $\lambda_{0}$ is $G(1,1)$. The DIC
and LPML values of these 8 candidate models are given in Table 7.

Table 7: DIC and LPML of Earthquake Data Model | DIC | LPML | Model | DIC | LPML
---|---|---|---|---|---
$\beta_{1},\beta_{2},\beta_{3}$ | -1268.687 | 634.329 | $\beta_{1},\beta_{2}$ | -1261.741 | 630.864
$\beta_{1},\beta_{3}$ | -1239.627 | 619.8027 | $\beta_{2},\beta_{3}$ | -1218.654 | 609.317
$\beta_{1}$ | -1179.762 | 589.878 | $\beta_{2}$ | -1219.484 | 609.7372
$\beta_{3}$ | -1206.92 | 603.452 | Homogenous | -1192.528 | 596.262

The DIC and LPML values in Table 7 suggest that the model with
$\beta_{1},\beta_{2},\beta_{3}$ is the best. The posterior means, the
posterior standard deviations and the 95% highest posterior density (HPD)
intervals Chen Shao (1999) under the best model are reported in Table 8.

Table 8: Best Model Parameter Estimation of Earthquake Data Parameter | Posterior Mean | Posterior Standard Deviation | HPD interval
---|---|---|---
$\beta_{1}$ | -1.975 | 0.352 | $(-2.699,-1.343)$
$\beta_{2}$ | 5.721 | 1.161 | $(3.621,7.884)$
$\beta_{3}$ | -2.744 | 1.213 | $(-4.909,-0.512)$
$\lambda_{0}$ | 37.850 | 5.491 | $(28.314,49.394)$

From Table 8, we see that the occurrence of earthquake has a higher intensity
around the New Madrid Seismic Zone and the southwestern area of U.S., which is
consistent with the findings in the literature Fuller (1988).

###  5.2 Forest of Barro Colorado Island (BCI)

BCI data have be analyzed in the literature Leininger . (2017); Thurman .
(2015); Yue Loh (2015); Thurman Zhu (2014). The data are obtained from a long-
term ecological monitoring program, which contain 400,000 individual trees
since the 1980s at the Barro Colorado Island (BCI) in central Panama Hubbell .
(1999); Condit (1998); Condit . (2012). We choose the same species like
Thurman Zhu (2014) in our data analysis. We model the intensity of the
B.pendula trees as a log-linear function of the incidences of six other tree
species (Eugenia nesiotica, Eugenia oerstediana, Piper cordulatum, Protium
panamense, Sorocea affinis, and Talisia croatii). The map of the locations of
B.pendula tress mentioned above is shown in Figure 5.

Figure 5: Locations Map of B.pendula

We choose three different pixel resolutions, $5m\times 5m$, $10m\times 10m$,
and $20m\times 20m$ for the incidences of six other tree species (6 spatial
covariates). The full spatial intensity model for the spatial point process is
given by

$\displaystyle\lambda(s)=\lambda_{0}\exp(\beta_{1}Z_{1}(s)+\cdots+\beta_{6}Z_{6}(s)),$
(21)

where $Z_{i}(s)$ is the $i$-th specie’s pixel value on location $s$. Including
the homogeneous spatial Poisson process model, we have 64 candidate models.
8,000 samples are kept for calculation after a burn-in of 10,000 samples using
NIMBLE in R for each model. The prior for $\beta$s is $N(0,10^{2})$ and the
prior for $\lambda_{0}$ is $G(1,1)$. The DIC and the LPML values of the best
models for the three pixel resolutions are given in Table 9.

Table 9: Model Comparison for Different Pixel Resolutions Pixel Resolution | Best Model | DIC | LPML
---|---|---|---
$5m\times 5m$ | $(\beta_{4},\beta_{5})$ | -57521.48 | 28760.74
$10m\times 10m$ | $(\beta_{1},\beta_{3},\beta_{4},\beta_{5})$ | -58255.15 | 29127.57
$20m\times 20m$ | $(\beta_{3},\beta_{4},\beta_{5})$ | -58188.46 | 29094.22

The combination of the covariate model and the resolution with the smallest
DIC value and the largest LPML value is the model with
$(\beta_{1},\beta_{3},\beta_{4},\beta_{5})$ plus the $10m\times 10m$
resolution. The posterior means, the posterior standard deviations and the 95%
HPD intervals under this best combined covariate model and resolution are
reported in Table 10. Figure 6 shows the plots of 6 different species in the
$10m\times 10m$ pixel resolution.

Figure 6: $10m\times 10m$ Pixel Resolution Covariates Map (Red indicate larger values) |  |
---|---|---
|  |
Table 10: Posterior Estimation Results of BCI Data Parameter | Posterior Mean | Posterior SD | 95%HPD interval
---|---|---|---
$\beta_{1}$ | 0.2986 | 0.0275 | (0.2415,0.3511)
$\beta_{3}$ | 0.1083 | 0.0037 | (0.1013,0.1156)
$\beta_{4}$ | 0.2305 | 0.0064 | (0.2180,0.2420)
$\beta_{5}$ | 0.2825 | 0.0122 | (0.2595,0.3053)
$\lambda_{0}$ | 1200.4259 | 25.3233 | (1157.8434,1254.6554)

From Table 10, we see that (i) the B.pendula trees seem to be attracted to
Eugenia oerstediana, Piper cordulatum, Protium panamense, and Sorocea affinis;
and (ii) Eugenia nesiotica and Talisia croatii have no effects on B.pendula
trees.

##  6 Discussion

In this paper, we have developed the Bayesian parameter estimation and model
selection in Poisson point process models in regression context. Our
simulation results indicate that our proposed method is more likely to select
the correct model and obtain the estimators that are closer to the true values
if we had fit the true model. Like most model selection methods, our proposed
method needs calculating DIC and LPML for all candidate models. Developing a
computational algorithm as in Chen . (2008) is an important future work. In
addition, extending our proposed criterion to a more complicated spatial point
process like Gibbs point process Moller Waagepetersen (2007) should be a
natural extension of our work. Furthermore, another future work is to
investigate the theoretical properties of the proposed procedure like Li .
(2017). Moreover, the development of the predictive distribution-based
criterion like WAIC Vehtari . (2017) or the L-measure Ibrahim . (2001) for
spatial point process is another interesting future work.

## Acknowledgement

Dr. Chen’s research was partially supported by NIH grants #GM70335 and
#P01CA142538. Dr. Hu’s research was supported by Dean’s office of College of
Liberal Arts and Sciences in University of Connecticut.

## References

  * Akaike (1974) akaike1974newAkaike, H.  1974\.  A new look at the statistical model identification A new look at the statistical model identification.  IEEE transactions on automatic control196716–723.
  * Banerjee . (2014) banerjee2014hierarchicalBanerjee, S., Carlin, BP. Gelfand, AE.  2014\.  Hierarchical modeling and analysis for spatial data Hierarchical modeling and analysis for spatial data.  Crc Press.
  * Beale . (2010) beale2010regressionBeale, CM., Lennon, JJ., Yearsley, JM., Brewer, MJ. Elston, DA.  2010\.  Regression analysis of spatial data Regression analysis of spatial data.  Ecology letters132246–264.
  * Chen . (2008) chen2008bayesianChen, MH., Huang, L., Ibrahim, JG. Kim, S.  2008\.  Bayesian variable selection and computation for generalized linear models with conjugate priors Bayesian variable selection and computation for generalized linear models with conjugate priors.  Bayesian analysis (Online)33585.
  * Chen Shao (1999) chen1999monteChen, MH. Shao, QM.  1999\.  Monte Carlo estimation of Bayesian credible and HPD intervals Monte carlo estimation of bayesian credible and hpd intervals.  Journal of Computational and Graphical Statistics8169–92.
  * Chen . (2012) chen2012monteChen, MH., Shao, QM. Ibrahim, JG.  2012\.  Monte Carlo methods in Bayesian computation Monte carlo methods in bayesian computation.  Springer Science & Business Media.
  * Condit (1998) condit1998tropicalCondit, R.  1998\.  Tropical forest census plots: methods and results from Barro Colorado Island, Panama and a comparison with other plots Tropical forest census plots: methods and results from barro colorado island, panama and a comparison with other plots.  Springer Science & Business Media.
  * Condit . (2012) condit2012datasetCondit, R., Lao, S., Pérez, R., Dolins, SB., Foster, R. Hubbell, S.  2012\.  [Dataset:] Barro Colorado Forest Census Plot Data (Version 2012) [dataset:] barro colorado forest census plot data (version 2012).
  * Cressie (1993) cressie1993statisticsCressie, NA.  1993\.  Statistics for spatial data Statistics for spatial data.  Wiley Online Library.
  * de Valpine . (2017) de2017programmingde Valpine, P., Turek, D., Paciorek, CJ., Anderson-Bergman, C., Lang, DT. Bodik, R.  2017\.  Programming with Models: Writing Statistical Algorithms for General Model Structures with NIMBLE Programming with models: Writing statistical algorithms for general model structures with NIMBLE.  Journal of Computational and Graphical Statistics262403–413.
  * Diggle (2013) diggle2013statisticalDiggle, PJ.  2013\.  Statistical analysis of spatial and spatio-temporal point patterns Statistical analysis of spatial and spatio-temporal point patterns.  CRC Press.
  * Fuller (1988) fuller1988newFuller, ML.  1988\.  The new Madrid earthquake The new madrid earthquake ( 494).  Central United States Earthquake Consortium.
  * Gelfand Dey (1994) gelfand1994bayesianGelfand, AE. Dey, DK.  1994\.  Bayesian model choice: asymptotics and exact calculations Bayesian model choice: asymptotics and exact calculations.  Journal of the Royal Statistical Society. Series B (Methodological)501–514.
  * Hoerl Kennard (1970) hoerl1970ridgeHoerl, AE. Kennard, RW.  1970\.  Ridge regression: Biased estimation for nonorthogonal problems Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics12155–67.
  * Hubbell . (1999) hubbell1999lightHubbell, SP., Foster, RB., O’Brien, ST., Harms, KE., Condit, R., Wechsler, B.De Lao, SL.  1999\.  Light-gap disturbances, recruitment limitation, and tree diversity in a neotropical forest Light-gap disturbances, recruitment limitation, and tree diversity in a neotropical forest.  Science2835401554–557.
  * Ibrahim . (2001) ibrahim2001criterionIbrahim, JG., Chen, MH. Sinha, D.  2001\.  Criterion-based methods for Bayesian model assessment Criterion-based methods for bayesian model assessment.  Statistica Sinica419–443.
  * Ishwaran . (2005) ishwaran2005spikeIshwaran, H., Rao, JS. .  2005\.  Spike and slab variable selection: frequentist and Bayesian strategies Spike and slab variable selection: frequentist and bayesian strategies.  The Annals of Statistics332730–773.
  * Leininger . (2017) leininger2017bayesianLeininger, TJ., Gelfand, AE. .  2017\.  Bayesian inference and model assessment for spatial point patterns using posterior predictive samples Bayesian inference and model assessment for spatial point patterns using posterior predictive samples.  Bayesian Analysis1211–30.
  * Li . (2017) li2017devianceLi, Y., Jun, Y. Zeng, T.  2017\.  Deviance information criterion for Bayesian model selection: Justification and variation Deviance information criterion for bayesian model selection: Justification and variation.
  * Moller Waagepetersen (2003) moller2003statisticalMoller, J. Waagepetersen, RP.  2003\.  Statistical inference and simulation for spatial point processes Statistical inference and simulation for spatial point processes.  CRC Press.
  * Moller Waagepetersen (2007) moller2007modernMoller, J. Waagepetersen, RP.  2007\.  Modern statistics for spatial point processes Modern statistics for spatial point processes.  Scandinavian Journal of Statistics344643–684.
  * Polson Scott (2010) polson2010shrinkPolson, NG. Scott, JG.  2010\.  Shrink globally, act locally: Sparse Bayesian regularization and prediction Shrink globally, act locally: Sparse bayesian regularization and prediction.  Bayesian statistics9501–538.
  * Reich . (2010) reich2010bayesianReich, BJ., Fuentes, M., Herring, AH. Evenson, KR.  2010\.  Bayesian variable selection for multivariate spatially varying coefficient regression Bayesian variable selection for multivariate spatially varying coefficient regression.  Biometrics663772–782.
  * Schwarz . (1978) schwarz1978estimatingSchwarz, G. .  1978\.  Estimating the dimension of a model Estimating the dimension of a model.  The annals of statistics62461–464.
  * Spiegelhalter . (2002) spiegelhalter2002bayesianSpiegelhalter, DJ., Best, NG., Carlin, BP. Van Der Linde, A.  2002\.  Bayesian measures of model complexity and fit Bayesian measures of model complexity and fit.  Journal of the Royal Statistical Society: Series B (Statistical Methodology)644583–639.
  * Thurman . (2015) thurman2015regularizedThurman, AL., Fu, R., Guan, Y. Zhu, J.  2015\.  Regularized estimating equations for model selection of clustered spatial point processes Regularized estimating equations for model selection of clustered spatial point processes.  Statistica Sinica173–188.
  * Thurman Zhu (2014) thurman2014variableThurman, AL. Zhu, J.  2014\.  Variable selection for spatial Poisson point processes via a regularization method Variable selection for spatial poisson point processes via a regularization method.  Statistical Methodology17113–125.
  * Tibshirani (1996) tibshirani1996regressionTibshirani, R.  1996\.  Regression shrinkage and selection via the lasso Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological)267–288.
  * Tuttle . (2005) tuttle2005evidenceTuttle, MP., Schweig, ES., Campbell, J., Thomas, PM., Sims, JD. Lafferty, RH.  2005\.  Evidence for New Madrid earthquakes in AD 300 and 2350 BC Evidence for new madrid earthquakes in ad 300 and 2350 bc.  Seismological Research Letters764489–501.
  * Vehtari . (2017) vehtari2017practicalVehtari, A., Gelman, A. Gabry, J.  2017\.  Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC Practical bayesian model evaluation using leave-one-out cross-validation and waic.  Statistics and Computing2751413–1432.
  * Waagepetersen (2004) waagepetersen2004convergenceWaagepetersen, R.  2004\.  Convergence of posteriors for discretized log Gaussian Cox processes Convergence of posteriors for discretized log gaussian cox processes.  Statistics & Probability Letters663229–235.
  * Wang Zhu (2009) wang2009variableWang, H. Zhu, J.  2009\.  Variable selection in spatial regression via penalized least squares Variable selection in spatial regression via penalized least squares.  Canadian Journal of Statistics374607–624.
  * Yue Loh (2015) yue2015variableYue, Y. Loh, JM.  2015\.  Variable selection for inhomogeneous spatial point process models Variable selection for inhomogeneous spatial point process models.  Canadian Journal of Statistics432288–305.
  * Zou Hastie (2005) zou2005regularizationZou, H. Hastie, T.  2005\.  Regularization and variable selection via the elastic net Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society: Series B (Statistical Methodology)672301–320.