Version:	3.4-0
Date:	2025-07-23
Title:	Parametric Statistical Modelling and Inference for the 'spatstat' Family
Maintainer:	Adrian Baddeley <Adrian.Baddeley@curtin.edu.au>
Depends:	R (≥ 3.5.0), spatstat.data (≥ 3.1-6), spatstat.univar (≥ 3.1-4), spatstat.geom (≥ 3.5-0), spatstat.random (≥ 3.4-1), spatstat.explore (≥ 3.5-2), stats, graphics, grDevices, utils, methods, nlme, rpart
Imports:	spatstat.utils (≥ 3.1-5), spatstat.sparse (≥ 3.1-0), mgcv, Matrix, abind, tensor, goftest (≥ 1.2-2)
Suggests:	sm, gsl, locfit, spatial, fftwtools (≥ 0.9-8), nleqslv, glmnet, spatstat.linnet (≥ 3.2-6), spatstat (≥ 3.3)
Description:	Functionality for parametric statistical modelling and inference for spatial data, mainly spatial point patterns, in the 'spatstat' family of packages. (Excludes analysis of spatial data on a linear network, which is covered by the separate package 'spatstat.linnet'.) Supports parametric modelling, formal statistical inference, and model validation. Parametric models include Poisson point processes, Cox point processes, Neyman-Scott cluster processes, Gibbs point processes and determinantal point processes. Models can be fitted to data using maximum likelihood, maximum pseudolikelihood, maximum composite likelihood and the method of minimum contrast. Fitted models can be simulated and predicted. Formal inference includes hypothesis tests (quadrat counting tests, Cressie-Read tests, Clark-Evans test, Berman test, Diggle-Cressie-Loosmore-Ford test, scan test, studentised permutation test, segregation test, ANOVA tests of fitted models, adjusted composite likelihood ratio test, envelope tests, Dao-Genton test, balanced independent two-stage test), confidence intervals for parameters, and prediction intervals for point counts. Model validation techniques include leverage, influence, partial residuals, added variable plots, diagnostic plots, pseudoscore residual plots, model compensators and Q-Q plots.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	http://spatstat.org/
NeedsCompilation:	yes
ByteCompile:	true
BugReports:	https://github.com/spatstat/spatstat.model/issues
Packaged:	2025-07-23 01:02:41 UTC; adrian
Author:	Adrian Baddeley [aut, cre, cph], Rolf Turner [aut, cph], Ege Rubak [aut, cph], Kasper Klitgaard Berthelsen [ctb], Achmad Choiruddin [ctb, cph], Jean-Francois Coeurjolly [ctb], Ottmar Cronie [ctb], Tilman Davies [ctb], Julian Gilbey [ctb], Yongtao Guan [ctb], Ute Hahn [ctb], Martin Hazelton [ctb], Kassel Hingee [ctb], Abdollah Jalilian [ctb], Frederic Lavancier [ctb], Marie-Colette van Lieshout [ctb], Bethany Macdonald [ctb], Greg McSwiggan [ctb], Tuomas Rajala [ctb], Suman Rakshit [ctb, cph], Dominic Schuhmacher [ctb], Rasmus Plenge Waagepetersen [ctb], Hangsheng Wang [ctb]
Repository:	CRAN
Date/Publication:	2025-07-23 05:10:33 UTC

The spatstat.model Package

Description

The spatstat.model package belongs to the spatstat family of packages. It contains the core functionality for parametric statistical modelling of spatial data.

Details

spatstat is a family of R packages for the statistical analysis of spatial data. Its main focus is the analysis of spatial patterns of points in two-dimensional space.

The original spatstat package has now been split into several sub-packages.

This sub-package spatstat.model contains all the main user-level functions that perform parametric statistical modelling of spatial data.

(The main exception is that functions for linear networks are in the separate sub-package spatstat.linnet.)

Structure of the spatstat family

The orginal spatstat package grew to be very large. It has now been divided into several sub-packages:

• spatstat.utils containing basic utilities

• spatstat.sparse containing linear algebra utilities

• spatstat.data containing datasets

• spatstat.univar containing functions for estimating probability distributions of random variables

• spatstat.geom containing geometrical objects and geometrical operations

• spatstat.explore containing the functionality for exploratory analysis and nonparametric modelling of spatial data

• spatstat.model containing the main functionality for parametric modelling, analysis and inference for spatial data

• spatstat.linnet containing functions for spatial data on a linear network

• spatstat, which simply loads the other sub-packages listed above, and provides documentation.

When you install spatstat, these sub-packages are also installed. Then if you load the spatstat package by typing library(spatstat), the other sub-packages listed above will automatically be loaded or imported.

For an overview of all the functions available in the sub-packages of spatstat, see the help file for "spatstat-package" in the spatstat package.

Additionally there are several extension packages:

• spatstat.gui for interactive graphics

• spatstat.local for local likelihood (including geographically weighted regression)

• spatstat.Knet for additional, computationally efficient code for linear networks

• spatstat.sphere (under development) for spatial data on a sphere, including spatial data on the earth's surface

The extension packages must be installed separately and loaded explicitly if needed. They also have separate documentation.

Overview of Functionality in spatstat.model

The spatstat family of packages is designed to support a complete statistical analysis of spatial data. It supports

• creation, manipulation and plotting of point patterns;

• exploratory data analysis;

• spatial random sampling;

• simulation of point process models;

• parametric model-fitting;

• non-parametric smoothing and regression;

• formal inference (hypothesis tests, confidence intervals);

• model diagnostics.

For an overview, see the help file for "spatstat-package" in the spatstat package.

Following is a list of the functionality provided in the spatstat.model package only.

To simulate a random point pattern:

Functions for generating random point patterns are now contained in the spatstat.random package.

Exploratory analysis

Exploratory graphics, smoothing, and exploratory analysis of spatial data are now provided in the spatstat.explore package.

Model fitting (Cox and cluster models)

Cluster process models (with homogeneous or inhomogeneous intensity) and Cox processes can be fitted by the function kppm. Its result is an object of class "kppm". The fitted model can be printed, plotted, predicted, simulated and updated.

For model selection, you can also use the generic functions step, drop1 and AIC on fitted point process models. For variable selection, see sdr.

The theoretical models can also be simulated, for any choice of parameter values, using rThomas, rMatClust, rCauchy, rVarGamma, and rLGCP.

Lower-level fitting functions include:

Model fitting (Poisson and Gibbs models)

Poisson point processes are the simplest models for point patterns. A Poisson model assumes that the points are stochastically independent. It may allow the points to have a non-uniform spatial density. The special case of a Poisson process with a uniform spatial density is often called Complete Spatial Randomness.

Poisson point processes are included in the more general class of Gibbs point process models. In a Gibbs model, there is interaction or dependence between points. Many different types of interaction can be specified.

For a detailed explanation of how to fit Poisson or Gibbs point process models to point pattern data using spatstat, see Baddeley and Turner (2005b) or Baddeley (2008).

To fit a Poisson or Gibbs point process model:

Model fitting in spatstat is performed mainly by the function ppm. Its result is an object of class "ppm".

Here are some examples, where X is a point pattern (class "ppp"):

It is also possible to fit models that depend on other covariates.

Manipulating the fitted model:

For model selection, you can also use the generic functions step, drop1 and AIC on fitted point process models. For variable selection, see sdr.

See spatstat.options to control plotting of fitted model.

To specify a point process model:

The first order “trend” of the model is determined by an R language formula. The formula specifies the form of the logarithm of the trend.

The higher order (“interaction”) components are described by an object of class "interact". Such objects are created by:

Note that it is also possible to combine several such interactions using Hybrid.

Simulation and goodness-of-fit for fitted models:

Model fitting (determinantal point process models)

Code for fitting determinantal point process models has recently been added to spatstat.

For information, see the help file for dppm.

Model fitting (spatial logistic regression)

Pixel-based spatial logistic regression is an alternative technique for analysing spatial point patterns that is widely used in Geographical Information Systems. It is approximately equivalent to fitting a Poisson point process model.

In pixel-based logistic regression, the spatial domain is divided into small pixels, the presence or absence of a data point in each pixel is recorded, and logistic regression is used to model the presence/absence indicators as a function of any covariates.

Facilities for performing spatial logistic regression are provided in spatstat for comparison purposes.

Fitting a spatial logistic regression

Spatial logistic regression is performed by the function slrm. Its result is an object of class "slrm". There are many methods for this class, including methods for print, fitted, predict, simulate, anova, coef, logLik, terms, update, formula and vcov.

For example, if X is a point pattern (class "ppp"):

Manipulating a fitted spatial logistic regression

There are many other undocumented methods for this class, including methods for print, update, formula and terms. Stepwise model selection is possible using step or stepAIC. For variable selection, see sdr.

Simulation

There are many ways to generate a random point pattern, line segment pattern, pixel image or tessellation in spatstat.

Random point patterns: Functions for random generation are now contained in the spatstat.random package.

See also varblock for estimating the variance of a summary statistic by block resampling, and lohboot for another bootstrap technique.

Fitted point process models:

If you have fitted a point process model to a point pattern dataset, the fitted model can be simulated.

Cluster process models are fitted by the function kppm yielding an object of class "kppm". To generate one or more simulated realisations of this fitted model, use simulate.kppm.

Gibbs point process models are fitted by the function ppm yielding an object of class "ppm". To generate a simulated realisation of this fitted model, use rmh.ppm. To generate one or more simulated realisations of the fitted model, use simulate.ppm.

Other random patterns: Functions for random generation are now contained in the spatstat.random package.

Simulation-based inference

Simulation-based inference including simulation envelopes and hypothesis tests is now supported by the package spatstat.explore.

Sensitivity diagnostics:

Classical measures of model sensitivity such as leverage and influence have been adapted to point process models.

Diagnostics for covariate effect:

Classical diagnostics for covariate effects have been adapted to point process models.

Residual diagnostics:

Residuals for a fitted point process model, and diagnostic plots based on the residuals, were introduced in Baddeley et al (2005) and Baddeley, Rubak and Moller (2011).

Type demo(diagnose) for a demonstration of the diagnostics features.

Resampling and randomisation procedures

You can build your own tests based on randomisation and resampling using the following capabilities:

Licence

This library and its documentation are usable under the terms of the "GNU General Public License", a copy of which is distributed with the package.

Acknowledgements

Kasper Klitgaard Berthelsen, Ottmar Cronie, Tilman Davies, Julian Gilbey, Yongtao Guan, Ute Hahn, Kassel Hingee, Abdollah Jalilian, Marie-Colette van Lieshout, Greg McSwiggan, Tuomas Rajala, Suman Rakshit, Dominic Schuhmacher, Rasmus Waagepetersen and Hangsheng Wang made substantial contributions of code.

For comments, corrections, bug alerts and suggestions, we thank Monsuru Adepeju, Corey Anderson, Ang Qi Wei, Ryan Arellano, Jens Astrom, Robert Aue, Marcel Austenfeld, Sandro Azaele, Malissa Baddeley, Guy Bayegnak, Colin Beale, Melanie Bell, Thomas Bendtsen, Ricardo Bernhardt, Andrew Bevan, Brad Biggerstaff, Anders Bilgrau, Leanne Bischof, Christophe Biscio, Roger Bivand, Jose M. Blanco Moreno, Florent Bonneu, Jordan Brown, Ian Buller, Julian Burgos, Simon Byers, Ya-Mei Chang, Jianbao Chen, Igor Chernayavsky, Y.C. Chin, Bjarke Christensen, Lucia Cobo Sanchez, Jean-Francois Coeurjolly, Kim Colyvas, Hadrien Commenges, Rochelle Constantine, Robin Corria Ainslie, Richard Cotton, Marcelino de la Cruz, Peter Dalgaard, Mario D'Antuono, Sourav Das, Peter Diggle, Patrick Donnelly, Ian Dryden, Stephen Eglen, Ahmed El-Gabbas, Belarmain Fandohan, Olivier Flores, David Ford, Peter Forbes, Shane Frank, Janet Franklin, Funwi-Gabga Neba, Oscar Garcia, Agnes Gault, Jonas Geldmann, Marc Genton, Shaaban Ghalandarayeshi, Jason Goldstick, Pavel Grabarnik, C. Graf, Ute Hahn, Andrew Hardegen, Martin Bogsted Hansen, Martin Hazelton, Juha Heikkinen, Mandy Hering, Markus Herrmann, Maximilian Hesselbarth, Paul Hewson, Hamidreza Heydarian, Kurt Hornik, Philipp Hunziker, Jack Hywood, Ross Ihaka, Cenk Icos, Aruna Jammalamadaka, Robert John-Chandran, Devin Johnson, Mahdieh Khanmohammadi, Bob Klaver, Lily Kozmian-Ledward, Peter Kovesi, Mike Kuhn, Jeff Laake, Robert Lamb, Frederic Lavancier, Tom Lawrence, Tomas Lazauskas, Jonathan Lee, George Leser, Angela Li, Li Haitao, George Limitsios, Andrew Lister, Nestor Luambua, Ben Madin, Martin Maechler, Kiran Marchikanti, Jeff Marcus, Robert Mark, Peter McCullagh, Monia Mahling, Jorge Mateu Mahiques, Ulf Mehlig, Frederico Mestre, Sebastian Wastl Meyer, Mi Xiangcheng, Lore De Middeleer, Robin Milne, Enrique Miranda, Jesper Moller, Annie Mollie, Ines Moncada, Mehdi Moradi, Virginia Morera Pujol, Erika Mudrak, Gopalan Nair, Nader Najari, Nicoletta Nava, Linda Stougaard Nielsen, Felipe Nunes, Jens Randel Nyengaard, Jens Oehlschlaegel, Thierry Onkelinx, Sean O'Riordan, Evgeni Parilov, Jeff Picka, Nicolas Picard, Tim Pollington, Mike Porter, Sergiy Protsiv, Adrian Raftery, Ben Ramage, Pablo Ramon, Xavier Raynaud, Nicholas Read, Matt Reiter, Ian Renner, Tom Richardson, Brian Ripley, Ted Rosenbaum, Barry Rowlingson, Jason Rudokas, Tyler Rudolph, John Rudge, Christopher Ryan, Farzaneh Safavimanesh, Aila Sarkka, Cody Schank, Katja Schladitz, Sebastian Schutte, Bryan Scott, Olivia Semboli, Francois Semecurbe, Vadim Shcherbakov, Shen Guochun, Shi Peijian, Harold-Jeffrey Ship, Tammy L Silva, Ida-Maria Sintorn, Yong Song, Malte Spiess, Mark Stevenson, Kaspar Stucki, Jan Sulavik, Michael Sumner, P. Surovy, Ben Taylor, Thordis Linda Thorarinsdottir, Leigh Torres, Berwin Turlach, Torben Tvedebrink, Kevin Ummer, Medha Uppala, Andrew van Burgel, Tobias Verbeke, Mikko Vihtakari, Alexendre Villers, Fabrice Vinatier, Maximilian Vogtland, Sasha Voss, Sven Wagner, Hao Wang, H. Wendrock, Jan Wild, Carl G. Witthoft, Selene Wong, Maxime Woringer, Luke Yates, Mike Zamboni and Achim Zeileis.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

The Area Interaction Point Process Model

Description

Creates an instance of the Area Interaction point process model (Widom-Rowlinson penetrable spheres model) which can then be fitted to point pattern data.

Usage

  AreaInter(r)

Arguments

Details

This function defines the interpoint interaction structure of a point process called the Widom-Rowlinson penetrable sphere model or area-interaction process. It can be used to fit this model to point pattern data.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the area interaction structure is yielded by the function AreaInter(). See the examples below.

In standard form, the area-interaction process (Widom and Rowlinson, 1970; Baddeley and Van Lieshout, 1995) with disc radius r, intensity parameter \kappa and interaction parameter \gamma is a point process with probability density

f(x_1,\ldots,x_n) = \alpha \kappa^{n(x)} \gamma^{-A(x)}

for a point pattern x, where x_1,\ldots,x_n represent the points of the pattern, n(x) is the number of points in the pattern, and A(x) is the area of the region formed by the union of discs of radius r centred at the points x_1,\ldots,x_n. Here \alpha is a normalising constant.

The interaction parameter \gamma can be any positive number. If \gamma = 1 then the model reduces to a Poisson process with intensity \kappa. If \gamma < 1 then the process is regular, while if \gamma > 1 the process is clustered. Thus, an area interaction process can be used to model either clustered or regular point patterns. Two points interact if the distance between them is less than 2r.

The standard form of the model, shown above, is a little complicated to interpret in practical applications. For example, each isolated point of the pattern x contributes a factor \kappa \gamma^{-\pi r^2} to the probability density.

In spatstat, the model is parametrised in a different form, which is easier to interpret. In canonical scale-free form, the probability density is rewritten as

f(x_1,\ldots,x_n) = \alpha \beta^{n(x)} \eta^{-C(x)}

where \beta is the new intensity parameter, \eta is the new interaction parameter, and C(x) = B(x) - n(x) is the interaction potential. Here

B(x) = \frac{A(x)}{\pi r^2}

is the normalised area (so that the discs have unit area). In this formulation, each isolated point of the pattern contributes a factor \beta to the probability density (so the first order trend is \beta). The quantity C(x) is a true interaction potential, in the sense that C(x) = 0 if the point pattern x does not contain any points that lie close together (closer than 2r units apart).

When a new point u is added to an existing point pattern x, the rescaled potential -C(x) increases by a value between 0 and 1. The increase is zero if u is not close to any point of x. The increase is 1 if the disc of radius r centred at u is completely contained in the union of discs of radius r centred at the data points x_i. Thus, the increase in potential is a measure of how close the new point u is to the existing pattern x. Addition of the point u contributes a factor \beta \eta^\delta to the probability density, where \delta is the increase in potential.

The old parameters \kappa,\gamma of the standard form are related to the new parameters \beta,\eta of the canonical scale-free form, by

\beta = \kappa \gamma^{-\pi r^2} = \kappa /\eta

and

\eta = \gamma^{\pi r^2}

provided \gamma and \kappa are positive and finite.

In the canonical scale-free form, the parameter \eta can take any nonnegative value. The value \eta = 1 again corresponds to a Poisson process, with intensity \beta. If \eta < 1 then the process is regular, while if \eta > 1 the process is clustered. The value \eta = 0 corresponds to a hard core process with hard core radius r (interaction distance 2r).

The nonstationary area interaction process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location, rather than a constant beta.

Note the only argument of AreaInter() is the disc radius r. When r is fixed, the model becomes an exponential family. The canonical parameters \log(\beta) and \log(\eta) are estimated by ppm(), not fixed in AreaInter().

Value

An object of class "interact" describing the interpoint interaction structure of the area-interaction process with disc radius r.

Warnings

The interaction distance of this process is equal to 2 * r. Two discs of radius r overlap if their centres are closer than 2 * r units apart.

The estimate of the interaction parameter \eta is unreliable if the interaction radius r is too small or too large. In these situations the model is approximately Poisson so that \eta is unidentifiable. As a rule of thumb, one can inspect the empty space function of the data, computed by Fest. The value F(r) of the empty space function at the interaction radius r should be between 0.2 and 0.8.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Baddeley, A.J. and Van Lieshout, M.N.M. (1995). Area-interaction point processes. Annals of the Institute of Statistical Mathematics 47 (1995) 601–619.

Widom, B. and Rowlinson, J.S. (1970). New model for the study of liquid-vapor phase transitions. The Journal of Chemical Physics 52 (1970) 1670–1684.

See Also

ppm, pairwise.family, ppm.object

ragsAreaInter and rmh for simulation of area-interaction models.

Examples

   

   # prints a sensible description of itself
   AreaInter(r=0.1)

   # Note the reach is twice the radius
   reach(AreaInter(r=1))

   # Fit the stationary area interaction process to Swedish Pines data
   ppm(swedishpines ~1, AreaInter(r=7))

   # Fit the stationary area interaction process to `cells'
   ppm(cells ~1, AreaInter(r=0.06))
   # eta=0 indicates hard core process.

   # Fit a nonstationary area interaction with log-cubic polynomial trend
   
     ppm(swedishpines ~polynom(x/10,y/10,3), AreaInter(r=7))
   

   

Hybrid Geyer Point Process Model

Description

Creates an instance of the Baddeley-Geyer point process model, defined as a hybrid of several Geyer interactions. The model can then be fitted to point pattern data.

Usage

  BadGey(r, sat)

Arguments

Details

This is Baddeley's generalisation of the Geyer saturation point process model, described in Geyer, to a process with multiple interaction distances.

The BadGey point process with interaction radii r_1,\ldots,r_k, saturation thresholds s_1,\ldots,s_k, intensity parameter \beta and interaction parameters \gamma_1,\ldots,gamma_k, is the point process in which each point x_i in the pattern X contributes a factor

\beta \gamma_1^{v_1(x_i, X)} \ldots gamma_k^{v_k(x_i,X)}

to the probability density of the point pattern, where

v_j(x_i, X) = \min( s_j, t_j(x_i,X) )

where t_j(x_i, X) denotes the number of points in the pattern X which lie within a distance r_j from the point x_i.

BadGey is used to fit this model to data. The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the piecewise constant Saturated pairwise interaction is yielded by the function BadGey(). See the examples below.

The argument r specifies the vector of interaction distances. The entries of r must be strictly increasing, positive numbers.

The argument sat specifies the vector of saturation parameters that are applied to the point counts t_j(x_i, X). It should be a vector of the same length as r, and its entries should be nonnegative numbers. Thus sat[1] is applied to the count of points within a distance r[1], and sat[2] to the count of points within a distance r[2], etc. Alternatively sat may be a single number, and this saturation value will be applied to every count.

