kppm
Fit Cluster or Cox Point Process Model
Fit a homogeneous or inhomogeneous cluster process or Cox point process model to a point pattern.
Usage
kppm(X, …) # S3 method for formula
kppm(X,
clusters = c("Thomas","MatClust","Cauchy","VarGamma","LGCP"),
…,
data=NULL)
# S3 method for ppp
kppm(X,
trend = ~1,
clusters = c("Thomas","MatClust","Cauchy","VarGamma","LGCP"),
data = NULL,
...,
covariates=data,
subset,
method = c("mincon", "clik2", "palm"),
improve.type = c("none", "clik1", "wclik1", "quasi"),
improve.args = list(),
weightfun=NULL,
control=list(),
algorithm="Nelder-Mead",
statistic="K",
statargs=list(),
rmax = NULL,
covfunargs=NULL,
use.gam=FALSE,
nd=NULL, eps=NULL)
# S3 method for quad
kppm(X,
trend = ~1,
clusters = c("Thomas","MatClust","Cauchy","VarGamma","LGCP"),
data = NULL,
...,
covariates=data,
subset,
method = c("mincon", "clik2", "palm"),
improve.type = c("none", "clik1", "wclik1", "quasi"),
improve.args = list(),
weightfun=NULL,
control=list(),
algorithm="Nelder-Mead",
statistic="K",
statargs=list(),
rmax = NULL,
covfunargs=NULL,
use.gam=FALSE,
nd=NULL, eps=NULL)
Arguments
- X
A point pattern dataset (object of class
"ppp"
or"quad"
) to which the model should be fitted, or aformula
in the R language defining the model. See Details.- trend
An R formula, with no left hand side, specifying the form of the log intensity.
- clusters
Character string determining the cluster model. Partially matched. Options are
"Thomas"
,"MatClust"
,"Cauchy"
,"VarGamma"
and"LGCP"
.- data,covariates
The values of spatial covariates (other than the Cartesian coordinates) required by the model. A named list of pixel images, functions, windows, tessellations or numeric constants.
- …
Additional arguments. See Details.
- subset
Optional. A subset of the spatial domain, to which the model-fitting should be restricted. A window (object of class
"owin"
) or a logical-valued pixel image (object of class"im"
), or an expression (possibly involving the names of entries indata
) which can be evaluated to yield a window or pixel image.- method
The fitting method. Either
"mincon"
for minimum contrast,"clik2"
for second order composite likelihood, or"palm"
for Palm likelihood. Partially matched.- improve.type
Method for updating the initial estimate of the trend. Initially the trend is estimated as if the process is an inhomogeneous Poisson process. The default,
improve.type = "none"
, is to use this initial estimate. Otherwise, the trend estimate is updated byimprove.kppm
, using information about the pair correlation function. Options are"clik1"
(first order composite likelihood, essentially equivalent to"none"
),"wclik1"
(weighted first order composite likelihood) and"quasi"
(quasi likelihood).- improve.args
Additional arguments passed to
improve.kppm
whenimprove.type != "none"
. See Details.- weightfun
Optional weighting function \(w\) in the composite likelihood or Palm likelihood. A
function
in the R language. See Details.- control
List of control parameters passed to the optimization function
optim
.- algorithm
Character string determining the mathematical optimisation algorithm to be used by
optim
. See the argumentmethod
ofoptim
.- statistic
Name of the summary statistic to be used for minimum contrast estimation: either
"K"
or"pcf"
.- statargs
Optional list of arguments to be used when calculating the
statistic
. See Details.- rmax
Maximum value of interpoint distance to use in the composite likelihood.
- covfunargs,use.gam,nd,eps
Arguments passed to
ppm
when fitting the intensity.
Details
This function fits a clustered point process model to the
point pattern dataset X
.
The model may be either a Neyman-Scott cluster process
or another Cox process.
The type of model is determined by the argument clusters
.
Currently the options
are clusters="Thomas"
for the Thomas process,
clusters="MatClust"
for the Matern cluster process,
clusters="Cauchy"
for the Neyman-Scott cluster process
with Cauchy kernel,
clusters="VarGamma"
for the Neyman-Scott cluster process
with Variance Gamma kernel (requires an additional argument nu
to be passed through the dots; see rVarGamma
for details),
and clusters="LGCP"
for the log-Gaussian Cox process (may
require additional arguments passed through …
; see
rLGCP
for details on argument names).
The first four models are Neyman-Scott cluster processes.
The algorithm first estimates the intensity function
of the point process using ppm
.
The argument X
may be a point pattern
(object of class "ppp"
) or a quadrature scheme
(object of class "quad"
). The intensity is specified by
the trend
argument.
If the trend formula is ~1
(the default)
then the model is homogeneous. The algorithm begins by
estimating the intensity as the number of points divided by
the area of the window.
Otherwise, the model is inhomogeneous.
The algorithm begins by fitting a Poisson process with log intensity
of the form specified by the formula trend
.
(See ppm
for further explanation).
The argument X
may also be a formula
in the
R language. The right hand side of the formula gives the
trend
as described above. The left hand side of the formula
gives the point pattern dataset to which the model should be fitted.
If improve.type="none"
this is the final estimate of the
intensity. Otherwise, the intensity estimate is updated, as explained in
improve.kppm
. Additional arguments to
improve.kppm
are passed as a named list in
improve.args
.
The clustering parameters of the model are then fitted either by minimum contrast estimation, or by maximum composite likelihood.
- Minimum contrast:
If
method = "mincon"
(the default) clustering parameters of the model will be fitted by minimum contrast estimation, that is, by matching the theoretical \(K\)-function of the model to the empirical \(K\)-function of the data, as explained inmincontrast
.For a homogeneous model (
trend = ~1
) the empirical \(K\)-function of the data is computed usingKest
, and the parameters of the cluster model are estimated by the method of minimum contrast.For an inhomogeneous model, the inhomogeneous \(K\) function is estimated by
Kinhom
using the fitted intensity. Then the parameters of the cluster model are estimated by the method of minimum contrast using the inhomogeneous \(K\) function. This two-step estimation procedure is due to Waagepetersen (2007).If
statistic="pcf"
then instead of using the \(K\)-function, the algorithm will use the pair correlation functionpcf
for homogeneous models and the inhomogeneous pair correlation functionpcfinhom
for inhomogeneous models. In this case, the smoothing parameters of the pair correlation can be controlled using the argumentstatargs
, as shown in the Examples.Additional arguments
…
will be passed tomincontrast
to control the minimum contrast fitting algorithm.- Composite likelihood:
If
method = "clik2"
the clustering parameters of the model will be fitted by maximising the second-order composite likelihood (Guan, 2006). The log composite likelihood is $$ \sum_{i,j} w(d_{ij}) \log\rho(d_{ij}; \theta) - \left( \sum_{i,j} w(d_{ij}) \right) \log \int_D \int_D w(\|u-v\|) \rho(\|u-v\|; \theta)\, du\, dv $$ where the sums are taken over all pairs of data points \(x_i, x_j\) separated by a distance \(d_{ij} = \| x_i - x_j\|\) less thanrmax
, and the double integral is taken over all pairs of locations \(u,v\) in the spatial window of the data. Here \(\rho(d;\theta)\) is the pair correlation function of the model with cluster parameters \(\theta\).The function \(w\) in the composite likelihood is a weighting function and may be chosen arbitrarily. It is specified by the argument
weightfun
. If this is missing orNULL
then the default is a threshold weight function, \(w(d) = 1(d \le R)\), where \(R\) isrmax/2
.- Palm likelihood:
If
method = "palm"
the clustering parameters of the model will be fitted by maximising the Palm loglikelihood (Tanaka et al, 2008) $$ \sum_{i,j} w(x_i, x_j) \log \lambda_P(x_j \mid x_i; \theta) - \int_D w(x_i, u) \lambda_P(u \mid x_i; \theta) {\rm d} u $$ with the same notation as above. Here \(\lambda_P(u|v;\theta\) is the Palm intensity of the model at location \(u\) given there is a point at \(v\).
In all three methods, the optimisation is performed by the generic
optimisation algorithm optim
.
The behaviour of this algorithm can be modified using the
argument control
.
Useful control arguments include
trace
, maxit
and abstol
(documented in the help for optim
).
Fitting the LGCP model requires the RandomFields package, except in the default case where the exponential covariance is assumed.
Value
An object of class "kppm"
representing the fitted model.
There are methods for printing, plotting, predicting, simulating
and updating objects of this class.
Log-Gaussian Cox Models
To fit a log-Gaussian Cox model with non-exponential covariance,
specify clusters="LGCP"
and use additional arguments
to specify the covariance structure. These additional arguments can
be given individually in the call to kppm
, or they can be
collected together in a list called covmodel
.
For example a Matern model with parameter \(\nu=0.5\) could be specified
either by kppm(X, clusters="LGCP", model="matern", nu=0.5)
or by
kppm(X, clusters="LGCP", covmodel=list(model="matern", nu=0.5))
.
The argument model
specifies the type of covariance
model: the default is model="exp"
for an exponential covariance.
Alternatives include "matern"
, "cauchy"
and "spheric"
.
Model names correspond to functions beginning with RM
in the
RandomFields package: for example model="matern"
corresponds to the function RMmatern
in the
RandomFields package.
Additional arguments are passed to the
relevant function in the RandomFields package:
for example if model="matern"
then the additional argument
nu
is required, and is passed to the function
RMmatern
in the RandomFields package.
Note that it is not possible to use anisotropic covariance models
because the kppm
technique assumes the pair correlation function
is isotropic.
Error and warning messages
See ppm.ppp
for a list of common error messages
and warnings originating from the first stage of model-fitting.
References
Guan, Y. (2006) A composite likelihood approach in fitting spatial point process models. Journal of the American Statistical Association 101, 1502--1512.
Jalilian, A., Guan, Y. and Waagepetersen, R. (2012) Decomposition of variance for spatial Cox processes. Scandinavian Journal of Statistics 40, 119--137.
Tanaka, U. and Ogata, Y. and Stoyan, D. (2008) Parameter estimation and model selection for Neyman-Scott point processes. Biometrical Journal 50, 43--57.
Waagepetersen, R. (2007) An estimating function approach to inference for inhomogeneous Neyman-Scott processes. Biometrics 63, 252--258.
See Also
Methods for kppm
objects:
plot.kppm
,
fitted.kppm
,
predict.kppm
,
simulate.kppm
,
update.kppm
,
vcov.kppm
,
methods.kppm
,
as.ppm.kppm
,
as.fv.kppm
,
Kmodel.kppm
,
pcfmodel.kppm
.
Minimum contrast fitting algorithm:
mincontrast
.
Alternative fitting algorithms:
thomas.estK
,
matclust.estK
,
lgcp.estK
,
cauchy.estK
,
vargamma.estK
,
thomas.estpcf
,
matclust.estpcf
,
lgcp.estpcf
,
cauchy.estpcf
,
vargamma.estpcf
,
Summary statistics:
Kest
,
Kinhom
,
pcf
,
pcfinhom
.
See also ppm
Examples
# NOT RUN {
# method for point patterns
kppm(redwood, ~1, "Thomas")
# method for formulas
kppm(redwood ~ 1, "Thomas")
kppm(redwood ~ 1, "Thomas", method="c")
kppm(redwood ~ 1, "Thomas", method="p")
kppm(redwood ~ x, "MatClust")
kppm(redwood ~ x, "MatClust", statistic="pcf", statargs=list(stoyan=0.2))
kppm(redwood ~ x, cluster="Cauchy", statistic="K")
kppm(redwood, cluster="VarGamma", nu = 0.5, statistic="pcf")
# LGCP models
kppm(redwood ~ 1, "LGCP", statistic="pcf")
if(require("RandomFields")) {
kppm(redwood ~ x, "LGCP", statistic="pcf",
model="matern", nu=0.3,
control=list(maxit=10))
}
# fit with composite likelihood method
kppm(redwood ~ x, "VarGamma", method="clik2", nu.ker=-3/8)
# fit intensity with quasi-likelihood method
kppm(redwood ~ x, "Thomas", improve.type = "quasi")
# }