dppm: Fit Determinantal Point Process Model
In spatstat.core: Core Functionality of the 'spatstat' Family

View source: R/dppm.R

dppm	R Documentation

Fit Determinantal Point Process Model

Description

Fit a determinantal point process model to a point pattern.

Usage

  dppm(formula, family, data=NULL,
       ...,
       startpar = NULL,
       method = c("mincon", "clik2", "palm", "adapcl"),
       weightfun = NULL,
       control =  list(),
       algorithm,
       statistic = "K",
       statargs = list(),
       rmax = NULL,
       epsilon = 0.01,
       covfunargs = NULL,
       use.gam = FALSE,
       nd = NULL, eps = NULL)

Arguments

`formula`	A `formula` in the R language specifying the data (on the left side) and the form of the model to be fitted (on the right side). For a stationary model it suffices to provide a point pattern without a formula. See Details.
`family`	Information specifying the family of point processes to be used in the model. Typically one of the family functions `dppGauss`, `dppMatern`, `dppCauchy`, `dppBessel` or `dppPowerExp`. Alternatively a character string giving the name of a family function, or the result of calling one of the family functions. See Details.
`data`	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.
`startpar`	Named vector of starting parameter values for the optimization.
`method`	The fitting method. Either `"mincon"` for minimum contrast, `"clik2"` for second order composite likelihood, `"adapcl"` for adaptive second order composite likelihood, or `"palm"` for Palm likelihood. Partially matched.
`weightfun`	Optional weighting function w in the composite likelihoods 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 algorithm to be used to solve the fitting problem. If `method="mincon", "clik2"` or `"palm"` this argument is passed to the generic optimization function `optim` (renamed as the argument `method`) with default `"Nelder-Mead"`. If `method="adapcl"`) the argument is passed to the equation solver `nleqslv`, with default `"Bryden"`.
`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.
`epsilon`	Tuning parameter for the adaptive composite likelihood.
`covfunargs,use.gam,nd,eps`	Arguments passed to `ppm` when fitting the intensity.

Details

This function fits a determinantal point process model to a point pattern dataset as described in Lavancier et al. (2015).

The model to be fitted is specified by the arguments formula and family.

The argument formula should normally be a formula in the R language. The left hand side of the formula specifies the point pattern dataset to which the model should be fitted. This should be a single argument which may be a point pattern (object of class "ppp") or a quadrature scheme (object of class "quad"). The right hand side of the formula is called the trend and specifies the form of the logarithm of the intensity of the process. Alternatively the argument formula may be a point pattern or quadrature scheme, and the trend formula is taken to be ~1.

The argument family specifies the family of point processes to be used in the model. It is typically one of the family functions dppGauss, dppMatern, dppCauchy, dppBessel or dppPowerExp. Alternatively it may be a character string giving the name of a family function, or the result of calling one of the family functions. A family function belongs to class "detpointprocfamilyfun". The result of calling a family function is a point process family, which belongs to class "detpointprocfamily".

The algorithm first estimates the intensity function of the point process using ppm. If the trend formula is ~1 (the default if a point pattern or quadrature scheme is given rather than a "formula") 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 interaction parameters of the model are then fitted either by minimum contrast estimation, or by a composite likelihood method (maximum composite likelihood, maximum Palm likelihood, or by solving the adaptive composite likelihood estimating equation).

Minimum contrast:

If method = "mincon" (the default) interaction 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 in mincontrast.

For a homogeneous model ( trend = ~1 ) the empirical K-function of the data is computed using Kest, and the interaction parameters of the 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 interaction parameters of the model are estimated by the method of minimum contrast using the inhomogeneous K function. This two-step estimation procedure is heavily inspired by Waagepetersen (2007).

If statistic="pcf" then instead of using the K-function, the algorithm will use the pair correlation function pcf for homogeneous models and the inhomogeneous pair correlation function pcfinhom for inhomogeneous models. In this case, the smoothing parameters of the pair correlation can be controlled using the argument statargs, as shown in the Examples.

Additional arguments ... will be passed to clusterfit to control the minimum contrast fitting algorithm.

Composite likelihood:

If method = "clik2" the interaction parameters of the model will be fitted by maximising the second-order composite likelihood (Guan, 2006). The log composite likelihood is

∑[i,j] w(d[i,j]) log(ρ(d[i,j]; θ)) - (∑[i,j] w(d[i,j])) log(integral[D,D] w(||u-v||) ρ(||u-v||; θ) du dv)

where the sums are taken over all pairs of data points x[i], x[j] separated by a distance d[i,j] = ||x[i] - x[j]|| less than rmax, and the double integral is taken over all pairs of locations u,v in the spatial window of the data. Here ρ(d;θ) is the pair correlation function of the model with interaction parameters θ.

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 or NULL then the default is a threshold weight function, w(d) = 1(d ≤ R), where R is rmax/2.

Palm likelihood:

If method = "palm" the interaction parameters of the model will be fitted by maximising the Palm loglikelihood (Tanaka et al, 2008)

∑[i,j] w(x[i], x[j]) log(λ[P](x[j] | x[i]; θ) - integral[D] w(x[i], u) λ[P](u | x[i]; θ) du

with the same notation as above. Here λ[P](u|v;θ) is the Palm intensity of the model at location u given there is a point at v.

Adaptive Composite likelihood:

If method = "cladap" the clustering parameters of the model will be fitted by solving the adaptive second order composite likelihood estimating equation (Lavancier et al, 2021). The estimating function is

∑[u,v] w(epsilon |g(0; theta)-1|/(g(||u-v||; theta)-1)) g'(||u-v||; theta) / g(||u-v||; theta) - integral[D,D] w(epsilon |g(0; theta)-1|/(g(||u-v||; theta)-1)) g'(||u-v||; theta) rho(u) rho(v) du dv

where the sum is taken over all distinct pairs of points. Here g(d; theta) is the pair correlation function with parameters theta. The partial derivative with respect to theta is g'(d; theta), and rho(u) denotes the fitted intensity function of the model.

The tuning parameter epsilon is independent of the data. It can be specified by the argument epsilon and has default value 0.01.

The function w in the estimating function is a weighting function of bounded support [-1,1]. It is specified by the argument weightfun. If this is missing or NULL then the default is w(d) = 1(||d|| ≤ 1) exp(1/(d^2-1)) . The estimating equation is solved using the nonlinear equation solver nleqslv from the package nleqslv. The package nleqslv must be installed in order to use this option.

It is also possible to fix any parameters desired before the optimisation by specifying them as name=value in the call to the family function. See Examples.

Value

An object of class "dppm" representing the fitted model. There are methods for printing, plotting, predicting and simulating objects of this class.

Optimization algorithm

The following details allow greater control over the fitting procedure.

For the first three fitting methods (method="mincon", "clik2" and "palm"), the optimisation is performed by the generic optimisation algorithm optim. The behaviour of this algorithm can be modified using the arguments control and algorithm. Useful control arguments include trace, maxit and abstol (documented in the help for optim).

For method="adapcl", the estimating equation is solved using the nonlinear equation solver nleqslv from the package nleqslv. Arguments available for controlling the solver are documented in the help for nleqslv; they include control, globStrat, startparm for the initial estimates and algorithm for the method. The package nleqslv must be installed in order to use this option.

Author(s)

\spatstatAuthors

. Adaptive composite likelihood method contributed by Chiara Fend and modified by Adrian Baddeley.

References

Guan, Y. (2006) A composite likelihood approach in fitting spatial point process models. Journal of the American Statistical Association 101, 1502–1512.

Lavancier, F., \Moller, J. and Rubak, E. (2015) Determinantal point process models and statistical inference. Journal of the Royal Statistical Society, Series B 77, 853–977.

Lavancier, F., Poinas, A., and Waagepetersen, R. (2021) Adaptive estimating function inference for nonstationary determinantal point processes. Scandinavian Journal of Statistics, 48 (1), 87–107.

Tanaka, U., 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.

Examples

  jpines <- residualspaper$Fig1
  

  dppm(jpines ~ 1, dppGauss)

  dppm(jpines ~ 1, dppGauss, method="c")
  dppm(jpines ~ 1, dppGauss, method="p")
  dppm(jpines ~ 1, dppGauss, method="a")

  if(interactive()) {
     # Fixing the intensity at lambda=2 rather than the Poisson MLE 2.04:
     dppm(jpines ~ 1, dppGauss(lambda=2))

     # The following is quite slow (using K-function)
     dppm(jpines ~ x, dppMatern)
  }

   # much faster using pair correlation function
  dppm(jpines ~ x, dppMatern, statistic="pcf", statargs=list(stoyan=0.2))

  # Fixing the Matern shape parameter at nu=2 rather than estimating it:
  dppm(jpines ~ x, dppMatern(nu=2))

spatstat.core documentation built on May 18, 2022, 9:05 a.m.