Infinite values of the saturation parameters are also permitted; in this case v_j(x_i,X) = t_j(x_i,X) and there is effectively no ‘saturation’ for the distance range in question. If all the saturation parameters are set to Inf then the model is effectively a pairwise interaction process, equivalent to PairPiece (however the interaction parameters \gamma obtained from BadGey have a complicated relationship to the interaction parameters \gamma obtained from PairPiece).

If r is a single number, this model is virtually equivalent to the Geyer process, see Geyer.

Value

An object of class "interact" describing the interpoint interaction structure of a point process.

Hybrids

A ‘hybrid’ interaction is one which is built by combining several different interactions (Baddeley et al, 2013). The BadGey interaction can be described as a hybrid of several Geyer interactions.

The Hybrid command can be used to build hybrids of any interactions. If the Hybrid operator is applied to several Geyer models, the result is equivalent to a BadGey model. This can be useful for incremental model selection.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net in collaboration with Hao Wang and Jeff Picka

References

Baddeley, A., Turner, R., Mateu, J. and Bevan, A. (2013) Hybrids of Gibbs point process models and their implementation. Journal of Statistical Software 55:11, 1–43. DOI: 10.18637/jss.v055.i11

See Also

ppm, pairsat.family, Geyer, PairPiece, SatPiece, Hybrid

Examples

   BadGey(c(0.1,0.2), c(1,1))
   # prints a sensible description of itself
   BadGey(c(0.1,0.2), 1)

   # fit a stationary Baddeley-Geyer model
   ppm(cells ~1, BadGey(c(0.07, 0.1, 0.13), 2))

   # nonstationary process with log-cubic polynomial trend
   
     ppm(cells ~polynom(x,y,3), BadGey(c(0.07, 0.1, 0.13), 2))
   

The Connected Component Process Model

Description

Creates an instance of the Connected Component point process model which can then be fitted to point pattern data.

Usage

  Concom(r)

Arguments

Details

This function defines the interpoint interaction structure of a point process called the connected component process. It can be used to fit this model to point pattern data.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the connected component interaction is yielded by the function Concom(). See the examples below.

In standard form, the connected component process (Baddeley and Moller, 1989) with disc radius r, intensity parameter \kappa and interaction parameter \gamma is a point process with probability density

f(x_1,\ldots,x_n) = \alpha \kappa^{n(x)} \gamma^{-C(x)}

for a point pattern x, where x_1,\ldots,x_n represent the points of the pattern, n(x) is the number of points in the pattern, and C(x) is defined below. Here \alpha is a normalising constant.

To define the term C(x), suppose that we construct a planar graph by drawing an edge between each pair of points x_i,x_j which are less than r units apart. Two points belong to the same connected component of this graph if they are joined by a path in the graph. Then C(x) is the number of connected components of the graph.

The interaction parameter \gamma can be any positive number. If \gamma = 1 then the model reduces to a Poisson process with intensity \kappa. If \gamma < 1 then the process is regular, while if \gamma > 1 the process is clustered. Thus, a connected-component interaction process can be used to model either clustered or regular point patterns.

In spatstat, the model is parametrised in a different form, which is easier to interpret. In canonical form, the probability density is rewritten as

f(x_1,\ldots,x_n) = \alpha \beta^{n(x)} \gamma^{-U(x)}

where \beta is the new intensity parameter and U(x) = C(x) - n(x) is the interaction potential. In this formulation, each isolated point of the pattern contributes a factor \beta to the probability density (so the first order trend is \beta). The quantity U(x) is a true interaction potential, in the sense that U(x) = 0 if the point pattern x does not contain any points that lie close together.

When a new point u is added to an existing point pattern x, the rescaled potential -U(x) increases by zero or a positive integer. The increase is zero if u is not close to any point of x. The increase is a positive integer k if there are k different connected components of x that lie close to u. Addition of the point u contributes a factor \beta \eta^\delta to the probability density, where \delta is the increase in potential.

If desired, the original parameter \kappa can be recovered from the canonical parameter by \kappa = \beta\gamma.

The nonstationary connected component process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location, rather than a constant beta.

Note the only argument of Concom() is the threshold distance r. When r is fixed, the model becomes an exponential family. The canonical parameters \log(\beta) and \log(\gamma) are estimated by ppm(), not fixed in Concom().

Value

An object of class "interact" describing the interpoint interaction structure of the connected component process with disc radius r.

Edge correction

The interaction distance of this process is infinite. There are no well-established procedures for edge correction for fitting such models, and accordingly the model-fitting function ppm will give an error message saying that the user must specify an edge correction. A reasonable solution is to use the border correction at the same distance r, as shown in the Examples.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

References

Baddeley, A.J. and Moller, J. (1989) Nearest-neighbour Markov point processes and random sets. International Statistical Review 57, 89–121.

See Also

ppm, pairwise.family, ppm.object

Examples

   # prints a sensible description of itself
   Concom(r=0.1)

   # Fit the stationary connected component process to redwood data
   ppm(redwood ~1, Concom(r=0.07), rbord=0.07)

   # Fit the stationary connected component process to `cells' data
   ppm(cells ~1, Concom(r=0.06), rbord=0.06)
   # eta=0 indicates hard core process.

   # Fit a nonstationary connected component model
   # with log-cubic polynomial trend
   
     ppm(swedishpines ~polynom(x/10,y/10,3), Concom(r=7), rbord=7)
   

Diggle-Gates-Stibbard Point Process Model

Description

Creates an instance of the Diggle-Gates-Stibbard point process model which can then be fitted to point pattern data.

Usage

  DiggleGatesStibbard(rho)

Arguments

Details

Diggle, Gates and Stibbard (1987) proposed a pairwise interaction point process in which each pair of points separated by a distance d contributes a factor e(d) to the probability density, where

e(d) = \sin^2\left(\frac{\pi d}{2\rho}\right)

for d < \rho, and e(d) is equal to 1 for d \ge \rho.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the Diggle-Gates-Stibbard pairwise interaction is yielded by the function DiggleGatesStibbard(). See the examples below.

Note that this model does not have any regular parameters (as explained in the section on Interaction Parameters in the help file for ppm). The parameter \rho is not estimated by ppm.

Value

An object of class "interact" describing the interpoint interaction structure of the Diggle-Gates-Stibbard process with interaction range rho.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

References

Baddeley, A. and Turner, R. (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42, 283–322.

Ripley, B.D. (1981) Spatial statistics. John Wiley and Sons.

Diggle, P.J., Gates, D.J., and Stibbard, A. (1987) A nonparametric estimator for pairwise-interaction point processes. Biometrika 74, 763 – 770. Scandinavian Journal of Statistics 21, 359–373.

See Also

ppm, pairwise.family, DiggleGratton, rDGS, ppm.object

Examples

   DiggleGatesStibbard(0.02)
   # prints a sensible description of itself

   ppm(cells ~1, DiggleGatesStibbard(0.05))
   # fit the stationary D-G-S process to `cells'

   
     ppm(cells ~ polynom(x,y,3), DiggleGatesStibbard(0.05))
     # fit a nonstationary D-G-S process
     # with log-cubic polynomial trend
   

Diggle-Gratton model

Description

Creates an instance of the Diggle-Gratton pairwise interaction point process model, which can then be fitted to point pattern data.

Usage

  DiggleGratton(delta=NA, rho)

Arguments

Details

Diggle and Gratton (1984, pages 208-210) introduced the pairwise interaction point process with pair potential h(t) of the form

h(t) = \left( \frac{t-\delta}{\rho-\delta} \right)^\kappa \quad\quad \mbox{ if } \delta \le t \le \rho

with h(t) = 0 for t < \delta and h(t) = 1 for t > \rho. Here \delta, \rho and \kappa are parameters.

Note that we use the symbol \kappa where Diggle and Gratton (1984) and Diggle, Gates and Stibbard (1987) use \beta, since in spatstat we reserve the symbol \beta for an intensity parameter.

The parameters must all be nonnegative, and must satisfy \delta \le \rho.

The potential is inhibitory, i.e.\ this model is only appropriate for regular point patterns. The strength of inhibition increases with \kappa. For \kappa=0 the model is a hard core process with hard core radius \delta. For \kappa=\infty the model is a hard core process with hard core radius \rho.

The irregular parameters \delta, \rho must be given in the call to DiggleGratton, while the regular parameter \kappa will be estimated.

If the lower threshold delta is missing or NA, it will be estimated from the data when ppm is called. The estimated value of delta is the minimum nearest neighbour distance multiplied by n/(n+1), where n is the number of data points.

Value

An object of class "interact" describing the interpoint interaction structure of a point process.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

References

Diggle, P.J., Gates, D.J. and Stibbard, A. (1987) A nonparametric estimator for pairwise-interaction point processes. Biometrika 74, 763 – 770.

Diggle, P.J. and Gratton, R.J. (1984) Monte Carlo methods of inference for implicit statistical models. Journal of the Royal Statistical Society, series B 46, 193 – 212.

See Also

ppm, ppm.object, Pairwise

Examples

   ppm(cells ~1, DiggleGratton(0.05, 0.1))

Extract Subset of Influence Object

Description

Extract a subset of an influence object, or extract the influence values at specified locations.

Usage

  ## S3 method for class 'influence.ppm'
x[i, ...]

Arguments

Details

An object of class "influence.ppm" contains the values of the likelihood influence for a point process model, computed by influence.ppm. This is effectively a marked point pattern obtained by marking each of the original data points with its likelihood influence.

This function extracts a designated subset of the influence values, either as another influence object, or as a vector of numeric values.

The function [.influence.ppm is a method for [ for the class "influence.ppm". The argument i should be an index applicable to a point pattern. It may be either a spatial window (object of class "owin") or a sequence index. The result will be another influence object (of class influence.ppm).

To extract the influence values as a numeric vector, use marks(as.ppp(x)).

Value

Another object of class "influence.ppm".

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

See Also

influence.ppm.

Examples

   fit <- ppm(cells, ~x)
   infl <- influence(fit)
   b <- owin(c(0.1, 0.3), c(0.2, 0.4))
   infl[b]
   infl[1:5]
   marks(as.ppp(infl))[1:3]

Extract Subset of Leverage Object

Description

Extract a subset of a leverage map, or extract the leverage values at specified locations.

Usage

  ## S3 method for class 'leverage.ppm'
x[i, ..., update=TRUE]

Arguments

Details

An object of class "leverage.ppm" contains the values of the leverage function for a point process model, computed by leverage.ppm.

This function extracts a designated subset of the leverage values, either as another leverage object, or as a vector of numeric values.

The function [.leverage.ppm is a method for [ for the class "leverage.ppm". The argument i should be either

• a spatial window (object of class "owin") determining a region where the leverage map is required. The result will typically be another leverage map (object of class leverage.ppm).

• a spatial point pattern (object of class "ppp") specifying locations at which the leverage values are required. The result will be a numeric vector.

The subset operator for images, [.im, is applied to the leverage map. If this yields a pixel image, then the result of [.leverage.ppm is another leverage object. Otherwise, a vector containing the numeric values of leverage is returned.

Value

Another object of class "leverage.ppm", or a vector of numeric values of leverage.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

See Also

leverage.ppm.

Examples

   fit <- ppm(cells ~x)
   lev <- leverage(fit)
   b <- owin(c(0.1, 0.3), c(0.2, 0.4))
   lev[b]
   lev[cells]

Extract Subset of Signed or Vector Measure

Description

Extract a subset of a signed measure or vector-valued measure.

Usage

## S3 method for class 'msr'
x[i, j, ...]

Arguments

Details

This operator extracts a subset of the data which determines the signed measure or vector-valued measure x. The result is another measure.

Value

An object of class "msr".

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

See Also

msr

Examples

   X <- rpoispp(function(x,y) { exp(3+3*x) })
   fit <- ppm(X ~x+y)
   rp <- residuals(fit, type="pearson")
   rs <- residuals(fit, type="score")

   rp[square(0.5)]
   rs[ , 2:3]

The Fiksel Interaction

Description

Creates an instance of Fiksel's double exponential pairwise interaction point process model, which can then be fitted to point pattern data.

Usage

  Fiksel(r, hc=NA, kappa)

Arguments

Details

Fiksel (1984) introduced a pairwise interaction point process with the following interaction function c. For two points u and v separated by a distance d=||u-v||, the interaction c(u,v) is equal to 0 if d < h, equal to 1 if d > r, and equal to

\exp(a \exp(-\kappa d))

if h \le d \le r, where h,r,\kappa,a are parameters.

A graph of this interaction function is shown in the Examples. The interpretation of the parameters is as follows.

• h is the hard core distance: distinct points are not permitted to come closer than a distance h apart.

• r is the interaction range: points further than this distance do not interact.

• \kappa is the rate or slope parameter, controlling the decay of the interaction as distance increases.

• a is the interaction strength parameter, controlling the strength and type of interaction. If a is zero, the process is Poisson. If a is positive, the process is clustered. If a is negative, the process is inhibited (regular).

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the Fiksel pairwise interaction is yielded by the function Fiksel(). See the examples below.

The parameters h, r and \kappa must be fixed and given in the call to Fiksel, while the canonical parameter a is estimated by ppm().

To estimate h, r and\kappa it is possible to use profilepl. The maximum likelihood estimator ofh is the minimum interpoint distance.

If the hard core distance argument hc is missing or NA, it will be estimated from the data when ppm is called. The estimated value of hc is the minimum nearest neighbour distance multiplied by n/(n+1), where n is the number of data points.

See also Stoyan, Kendall and Mecke (1987) page 161.

Value

An object of class "interact" describing the interpoint interaction structure of the Fiksel process with interaction radius r, hard core distance hc and rate parameter kappa.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Baddeley, A. and Turner, R. (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42, 283–322.

Fiksel, T. (1984) Estimation of parameterized pair potentials of marked and non-marked Gibbsian point processes. Electronische Informationsverabeitung und Kybernetika 20, 270–278.

Stoyan, D, Kendall, W.S. and Mecke, J. (1987) Stochastic geometry and its applications. Wiley.

See Also

ppm, pairwise.family, ppm.object, StraussHard

Examples

   Fiksel(r=1,hc=0.02, kappa=2)
   # prints a sensible description of itself

   X <- unmark(spruces)

   fit <- ppm(X ~ 1, Fiksel(r=3.5, kappa=1))
   plot(fitin(fit))

Model Compensator of Nearest Neighbour Function

Description

Given a point process model fitted to a point pattern dataset, this function computes the compensator of the nearest neighbour distance distribution function G based on the fitted model (as well as the usual nonparametric estimates of G based on the data alone). Comparison between the nonparametric and model-compensated G functions serves as a diagnostic for the model.

Usage

Gcom(object, r = NULL, breaks = NULL, ...,
     correction = c("border", "Hanisch"),
     conditional = !is.poisson(object),
     restrict=FALSE,
     model=NULL,
     trend = ~1, interaction = Poisson(),
     rbord = reach(interaction),
     ppmcorrection="border",
     truecoef = NULL, hi.res = NULL)

Arguments

Details

This command provides a diagnostic for the goodness-of-fit of a point process model fitted to a point pattern dataset. It computes different estimates of the nearest neighbour distance distribution function G of the dataset, which should be approximately equal if the model is a good fit to the data.

The first argument, object, is usually a fitted point process model (object of class "ppm"), obtained from the model-fitting function ppm.

For convenience, object can also be a point pattern (object of class "ppp"). In that case, a point process model will be fitted to it, by calling ppm using the arguments trend (for the first order trend), interaction (for the interpoint interaction) and rbord (for the erosion distance in the border correction for the pseudolikelihood). See ppm for details of these arguments.

The algorithm first extracts the original point pattern dataset (to which the model was fitted) and computes the standard nonparametric estimates of the G function. It then also computes the model-compensated G function. The different functions are returned as columns in a data frame (of class "fv"). The interpretation of the columns is as follows (ignoring edge corrections):

bord:

the nonparametric border-correction estimate of G(r),

\hat G(r) = \frac{\sum_i I\{ d_i \le r\} I\{ b_i > r \}}{\sum_i I\{ b_i > r\}}

where d_i is the distance from the i-th data point to its nearest neighbour, and b_i is the distance from the i-th data point to the boundary of the window W.

bcom:

the model compensator of the border-correction estimate

{\bf C}\, \hat G(r) = \frac{\int \lambda(u,x) I\{ b(u) > r\} I\{ d(u,x) \le r\}}{ 1 + \sum_i I\{ b_i > r\} }

where \lambda(u,x) denotes the conditional intensity of the model at the location u, and d(u,x) denotes the distance from u to the nearest point in x, while b(u) denotes the distance from u to the boundary of the windowW.

han:

the nonparametric Hanisch estimate of G(r)

\hat G(r) = \frac{D(r)}{D(\infty)}

where

D(r) = \sum_i \frac{ I\{x_i \in W_{\ominus d_i}\} I\{d_i \le r\} }{ \mbox{area}(W_{\ominus d_i}) }

in which W_{\ominus r} denotes the erosion of the window W by a distance r.

hcom:

the corresponding model-compensated function

{\bf C} \, G(r) = \int_W \frac{ \lambda(u,x) I(u \in W_{\ominus d(u)}) I(d(u) \le r) }{ \hat D(\infty) \mbox{area}(W_{\ominus d(u)}) + 1 }

where d(u) = d(u, x) is the (‘empty space’) distance from location u to the nearest point of x.

If the fitted model is a Poisson point process, then the formulae above are exactly what is computed. If the fitted model is not Poisson, the formulae above are modified slightly to handle edge effects.

The modification is determined by the arguments conditional and restrict. The value of conditional defaults to FALSE for Poisson models and TRUE for non-Poisson models. If conditional=FALSE then the formulae above are not modified. If conditional=TRUE, then the algorithm calculates the restriction estimator if restrict=TRUE, and calculates the reweighting estimator if restrict=FALSE. See Appendix E of Baddeley, Rubak and Moller (2011). See also spatstat.options('eroded.intensity'). Thus, by default, the reweighting estimator is computed for non-Poisson models.

The border-corrected and Hanisch-corrected estimates of G(r) are approximately unbiased estimates of the G-function, assuming the point process is stationary. The model-compensated functions are unbiased estimates of the mean value of the corresponding nonparametric estimate, assuming the model is true. Thus, if the model is a good fit, the mean value of the difference between the nonparametric and model-compensated estimates is approximately zero.

To compute the difference between the nonparametric and model-compensated functions, use Gres.

Value

A function value table (object of class "fv"), essentially a data frame of function values. There is a plot method for this class. See fv.object.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ege Rubak rubak@math.aau.dk and Jesper Moller.

References

Baddeley, A., Rubak, E. and Moller, J. (2011) Score, pseudo-score and residual diagnostics for spatial point process models. Statistical Science 26, 613–646.

See Also

Related functions: Gest, Gres.

Alternative functions: Kcom, psstA, psstG, psst.

Model fitting: ppm.

Examples

    fit0 <- ppm(cells ~1) # uniform Poisson
    G0 <- Gcom(fit0)
    G0
    plot(G0)
# uniform Poisson is clearly not correct

# Hanisch estimates only
    plot(Gcom(fit0), cbind(han, hcom) ~ r)

    fit1 <- ppm(cells, ~1, Strauss(0.08))
    plot(Gcom(fit1), cbind(han, hcom) ~ r)

# Try adjusting interaction distance

    fit2 <- update(fit1, Strauss(0.10))
    plot(Gcom(fit2), cbind(han, hcom) ~ r)

    G3 <- Gcom(cells, interaction=Strauss(0.12))
    plot(G3, cbind(han, hcom) ~ r)

Geyer's Saturation Point Process Model

Description

Creates an instance of Geyer's saturation point process model which can then be fitted to point pattern data.

Usage

  Geyer(r,sat)

Arguments

Details

Geyer (1999) introduced the “saturation process”, a modification of the Strauss process (see Strauss) in which the total contribution to the potential from each point (from its pairwise interaction with all other points) is trimmed to a maximum value s. The interaction structure of this model is implemented in the function Geyer().

The saturation point process with interaction radius r, saturation threshold s, and parameters \beta and \gamma, is the point process in which each point x_i in the pattern X contributes a factor

\beta \gamma^{\min(s, t(x_i, X))}

to the probability density of the point pattern, where t(x_i, X) denotes the number of ‘close neighbours’ of x_i in the pattern X. A close neighbour of x_i is a point x_j with j \neq i such that the distance between x_i and x_j is less than or equal to r.

If the saturation threshold s is set to infinity, this model reduces to the Strauss process (see Strauss) with interaction parameter \gamma^2. If s = 0, the model reduces to the Poisson point process. If s is a finite positive number, then the interaction parameter \gamma may take any positive value (unlike the case of the Strauss process), with values \gamma < 1 describing an ‘ordered’ or ‘inhibitive’ pattern, and values \gamma > 1 describing a ‘clustered’ or ‘attractive’ pattern.

The nonstationary saturation process is similar except that the value \beta is replaced by a function \beta(x_i) of location.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the saturation process interaction is yielded by Geyer(r, sat) where the arguments r and sat specify the Strauss interaction radius r and the saturation threshold s, respectively. See the examples below.

Note the only arguments are the interaction radius r and the saturation threshold sat. When r and sat are fixed, the model becomes an exponential family. The canonical parameters \log(\beta) and \log(\gamma) are estimated by ppm(), not fixed in Geyer().

Value

An object of class "interact" describing the interpoint interaction structure of Geyer's saturation point process with interaction radius r and saturation threshold sat.

Zero saturation

The value sat=0 is permitted by Geyer, but this is not very useful. For technical reasons, when ppm fits a Geyer model with sat=0, the default behaviour is to return an “invalid” fitted model in which the estimate of \gamma is NA. In order to get a Poisson process model returned when sat=0, you would need to set emend=TRUE in the call to ppm.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Geyer, C.J. (1999) Likelihood Inference for Spatial Point Processes. Chapter 3 in O.E. Barndorff-Nielsen, W.S. Kendall and M.N.M. Van Lieshout (eds) Stochastic Geometry: Likelihood and Computation, Chapman and Hall / CRC, Monographs on Statistics and Applied Probability, number 80. Pages 79–140.

See Also

ppm, pairwise.family, ppm.object, Strauss.

To make an interaction object like Geyer but having multiple interaction radii, see BadGey or Hybrid.

Examples

   ppm(cells, ~1, Geyer(r=0.07, sat=2))
   # fit the stationary saturation process to `cells'

Residual G Function

Description

Given a point process model fitted to a point pattern dataset, this function computes the residual G function, which serves as a diagnostic for goodness-of-fit of the model.

Usage

   Gres(object, ...)

Arguments

Details

This command provides a diagnostic for the goodness-of-fit of a point process model fitted to a point pattern dataset. It computes a residual version of the G function of the dataset, which should be approximately zero if the model is a good fit to the data.

In normal use, object is a fitted point process model or a point pattern. Then Gres first calls Gcom to compute both the nonparametric estimate of the G function and its model compensator. Then Gres computes the difference between them, which is the residual G-function.

Alternatively, object may be a function value table (object of class "fv") that was returned by a previous call to Gcom. Then Gres computes the residual from this object.

Value

A function value table (object of class "fv"), essentially a data frame of function values. There is a plot method for this class. See fv.object.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ege Rubak rubak@math.aau.dk and Jesper Moller.

References

Baddeley, A., Rubak, E. and Moller, J. (2011) Score, pseudo-score and residual diagnostics for spatial point process models. Statistical Science 26, 613–646.

See Also

Related functions: Gcom, Gest.

Alternative functions: Kres, psstA, psstG, psst.

Model-fitting: ppm.

Examples

    fit0 <- ppm(cells, ~1) # uniform Poisson
    G0 <- Gres(fit0)
    plot(G0)
# Hanisch correction estimate
    plot(G0, hres ~ r)
# uniform Poisson is clearly not correct

    fit1 <- ppm(cells, ~1, Strauss(0.08))
    plot(Gres(fit1), hres ~ r)
# fit looks approximately OK; try adjusting interaction distance

    plot(Gres(cells, interaction=Strauss(0.12)))

# How to make envelopes
    if(interactive()) {
      E <- envelope(fit1, Gres, model=fit1, nsim=39)
      plot(E)
    }
# For computational efficiency
    Gc <- Gcom(fit1)
    G1 <- Gres(Gc)

The Hard Core Point Process Model

Description

Creates an instance of the hard core point process model which can then be fitted to point pattern data.

Usage

  Hardcore(hc=NA)

Arguments

Details

A hard core process with hard core distance h and abundance parameter \beta is a pairwise interaction point process in which distinct points are not allowed to come closer than a distance h apart.

The probability density is zero if any pair of points is closer than h units apart, and otherwise equals

f(x_1,\ldots,x_n) = \alpha \beta^{n(x)}

where x_1,\ldots,x_n represent the points of the pattern, n(x) is the number of points in the pattern, and \alpha is the normalising constant.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the hard core process pairwise interaction is yielded by the function Hardcore(). See the examples below.

If the hard core distance argument hc is missing or NA, it will be estimated from the data when ppm is called. The estimated value of hc is the minimum nearest neighbour distance multiplied by n/(n+1), where n is the number of data points.

Value

An object of class "interact" describing the interpoint interaction structure of the hard core process with hard core distance hc.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Baddeley, A. and Turner, R. (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42, 283–322.

Ripley, B.D. (1981) Spatial statistics. John Wiley and Sons.

See Also

Strauss, StraussHard, MultiHard, ppm, pairwise.family, ppm.object

Examples

   Hardcore(0.02)
   # prints a sensible description of itself

   ppm(cells ~1, Hardcore(0.05))
   # fit the stationary hard core process to `cells'

   # estimate hard core radius from data
   ppm(cells ~1, Hardcore())
   
     # equivalent:
     ppm(cells ~1, Hardcore)
   

   
     # fit a nonstationary hard core process
     # with log-cubic polynomial trend
     ppm(cells ~ polynom(x,y,3), Hardcore(0.05))
   

The Hierarchical Hard Core Point Process Model

Description

Creates an instance of the hierarchical hard core point process model which can then be fitted to point pattern data.

Usage

  HierHard(hradii=NULL, types=NULL, archy=NULL)

Arguments

Details

This is a hierarchical point process model for a multitype point pattern (Hogmander and Sarkka, 1999; Grabarnik and Sarkka, 2009). It is appropriate for analysing multitype point pattern data in which the types are ordered so that the points of type j depend on the points of type 1,2,\ldots,j-1.

The hierarchical version of the (stationary) hard core process with m types, with hard core distances h_{ij} and parameters \beta_j, is a point process in which each point of type j contributes a factor \beta_j to the probability density of the point pattern. If any pair of points of types i and j lies closer than h_{ij} units apart, the configuration of points is impossible (probability density zero).

The nonstationary hierarchical hard core process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location and type, rather than a constant beta.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the hierarchical hard core process pairwise interaction is yielded by the function HierHard(). See the examples below.

The argument types need not be specified in normal use. It will be determined automatically from the point pattern data set to which the HierHard interaction is applied, when the user calls ppm. However, the user should be confident that the ordering of types in the dataset corresponds to the ordering of rows and columns in the matrix radii.

The argument archy can be used to specify a hierarchical ordering of the types. It can be either a vector of integers or a character vector matching the possible types. The default is the sequence 1,2, \ldots, m meaning that type j depends on types 1,2, \ldots, j-1.

The matrix iradii must be square, with entries which are either positive numbers, or zero or NA. A value of zero or NA indicates that no hard core interaction term should be included for this combination of types.

Note that only the hard core distances are specified in HierHard. The canonical parameters \log(\beta_j) are estimated by ppm(), not fixed in HierHard().

Value

An object of class "interact" describing the interpoint interaction structure of the hierarchical hard core process with hard core distances hradii[i,j].

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

, Rolf Turner rolfturner@posteo.net

and Ege Rubak rubak@math.aau.dk.

References

Grabarnik, P. and Sarkka, A. (2009) Modelling the spatial structure of forest stands by multivariate point processes with hierarchical interactions. Ecological Modelling 220, 1232–1240.

Hogmander, H. and Sarkka, A. (1999) Multitype spatial point patterns with hierarchical interactions. Biometrics 55, 1051–1058.

See Also

MultiHard for the corresponding symmetrical interaction.

HierStrauss, HierStraussHard.

Examples

   h <- matrix(c(4, NA, 10, 15), 2, 2)
   HierHard(h)
   # prints a sensible description of itself
   ppm(ants ~1, HierHard(h))
   # fit the stationary hierarchical hard core process to ants data

The Hierarchical Strauss Point Process Model

Description

Creates an instance of the hierarchical Strauss point process model which can then be fitted to point pattern data.

Usage

  HierStrauss(radii, types=NULL, archy=NULL)

Arguments

Details

This is a hierarchical point process model for a multitype point pattern (Hogmander and Sarkka, 1999; Grabarnik and Sarkka, 2009). It is appropriate for analysing multitype point pattern data in which the types are ordered so that the points of type j depend on the points of type 1,2,\ldots,j-1.

The hierarchical version of the (stationary) Strauss process with m types, with interaction radii r_{ij} and parameters \beta_j and \gamma_{ij} is a point process in which each point of type j contributes a factor \beta_j to the probability density of the point pattern, and a pair of points of types i and j closer than r_{ij} units apart contributes a factor \gamma_{ij} to the density provided i \le j.

The nonstationary hierarchical Strauss process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location and type, rather than a constant beta.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the hierarchical Strauss process pairwise interaction is yielded by the function HierStrauss(). See the examples below.

The argument types need not be specified in normal use. It will be determined automatically from the point pattern data set to which the HierStrauss interaction is applied, when the user calls ppm. However, the user should be confident that the ordering of types in the dataset corresponds to the ordering of rows and columns in the matrix radii.

The argument archy can be used to specify a hierarchical ordering of the types. It can be either a vector of integers or a character vector matching the possible types. The default is the sequence 1,2, \ldots, m meaning that type j depends on types 1,2, \ldots, j-1.

The matrix radii must be symmetric, with entries which are either positive numbers or NA. A value of NA indicates that no interaction term should be included for this combination of types.

Note that only the interaction radii are specified in HierStrauss. The canonical parameters \log(\beta_j) and \log(\gamma_{ij}) are estimated by ppm(), not fixed in HierStrauss().

Value

An object of class "interact" describing the interpoint interaction structure of the hierarchical Strauss process with interaction radii radii[i,j].

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

, Rolf Turner rolfturner@posteo.net

and Ege Rubak rubak@math.aau.dk.

References

Grabarnik, P. and Sarkka, A. (2009) Modelling the spatial structure of forest stands by multivariate point processes with hierarchical interactions. Ecological Modelling 220, 1232–1240.

Hogmander, H. and Sarkka, A. (1999) Multitype spatial point patterns with hierarchical interactions. Biometrics 55, 1051–1058.

See Also

MultiStrauss for the corresponding symmetrical interaction.

HierHard, HierStraussHard.

Examples

   r <- matrix(10 * c(3,4,4,3), nrow=2,ncol=2)
   HierStrauss(r)
   # prints a sensible description of itself
   ppm(ants ~1, HierStrauss(r, , c("Messor", "Cataglyphis")))
   # fit the stationary hierarchical Strauss process to ants data

The Hierarchical Strauss Hard Core Point Process Model

Description

Creates an instance of the hierarchical Strauss-hard core point process model which can then be fitted to point pattern data.

Usage

  HierStraussHard(iradii, hradii=NULL, types=NULL, archy=NULL)

Arguments

Details

This is a hierarchical point process model for a multitype point pattern (Hogmander and Sarkka, 1999; Grabarnik and Sarkka, 2009). It is appropriate for analysing multitype point pattern data in which the types are ordered so that the points of type j depend on the points of type 1,2,\ldots,j-1.

The hierarchical version of the (stationary) Strauss hard core process with m types, with interaction radii r_{ij}, hard core distances h_{ij} and parameters \beta_j and \gamma_{ij} is a point process in which each point of type j contributes a factor \beta_j to the probability density of the point pattern, and a pair of points of types i and j closer than r_{ij} units apart contributes a factor \gamma_{ij} to the density provided i \le j. If any pair of points of types i and j lies closer than h_{ij} units apart, the configuration of points is impossible (probability density zero).

The nonstationary hierarchical Strauss hard core process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location and type, rather than a constant beta.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the hierarchical Strauss hard core process pairwise interaction is yielded by the function HierStraussHard(). See the examples below.

The argument types need not be specified in normal use. It will be determined automatically from the point pattern data set to which the HierStraussHard interaction is applied, when the user calls ppm. However, the user should be confident that the ordering of types in the dataset corresponds to the ordering of rows and columns in the matrix radii.

The argument archy can be used to specify a hierarchical ordering of the types. It can be either a vector of integers or a character vector matching the possible types. The default is the sequence 1,2, \ldots, m meaning that type j depends on types 1,2, \ldots, j-1.

The matrices iradii and hradii must be square, with entries which are either positive numbers or zero or NA. A value of zero or NA indicates that no interaction term should be included for this combination of types.

Note that only the interaction radii and hard core distances are specified in HierStraussHard. The canonical parameters \log(\beta_j) and \log(\gamma_{ij}) are estimated by ppm(), not fixed in HierStraussHard().

Value

An object of class "interact" describing the interpoint interaction structure of the hierarchical Strauss-hard core process with interaction radii iradii[i,j] and hard core distances hradii[i,j].

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

, Rolf Turner rolfturner@posteo.net

and Ege Rubak rubak@math.aau.dk.

References

Grabarnik, P. and Sarkka, A. (2009) Modelling the spatial structure of forest stands by multivariate point processes with hierarchical interactions. Ecological Modelling 220, 1232–1240.

Hogmander, H. and Sarkka, A. (1999) Multitype spatial point patterns with hierarchical interactions. Biometrics 55, 1051–1058.

See Also

MultiStraussHard for the corresponding symmetrical interaction.

HierHard, HierStrauss.

Examples

   r <- matrix(c(30, NA, 40, 30), nrow=2,ncol=2)
   h <- matrix(c(4, NA, 10, 15), 2, 2)
   HierStraussHard(r, h)
   # prints a sensible description of itself
   ppm(ants ~1, HierStraussHard(r, h))
   # fit the stationary hierarchical Strauss-hard core process to ants data

Hybrid Interaction Point Process Model

Description

Creates an instance of a hybrid point process model which can then be fitted to point pattern data.

Usage

Hybrid(...)

Arguments

Details

A hybrid (Baddeley, Turner, Mateu and Bevan, 2013) is a point process model created by combining two or more point process models, or an interpoint interaction created by combining two or more interpoint interactions.

The hybrid of two point processes, with probability densities f(x) and g(x) respectively, is the point process with probability density

h(x) = c \, f(x) \, g(x)

where c is a normalising constant.

Equivalently, the hybrid of two point processes with conditional intensities \lambda(u,x) and \kappa(u,x) is the point process with conditional intensity

\phi(u,x) = \lambda(u,x) \, \kappa(u,x).

The hybrid of m > 3 point processes is defined in a similar way.

The function ppm, which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of a hybrid interaction is yielded by the function Hybrid().

The arguments ... will be interpreted as interpoint interactions (objects of class "interact") and the result will be the hybrid of these interactions. Each argument must either be an interpoint interaction (object of class "interact"), or a point process model (object of class "ppm") from which the interpoint interaction will be extracted.

The arguments ... may also be given in the form name=value. This is purely cosmetic: it can be used to attach simple mnemonic names to the component interactions, and makes the printed output from print.ppm neater.

Value

An object of class "interact" describing an interpoint interaction structure.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Baddeley, A., Turner, R., Mateu, J. and Bevan, A. (2013) Hybrids of Gibbs point process models and their implementation. Journal of Statistical Software 55:11, 1–43. DOI: 10.18637/jss.v055.i11

See Also

ppm

Examples

  Hybrid(Strauss(0.1), Geyer(0.2, 3))

  Hybrid(Ha=Hardcore(0.05), St=Strauss(0.1), Ge=Geyer(0.2, 3))

  fit <- ppm(redwood, ~1, Hybrid(A=Strauss(0.02), B=Geyer(0.1, 2)))
  fit

  ctr <- rmhcontrol(nrep=5e4, expand=1)
  plot(simulate(fit, control=ctr))

  # hybrid components can be models (including hybrid models)
  Hybrid(fit, S=Softcore(0.5))

  # plot.fii only works if every component is a pairwise interaction
  fit2 <- ppm(swedishpines, ~1, Hybrid(DG=DiggleGratton(2,10), S=Strauss(5)))
  plot(fitin(fit2))
  plot(fitin(fit2), separate=TRUE, mar.panel=rep(4,4))

Model Compensator of K Function

Description

Given a point process model fitted to a point pattern dataset, this function computes the compensator of the K function based on the fitted model (as well as the usual nonparametric estimates of K based on the data alone). Comparison between the nonparametric and model-compensated K functions serves as a diagnostic for the model.

Usage

Kcom(object, r = NULL, breaks = NULL, ...,
     correction = c("border", "isotropic", "translate"),
     conditional = !is.poisson(object),
     restrict = FALSE,
     model = NULL,
     trend = ~1, interaction = Poisson(), rbord = reach(interaction),
     compute.var = TRUE,
     truecoef = NULL, hi.res = NULL)

Arguments

Details

This command provides a diagnostic for the goodness-of-fit of a point process model fitted to a point pattern dataset. It computes an estimate of the K function of the dataset, together with a model compensator of the K function, which should be approximately equal if the model is a good fit to the data.

The first argument, object, is usually a fitted point process model (object of class "ppm"), obtained from the model-fitting function ppm.

For convenience, object can also be a point pattern (object of class "ppp"). In that case, a point process model will be fitted to it, by calling ppm using the arguments trend (for the first order trend), interaction (for the interpoint interaction) and rbord (for the erosion distance in the border correction for the pseudolikelihood). See ppm for details of these arguments.

The algorithm first extracts the original point pattern dataset (to which the model was fitted) and computes the standard nonparametric estimates of the K function. It then also computes the model compensator of the K function. The different function estimates are returned as columns in a data frame (of class "fv").

The argument correction determines the edge correction(s) to be applied. See Kest for explanation of the principle of edge corrections. The following table gives the options for the correction argument, and the corresponding column names in the result:

The nonparametric estimates can all be expressed in the form

\hat K(r) = \sum_i \sum_{j < i} e(x_i,x_j,r,x) I\{ d(x_i,x_j) \le r \}

where x_i is the i-th data point, d(x_i,x_j) is the distance between x_i and x_j, and e(x_i,x_j,r,x) is a term that serves to correct edge effects and to re-normalise the sum. The corresponding model compensator is

{\bf C} \, \tilde K(r) = \int_W \lambda(u,x) \sum_j e(u,x_j,r,x \cup u) I\{ d(u,x_j) \le r\}

where the integral is over all locations u in the observation window, \lambda(u,x) denotes the conditional intensity of the model at the location u, and x \cup u denotes the data point pattern x augmented by adding the extra point u.

If the fitted model is a Poisson point process, then the formulae above are exactly what is computed. If the fitted model is not Poisson, the formulae above are modified slightly to handle edge effects.

The modification is determined by the arguments conditional and restrict. The value of conditional defaults to FALSE for Poisson models and TRUE for non-Poisson models. If conditional=FALSE then the formulae above are not modified. If conditional=TRUE, then the algorithm calculates the restriction estimator if restrict=TRUE, and calculates the reweighting estimator if restrict=FALSE. See Appendix D of Baddeley, Rubak and Moller (2011). Thus, by default, the reweighting estimator is computed for non-Poisson models.

The nonparametric estimates of K(r) are approximately unbiased estimates of the K-function, assuming the point process is stationary. The model compensators are unbiased estimates of the mean values of the corresponding nonparametric estimates, assuming the model is true. Thus, if the model is a good fit, the mean value of the difference between the nonparametric estimates and model compensators is approximately zero.

Value

A function value table (object of class "fv"), essentially a data frame of function values. There is a plot method for this class. See fv.object.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ege Rubak rubak@math.aau.dk and Jesper Moller.

References

Baddeley, A., Rubak, E. and Moller, J. (2011) Score, pseudo-score and residual diagnostics for spatial point process models. Statistical Science 26, 613–646.

See Also

Related functions: Kres, Kest.

Alternative functions: Gcom, psstG, psstA, psst.

Point process models: ppm.

Examples

    fit0 <- ppm(cells, ~1) # uniform Poisson
    

    if(interactive()) {
      plot(Kcom(fit0))
# compare the isotropic-correction estimates
      plot(Kcom(fit0), cbind(iso, icom) ~ r)
# uniform Poisson is clearly not correct
    }

    fit1 <- ppm(cells, ~1, Strauss(0.08))
    
    K1 <- Kcom(fit1)
    K1
    if(interactive()) {
      plot(K1)
      plot(K1, cbind(iso, icom) ~ r)
      plot(K1, cbind(trans, tcom) ~ r)
# how to plot the difference between nonparametric estimates and compensators
      plot(K1, iso - icom ~ r)
# fit looks approximately OK; try adjusting interaction distance
    }
    fit2 <- ppm(cells, ~1, Strauss(0.12))
    
    K2 <- Kcom(fit2)
    if(interactive()) {
      plot(K2)
      plot(K2, cbind(iso, icom) ~ r)
      plot(K2, iso - icom ~ r)
    }

K Function or Pair Correlation Function of a Point Process Model

Description

Returns the theoretical K function or the pair correlation function of a point process model.

Usage

   Kmodel(model, ...)

   pcfmodel(model, ...)

Arguments

Details

For certain types of point process models, it is possible to write down a mathematical expression for the K function or the pair correlation function of the model.

The functions Kmodel and pcfmodel give the theoretical K-function and the theoretical pair correlation function for a point process model that has been fitted to data.

The functions Kmodel and pcfmodel are generic, with methods for the classes "kppm" (cluster processes and Cox processes) and "ppm" (Gibbs processes).

The return value is a function in the R language, which takes one argument r. Evaluation of this function, on a numeric vector r, yields values of the desired K function or pair correlation function at these distance values.

Value

A function in the R language, which takes one argument r.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

Kest or pcf to estimate the K function or pair correlation function nonparametrically from data.

Kmodel.kppm for the method for cluster processes and Cox processes.

Kmodel.ppm for the method for Gibbs processes.

K-function or Pair Correlation Function of a Determinantal Point Process Model

Description

Returns the theoretical K-function or theoretical pair correlation function of a determinantal point process model as a function of one argument r.

Usage

   ## S3 method for class 'dppm'
Kmodel(model, ...)

   ## S3 method for class 'dppm'
pcfmodel(model, ...)

   ## S3 method for class 'detpointprocfamily'
Kmodel(model, ...)

   ## S3 method for class 'detpointprocfamily'
pcfmodel(model, ...)

Arguments

Value

A function in the R language, with one numeric argument r, that can be used to evaluate the theoretical K-function or pair correlation function of the model at distances r.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

Examples

  model <- dppMatern(lambda=100, alpha=.01, nu=1, d=2)
  KMatern <- Kmodel(model)
  pcfMatern <- pcfmodel(model)
  plot(KMatern, xlim = c(0,0.05))
  plot(pcfMatern, xlim = c(0,0.05))

K Function or Pair Correlation Function of Cluster Model or Cox model

Description

Returns the theoretical K function or the pair correlation function of a cluster point process model or Cox point process model.

Usage

   ## S3 method for class 'kppm'
Kmodel(model, ...)

   ## S3 method for class 'kppm'
pcfmodel(model, ...)

Arguments

Details

For certain types of point process models, it is possible to write down a mathematical expression for the K function or the pair correlation function of the model. In particular this is possible for a fitted cluster point process model (object of class "kppm" obtained from kppm).

The functions Kmodel and pcfmodel are generic. The functions documented here are the methods for the class "kppm".

The return value is a function in the R language, which takes one argument r. Evaluation of this function, on a numeric vector r, yields values of the desired K function or pair correlation function at these distance values.

Value

A function in the R language, which takes one argument r.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

See Also

Kest or pcf to estimate the K function or pair correlation function nonparametrically from data.

kppm to fit cluster models.

Kmodel for the generic functions.

Kmodel.ppm for the method for Gibbs processes.

Examples

  fit <- kppm(redwood, ~x, "MatClust")
  K <- Kmodel(fit)
  K(c(0.1, 0.2))
  curve(K(x), from=0, to=0.25)

K Function or Pair Correlation Function of Gibbs Point Process model

Description

Returns the theoretical K function or the pair correlation function of a fitted Gibbs point process model.

Usage

   ## S3 method for class 'ppm'
Kmodel(model, ...)

   ## S3 method for class 'ppm'
pcfmodel(model, ...)

Arguments

Details

This function computes an approximation to the K function or the pair correlation function of a Gibbs point process.

The functions Kmodel and pcfmodel are generic. The functions documented here are the methods for the class "ppm".

The approximation is only available for stationary pairwise-interaction models. It uses the second order Poisson-saddlepoint approximation (Baddeley and Nair, 2012b) which is a combination of the Poisson-Boltzmann-Emden and Percus-Yevick approximations.

The return value is a function in the R language, which takes one argument r. Evaluation of this function, on a numeric vector r, yields values of the desired K function or pair correlation function at these distance values.

Value

A function in the R language, which takes one argument r.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Gopalan Nair.

References

Baddeley, A. and Nair, G. (2012a) Fast approximation of the intensity of Gibbs point processes. Electronic Journal of Statistics 6 1155–1169.

Baddeley, A. and Nair, G. (2012b) Approximating the moments of a spatial point process. Stat 1, 1, 18–30. DOI: 10.1002/sta4.5

See Also

Kest or pcf to estimate the K function or pair correlation function nonparametrically from data.

ppm to fit Gibbs models.

Kmodel for the generic functions.

Kmodel.kppm for the method for cluster/Cox processes.

Examples

  fit <- ppm(swedishpines, ~1, Strauss(8))
  p <- pcfmodel(fit)
  K <- Kmodel(fit)
  p(6)
  K(8)
  curve(K(x), from=0, to=15)

Residual K Function

Description

Given a point process model fitted to a point pattern dataset, this function computes the residual K function, which serves as a diagnostic for goodness-of-fit of the model.

Usage

   Kres(object, ...)

Arguments

Details

This command provides a diagnostic for the goodness-of-fit of a point process model fitted to a point pattern dataset. It computes a residual version of the K function of the dataset, which should be approximately zero if the model is a good fit to the data.

In normal use, object is a fitted point process model or a point pattern. Then Kres first calls Kcom to compute both the nonparametric estimate of the K function and its model compensator. Then Kres computes the difference between them, which is the residual K-function.

Alternatively, object may be a function value table (object of class "fv") that was returned by a previous call to Kcom. Then Kres computes the residual from this object.

Value

A function value table (object of class "fv"), essentially a data frame of function values. There is a plot method for this class. See fv.object.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ege Rubak rubak@math.aau.dk and Jesper Moller.

References

Baddeley, A., Rubak, E. and Moller, J. (2011) Score, pseudo-score and residual diagnostics for spatial point process models. Statistical Science 26, 613–646.

See Also

Related functions: Kcom, Kest.

Alternative functions: Gres, psstG, psstA, psst.

Point process models: ppm.

Examples

    fit0 <- ppm(cells ~1) # uniform Poisson
    
    K0 <- Kres(fit0)
    K0
    plot(K0)
# isotropic-correction estimate
    plot(K0, ires ~ r)
# uniform Poisson is clearly not correct

    fit1 <- ppm(cells ~1, Strauss(0.08))
    
    K1 <- Kres(fit1)

    if(interactive()) {
      plot(K1, ires ~ r)
   # fit looks approximately OK; try adjusting interaction distance
      plot(Kres(cells, interaction=Strauss(0.12)))
    }

# How to make envelopes
    
      E <- envelope(fit1, Kres, model=fit1, nsim=19)
      plot(E)
    

# For computational efficiency
    Kc <- Kcom(fit1)
    K1 <- Kres(Kc)

Lambert's W Function

Description

Computes Lambert's W-function.

Usage

LambertW(x)

Arguments

Details

Lambert's W-function is the inverse function of f(y) = y e^y. That is, W is the function such that

W(x) e^{W(x)} = x

This command LambertW computes W(x) for each entry in the argument x. If the library gsl has been installed, then the function lambert_W0 in that library is invoked. Otherwise, values of the W-function are computed by root-finding, using the function uniroot.

Computation using gsl is about 100 times faster.

If any entries of x are infinite or NA, the corresponding results are NA.

Value

Numeric vector.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

and Rolf Turner rolfturner@posteo.net

References

Corless, R, Gonnet, G, Hare, D, Jeffrey, D and Knuth, D (1996), On the Lambert W function. Computational Mathematics, 5, 325–359.

Roy, R and Olver, F (2010), Lambert W function. In Olver, F, Lozier, D and Boisvert, R (eds.), NIST Handbook of Mathematical Functions, Cambridge University Press.

Examples

   LambertW(exp(1))

The Lennard-Jones Potential

Description

Creates the Lennard-Jones pairwise interaction structure which can then be fitted to point pattern data.

Usage

  LennardJones(sigma0=NA)

Arguments

Details

In a pairwise interaction point process with the Lennard-Jones pair potential (Lennard-Jones, 1924) each pair of points in the point pattern, a distance d apart, contributes a factor

v(d) = \exp \left\{ - 4\epsilon \left[ \left( \frac{\sigma}{d} \right)^{12} - \left( \frac{\sigma}{d} \right)^6 \right] \right\}

to the probability density, where \sigma and \epsilon are positive parameters to be estimated.

See Examples for a plot of this expression.

This potential causes very strong inhibition between points at short range, and attraction between points at medium range. The parameter \sigma is called the characteristic diameter and controls the scale of interaction. The parameter \epsilon is called the well depth and determines the strength of attraction. The potential switches from inhibition to attraction at d=\sigma. The maximum value of the pair potential is \exp(\epsilon) occuring at distance d = 2^{1/6} \sigma. Interaction is usually considered to be negligible for distances d > 2.5 \sigma \max\{1,\epsilon^{1/6}\}.

This potential is used to model interactions between uncharged molecules in statistical physics.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the Lennard-Jones pairwise interaction is yielded by the function LennardJones(). See the examples below.

Value

An object of class "interact" describing the Lennard-Jones interpoint interaction structure.

Rescaling

To avoid numerical instability, the interpoint distances d are rescaled when fitting the model.

Distances are rescaled by dividing by sigma0. In the formula for v(d) above, the interpoint distance d will be replaced by d/sigma0.

The rescaling happens automatically by default. If the argument sigma0 is missing or NA (the default), then sigma0 is taken to be the minimum nearest-neighbour distance in the data point pattern (in the call to ppm).

If the argument sigma0 is given, it should be a positive number, and it should be a rough estimate of the parameter \sigma.

The “canonical regular parameters” estimated by ppm are \theta_1 = 4 \epsilon (\sigma/\sigma_0)^{12} and \theta_2 = 4 \epsilon (\sigma/\sigma_0)^6.

Warnings and Errors

Fitting the Lennard-Jones model is extremely unstable, because of the strong dependence between the functions d^{-12} and d^{-6}. The fitting algorithm often fails to converge. Try increasing the number of iterations of the GLM fitting algorithm, by setting gcontrol=list(maxit=1e3) in the call to ppm.

Errors are likely to occur if this model is fitted to a point pattern dataset which does not exhibit both short-range inhibition and medium-range attraction between points. The values of the parameters \sigma and \epsilon may be NA (because the fitted canonical parameters have opposite sign, which usually occurs when the pattern is completely random).

An absence of warnings does not mean that the fitted model is sensible. A negative value of \epsilon may be obtained (usually when the pattern is strongly clustered); this does not correspond to a valid point process model, but the software does not issue a warning.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

and Rolf Turner rolfturner@posteo.net

References

Lennard-Jones, J.E. (1924) On the determination of molecular fields. Proc Royal Soc London A 106, 463–477.

See Also

ppm, pairwise.family, ppm.object

Examples

   badfit <- ppm(cells ~1, LennardJones(), rbord=0.1)
   badfit

   fit <- ppm(unmark(longleaf) ~1, LennardJones(), rbord=1)
   fit
   plot(fitin(fit))
   # Note the Longleaf Pines coordinates are rounded to the nearest decimetre
   # (multiple of 0.1 metres) so the apparent inhibition may be an artefact

The Multitype Hard Core Point Process Model

Description

Creates an instance of the multitype hard core point process model which can then be fitted to point pattern data.

Usage

  MultiHard(hradii, types=NULL)

Arguments

Details

This is a multitype version of the hard core process. A pair of points of types i and j must not lie closer than h_{ij} units apart.

The argument types need not be specified in normal use. It will be determined automatically from the point pattern data set to which the MultiStrauss interaction is applied, when the user calls ppm. However, the user should be confident that the ordering of types in the dataset corresponds to the ordering of rows and columns in the matrix hradii.

The matrix hradii must be symmetric, with entries which are either positive numbers or NA. A value of NA indicates that no distance constraint should be applied for this combination of types.

Note that only the hardcore radii are specified in MultiHard. The canonical parameters \log(\beta_j) are estimated by ppm(), not fixed in MultiHard().

Value

An object of class "interact" describing the interpoint interaction structure of the multitype hard core process with hard core radii hradii[i,j].

Warnings

In order that ppm can fit the multitype hard core model correctly to a point pattern X, this pattern must be marked, with markformat equal to vector and the mark vector marks(X) must be a factor. If the argument types is specified it is interpreted as a set of factor levels and this set must equal levels(marks(X)).

Changed Syntax

Before spatstat version 1.37-0, the syntax of this function was different: MultiHard(types=NULL, hradii). The new code attempts to handle the old syntax as well.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

ppm, pairwise.family, ppm.object, MultiStrauss, MultiStraussHard, Strauss.

See ragsMultiHard and rmh for simulation.

Examples

   h <- matrix(c(1,2,2,1), nrow=2,ncol=2)

   # prints a sensible description of itself
   MultiHard(h)

   # Fit the stationary multitype hardcore process to `amacrine'
   # with hard core operating only between cells of the same type.
   h <- 0.02 * matrix(c(1, NA, NA, 1), nrow=2,ncol=2)
   ppm(amacrine ~1, MultiHard(h))

The Multitype Strauss Point Process Model

Description

Creates an instance of the multitype Strauss point process model which can then be fitted to point pattern data.

Usage

  MultiStrauss(radii, types=NULL)

Arguments

Details

The (stationary) multitype Strauss process with m types, with interaction radii r_{ij} and parameters \beta_j and \gamma_{ij} is the pairwise interaction point process in which each point of type j contributes a factor \beta_j to the probability density of the point pattern, and a pair of points of types i and j closer than r_{ij} units apart contributes a factor \gamma_{ij} to the density.

The nonstationary multitype Strauss process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location and type, rather than a constant beta.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the multitype Strauss process pairwise interaction is yielded by the function MultiStrauss(). See the examples below.

The argument types need not be specified in normal use. It will be determined automatically from the point pattern data set to which the MultiStrauss interaction is applied, when the user calls ppm. However, the user should be confident that the ordering of types in the dataset corresponds to the ordering of rows and columns in the matrix radii.

The matrix radii must be symmetric, with entries which are either positive numbers or NA. A value of NA indicates that no interaction term should be included for this combination of types.

Note that only the interaction radii are specified in MultiStrauss. The canonical parameters \log(\beta_j) and \log(\gamma_{ij}) are estimated by ppm(), not fixed in MultiStrauss().

Value

An object of class "interact" describing the interpoint interaction structure of the multitype Strauss process with interaction radii radii[i,j].

Warnings

In order that ppm can fit the multitype Strauss model correctly to a point pattern X, this pattern must be marked, with markformat equal to vector and the mark vector marks(X) must be a factor. If the argument types is specified it is interpreted as a set of factor levels and this set must equal levels(marks(X)).

Changed Syntax

Before spatstat version 1.37-0, the syntax of this function was different: MultiStrauss(types=NULL, radii). The new code attempts to handle the old syntax as well.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

ppm, pairwise.family, ppm.object, Strauss, MultiHard

Examples

   r <- matrix(c(1,2,2,1), nrow=2,ncol=2)
   MultiStrauss(r)
   # prints a sensible description of itself
   r <- 0.03 * matrix(c(1,2,2,1), nrow=2,ncol=2)
   X <- amacrine
   
   ppm(X ~1, MultiStrauss(r))
   # fit the stationary multitype Strauss process to `amacrine'

   
     ppm(X ~polynom(x,y,3), MultiStrauss(r, c("off","on")))
     # fit a nonstationary multitype Strauss process with log-cubic trend
   

The Multitype/Hard Core Strauss Point Process Model

Description

Creates an instance of the multitype/hard core Strauss point process model which can then be fitted to point pattern data.

Usage

  MultiStraussHard(iradii, hradii, types=NULL)

Arguments

Details

This is a hybrid of the multitype Strauss process (see MultiStrauss) and the hard core process (case \gamma=0 of the Strauss process). A pair of points of types i and j must not lie closer than h_{ij} units apart; if the pair lies more than h_{ij} and less than r_{ij} units apart, it contributes a factor \gamma_{ij} to the probability density.

The argument types need not be specified in normal use. It will be determined automatically from the point pattern data set to which the MultiStraussHard interaction is applied, when the user calls ppm. However, the user should be confident that the ordering of types in the dataset corresponds to the ordering of rows and columns in the matrices iradii and hradii.

The matrices iradii and hradii must be symmetric, with entries which are either positive numbers or NA. A value of NA indicates that no interaction term should be included for this combination of types.

Note that only the interaction radii and hardcore radii are specified in MultiStraussHard. The canonical parameters \log(\beta_j) and \log(\gamma_{ij}) are estimated by ppm(), not fixed in MultiStraussHard().

Value

An object of class "interact" describing the interpoint interaction structure of the multitype/hard core Strauss process with interaction radii iradii[i,j] and hard core radii hradii[i,j].

Warnings

In order that ppm can fit the multitype/hard core Strauss model correctly to a point pattern X, this pattern must be marked, with markformat equal to vector and the mark vector marks(X) must be a factor. If the argument types is specified it is interpreted as a set of factor levels and this set must equal levels(marks(X)).

Changed Syntax

Before spatstat version 1.37-0, the syntax of this function was different: MultiStraussHard(types=NULL, iradii, hradii). The new code attempts to handle the old syntax as well.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

, Rolf Turner rolfturner@posteo.net

and Ege Rubak rubak@math.aau.dk

See Also

ppm, pairwise.family, ppm.object, MultiStrauss, MultiHard, Strauss

Examples

   r <- matrix(3, nrow=2,ncol=2)
   h <- matrix(c(1,2,2,1), nrow=2,ncol=2)
   MultiStraussHard(r,h)
   # prints a sensible description of itself
   r <- 0.04 * matrix(c(1,2,2,1), nrow=2,ncol=2)
   h <- 0.02 * matrix(c(1,NA,NA,1), nrow=2,ncol=2)
   X <- amacrine
   
   fit <- ppm(X ~1, MultiStraussHard(r,h))
   # fit stationary multitype hardcore Strauss process to `amacrine'

Arithmetic Operations on Measures

Description

These group generic methods for the class "msr" allow the arithmetic operators +, -, * and / to be applied directly to measures.

Usage

## S3 methods for group generics have prototypes:
Ops(e1, e2)

Arguments

Details

Arithmetic operators on a measure A are only defined in some cases. The arithmetic operator is effectively applied to the value of A(W) for every spatial domain W. If the result is a measure, then this operation is valid.

If A is a measure (object of class "msr") then the operations -A and +A are defined.

If A and B are measures with the same dimension (i.e. both are scalar-valued, or both are k-dimensional vector-valued) then A + B and A - B are defined.

If A is a measure and z is a numeric value, then A * z and A / z are defined, and z * A is defined.

Value

Another measure (object of class "msr").

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

with.msr

Examples

   X <- rpoispp(function(x,y) { exp(3+3*x) })
   fit <- ppm(X, ~x+y)
   rp <- residuals(fit, type="pearson")
   rp

   -rp
   2 * rp
   rp /2

   rp - rp

   rr <- residuals(fit, type="raw")
   rp - rr

Generic Ord Interaction model

Description

Creates an instance of an Ord-type interaction point process model which can then be fitted to point pattern data.

Usage

  Ord(pot, name)

Arguments

Details

Ord's point process model (Ord, 1977) is a Gibbs point process of infinite order. Each point x_i in the point pattern x contributes a factor g(a_i) where a_i = a(x_i, x) is the area of the tile associated with x_i in the Dirichlet tessellation of x.

Ord (1977) proposed fitting this model to forestry data when g(a) has a simple “threshold” form. That model is implemented in our function OrdThresh. The present function Ord implements the case of a completely general Ord potential g(a) specified as an S language function pot.

This is experimental.

Value

An object of class "interact" describing the interpoint interaction structure of a point process.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

and Rolf Turner rolfturner@posteo.net

References

Baddeley, A. and Turner, R. (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42, 283–322.

Ord, J.K. (1977) Contribution to the discussion of Ripley (1977).

Ord, J.K. (1978) How many trees in a forest? Mathematical Scientist 3, 23–33.

Ripley, B.D. (1977) Modelling spatial patterns (with discussion). Journal of the Royal Statistical Society, Series B, 39, 172 – 212.

See Also

ppm, ppm.object, OrdThresh

Ord's Interaction model

Description

Creates an instance of Ord's point process model which can then be fitted to point pattern data.

Usage

  OrdThresh(r)

Arguments

Details

Ord's point process model (Ord, 1977) is a Gibbs point process of infinite order. Each point x_i in the point pattern x contributes a factor g(a_i) where a_i = a(x_i, x) is the area of the tile associated with x_i in the Dirichlet tessellation of x. The function g is simply g(a) = 1 if a \ge r and g(a) = \gamma < 1 if a < r, where r is called the threshold value.

This function creates an instance of Ord's model with a given value of r. It can then be fitted to point process data using ppm.

Value

An object of class "interact" describing the interpoint interaction structure of a point process.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

and Rolf Turner rolfturner@posteo.net

References

Baddeley, A. and Turner, R. (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42, 283–322.

Ord, J.K. (1977) Contribution to the discussion of Ripley (1977).

Ord, J.K. (1978) How many trees in a forest? Mathematical Scientist 3, 23–33.

Ripley, B.D. (1977) Modelling spatial patterns (with discussion). Journal of the Royal Statistical Society, Series B, 39, 172 – 212.

See Also

ppm, ppm.object

The Piecewise Constant Pairwise Interaction Point Process Model

Description

Creates an instance of a pairwise interaction point process model with piecewise constant potential function. The model can then be fitted to point pattern data.

Usage

  PairPiece(r)

Arguments

Details

A pairwise interaction point process in a bounded region is a stochastic point process with probability density of the form

f(x_1,\ldots,x_n) = \alpha \prod_i b(x_i) \prod_{i < j} h(x_i, x_j)

where x_1,\ldots,x_n represent the points of the pattern. The first product on the right hand side is over all points of the pattern; the second product is over all unordered pairs of points of the pattern.

Thus each point x_i of the pattern contributes a factor b(x_i) to the probability density, and each pair of points x_i, x_j contributes a factor h(x_i,x_j) to the density.

The pairwise interaction term h(u, v) is called piecewise constant if it depends only on the distance between u and v, say h(u,v) = H(||u-v||), and H is a piecewise constant function (a function which is constant except for jumps at a finite number of places). The use of piecewise constant interaction terms was first suggested by Takacs (1986).

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the piecewise constant pairwise interaction is yielded by the function PairPiece(). See the examples below.

The entries of r must be strictly increasing, positive numbers. They are interpreted as the points of discontinuity of H. It is assumed that H(s) =1 for all s > r_{max} where r_{max} is the maximum value in r. Thus the model has as many regular parameters (see ppm) as there are entries in r. The i-th regular parameter \theta_i is the logarithm of the value of the interaction function H on the interval [r_{i-1},r_i).

If r is a single number, this model is similar to the Strauss process, see Strauss. The difference is that in PairPiece the interaction function is continuous on the right, while in Strauss it is continuous on the left.

The analogue of this model for multitype point processes has not yet been implemented.

Value

An object of class "interact" describing the interpoint interaction structure of a point process. The process is a pairwise interaction process, whose interaction potential is piecewise constant, with jumps at the distances given in the vector r.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Takacs, R. (1986) Estimator for the pair potential of a Gibbsian point process. Statistics 17, 429–433.

See Also

ppm, pairwise.family, ppm.object, Strauss rmh.ppm

Examples

   PairPiece(c(0.1,0.2))
   # prints a sensible description of itself

   ppm(cells ~1, PairPiece(r = c(0.05, 0.1, 0.2)))
   # fit a stationary piecewise constant pairwise interaction process

   
     ppm(cells ~polynom(x,y,3), PairPiece(c(0.05, 0.1)))
     # nonstationary process with log-cubic polynomial trend
   

Generic Pairwise Interaction model

Description

Creates an instance of a pairwise interaction point process model which can then be fitted to point pattern data.

Usage

  Pairwise(pot, name, par, parnames, printfun)

Arguments

Details

This code constructs a member of the pairwise interaction family pairwise.family with arbitrary pairwise interaction potential given by the user.

Each pair of points in the point pattern contributes a factor h(d) to the probability density, where d is the distance between the two points. The factor term h(d) is

h(d) = \exp(-\theta \mbox{pot}(d))

provided \mbox{pot}(d) is finite, where \theta is the coefficient vector in the model.

The function pot must take as its first argument a matrix of interpoint distances, and evaluate the potential for each of these distances. The result must be either a matrix with the same dimensions as its input, or an array with its first two dimensions the same as its input (the latter case corresponds to a vector-valued potential).

If irregular parameters are present, then the second argument to pot should be a vector of the same type as par giving those parameter values.

The values returned by pot may be finite numeric values, or -Inf indicating a hard core (that is, the corresponding interpoint distance is forbidden). We define h(d) = 0 if \mbox{pot}(d) = -\infty. Thus, a potential value of minus infinity is always interpreted as corresponding to h(d) = 0, regardless of the sign and magnitude of \theta.

Value

An object of class "interact" describing the interpoint interaction structure of a point process.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

See Also

ppm, pairwise.family, ppm.object

Examples

   #This is the same as StraussHard(r=0.7,h=0.05)
   strpot <- function(d,par) {
         r <- par$r
         h <- par$h
         value <- (d <= r)
         value[d < h] <- -Inf
         value
   }
   mySH <- Pairwise(strpot, "StraussHard process", list(r=0.7,h=0.05),
           c("interaction distance r", "hard core distance h"))
   ppm(cells ~ 1, mySH, correction="isotropic")

   # Fiksel (1984) double exponential interaction
   # see Stoyan, Kendall, Mecke 1987 p 161

   fikspot <- function(d, par) {
      r <- par$r
      h <- par$h
      zeta <- par$zeta
      value <- exp(-zeta * d)
      value[d < h] <- -Inf
      value[d > r] <- 0
      value
   }
   Fiksel <- Pairwise(fikspot, "Fiksel double exponential process",
                      list(r=3.5, h=1, zeta=1),
                      c("interaction distance r",
                        "hard core distance h",
                        "exponential coefficient zeta"))
   fit <- ppm(unmark(spruces) ~1, Fiksel, rbord=3.5)
   fit
   plot(fitin(fit), xlim=c(0,4))
   coef(fit)
   # corresponding values obtained by Fiksel (1984) were -1.9 and -6.0

Penttinen Interaction

Description

Creates an instance of the Penttinen pairwise interaction point process model, which can then be fitted to point pattern data.

Usage

  Penttinen(r)

Arguments

Details

Penttinen (1984, Example 2.1, page 18), citing Cormack (1979), described the pairwise interaction point process with interaction factor

h(d) = e^{\theta A(d)} = \gamma^{A(d)}

between each pair of points separated by a distance $d$. Here A(d) is the area of intersection between two discs of radius r separated by a distance d, normalised so that A(0) = 1.

The scale of interaction is controlled by the disc radius r: two points interact if they are closer than 2 r apart. The strength of interaction is controlled by the canonical parameter \theta, which must be less than or equal to zero, or equivalently by the parameter \gamma = e^\theta, which must lie between 0 and 1.

The potential is inhibitory, i.e.\ this model is only appropriate for regular point patterns. For \gamma=0 the model is a hard core process with hard core diameter 2 r. For \gamma=1 the model is a Poisson process.

The irregular parameter r must be given in the call to Penttinen, while the regular parameter \theta will be estimated.

This model can be considered as a pairwise approximation to the area-interaction model AreaInter.

Value

An object of class "interact" describing the interpoint interaction structure of a point process.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

References

Cormack, R.M. (1979) Spatial aspects of competition between individuals. Pages 151–212 in Spatial and Temporal Analysis in Ecology, eds. R.M. Cormack and J.K. Ord, International Co-operative Publishing House, Fairland, MD, USA.

Penttinen, A. (1984) Modelling Interaction in Spatial Point Patterns: Parameter Estimation by the Maximum Likelihood Method. Jyvaskyla Studies in Computer Science, Economics and Statistics 7, University of Jyvaskyla, Finland.

See Also

ppm, ppm.object, Pairwise, AreaInter.

Examples

   fit <- ppm(cells ~ 1, Penttinen(0.07))
   fit
   reach(fit) # interaction range is circle DIAMETER

Poisson Point Process Model

Description

Creates an instance of the Poisson point process model which can then be fitted to point pattern data.

Usage

 Poisson()

Details

The function ppm, which fits point process models to point pattern data, requires an argument interaction of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the Poisson process is provided by the value of the function Poisson.

This works for all types of Poisson processes including multitype and nonstationary Poisson processes.

Value

An object of class "interact" describing the interpoint interaction structure of the Poisson point process (namely, there are no interactions).

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

See Also

ppm, Strauss

Examples

 ppm(nztrees ~1, Poisson())
 # fit the stationary Poisson process to 'nztrees'
 # no edge correction needed

 lon <- longleaf
 
 longadult <- unmark(subset(lon, marks >= 30))
 ppm(longadult ~ x, Poisson())
 # fit the nonstationary Poisson process 
 # with intensity lambda(x,y) = exp( a + bx)

 # trees marked by species
 lans <- lansing
 
 ppm(lans ~ marks, Poisson())
 # fit stationary marked Poisson process
 # with different intensity for each species

 
   ppm(lansing ~ marks * polynom(x,y,3), Poisson())
   # fit nonstationary marked Poisson process
   # with different log-cubic trend for each species
 
 

Piecewise Constant Saturated Pairwise Interaction Point Process Model

Description

Creates an instance of a saturated pairwise interaction point process model with piecewise constant potential function. The model can then be fitted to point pattern data.

Usage

  SatPiece(r, sat)

Arguments

Details

This is a generalisation of the Geyer saturation point process model, described in Geyer, to the case of multiple interaction distances. It can also be described as the saturated analogue of a pairwise interaction process with piecewise-constant pair potential, described in PairPiece.

The saturated point process with interaction radii r_1,\ldots,r_k, saturation thresholds s_1,\ldots,s_k, intensity parameter \beta and interaction parameters \gamma_1,\ldots,gamma_k, is the point process in which each point x_i in the pattern X contributes a factor

\beta \gamma_1^{v_1(x_i, X)} \ldots gamma_k^{v_k(x_i,X)}

to the probability density of the point pattern, where

v_j(x_i, X) = \min( s_j, t_j(x_i,X) )

where t_j(x_i, X) denotes the number of points in the pattern X which lie at a distance between r_{j-1} and r_j from the point x_i. We take r_0 = 0 so that t_1(x_i,X) is the number of points of X that lie within a distance r_1 of the point x_i.

SatPiece is used to fit this model to data. The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the piecewise constant Saturated pairwise interaction is yielded by the function SatPiece(). See the examples below.

Simulation of this point process model is not yet implemented. This model is not locally stable (the conditional intensity is unbounded).

The argument r specifies the vector of interaction distances. The entries of r must be strictly increasing, positive numbers.

The argument sat specifies the vector of saturation parameters. It should be a vector of the same length as r, and its entries should be nonnegative numbers. Thus sat[1] corresponds to the distance range from 0 to r[1], and sat[2] to the distance range from r[1] to r[2], etc. Alternatively sat may be a single number, and this saturation value will be applied to every distance range.

Infinite values of the saturation parameters are also permitted; in this case v_j(x_i,X) = t_j(x_i,X) and there is effectively no ‘saturation’ for the distance range in question. If all the saturation parameters are set to Inf then the model is effectively a pairwise interaction process, equivalent to PairPiece (however the interaction parameters \gamma obtained from SatPiece are the square roots of the parameters \gamma obtained from PairPiece).

If r is a single number, this model is virtually equivalent to the Geyer process, see Geyer.

Value

An object of class "interact" describing the interpoint interaction structure of a point process.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net in collaboration with Hao Wang and Jeff Picka

See Also

ppm, pairsat.family, Geyer, PairPiece, BadGey.

Examples

   SatPiece(c(0.1,0.2), c(1,1))
   # prints a sensible description of itself
   SatPiece(c(0.1,0.2), 1)

   ppm(cells ~1, SatPiece(c(0.07, 0.1, 0.13), 2))
   # fit a stationary piecewise constant Saturated pairwise interaction process

   
     ppm(cells ~polynom(x,y,3), SatPiece(c(0.07, 0.1, 0.13), 2))
     # nonstationary process with log-cubic polynomial trend
   

Saturated Pairwise Interaction model

Description

Experimental.

Usage

  Saturated(pot, name)

Arguments

Details

This is experimental. It constructs a member of the “saturated pairwise” family pairsat.family.

Value

An object of class "interact" describing the interpoint interaction structure of a point process.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

and Rolf Turner rolfturner@posteo.net

See Also

ppm, pairsat.family, Geyer, SatPiece, ppm.object

Smooth a Signed or Vector-Valued Measure

Description

Apply kernel smoothing to a signed measure or vector-valued measure.

Usage

 ## S3 method for class 'msr'
Smooth(X, ..., drop=TRUE)

Arguments

Details

This function applies kernel smoothing to a signed measure or vector-valued measure X. The Gaussian kernel is used.

The object X would typically have been created by residuals.ppm or msr.

Value

A pixel image or a list of pixel images. For scalar-valued measures, a pixel image (object of class "im") provided drop=TRUE. For vector-valued measures (or if drop=FALSE), a list of pixel images; the list also belongs to the class "solist" so that it can be printed and plotted.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au

References

Baddeley, A., Turner, R., Moller, J. and Hazelton, M. (2005) Residual analysis for spatial point processes. Journal of the Royal Statistical Society, Series B 67, 617–666.

Baddeley, A., Moller, J. and Pakes, A.G. (2008) Properties of residuals for spatial point processes. Annals of the Institute of Statistical Mathematics 60, 627–649.

See Also

Smooth, msr, plot.msr

Examples

   X <- rpoispp(function(x,y) { exp(3+3*x) })
   fit <- ppm(X, ~x+y)
   rp <- residuals(fit, type="pearson")
   rs <- residuals(fit, type="score")

   plot(Smooth(rp))
   plot(Smooth(rs))

The Soft Core Point Process Model

Description

Creates an instance of the Soft Core point process model which can then be fitted to point pattern data.

Usage

  Softcore(kappa, sigma0=NA)

Arguments

Details

The (stationary) Soft Core point process with parameters \beta and \sigma and exponent \kappa is the pairwise interaction point process in which each point contributes a factor \beta to the probability density of the point pattern, and each pair of points contributes a factor

\exp \left\{ - \left( \frac{\sigma}{d} \right)^{2/\kappa} \right\}

to the density, where d is the distance between the two points. See the Examples for a plot of this interaction curve.

Thus the process has probability density

f(x_1,\ldots,x_n) = \alpha \beta^{n(x)} \exp \left\{ - \sum_{i < j} \left( \frac{\sigma}{||x_i-x_j||} \right)^{2/\kappa} \right\}

where x_1,\ldots,x_n represent the points of the pattern, n(x) is the number of points in the pattern, \alpha is the normalising constant, and the sum on the right hand side is over all unordered pairs of points of the pattern.

This model describes an “ordered” or “inhibitive” process, with the strength of inhibition decreasing smoothly with distance. The interaction is controlled by the parameters \sigma and \kappa.

• The spatial scale of interaction is controlled by the parameter \sigma, which is a positive real number interpreted as a distance, expressed in the same units of distance as the spatial data. The parameter \sigma is the distance at which the pair potential reaches the threshold value 0.37.

• The shape of the interaction function is controlled by the exponent \kappa which is a dimensionless number in the range (0,1), with larger values corresponding to a flatter shape (or a more gradual decay rate). The process is well-defined only for \kappa in (0,1). The limit of the model as \kappa \to 0 is the hard core process with hard core distance h=\sigma.

• The “strength” of the interaction is determined by both of the parameters \sigma and \kappa. The larger the value of \kappa, the wider the range of distances over which the interaction has an effect. If \sigma is very small, the interaction is very weak for all practical purposes (theoretically if \sigma = 0 the model reduces to the Poisson point process).

The nonstationary Soft Core process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location, rather than a constant beta.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the Soft Core process pairwise interaction is yielded by the function Softcore(). See the examples below.

The main argument is the exponent kappa. When kappa is fixed, the model becomes an exponential family with canonical parameters \log \beta and

\log \gamma = \frac{2}{\kappa} \log\sigma

The canonical parameters are estimated by ppm(), not fixed in Softcore().

The optional argument sigma0 can be used to improve numerical stability. If sigma0 is given, it should be a positive number, and it should be a rough estimate of the parameter \sigma.

Value

An object of class "interact" describing the interpoint interaction structure of the Soft Core process with exponent \kappa.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

References

Ogata, Y, and Tanemura, M. (1981). Estimation of interaction potentials of spatial point patterns through the maximum likelihood procedure. Annals of the Institute of Statistical Mathematics, B 33, 315–338.

Ogata, Y, and Tanemura, M. (1984). Likelihood analysis of spatial point patterns. Journal of the Royal Statistical Society, series B 46, 496–518.

See Also

ppm, pairwise.family, ppm.object

Examples

   # fit the stationary Soft Core process to `cells' 
   fit5 <- ppm(cells ~1, Softcore(kappa=0.5), correction="isotropic")

   # study shape of interaction and explore effect of parameters
   fit2 <- update(fit5, Softcore(kappa=0.2))
   fit8 <- update(fit5, Softcore(kappa=0.8))
   plot(fitin(fit2), xlim=c(0, 0.4),
        main="Pair potential (sigma = 0.1)", 
        xlab=expression(d), ylab=expression(h(d)), legend=FALSE)
   plot(fitin(fit5), add=TRUE, col=4)
   plot(fitin(fit8), add=TRUE, col=3)
   legend("bottomright", col=c(1,4,3), lty=1,
          legend=expression(kappa==0.2, kappa==0.5, kappa==0.8))

The Strauss Point Process Model

Description

Creates an instance of the Strauss point process model which can then be fitted to point pattern data.

Usage

  Strauss(r)

Arguments

Details

The (stationary) Strauss process with interaction radius r and parameters \beta and \gamma is the pairwise interaction point process in which each point contributes a factor \beta to the probability density of the point pattern, and each pair of points closer than r units apart contributes a factor \gamma to the density.

Thus the probability density is

f(x_1,\ldots,x_n) = \alpha \beta^{n(x)} \gamma^{s(x)}

where x_1,\ldots,x_n represent the points of the pattern, n(x) is the number of points in the pattern, s(x) is the number of distinct unordered pairs of points that are closer than r units apart, and \alpha is the normalising constant.

The interaction parameter \gamma must be less than or equal to 1 so that this model describes an “ordered” or “inhibitive” pattern.

The nonstationary Strauss process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location, rather than a constant beta.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the Strauss process pairwise interaction is yielded by the function Strauss(). See the examples below.

Note the only argument is the interaction radius r. When r is fixed, the model becomes an exponential family. The canonical parameters \log(\beta) and \log(\gamma) are estimated by ppm(), not fixed in Strauss().

Value

An object of class "interact" describing the interpoint interaction structure of the Strauss process with interaction radius r.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net.

References

Kelly, F.P. and Ripley, B.D. (1976) On Strauss's model for clustering. Biometrika 63, 357–360.

Strauss, D.J. (1975) A model for clustering. Biometrika 62, 467–475.

See Also

ppm, pairwise.family, ppm.object

Examples

   Strauss(r=0.1)
   # prints a sensible description of itself

   
     ppm(cells ~1, Strauss(r=0.07))
     # fit the stationary Strauss process to `cells'
   

   ppm(cells ~polynom(x,y,3), Strauss(r=0.07))
   # fit a nonstationary Strauss process with log-cubic polynomial trend

The Strauss / Hard Core Point Process Model

Description

Creates an instance of the “Strauss/ hard core” point process model which can then be fitted to point pattern data.

Usage

  StraussHard(r, hc=NA)

Arguments

Details

A Strauss/hard core process with interaction radius r, hard core distance h < r, and parameters \beta and \gamma, is a pairwise interaction point process in which

• distinct points are not allowed to come closer than a distance h apart

• each pair of points closer than r units apart contributes a factor \gamma to the probability density.

This is a hybrid of the Strauss process and the hard core process.

The probability density is zero if any pair of points is closer than h units apart, and otherwise equals

f(x_1,\ldots,x_n) = \alpha \beta^{n(x)} \gamma^{s(x)}

where x_1,\ldots,x_n represent the points of the pattern, n(x) is the number of points in the pattern, s(x) is the number of distinct unordered pairs of points that are closer than r units apart, and \alpha is the normalising constant.

The interaction parameter \gamma may take any positive value (unlike the case for the Strauss process). If \gamma < 1, the model describes an “ordered” or “inhibitive” pattern. If \gamma > 1, the model is “ordered” or “inhibitive” up to the distance h, but has an “attraction” between points lying at distances in the range between h and r.

If \gamma = 1, the process reduces to a classical hard core process with hard core distance h. If \gamma = 0, the process reduces to a classical hard core process with hard core distance r.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the Strauss/hard core process pairwise interaction is yielded by the function StraussHard(). See the examples below.

The canonical parameter \log(\gamma) is estimated by ppm(), not fixed in StraussHard().

If the hard core distance argument hc is missing or NA, it will be estimated from the data when ppm is called. The estimated value of hc is the minimum nearest neighbour distance multiplied by n/(n+1), where n is the number of data points.

Value

An object of class "interact" describing the interpoint interaction structure of the “Strauss/hard core” process with Strauss interaction radius r and hard core distance hc.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Baddeley, A. and Turner, R. (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42, 283–322.

Ripley, B.D. (1981) Spatial statistics. John Wiley and Sons.

Strauss, D.J. (1975) A model for clustering. Biometrika 62, 467–475.

See Also

ppm, pairwise.family, ppm.object

Examples

   StraussHard(r=1,hc=0.02)
   # prints a sensible description of itself

   # ppm(cells ~1, StraussHard(r=0.1, hc=0.05))
   # fit the stationary Strauss/hard core  process to `cells'

   ppm(cells ~ polynom(x,y,3), StraussHard(r=0.1, hc=0.05))
   # fit a nonstationary Strauss/hard core process
   # with log-cubic polynomial trend

The Triplet Point Process Model

Description

Creates an instance of Geyer's triplet interaction point process model which can then be fitted to point pattern data.

Usage

  Triplets(r)

Arguments

Details

The (stationary) Geyer triplet process (Geyer, 1999) with interaction radius r and parameters \beta and \gamma is the point process in which each point contributes a factor \beta to the probability density of the point pattern, and each triplet of close points contributes a factor \gamma to the density. A triplet of close points is a group of 3 points, each pair of which is closer than r units apart.

Thus the probability density is

f(x_1,\ldots,x_n) = \alpha \beta^{n(x)} \gamma^{s(x)}

where x_1,\ldots,x_n represent the points of the pattern, n(x) is the number of points in the pattern, s(x) is the number of unordered triples of points that are closer than r units apart, and \alpha is the normalising constant.

The interaction parameter \gamma must be less than or equal to 1 so that this model describes an “ordered” or “inhibitive” pattern.

The nonstationary Triplets process is similar except that the contribution of each individual point x_i is a function \beta(x_i) of location, rather than a constant beta.

The function ppm(), which fits point process models to point pattern data, requires an argument of class "interact" describing the interpoint interaction structure of the model to be fitted. The appropriate description of the Triplets process pairwise interaction is yielded by the function Triplets(). See the examples below.

Note the only argument is the interaction radius r. When r is fixed, the model becomes an exponential family. The canonical parameters \log(\beta) and \log(\gamma) are estimated by ppm(), not fixed in Triplets().

Value

An object of class "interact" describing the interpoint interaction structure of the Triplets process with interaction radius r.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Geyer, C.J. (1999) Likelihood Inference for Spatial Point Processes. Chapter 3 in O.E. Barndorff-Nielsen, W.S. Kendall and M.N.M. Van Lieshout (eds) Stochastic Geometry: Likelihood and Computation, Chapman and Hall / CRC, Monographs on Statistics and Applied Probability, number 80. Pages 79–140.

See Also

ppm, triplet.family, ppm.object

Examples

   Triplets(r=0.1)
   # prints a sensible description of itself

   ppm(cells ~1, Triplets(r=0.2))
   # fit the stationary Triplets process to `cells'

   
     ppm(cells ~polynom(x,y,3), Triplets(r=0.2))
     # fit a nonstationary Triplets process with log-cubic polynomial trend
   

Extract Window of Spatial Object

Description

Given a spatial object (such as a point pattern or pixel image) in two dimensions, these functions extract the window in which the object is defined.

Usage

 ## S3 method for class 'ppm'
Window(X, ..., from=c("points", "covariates"))

 ## S3 method for class 'kppm'
Window(X, ..., from=c("points", "covariates"))

 ## S3 method for class 'dppm'
Window(X, ..., from=c("points", "covariates"))

 ## S3 method for class 'slrm'
Window(X, ..., from=c("points", "covariates"))

 ## S3 method for class 'msr'
Window(X, ...)

Arguments

Details

These are methods for the generic function Window which extract the spatial window in which the object X is defined. The argument from applies when X is a fitted two-dimensional point process model (object of class "ppm", "kppm", "slrm" or "dppm"). If from="data" (the default), Window extracts the window of the original point pattern data to which the model was fitted. If from="covariates" then Window returns the window in which the spatial covariates of the model were provided.

Value

An object of class "owin" (see owin.object) specifying an observation window.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

Window, Window.ppp, Window.psp.

owin.object

Examples

   A <- ppm(cells ~ 1)
   Window(A)

ROC Curves for Single Term Additions to a Model

Description

Given a fitted point process model, consider adding new explanatory variables, and compute the ROC curve for each new variable.

Usage

addROC(object, scope, high=TRUE, ...)

Arguments

Details

This function is like add1 in that it considers each possible term that could be added to the model object (or only the terms listed in the scope argument), adds each such term to the model, and measures the change in the model. In this case the change is measured by computing the ROC curve for the added covariate, using the original model object as a baseline.

Either object or scope should be a fitted point process model, and the other argument may be a fitted point process model or a formula. If object is a fitted model then scope may be a character vector of the names of variables to be added.

Value

A named list containing the ROC curves for each new explanatory variable. The individual entries belong to class "fv", so they can be plotted. The list belongs to the class "anylist" so it can be plotted in its entirety.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ege Rubak rubak@math.aau.dk and Suman Rakshit Suman.Rakshit@curtin.edu.au.

References

Baddeley, A., Rubak, E., Rakshit, S. and Nair, G. (2025) ROC curves for spatial point patterns and presence-absence data. doi:10.48550/arXiv.2506.03414.

See Also

addapply, roc.ppm.

Examples

  dimyx <- if(interactive()) NULL else 32
  fit0 <- ppm(bei ~ 1, data=bei.extra)

  z <- addROC(fit0, . ~ grad + elev, dimyx=dimyx)

  plot(z)

  ## how to compute AUC for each curve
  sapply(z, auc)

Significance Tests or Effect Size for Single Term Additions to a Model

Description

Given a fitted point process model, consider adding new explanatory variables, and apply a significance test (or effect size calculation) for the effect of each new variable.

Usage

addapply(object, 
         action = c("berman.test", "cdf.test", "rhohat", "roc"),
         scope, 
         ..., high = TRUE)

Arguments

Details

This function is like add1 in that it considers adding new terms to the model object and measures the change in the model. In this case the change is measured by performing the action. Options are:

action="roc":

the ROC curve for the added covariate is computed using the original model object as a baseline.

action="berman.test":

One of Berman's tests is applied (see berman.test), using the original model object as the null hypothesis, and the extended model as the alternative.

action="cdf.test":

One of the CDF tests is applied (see cdf.test), using the original model object as the null hypothesis, and the extended model as the alternative.

action="rhohat":

taking the original model object as a baseline, the true intensity (ratio of true intensity to baseline intensity) is estimated as a function of the added explanatory variable, using the function rhohat.

Note that addapply(object, "roc", scope) is equivalent to addROC(object, scope).

Either object or scope should be a fitted point process model, and the other argument may be a fitted point process model or a formula. If object is a fitted model then scope may be a character vector of the names of variables to be added.

Value

A named list containing the results for each added variable.

The list belongs to the class "anylist" so it can be printed and plotted in its entirety.

If action="roc" the individual entries are ROC curves belonging to class "fv". If action="rhohat" the individual entries are curves belonging to class "fv" and class "rhohat". If action="berman.test" the individual entries are hypothesis tests of class "htest" and "bermantest". If action="cdf.test" the individual entries are hypothesis tests of class "htest" and "cdftest".

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ege Rubak rubak@math.aau.dk and Suman Rakshit Suman.Rakshit@curtin.edu.au.

References

Baddeley, A., Rubak, E., Rakshit, S. and Nair, G. (2025) ROC curves for spatial point patterns and presence-absence data. doi:10.48550/arXiv.2506.03414.

See Also

addROC, dropROC, roc.ppm, rhohat, berman.test, cdf.test

Examples

  dimyx <- if(interactive()) NULL else 32
  fit0 <- ppm(bei ~ 1, data=bei.extra)
  z <- addapply(fit0, "ber", . ~ grad+elev, dimyx=dimyx)
  z
  plot(z, mar.panel=5)

  ## how to extract p-values from each test
  sapply(z, getElement, name="p.value")

Added Variable Plot for Point Process Model

Description

Computes the coordinates for an Added Variable Plot for a fitted point process model.

Usage

addvar(model, covariate, ...,
                   subregion=NULL,
                   bw="nrd0", adjust=1,
                   from=NULL, to=NULL, n=512,
                   bw.input = c("points", "quad"),
                   bw.restrict = FALSE,
                   covname, crosscheck=FALSE)

Arguments

Details

This command generates the plot coordinates for an Added Variable Plot for a spatial point process model.

Added Variable Plots (Cox, 1958, sec 4.5; Wang, 1985) are commonly used in linear models and generalized linear models, to decide whether a model with response y and predictors x would be improved by including another predictor z.

In a (generalised) linear model with response y and predictors x, the Added Variable Plot for a new covariate z is a plot of the smoothed Pearson residuals from the original model against the scaled residuals from a weighted linear regression of z on x. If this plot has nonzero slope, then the new covariate z is needed. For general advice see Cook and Weisberg(1999); Harrell (2001).

Essentially the same technique can be used for a spatial point process model (Baddeley et al, 2012).

The argument model should be a fitted spatial point process model (object of class "ppm").

The argument covariate identifies the covariate that is to be considered for addition to the model. It should be either a pixel image (object of class "im") or a function(x,y) giving the values of the covariate at any spatial location. Alternatively covariate may be a character string, giving the name of a covariate that was supplied (in the covariates argument to ppm) when the model was fitted, but was not used in the model.

The result of addvar(model, covariate) is an object belonging to the classes "addvar" and "fv". Plot this object to generate the added variable plot.

Note that the plot method shows the pointwise significance bands for a test of the null model, i.e. the null hypothesis that the new covariate has no effect.

The smoothing bandwidth is controlled by the arguments bw, adjust, bw.input and bw.restrict. If bw is a numeric value, then the bandwidth is taken to be adjust * bw. If bw is a string representing a bandwidth selection rule (recognised by density.default) then the bandwidth is selected by this rule.

The data used for automatic bandwidth selection are specified by bw.input and bw.restrict. If bw.input="points" (the default) then bandwidth selection is based on the covariate values at the points of the original point pattern dataset to which the model was fitted. If bw.input="quad" then bandwidth selection is based on the covariate values at every quadrature point used to fit the model. If bw.restrict=TRUE then the bandwidth selection is performed using only data from inside the subregion.

Value

An object of class "addvar" containing the coordinates for the added variable plot. There is a plot method.

Slow computation

In a large dataset, computation can be very slow if the default settings are used, because the smoothing bandwidth is selected automatically. To avoid this, specify a numerical value for the bandwidth bw. One strategy is to use a coarser subset of the data to select bw automatically. The selected bandwidth can be read off the print output for addvar.

Internal data

The return value has an attribute "spatial" which contains the internal data: the computed values of the residuals, and of all relevant covariates, at each quadrature point of the model. It is an object of class "ppp" with a data frame of marks.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net, Ya-Mei Chang and Yong Song.

References

Baddeley, A., Chang, Y.-M., Song, Y. and Turner, R. (2013) Residual diagnostics for covariate effects in spatial point process models. Journal of Computational and Graphical Statistics, 22, 886–905.

Cook, R.D. and Weisberg, S. (1999) Applied regression, including computing and graphics. New York: Wiley.

Cox, D.R. (1958) Planning of Experiments. New York: Wiley.

Harrell, F. (2001) Regression Modeling Strategies. New York: Springer.

Wang, P. (1985) Adding a variable in generalized linear models. Technometrics 27, 273–276.

See Also

parres, rhohat, rho2hat.

Examples

  X <-  rpoispp(function(x,y){exp(3+3*x)})
  model <- ppm(X, ~y)
  adv <- addvar(model, "x")
  plot(adv)
  adv <- addvar(model, "x", subregion=square(0.5))

ANOVA for Fitted Point Process Models for Replicated Patterns

Description

Performs analysis of deviance for one or more point process models fitted to replicated point pattern data.

Usage

  ## S3 method for class 'mppm'
anova(object, ...,
                  test=NULL, adjust=TRUE,
                  fine=FALSE, warn=TRUE)

Arguments

Details

This is a method for anova for comparing several fitted point process models of class "mppm", usually generated by the model-fitting function mppm).

If the fitted models are all Poisson point processes, then this function performs an Analysis of Deviance of the fitted models. The output shows the deviance differences (i.e. 2 times log likelihood ratio), the difference in degrees of freedom, and (if test="Chi") the two-sided p-values for the chi-squared tests. Their interpretation is very similar to that in anova.glm.

If some of the fitted models are not Poisson point processes, the ‘deviance’ differences in this table are 'pseudo-deviances' equal to 2 times the differences in the maximised values of the log pseudolikelihood (see ppm). It is not valid to compare these values to the chi-squared distribution. In this case, if adjust=TRUE (the default), the pseudo-deviances will be adjusted using the method of Pace et al (2011) and Baddeley, Turner and Rubak (2015) so that the chi-squared test is valid. It is strongly advisable to perform this adjustment.

The argument test determines which hypothesis test, if any, will be performed to compare the models. The argument test should be a character string, partially matching one of "Chisq", "F" or "Cp", or NULL. The first option "Chisq" gives the likelihood ratio test based on the asymptotic chi-squared distribution of the deviance difference. The meaning of the other options is explained in anova.glm.

Value

An object of class "anova", or NULL.

Random effects models are currently not supported

For models with random effects (i.e. where the call to mppm included the argument random), analysis of deviance is currently not supported, due to changes in the nlme package. We will try to find a solution.

Error messages

An error message that reports system is computationally singular indicates that the determinant of the Fisher information matrix of one of the models was either too large or too small for reliable numerical calculation. See vcov.ppm for suggestions on how to handle this.

Author(s)

Adrian Baddeley, Ida-Maria Sintorn and Leanne Bischoff. Implemented by Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

References

Baddeley, A., Rubak, E. and Turner, R. (2015) Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC Press.

Baddeley, A., Turner, R. and Rubak, E. (2015) Adjusted composite likelihood ratio test for Gibbs point processes. Journal of Statistical Computation and Simulation 86 (5) 922–941. DOI: 10.1080/00949655.2015.1044530.

Pace, L., Salvan, A. and Sartori, N. (2011) Adjusting composite likelihood ratio statistics. Statistica Sinica 21, 129–148.

See Also

mppm

Examples

 H <- hyperframe(X=waterstriders)
 #' test for loglinear trend in x coordinate
 mod0 <- mppm(X~1, data=H, Poisson())
 modx <- mppm(X~x, data=H, Poisson())
 anova(mod0, modx, test="Chi")
 # not significant
 anova(modx, test="Chi")
 # not significant

 #' test for inhibition
 mod0S <- mppm(X~1, data=H, Strauss(2))
 anova(mod0, mod0S, test="Chi")
 # significant! 

 #' test for trend after accounting for inhibition
 modxS <- mppm(X~x, data=H, Strauss(2))
 anova(mod0S, modxS, test="Chi")
 # not significant

ANOVA for Fitted Point Process Models

Description

Performs analysis of deviance for one or more fitted point process models.

Usage

  ## S3 method for class 'ppm'
anova(object, ..., test=NULL,
                      adjust=TRUE, warn=TRUE, fine=FALSE)

Arguments

Details

This is a method for anova for fitted point process models (objects of class "ppm", usually generated by the model-fitting function ppm).

If the fitted models are all Poisson point processes, then by default, this function performs an Analysis of Deviance of the fitted models. The output shows the deviance differences (i.e. 2 times log likelihood ratio), the difference in degrees of freedom, and (if test="Chi" or test="LRT") the two-sided p-values for the chi-squared tests. Their interpretation is very similar to that in anova.glm. If test="Rao" or test="score", the score test (Rao, 1948) is performed instead.

If some of the fitted models are not Poisson point processes, the ‘deviance’ differences in this table are 'pseudo-deviances' equal to 2 times the differences in the maximised values of the log pseudolikelihood (see ppm). It is not valid to compare these values to the chi-squared distribution. In this case, if adjust=TRUE (the default), the pseudo-deviances will be adjusted using the method of Pace et al (2011) and Baddeley et al (2015) so that the chi-squared test is valid. It is strongly advisable to perform this adjustment.

Value

An object of class "anova", or NULL.

Errors and warnings

models not nested:

There may be an error message that the models are not “nested”. For an Analysis of Deviance the models must be nested, i.e. one model must be a special case of the other. For example the point process model with formula ~x is a special case of the model with formula ~x+y, so these models are nested. However the two point process models with formulae ~x and ~y are not nested.

If you get this error message and you believe that the models should be nested, the problem may be the inability of R to recognise that the two formulae are nested. Try modifying the formulae to make their relationship more obvious.

different sizes of dataset:

There may be an error message from anova.glmlist that “models were not all fitted to the same size of dataset”. This implies that the models were fitted using different quadrature schemes (see quadscheme) and/or with different edge corrections or different values of the border edge correction distance rbord.

To ensure that models are comparable, check the following:

• the models must all have been fitted to the same point pattern dataset, in the same window.

• all models must have been fitted by the same fitting method as specified by the argument method in ppm.

• If some of the models depend on covariates, then they should all have been fitted using the same list of covariates, and using allcovar=TRUE to ensure that the same quadrature scheme is used.

• all models must have been fitted using the same edge correction as specified by the arguments correction and rbord. If you did not specify the value of rbord, then it may have taken a different value for different models. The default value of rbord is equal to zero for a Poisson model, and otherwise equals the reach (interaction distance) of the interaction term (see reach). To ensure that the models are comparable, set rbord to equal the maximum reach of the interactions that you are fitting.

Error messages

An error message that reports system is computationally singular indicates that the determinant of the Fisher information matrix of one of the models was either too large or too small for reliable numerical calculation. See vcov.ppm for suggestions on how to handle this.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

References

Baddeley, A., Turner, R. and Rubak, E. (2015) Adjusted composite likelihood ratio test for Gibbs point processes. Journal of Statistical Computation and Simulation 86 (5) 922–941. DOI: 10.1080/00949655.2015.1044530.

Pace, L., Salvan, A. and Sartori, N. (2011) Adjusting composite likelihood ratio statistics. Statistica Sinica 21, 129–148.

Rao, C.R. (1948) Large sample tests of statistical hypotheses concerning several parameters with applications to problems of estimation. Proceedings of the Cambridge Philosophical Society 44, 50–57.

See Also

ppm, vcov.ppm

Examples

 mod0 <- ppm(swedishpines ~1)
 modx <- ppm(swedishpines ~x)
 # Likelihood ratio test
 anova(mod0, modx, test="Chi")
 # Score test
 anova(mod0, modx, test="Rao")

 # Single argument
 modxy <- ppm(swedishpines ~x + y)
 anova(modxy, test="Chi")

 # Adjusted composite likelihood ratio test
 modP <- ppm(swedishpines ~1, rbord=9)
 modS <- ppm(swedishpines ~1, Strauss(9))
 anova(modP, modS, test="Chi")

Analysis of Deviance for Spatial Logistic Regression Models

Description

Performs Analysis of Deviance for two or more fitted Spatial Logistic Regression models.

Usage

  ## S3 method for class 'slrm'
anova(object, ..., test = NULL)

Arguments

Details

This is a method for anova for fitted spatial logistic regression models (objects of class "slrm", usually obtained from the function slrm).

The output shows the deviance differences (i.e. 2 times log likelihood ratio), the difference in degrees of freedom, and (if test="Chi") the two-sided p-values for the chi-squared tests. Their interpretation is very similar to that in anova.glm.

Value

An object of class "anova", inheriting from class "data.frame", representing the analysis of deviance table.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

See Also

slrm

Examples

  X <- rpoispp(42)
  fit0 <- slrm(X ~ 1)
  fit1 <- slrm(X ~ x+y)
  anova(fit0, fit1, test="Chi")

Convert Leverage Object to Function of Coordinates

Description

Converts an object of class "leverage.ppm" to a function of the x and y coordinates.

Usage

 ## S3 method for class 'leverage.ppm'
as.function(x, ...)

Arguments

Details

An object of class "leverage.ppm" represents the leverage function of a fitted point process model. This command converts the object to a function(x,y) where the arguments x and y are (vectors of) spatial coordinates. This function returns the leverage values at the specified locations (calculated by referring to the nearest location where the leverage has been computed).

Value

A function in the R language, also belonging to the class "funxy".

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

as.im.leverage.ppm

Examples

  X <- rpoispp(function(x,y) { exp(3+3*x) })
  fit <- ppm(X ~x+y)
  lev <- leverage(fit)
  f <- as.function(lev)
  
  f(0.2, 0.3)  # evaluate at (x,y) coordinates
  y <- f(X)    # evaluate at a point pattern

Convert Fitted Model To Class fv

Description

Converts fitted model into a function table (an object of class "fv").

Usage

  ## S3 method for class 'kppm'
as.fv(x)

  ## S3 method for class 'dppm'
as.fv(x)

  ## S3 method for class 'minconfit'
as.fv(x)

Arguments

Details

The generic command as.fv converts data x, that could be interpreted as the values of a function, into a function value table (object of the class "fv" as described in fv.object). This object can then be plotted easily using plot.fv.

Objects of class "kppm" (and related classes) represent a model that has been fitted to a dataset by computing a summary function of the dataset and matching it to the corresponding summary function of the model. The methods for as.fv for classes "kppm", "dppm" and "minconfit" extract this information: the result is a function table containing the observed summary function and the best fit summary function.

Value

An object of class "fv" (see fv.object).

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk

Examples

  as.fv(kppm(redwood))

Extract Interaction Structure

Description

Extracts the interpoint interaction structure from a point pattern model.

Usage

as.interact(object)
## S3 method for class 'fii'
as.interact(object)
## S3 method for class 'interact'
as.interact(object)
## S3 method for class 'ppm'
as.interact(object)

Arguments

Details

The function as.interact extracts the interpoint interaction structure from a suitable object.

An object of class "interact" describes an interpoint interaction structure, before it has been fitted to point pattern data. The irregular parameters of the interaction (such as the interaction range) are fixed, but the regular parameters (such as interaction strength) are undetermined. Objects of this class are created by the functions Poisson, Strauss and so on. The main use of such objects is in a call to ppm.

The function as.interact is generic, with methods for the classes "ppm", "fii" and "interact". The result is an object of class "interact" which can be printed.

Value

An object of class "interact" representing the interpoint interaction. This object can be printed and plotted.

Note on parameters

This function does not extract the fitted coefficients of the interaction. To extract the fitted interaction including the fitted coefficients, use fitin.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

See Also

fitin, ppm.

Examples

   model <- ppm(cells ~1, Strauss(0.07))
   f <- as.interact(model)
   f

Convert Measure To Layered Object

Description

Converts a measure into a layered object.

Usage

 ## S3 method for class 'msr'
as.layered(X)

Arguments

Details

This function converts the object X into an object of class "layered".

It is a method for the generic as.layered for the class of measures.

If X is a vector-valued measure, then as.layered(X) consists of several layers, each containing a scalar-valued measure.

Value

An object of class "layered" (see layered).

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

as.layered, msr.

Examples

   P <- rpoispp(100)
   fit <- ppm(P ~ x+y)
   rs <- residuals(fit, type="score")
   as.layered(rs)

Convert Data To Class owin

Description

Converts data specifying an observation window in any of several formats, into an object of class "owin".

Usage

 ## S3 method for class 'ppm'
as.owin(W, ..., from=c("points", "covariates"), fatal=TRUE)

 ## S3 method for class 'kppm'
as.owin(W, ..., from=c("points", "covariates"), fatal=TRUE)

 ## S3 method for class 'dppm'
as.owin(W, ..., from=c("points", "covariates"), fatal=TRUE)

 ## S3 method for class 'slrm'
as.owin(W, ..., from=c("points", "covariates"))

 ## S3 method for class 'msr'
as.owin(W, ..., fatal=TRUE)

Arguments

Details

The class "owin" is a way of specifying the observation window for a point pattern. See owin.object for an overview.

The generic function as.owin converts data in any of several formats into an object of class "owin" for use by the spatstat package. The function as.owin is generic, with methods for different classes of objects, and a default method.

The argument W may be

• an object of class "owin"

• a structure with entries xrange, yrange specifying the x and y dimensions of a rectangle

• a structure with entries named xmin, xmax, ymin, ymax (in any order) specifying the x and y dimensions of a rectangle. This will accept objects of class bbox in the sf package.

• a numeric vector of length 4 (interpreted as (xmin, xmax, ymin, ymax) in that order) specifying the x and y dimensions of a rectangle

• a structure with entries named xl, xu, yl, yu (in any order) specifying the x and y dimensions of a rectangle as (xmin, xmax) = (xl, xu) and (ymin, ymax) = (yl, yu). This will accept objects of class spp used in the Venables and Ripley spatial package.

• an object of class "ppp" representing a point pattern. In this case, the object's window structure will be extracted.

• an object of class "psp" representing a line segment pattern. In this case, the object's window structure will be extracted.

• an object of class "tess" representing a tessellation. In this case, the object's window structure will be extracted.

• an object of class "quad" representing a quadrature scheme. In this case, the window of the data component will be extracted.

• an object of class "im" representing a pixel image. In this case, a window of type "mask" will be returned, with the same pixel raster coordinates as the image. An image pixel value of NA, signifying that the pixel lies outside the window, is transformed into the logical value FALSE, which is the corresponding convention for window masks.

• an object of class "ppm", "kppm", "slrm" or "dppm" representing a fitted point process model. In this case, if from="data" (the default), as.owin extracts the original point pattern data to which the model was fitted, and returns the observation window of this point pattern. If from="covariates" then as.owin extracts the covariate images to which the model was fitted, and returns a binary mask window that specifies the pixel locations.

• an object of class "lpp" representing a point pattern on a linear network. In this case, as.owin extracts the linear network and returns a window containing this network.

• an object of class "lppm" representing a fitted point process model on a linear network. In this case, as.owin extracts the linear network and returns a window containing this network.

• A data.frame with exactly three columns. Each row of the data frame corresponds to one pixel. Each row contains the x and y coordinates of a pixel, and a logical value indicating whether the pixel lies inside the window.

• A data.frame with exactly two columns. Each row of the data frame contains the x and y coordinates of a pixel that lies inside the window.

• an object of class "distfun", "nnfun" or "funxy" representing a function of spatial location, defined on a spatial domain. The spatial domain of the function will be extracted.

• an object of class "rmhmodel" representing a point process model that can be simulated using rmh. The window (spatial domain) of the model will be extracted. The window may be NULL in some circumstances (indicating that the simulation window has not yet been determined). This is not treated as an error, because the argument fatal defaults to FALSE for this method.

• an object of class "layered" representing a list of spatial objects. See layered. In this case, as.owin will be applied to each of the objects in the list, and the union of these windows will be returned.

• an object of some other suitable class from another package. For full details, see vignette('shapefiles').

If the argument W is not in one of these formats and cannot be converted to a window, then an error will be generated (if fatal=TRUE) or a value of NULL will be returned (if fatal=FALSE).

When W is a data frame, the argument step can be used to specify the pixel grid spacing; otherwise, the spacing will be guessed from the data.

Value

An object of class "owin" (see owin.object) specifying an observation window.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

as.owin, as.owin.rmhmodel, as.owin.lpp.

owin.object, owin.

Additional methods for as.owin may be provided by other packages outside the spatstat family.

Examples

  fit <- ppm(cells ~ 1)
  as.owin(fit)

Extract Fitted Point Process Model

Description

Extracts the fitted point process model from some kind of fitted model.

Usage

as.ppm(object)

## S3 method for class 'ppm'
as.ppm(object)

## S3 method for class 'profilepl'
as.ppm(object)

## S3 method for class 'kppm'
as.ppm(object)

## S3 method for class 'dppm'
as.ppm(object)

## S3 method for class 'rppm'
as.ppm(object)

Arguments

Details

The function as.ppm extracts the fitted point process model (of class "ppm") from a suitable object.

The function as.ppm is generic, with methods for the classes "ppm", "profilepl", "kppm", "dppm" and "rppm", and possibly for other classes.

For the class "profilepl" of models fitted by maximum profile pseudolikelihood, the method as.ppm.profilepl extracts the fitted point process model (with the optimal values of the irregular parameters).

For the class "kppm" of models fitted by minimum contrast (or Palm or composite likelihood) using Waagepetersen's two-step estimation procedure (see kppm), the method as.ppm.kppm extracts the Poisson point process model that is fitted in the first stage of the procedure.

The behaviour for the class "dppm" is analogous to the "kppm" case above.

For the class "rppm" of models fitted by recursive partitioning (regression trees), the method as.ppm.rppm extracts the corresponding loglinear model that is fitted in the first stage of the procedure (whose purpose is merely to identify and evaluate the explanatory variables).

Value

An object of class "ppm".

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

ppm, profilepl.

Examples

   # fit a model by profile maximum pseudolikelihood
   rvals <- data.frame(r=(1:10)/100)
   pfit <- profilepl(rvals, Strauss, cells, ~1)
   # extract the fitted model
   fit <- as.ppm(pfit)

Area Under ROC Curve

Description

Compute the AUC (area under the Receiver Operating Characteristic curve) for a fitted point process model.

Usage

## S3 method for class 'ppm'
auc(X, ..., subset=NULL)

## S3 method for class 'kppm'
auc(X, ..., subset=NULL)

## S3 method for class 'slrm'
auc(X, ..., subset=NULL)

Arguments

Details

This command computes the AUC, the area under the Receiver Operating Characteristic curve. The ROC itself is computed by roc.

For a fitted point process model X, the AUC measures the ability of the fitted model intensity to separate the spatial domain into areas of high and low density of points. Suppose \lambda(u) is the intensity function of the model. The AUC is the probability that \lambda(x_i) > \lambda(U). That is, AUC is the probability that a randomly-selected data point has higher predicted intensity than does a randomly-selected spatial location. The AUC is not a measure of the goodness-of-fit of the model (Lobo et al, 2007).

(For spatial logistic regression models (class "slrm") replace “intensity” by “probability of presence” in the text above.)

The algorithm also calculates the theoretically expected AUC value for this model, as described in Baddeley et al (2025).

Value

Numeric vector of length 2 giving the AUC value and the theoretically expected AUC value for this model.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ege Rubak rubak@math.aau.dk and Suman Rakshit Suman.Rakshit@curtin.edu.au.

References

Baddeley, A., Rubak, E., Rakshit, S. and Nair, G. (2025) ROC curves for spatial point patterns and presence-absence data. doi:10.48550/arXiv.2506.03414..

Lobo, J.M., Jimenez-Valverde, A. and Real, R. (2007) AUC: a misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography 17(2) 145–151.

Nam, B.-H. and D'Agostino, R. (2002) Discrimination index, the area under the ROC curve. Pages 267–279 in Huber-Carol, C., Balakrishnan, N., Nikulin, M.S. and Mesbah, M., Goodness-of-fit tests and model validity, Birkhauser, Basel.

See Also

roc, roc.ppm, youden.

Examples

  fit <- ppm(swedishpines ~ x+y)
  auc(fit)
  auc(swedishpines, "x")

Bias Correction for Fitted Model

Description

Applies a first-order bias correction to a fitted model.

Usage

  bc(fit, ...)

  ## S3 method for class 'ppm'
bc(fit, ..., nfine = 256)

Arguments

Details

This command applies the first order Newton-Raphson bias correction method of Baddeley and Turner (2014, sec 4.2) to a fitted model. The function bc is generic, with a method for fitted point process models of class "ppm".

A fine grid of locations, of dimensions nfine * nfine or nfine[2] * nfine[1], is created over the original window of the data, and the intensity or conditional intensity of the fitted model is calculated on this grid. The result is used to update the fitted model parameters once by a Newton-Raphson update.

This is only useful if the quadrature points used to fit the original model fit are coarser than the grid of points specified by nfine.

Value

A numeric vector, of the same length as coef(fit), giving updated values for the fitted model coefficients.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net.

References

Baddeley, A. and Turner, R. (2014) Bias correction for parameter estimates of spatial point process models. Journal of Statistical Computation and Simulation 84, 1621–1643. DOI: 10.1080/00949655.2012.755976

See Also

rex

Examples

  fit <- ppm(cells ~ x, Strauss(0.07))
  coef(fit)
  if(!interactive()) {
    bc(fit, nfine=64)
  } else {
    bc(fit)
  }

Berman's Tests for Point Process Model

Description

Tests the goodness-of-fit of a Poisson point process model using methods of Berman (1986).

Usage

## S3 method for class 'ppm'
berman.test(model, covariate,
                         which = c("Z1", "Z2"),
               alternative = c("two.sided", "less", "greater"), ...)

Arguments

Details

These functions perform a goodness-of-fit test of a Poisson point process model fitted to point pattern data. The observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same values under the model, are compared using either of two test statistics Z_1 and Z_2 proposed by Berman (1986). The Z_1 test is also known as the Lawson-Waller test.

The function berman.test is generic, with methods for point patterns ("ppp" or "lpp") and point process models ("ppm" or "lppm").

• If X is a point pattern dataset (object of class "ppp" or "lpp"), then berman.test(X, ...) performs a goodness-of-fit test of the uniform Poisson point process (Complete Spatial Randomness, CSR) for this dataset.

• If model is a fitted point process model (object of class "ppm" or "lppm") then berman.test(model, ...) performs a test of goodness-of-fit for this fitted model. In this case, model should be a Poisson point process.

The test is performed by comparing the observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same covariate under the model. Thus, you must nominate a spatial covariate for this test.

The argument covariate should be either a function(x,y) or a pixel image (object of class "im" containing the values of a spatial function. If covariate is an image, it should have numeric values, and its domain should cover the observation window of the model. If covariate is a function, it should expect two arguments x and y which are vectors of coordinates, and it should return a numeric vector of the same length as x and y.

First the original data point pattern is extracted from model. The values of the covariate at these data points are collected.

Next the values of the covariate at all locations in the observation window are evaluated. The point process intensity of the fitted model is also evaluated at all locations in the window.

• If which="Z1", the test statistic Z_1 is computed as follows. The sum S of the covariate values at all data points is evaluated. The predicted mean \mu and variance \sigma^2 of S are computed from the values of the covariate at all locations in the window. Then we compute Z_1 = (S-\mu)/\sigma. Closely-related tests were proposed independently by Waller et al (1993) and Lawson (1993) so this test is often termed the Lawson-Waller test in epidemiological literature.

• If which="Z2", the test statistic Z_2 is computed as follows. The values of the covariate at all locations in the observation window, weighted by the point process intensity, are compiled into a cumulative distribution function F. The probability integral transformation is then applied: the values of the covariate at the original data points are transformed by the predicted cumulative distribution function F into numbers between 0 and 1. If the model is correct, these numbers are i.i.d. uniform random numbers. The standardised sample mean of these numbers is the statistic Z_2.

In both cases the null distribution of the test statistic is the standard normal distribution, approximately.

The return value is an object of class "htest" containing the results of the hypothesis test. The print method for this class gives an informative summary of the test outcome.

Value

An object of class "htest" (hypothesis test) and also of class "bermantest", containing the results of the test. The return value can be plotted (by plot.bermantest) or printed to give an informative summary of the test.

Warning

The meaning of a one-sided test must be carefully scrutinised: see the printed output.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

References

Berman, M. (1986) Testing for spatial association between a point process and another stochastic process. Applied Statistics 35, 54–62.

Lawson, A.B. (1993) On the analysis of mortality events around a prespecified fixed point. Journal of the Royal Statistical Society, Series A 156 (3) 363–377.

Waller, L., Turnbull, B., Clark, L.C. and Nasca, P. (1992) Chronic Disease Surveillance and testing of clustering of disease and exposure: Application to leukaemia incidence and TCE-contaminated dumpsites in upstate New York. Environmetrics 3, 281–300.

See Also

cdf.test, quadrat.test, ppm

Examples

   # Berman's data
   X <- copper$SouthPoints
   L <- copper$SouthLines
   D <- distmap(L, eps=1)
   # test of fitted model
   fit <- ppm(X ~ x+y)
   berman.test(fit, D)

Fit the Neyman-Scott cluster process with Cauchy kernel

Description

Fits the Neyman-Scott Cluster point process with Cauchy kernel to a point pattern dataset by the Method of Minimum Contrast.

Usage

cauchy.estK(X, startpar=c(kappa=1,scale=1), lambda=NULL,
            q = 1/4, p = 2, rmin = NULL, rmax = NULL, ...)

Arguments

Details

This algorithm fits the Neyman-Scott cluster point process model with Cauchy kernel to a point pattern dataset by the Method of Minimum Contrast, using the K function.

The argument X can be either

a point pattern:

An object of class "ppp" representing a point pattern dataset. The K function of the point pattern will be computed using Kest, and the method of minimum contrast will be applied to this.

a summary statistic:

An object of class "fv" containing the values of a summary statistic, computed for a point pattern dataset. The summary statistic should be the K function, and this object should have been obtained by a call to Kest or one of its relatives.

The algorithm fits the Neyman-Scott cluster point process with Cauchy kernel to X, by finding the parameters of the Matern Cluster model which give the closest match between the theoretical K function of the Matern Cluster process and the observed K function. For a more detailed explanation of the Method of Minimum Contrast, see mincontrast.

The model is described in Jalilian et al (2013). It is a cluster process formed by taking a pattern of parent points, generated according to a Poisson process with intensity \kappa, and around each parent point, generating a random number of offspring points, such that the number of offspring of each parent is a Poisson random variable with mean \mu, and the locations of the offspring points of one parent follow a common distribution described in Jalilian et al (2013).

If the argument lambda is provided, then this is used as the value of the point process intensity \lambda. Otherwise, if X is a point pattern, then \lambda will be estimated from X. If X is a summary statistic and lambda is missing, then the intensity \lambda cannot be estimated, and the parameter \mu will be returned as NA.

The remaining arguments rmin,rmax,q,p control the method of minimum contrast; see mincontrast.

The corresponding model can be simulated using rCauchy.

For computational reasons, the optimisation procedure uses the parameter eta2, which is equivalent to 4 * scale^2 where scale is the scale parameter for the model as used in rCauchy.

Homogeneous or inhomogeneous Neyman-Scott/Cauchy models can also be fitted using the function kppm and the fitted models can be simulated using simulate.kppm.

The optimisation algorithm can be controlled through the additional arguments "..." which are passed to the optimisation function optim. For example, to constrain the parameter values to a certain range, use the argument method="L-BFGS-B" to select an optimisation algorithm that respects box constraints, and use the arguments lower and upper to specify (vectors of) minimum and maximum values for each parameter.

Value

An object of class "minconfit". There are methods for printing and plotting this object. It contains the following main components:

Author(s)

Abdollah Jalilian and Rasmus Waagepetersen. Adapted for spatstat by Adrian Baddeley Adrian.Baddeley@curtin.edu.au

References

Ghorbani, M. (2013) Cauchy cluster process. Metrika 76, 697–706.

Jalilian, A., Guan, Y. and Waagepetersen, R. (2013) Decomposition of variance for spatial Cox processes. Scandinavian Journal of Statistics 40, 119-137.

Waagepetersen, R. (2007) An estimating function approach to inference for inhomogeneous Neyman-Scott processes. Biometrics 63, 252–258.

See Also

kppm, cauchy.estpcf, lgcp.estK, thomas.estK, vargamma.estK, mincontrast, Kest, Kmodel.

rCauchy to simulate the model.

Examples

    u <- cauchy.estK(redwood)
    u
    plot(u)

Fit the Neyman-Scott cluster process with Cauchy kernel

Description

Fits the Neyman-Scott Cluster point process with Cauchy kernel to a point pattern dataset by the Method of Minimum Contrast, using the pair correlation function.

Usage

cauchy.estpcf(X, startpar=c(kappa=1,scale=1), lambda=NULL,
            q = 1/4, p = 2, rmin = NULL, rmax = NULL, ...,
            pcfargs = list())

Arguments

Details

This algorithm fits the Neyman-Scott cluster point process model with Cauchy kernel to a point pattern dataset by the Method of Minimum Contrast, using the pair correlation function.

The argument X can be either

a point pattern:

An object of class "ppp" representing a point pattern dataset. The pair correlation function of the point pattern will be computed using pcf, and the method of minimum contrast will be applied to this.

a summary statistic:

An object of class "fv" containing the values of a summary statistic, computed for a point pattern dataset. The summary statistic should be the pair correlation function, and this object should have been obtained by a call to pcf or one of its relatives.

The algorithm fits the Neyman-Scott cluster point process with Cauchy kernel to X, by finding the parameters of the Matern Cluster model which give the closest match between the theoretical pair correlation function of the Matern Cluster process and the observed pair correlation function. For a more detailed explanation of the Method of Minimum Contrast, see mincontrast.

The model is described in Jalilian et al (2013). It is a cluster process formed by taking a pattern of parent points, generated according to a Poisson process with intensity \kappa, and around each parent point, generating a random number of offspring points, such that the number of offspring of each parent is a Poisson random variable with mean \mu, and the locations of the offspring points of one parent follow a common distribution described in Jalilian et al (2013).

If the argument lambda is provided, then this is used as the value of the point process intensity \lambda. Otherwise, if X is a point pattern, then \lambda will be estimated from X. If X is a summary statistic and lambda is missing, then the intensity \lambda cannot be estimated, and the parameter \mu will be returned as NA.

The remaining arguments rmin,rmax,q,p control the method of minimum contrast; see mincontrast.

The corresponding model can be simulated using rCauchy.

For computational reasons, the optimisation procedure internally uses the parameter eta2, which is equivalent to 4 * scale^2 where scale is the scale parameter for the model as used in rCauchy.

Homogeneous or inhomogeneous Neyman-Scott/Cauchy models can also be fitted using the function kppm and the fitted models can be simulated using simulate.kppm.

The optimisation algorithm can be controlled through the additional arguments "..." which are passed to the optimisation function optim. For example, to constrain the parameter values to a certain range, use the argument method="L-BFGS-B" to select an optimisation algorithm that respects box constraints, and use the arguments lower and upper to specify (vectors of) minimum and maximum values for each parameter.

Value

An object of class "minconfit". There are methods for printing and plotting this object. It contains the following main components:

Author(s)

Abdollah Jalilian and Rasmus Waagepetersen. Adapted for spatstat by Adrian Baddeley Adrian.Baddeley@curtin.edu.au

References

Ghorbani, M. (2013) Cauchy cluster process. Metrika 76, 697–706.

Jalilian, A., Guan, Y. and Waagepetersen, R. (2013) Decomposition of variance for spatial Cox processes. Scandinavian Journal of Statistics 40, 119-137.

Waagepetersen, R. (2007) An estimating function approach to inference for inhomogeneous Neyman-Scott processes. Biometrics 63, 252–258.

See Also

kppm, cauchy.estK, lgcp.estpcf, thomas.estpcf, vargamma.estpcf, mincontrast, pcf, pcfmodel.

rCauchy to simulate the model.

Examples

    u <- cauchy.estpcf(redwood)
    u
    plot(u, legendpos="topright")

Spatial Distribution Test for Multiple Point Process Model

Description

Performs a spatial distribution test of a point process model fitted to multiple spatial point patterns. The test compares the observed and predicted distributions of the values of a spatial covariate, using either the Kolmogorov-Smirnov, Cramer-von Mises or Anderson-Darling test of goodness-of-fit.

Usage

## S3 method for class 'mppm'
cdf.test(model, covariate, test=c("ks", "cvm", "ad"), ...,
            nsim=19, verbose=TRUE, interpolate=FALSE, fast=TRUE, jitter=TRUE)

Arguments

Details

This function is a method for the generic function cdf.test for the class mppm.

This function performs a goodness-of-fit test of a point process model that has been fitted to multiple point patterns. The observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same values under the model, are compared using the Kolmogorov-Smirnov, Cramer-von Mises or Anderson-Darling test of goodness-of-fit. These are exact tests if the model is Poisson; otherwise, for a Gibbs model, a Monte Carlo p-value is computed by generating simulated realisations of the model and applying the selected goodness-of-fit test to each simulation.

The argument model should be a fitted point process model fitted to multiple point patterns (object of class "mppm").

The argument covariate contains the values of a spatial function. It can be

• a function(x,y)

• a pixel image (object of class "im"

• a list of function(x,y), one for each point pattern

• a list of pixel images, one for each point pattern

• a hyperframe (see hyperframe) of which the first column will be taken as containing the covariate

• a character string giving the name of one of the covariates in model

• one of the character strings "x" or "y", indicating the spatial coordinates.

If covariate is an image, it should have numeric values, and its domain should cover the observation window of the model. If covariate is a function, it should expect two arguments x and y which are vectors of coordinates, and it should return a numeric vector of the same length as x and y.

First the original data point pattern is extracted from model. The values of the covariate at these data points are collected.

The predicted distribution of the values of the covariate under the fitted model is computed as follows. The values of the covariate at all locations in the observation window are evaluated, weighted according to the point process intensity of the fitted model, and compiled into a cumulative distribution function F using ewcdf.

The probability integral transformation is then applied: the values of the covariate at the original data points are transformed by the predicted cumulative distribution function F into numbers between 0 and 1. If the model is correct, these numbers are i.i.d. uniform random numbers. A goodness-of-fit test of the uniform distribution is applied to these numbers using ks.test, cvm.test or ad.test.

The argument interpolate determines how pixel values will be handled when covariate is a pixel image. The value of the covariate at a data point is obtained by looking up the value of the nearest pixel if interpolate=FALSE, or by linearly interpolating between the values of the four nearest pixels if interpolate=TRUE. Linear interpolation is slower, but is sometimes necessary to avoid tied values of the covariate arising when the pixel grid is coarse.

If model is a Poisson point process, then the Kolmogorov-Smirnov, Cramer-von Mises and Anderson-Darling tests are theoretically exact. This test was apparently first described (in the context of spatial data, and for Kolmogorov-Smirnov) by Berman (1986). See also Baddeley et al (2005).

If model is not a Poisson point process, then the Kolmogorov-Smirnov, Cramer-von Mises and Anderson-Darling tests are biased. Instead they are used as the basis of a Monte Carlo test. First nsim simulated realisations of the model will be generated. Each simulated realisation consists of a list of simulated point patterns, one for each of the original data patterns. This can take a very long time. The model is then re-fitted to each simulation, and the refitted model is subjected to the goodness-of-fit test described above. A Monte Carlo p-value is then computed by comparing the p-value of the original test with the p-values obtained from the simulations.

Value

An object of class "cdftest" and "htest" containing the results of the test. See cdf.test for details.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ida-Maria Sintorn and Leanne Bischoff. Implemented by Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

References

Baddeley, A., Rubak, E. and Turner, R. (2015) Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC Press.

Baddeley, A., Turner, R., Moller, J. and Hazelton, M. (2005) Residual analysis for spatial point processes. Journal of the Royal Statistical Society, Series B 67, 617–666.

Berman, M. (1986) Testing for spatial association between a point process and another stochastic process. Applied Statistics 35, 54–62.

See Also

cdf.test, quadrat.test, mppm

Examples

   # three i.i.d. realisations of nonuniform Poisson process
   lambda <- as.im(function(x,y) { 200 * exp(x) }, square(1))
   dat <- hyperframe(X=list(rpoispp(lambda), rpoispp(lambda), rpoispp(lambda)))

   # fit uniform Poisson process
   fit0 <- mppm(X~1, dat)
   # fit correct nonuniform Poisson process
   fit1 <- mppm(X~x, dat)

   # test wrong model
   cdf.test(fit0, "x")
   # test right model
   cdf.test(fit1, "x")

   # Gibbs model
   fitGibbs <- update(fit0, interaction=Strauss(0.05))
   ns <- if(interactive()) 19 else 2
   cdf.test(fitGibbs, "x", nsim=ns)

Spatial Distribution Test for Point Pattern or Point Process Model

Description

Performs a test of goodness-of-fit of a point process model. The observed and predicted distributions of the values of a spatial covariate are compared using either the Kolmogorov-Smirnov test, Cramer-von Mises test or Anderson-Darling test. For non-Poisson models, a Monte Carlo test is used.

Usage

## S3 method for class 'ppm'
cdf.test(model, covariate,  test=c("ks", "cvm", "ad"), ...,
          interpolate=TRUE, jitter=TRUE, nsim=99, verbose=TRUE)

## S3 method for class 'slrm'
cdf.test(model, covariate,  test=c("ks", "cvm", "ad"), ...,
          modelname=NULL, covname=NULL)

Arguments

Details

These functions perform a goodness-of-fit test of a Poisson or Gibbs point process model fitted to point pattern data. The observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same values under the model, are compared using the Kolmogorov-Smirnov test, the Cramer-von Mises test or the Anderson-Darling test. For Gibbs models, a Monte Carlo test is performed using these test statistics.

The function cdf.test is generic, with methods for point patterns ("ppp" or "lpp"), point process models ("ppm" or "lppm") and spatial logistic regression models ("slrm").

• If X is a point pattern dataset (object of class "ppp"), then cdf.test(X, ...) performs a goodness-of-fit test of the uniform Poisson point process (Complete Spatial Randomness, CSR) for this dataset. For a multitype point pattern, the uniform intensity is assumed to depend on the type of point (sometimes called Complete Spatial Randomness and Independence, CSRI).

• If model is a fitted point process model (object of class "ppm" or "lppm") then cdf.test(model, ...) performs a test of goodness-of-fit for this fitted model.

• If model is a fitted spatial logistic regression (object of class "slrm") then cdf.test(model, ...) performs a test of goodness-of-fit for this fitted model.

The test is performed by comparing the observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same covariate under the model, using a classical goodness-of-fit test. Thus, you must nominate a spatial covariate for this test.

If X is a point pattern that does not have marks, the argument covariate should be either a function(x,y) or a pixel image (object of class "im" containing the values of a spatial function, or one of the characters "x" or "y" indicating the Cartesian coordinates. If covariate is an image, it should have numeric values, and its domain should cover the observation window of the model. If covariate is a function, it should expect two arguments x and y which are vectors of coordinates, and it should return a numeric vector of the same length as x and y.

If X is a multitype point pattern, the argument covariate can be either a function(x,y,marks), or a pixel image, or a list of pixel images corresponding to each possible mark value, or one of the characters "x" or "y" indicating the Cartesian coordinates.

First the original data point pattern is extracted from model. The values of the covariate at these data points are collected.

The predicted distribution of the values of the covariate under the fitted model is computed as follows. The values of the covariate at all locations in the observation window are evaluated, weighted according to the point process intensity of the fitted model, and compiled into a cumulative distribution function F using ewcdf.

The probability integral transformation is then applied: the values of the covariate at the original data points are transformed by the predicted cumulative distribution function F into numbers between 0 and 1. If the model is correct, these numbers are i.i.d. uniform random numbers. The A goodness-of-fit test of the uniform distribution is applied to these numbers using stats::ks.test, goftest::cvm.test or goftest::ad.test.

This test was apparently first described (in the context of spatial data, and using Kolmogorov-Smirnov) by Berman (1986). See also Baddeley et al (2005).

If model is not a Poisson process, then a Monte Carlo test is performed, by generating nsim point patterns which are simulated realisations of the model, re-fitting the model to each simulated point pattern, and calculating the test statistic for each fitted model. The Monte Carlo p value is determined by comparing the simulated values of the test statistic with the value for the original data.

The return value is an object of class "htest" containing the results of the hypothesis test. The print method for this class gives an informative summary of the test outcome.

The return value also belongs to the class "cdftest" for which there is a plot method plot.cdftest. The plot method displays the empirical cumulative distribution function of the covariate at the data points, and the predicted cumulative distribution function of the covariate under the model, plotted against the value of the covariate.

The argument jitter controls whether covariate values are randomly perturbed, in order to avoid ties. If the original data contains any ties in the covariate (i.e. points with equal values of the covariate), and if jitter=FALSE, then the Kolmogorov-Smirnov test implemented in ks.test will issue a warning that it cannot calculate the exact p-value. To avoid this, if jitter=TRUE each value of the covariate will be perturbed by adding a small random value. The perturbations are normally distributed with standard deviation equal to one hundredth of the range of values of the covariate. This prevents ties, and the p-value is still correct. There is a very slight loss of power.

Value

An object of class "htest" containing the results of the test. See ks.test for details. The return value can be printed to give an informative summary of the test.

The value also belongs to the class "cdftest" for which there is a plot method.

Warning

The outcome of the test involves a small amount of random variability, because (by default) the coordinates are randomly perturbed to avoid tied values. Hence, if cdf.test is executed twice, the p-values will not be exactly the same. To avoid this behaviour, set jitter=FALSE.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

References

Baddeley, A., Turner, R., Moller, J. and Hazelton, M. (2005) Residual analysis for spatial point processes. Journal of the Royal Statistical Society, Series B 67, 617–666.

Berman, M. (1986) Testing for spatial association between a point process and another stochastic process. Applied Statistics 35, 54–62.

See Also

plot.cdftest, quadrat.test, berman.test, ks.test, cvm.test, ad.test, ppm

Examples

   op <- options(useFancyQuotes=FALSE)


   # fit inhomogeneous Poisson model and test
   model <- ppm(nztrees ~x)
   cdf.test(model, "x")

   if(interactive()) {
     # synthetic data: nonuniform Poisson process
     X <- rpoispp(function(x,y) { 100 * exp(x) }, win=square(1))

     # fit uniform Poisson process
     fit0 <- ppm(X ~1)
     # fit correct nonuniform Poisson process
     fit1 <- ppm(X ~x)

     # test wrong model
     cdf.test(fit0, "x")
     # test right model
     cdf.test(fit1, "x")
   }

   # multitype point pattern
   yimage <- as.im(function(x,y){y}, W=Window(amacrine))
   cdf.test(ppm(amacrine ~marks+y), yimage)

   options(op)

Count Close Pairs of Points

Description

Low-level functions to count the number of close pairs of points.

Usage

closepaircounts(X, r)

crosspaircounts(X, Y, r)

Arguments

Details

These are the efficient low-level functions used by spatstat to count close pairs of points in a point pattern or between two point patterns.

closepaircounts(X,r) counts the number of neighbours for each point in the pattern X. That is, for each point X[i], it counts the number of other points X[j] with j != i such that d(X[i],X[j]) <= r where d denotes Euclidean distance. The result is an integer vector v such that v[i] is the number of neighbours of X[i].

crosspaircounts(X,Y,r) counts, for each point in the pattern X, the number of neighbours in the pattern Y. That is, for each point X[i], it counts the number of points Y[j] such that d(X[i],X[j]) <= r. The result is an integer vector v such that v[i] is the number of neighbours of X[i] in the pattern Y.

Value

An integer vector of length equal to the number of points in X.

Warning about accuracy

The results of these functions may not agree exactly with the correct answer (as calculated by a human) and may not be consistent between different computers and different installations of R. The discrepancies arise in marginal cases where the interpoint distance is equal to, or very close to, the threshold rmax.

Floating-point numbers in a computer are not mathematical Real Numbers: they are approximations using finite-precision binary arithmetic. The approximation is accurate to a tolerance of about .Machine$double.eps.

If the true interpoint distance d and the threshold rmax are equal, or if their difference is no more than .Machine$double.eps, the result may be incorrect.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

See Also

closepairs to identify all close pairs of points.

Examples

   a <- closepaircounts(cells, 0.1)
   sum(a)
   Y <- split(amacrine)
   b <- crosspaircounts(Y$on, Y$off, 0.1)

Field of clusters

Description

Calculate the superposition of cluster kernels at the location of a point pattern.

Usage

  ## S3 method for class 'kppm'
clusterfield(model, locations = NULL, ...)

Arguments

Details

The function clusterfield is generic, with a method for "kppm" (described here) and methods for "character" and "function".

The method clusterfield.kppm extracts the relevant information from the fitted model and calls clusterfield.function.

The calculations are performed by density.ppp and ... arguments are passed thereto for control over the pixel resolution etc. (These arguments are then passed on to pixellate.ppp and as.mask.)

Value

A pixel image (object of class "im").

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

kppm,

clusterfield.

Examples

  fit <- kppm(redwood~1, "Thomas")
  Z <- clusterfield(fit, eps = 0.01)

Fit Cluster or Cox Point Process Model via Minimum Contrast

Description

Fit a homogeneous or inhomogeneous cluster process or Cox point process model to a point pattern by the Method of Minimum Contrast.

Usage

clusterfit(X, clusters, lambda = NULL, startpar = NULL, ...,
           q = 1/4, p = 2, rmin = NULL, rmax = NULL,
           ctrl=list(q=q, p=p, rmin=rmin, rmax=rmax),
           statistic = NULL, statargs = NULL, algorithm="Nelder-Mead",
           verbose=FALSE, pspace=NULL)

Arguments

Details

This function fits the clustering parameters of a cluster or Cox point process model by the Method of Minimum Contrast, that is, by matching the theoretical K-function of the model to the empirical K-function of the data, as explained in mincontrast.

If statistic="pcf" (or X appears to be an estimated pair correlation function) then instead of using the K-function, the algorithm will use the pair correlation function.

If X is a point pattern of class "ppp" an estimate of the summary statistic specfied by statistic (defaults to "K") is first computed before minimum contrast estimation is carried out as described above. In this case the argument statargs can be used for controlling the summary statistic estimation. The precise algorithm for computing the summary statistic depends on whether the intensity specification (lambda) is:

homogeneous:

If lambda is NUll or a single numeric the pattern is considered homogeneous and either Kest or pcf is invoked. In this case lambda is not used for anything when estimating the summary statistic.

inhomogeneous:

If lambda is a pixel image (object of class "im"), a fitted point process model (object of class "ppm" or "kppm") or a function(x,y) the pattern is considered inhomogeneous. In this case either Kinhom or pcfinhom is invoked with lambda as an argument.

After the clustering parameters of the model have been estimated by minimum contrast lambda (if non-null) is used to compute the additional model parameter \mu.

The algorithm parameters q,p,rmax,rmin are described in the help for mincontrast. They may be provided either as individually-named arguments, or as entries in the list ctrl. The individually-named arguments q,p,rmax,rmin override the entries in the list ctrl.

Value

An object of class "minconfit". There are methods for printing and plotting this object. See mincontrast.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

References

Diggle, P.J. and Gratton, R.J. (1984) Monte Carlo methods of inference for implicit statistical models. Journal of the Royal Statistical Society, series B 46, 193 – 212.

Moller, J. and Waagepetersen, R. (2003). Statistical Inference and Simulation for Spatial Point Processes. Chapman and Hall/CRC, Boca Raton.

Waagepetersen, R. (2007). An estimating function approach to inference for inhomogeneous Neyman-Scott processes. Biometrics 63 (2007) 252–258.

See Also

kppm

Examples

  fit <- clusterfit(redwood, "Thomas")
  fit
  if(interactive()){
    plot(fit)
  }
  K <- Kest(redwood)
  fit2 <- clusterfit(K, "MatClust")

Extract Cluster Offspring Kernel

Description

Given a fitted cluster point process model, this command returns the probability density of the cluster offspring.

Usage

## S3 method for class 'kppm'
clusterkernel(model, ...)

Arguments

Details

Given a cluster point process model, this command returns a function(x,y) giving the two-dimensional probability density of the cluster offspring points assuming a cluster parent located at the origin.

The function clusterkernel is generic, with methods for class "kppm" (described here) and "character" (described in clusterkernel.character).

Value

A function in the R language with arguments x,y,....

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

clusterkernel.character, clusterfield, kppm

Examples

  fit <- kppm(redwood ~ x, "MatClust")
  f <- clusterkernel(fit)
  f(0.05, 0.02)

Compute or Extract Effective Range of Cluster Kernel

Description

Given a cluster point process model, this command returns a value beyond which the the probability density of the cluster offspring is neglible.

Usage

## S3 method for class 'kppm'
clusterradius(model, ..., thresh = NULL, precision = FALSE)

Arguments

Details

Given a cluster model this function by default returns the effective range of the model with the given parameters as used in spatstat. For the Matern cluster model (see e.g. rMatClust) this is simply the finite radius of the offsring density given by the paramter scale irrespective of other options given to this function. The remaining models in spatstat have infinite theoretical range, and an effective finite value is given as follows: For the Thomas model (see e.g. rThomas the default is 4*scale where scale is the scale or standard deviation parameter of the model. If thresh is given the value is instead found as described for the other models below.

For the Cauchy model (see e.g. rCauchy) and the Variance Gamma (Bessel) model (see e.g. rVarGamma) the value of thresh defaults to 0.001, and then this is used to compute the range numerically as follows. If k(x,y)=k_0(r) with r=\sqrt(x^2+y^2) denotes the isotropic cluster kernel then f(r) = 2 \pi r k_0(r) is the density function of the offspring distance from the parent. The range is determined as the value of r where f(r) falls below thresh times k_0(r).

If precision=TRUE the precision related to the chosen range is returned as an attribute. Here the precision is defined as the polar integral of the kernel from distance 0 to the calculated range. Ideally this should be close to the value 1 which would be obtained for the true theretical infinite range.

Value

A positive numeric.

Additionally, the precision related to this range value is returned as an attribute "prec", if precision=TRUE.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

See Also

clusterkernel, kppm, rMatClust, rThomas, rCauchy, rVarGamma, rNeymanScott.

Examples

  fit <- kppm(redwood ~ x, "MatClust")
  clusterradius(fit)

Coefficients of Point Process Model Fitted to Multiple Point Patterns

Description

Given a point process model fitted to a list of point patterns, extract the coefficients of the fitted model. A method for coef.

Usage

  ## S3 method for class 'mppm'
coef(object, ...)

Arguments

Details

This function is a method for the generic function coef.

The argument object must be a fitted point process model (object of class "mppm") produced by the fitting algorithm mppm). This represents a point process model that has been fitted to a list of several point pattern datasets. See mppm for information.

This function extracts the vector of coefficients of the fitted model. This is the estimate of the parameter vector \theta such that the conditional intensity of the model is of the form

\lambda(u,x) = \exp(\theta S(u,x))

where S(u,x) is a (vector-valued) statistic.

For example, if the model object is the uniform Poisson process, then coef(object) will yield a single value (named "(Intercept)") which is the logarithm of the fitted intensity of the Poisson process.

If the fitted model includes random effects (i.e. if the argument random was specified in the call to mppm), then the fitted coefficients are different for each point pattern in the original data, so coef(object) is a data frame with one row for each point pattern, and one column for each parameter. Use fixef.mppm to extract the vector of fixed effect coefficients, and ranef.mppm to extract the random effect coefficients at each level.

Use print.mppm to print a more useful description of the fitted model.

Value

Either a vector containing the fitted coefficients, or a data frame containing the fitted coefficients for each point pattern.

Author(s)

Adrian Baddeley, Ida-Maria Sintorn and Leanne Bischoff. Implemented in spatstat by Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Rolf Turner rolfturner@posteo.net and Ege Rubak rubak@math.aau.dk.

References

Baddeley, A., Rubak, E. and Turner, R. (2015) Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC Press.

See Also

fixef.mppm and ranef.mppm for the fixed and random effect coefficients in a model that includes random effects.

print.mppm, mppm

Examples

    H <- hyperframe(X=waterstriders)

    fit.Poisson <- mppm(X ~ 1, H)
    coef(fit.Poisson)

    # The single entry "(Intercept)" 
    # is the log of the fitted intensity of the Poisson process

    fit.Strauss <- mppm(X~1, H, Strauss(7))
    coef(fit.Strauss)

    # The two entries "(Intercept)" and "Interaction"
    # are respectively log(beta) and log(gamma)
    # in the usual notation for Strauss(beta, gamma, r)

    # Tweak data to exaggerate differences
    H$X[[1]] <- rthin(H$X[[1]], 0.3)
    # Model with random effects
    fitran <- mppm(X ~ 1, H, random=~1|id)
    coef(fitran)

Coefficients of Fitted Point Process Model

Description

Given a point process model fitted to a point pattern, extract the coefficients of the fitted model. A method for coef.

Usage

  ## S3 method for class 'ppm'
coef(object, ...)

Arguments

Details

This function is a method for the generic function coef.

The argument object must be a fitted point process model (object of class "ppm"). Such objects are produced by the maximum pseudolikelihood fitting algorithm ppm).

This function extracts the vector of coefficients of the fitted model. This is the estimate of the parameter vector \theta such that the conditional intensity of the model is of the form

\lambda(u,x) = \exp(\theta S(u,x))

where S(u,x) is a (vector-valued) statistic.

For example, if the model object is the uniform Poisson process, then coef(object) will yield a single value (named "(Intercept)") which is the logarithm of the fitted intensity of the Poisson process.

Use print.ppm to print a more useful description of the fitted model.

Value

A vector containing the fitted coefficients.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

See Also

print.ppm, ppm.object, ppm

Examples

    poi <- ppm(cells, ~1, Poisson())
    coef(poi)
    # This is the log of the fitted intensity of the Poisson process

    stra <- ppm(cells, ~1, Strauss(r=0.07))
    coef(stra)

    # The two entries "(Intercept)" and "Interaction"
    # are respectively log(beta) and log(gamma)
    # in the usual notation for Strauss(beta, gamma, r)

Coefficients of Fitted Spatial Logistic Regression Model

Description

Extracts the coefficients (parameters) from a fitted Spatial Logistic Regression model.

Usage

  ## S3 method for class 'slrm'
coef(object, ...)

Arguments

Details

This is a method for coef for fitted spatial logistic regression models (objects of class "slrm", usually obtained from the function slrm).

It extracts the fitted canonical parameters, i.e.\ the coefficients in the linear predictor of the spatial logistic regression.

Value

Numeric vector of coefficients.

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net

See Also

slrm

Examples

  X <- rpoispp(42)
  fit <- slrm(X ~ x+y)
  coef(fit)

Residual Diagnostics for Multiple Fitted Models

Description

Compares several fitted point process models using the same residual diagnostic.

Usage

compareFit(object, Fun, r = NULL, breaks = NULL, ...,
         trend = ~1, interaction = Poisson(), rbord = NULL,
         modelnames = NULL, same = NULL, different = NULL)

Arguments

Details

This is a convenient way to collect diagnostic information for several different point process models fitted to the same point pattern dataset, or for point process models of the same form fitted to several different datasets, etc.

The first argument, object, is usually a list of fitted point process models (objects of class "ppm"), obtained from the model-fitting function ppm.

For convenience, object can also be a list of point patterns (objects of class "ppp"). In that case, point process models will be fitted to each of the point pattern datasets, by calling ppm using the arguments trend (for the first order trend), interaction (for the interpoint interaction) and rbord (for the erosion distance in the border correction for the pseudolikelihood). See ppm for details of these arguments.

Alternatively object can be a single point pattern (object of class "ppp") and one or more of the arguments trend, interaction or rbord can be a list. In this case, point process models will be fitted to the same point pattern dataset, using each of the model specifications listed.

The diagnostic function Fun will be applied to each of the point process models. The results will be collected into a single function value table. The modelnames are used to label the results from each fitted model.

Value

Function value table (object of class "fv").

Author(s)

Ege Rubak rubak@math.aau.dk, Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Jesper Moller.

See Also

ppm, Kcom, Kres, Gcom, Gres, psst, psstA, psstG, collapse.fv

Examples

   nd <- 40
   
   ilist <- list(Poisson(), Geyer(7, 2), Strauss(7))
   iname <- c("Poisson", "Geyer", "Strauss")
   
   K <- compareFit(swedishpines, Kcom, interaction=ilist, rbord=9,
            correction="translate",
            same="trans", different="tcom", modelnames=iname, nd=nd)
   K

Extract Original Data from a Fitted Point Process Model

Description

Given a fitted point process model, this function extracts the original point pattern dataset to which the model was fitted.

Usage

  data.ppm(object)

Arguments

Details

An object of class "ppm" represents a point process model that has been fitted to data. It is typically produced by the model-fitting algorithm ppm. The object contains complete information about the original data point pattern to which the model was fitted. This function extracts the original data pattern.

See ppm.object for a list of all operations that can be performed on objects of class "ppm".

Value

A point pattern (object of class "ppp").

Author(s)

Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Rolf Turner rolfturner@posteo.net.

See Also

ppm.object, ppp.object

Examples

 fit <- ppm(cells, ~1, Strauss(r=0.1))
 X <- data.ppm(fit)
 # 'X' is identical to 'cells'

Construct a New Determinantal Point Process Model Family Function

Description

Function to ease the implementation of a new determinantal point process model family.

Usage

detpointprocfamilyfun(kernel = NULL,
    specden = NULL, basis = "fourierbasis", 
    convkernel = NULL, Kfun = NULL, valid = NULL, intensity = NULL, 
    dim = 2, name = "User-defined", isotropic = TRUE, range = NULL, 
    parbounds = NULL, specdenrange = NULL, startpar = NULL, ...)

Arguments