ppm
Fit Point Process Model to Data
Fits a point process model to an observed point pattern.
Usage
ppm(Q, ...) ## S3 method for class 'formula':
ppm(Q, interaction=NULL, \dots, data=NULL)
Arguments
- Q
- A
formula
in the Rlanguage describing the model to be fitted. - interaction
- An object of class
"interact"
describing the point process interaction structure, orNULL
indicating that a Poisson process (stationary or nonstationary) should be fitted. - ...
- Arguments passed to
ppm.ppp
orppm.quad
to control the model-fitting process. - data
- Optional. The values of spatial covariates (other than the Cartesian coordinates) required by the model. Either a data frame, or a list whose entries are images, functions, windows, tessellations or single numbers. See Details.
Details
This function fits a point process model to an observed point pattern. The model may include spatial trend, interpoint interaction, and dependence on covariates.
The function ppm
is generic, with methods for
the classes formula
, ppp
and quad
.
This page describes the method for a formula
.
The first argument is a formula
in the Rlanguage
describing the spatial trend model to be fitted. It has the general form
pattern ~ trend
where the left hand side pattern
is usually
the name of a spatial point pattern (object of class "ppp"
)
to which the model should be fitted, or an expression which evaluates
to a point pattern;
and the right hand side trend
is an expression specifying the
spatial trend of the model.
Systematic effects (spatial trend and/or dependence on
spatial covariates) are specified by the
trend
expression on the right hand side of the formula.
The trend may involve
the Cartesian coordinates x
, y
,
the marks marks
,
the names of entries in the argument data
(if supplied),
or the names of objects that exist in the Rsession.
The trend formula specifies the logarithm of the
intensity of a Poisson process, or in general, the logarithm of
the first order potential of the Gibbs process.
The formula should not use any names beginning with .mpl
as these are reserved for internal use.
If the formula is pattern~1
, then
the model to be fitted is stationary (or at least, its first order
potential is constant).
Stochastic interactions between random points of the point process
are defined by the argument interaction
. This is an object of
class "interact"
which is initialised in a very similar way to the
usage of family objects in glm
and gam
.
The models currently available are:
Poisson
,
AreaInter
,
BadGey
,
Concom
,
DiggleGatesStibbard
,
DiggleGratton
,
Fiksel
,
Geyer
,
Hardcore
,
LennardJones
,
MultiStrauss
,
MultiStraussHard
,
OrdThresh
,
Ord
,
Pairwise
,
PairPiece
,
Saturated
,
SatPiece
,
Softcore
,
Strauss
and
StraussHard
.
See the examples below.
It is also possible to combine several interactions
using Hybrid
.
If interaction
is missing or NULL
,
then the model to be fitted
has no interpoint interactions, that is, it is a Poisson process
(stationary or nonstationary according to trend
). In this case
the methods of maximum pseudolikelihood and maximum logistic likelihood
coincide with maximum likelihood.
The fitted point process model returned by this function can be printed
(by the print method print.ppm
)
to inspect the fitted parameter values.
If a nonparametric spatial trend was fitted, this can be extracted using
the predict method predict.ppm
.
To fit a model involving spatial covariates
other than the Cartesian coordinates $x$ and $y$,
the values of the covariates should either be supplied in the
argument data
, or should be stored in objects that exist
in the Rsession.
Note that it is not sufficient to have observed
the covariate only at the points of the data point pattern;
the covariate must also have been observed at other
locations in the window.
If it is given, the argument data
is typically
a list, with names corresponding to variables in the trend
formula.
Each entry in the list is either
[object Object],[object Object],[object Object],[object Object],[object Object]
The software will look up
the values of each covariate at the required locations
(quadrature points).
Note that, for covariate functions, only the name of the function appears in the trend formula. A covariate function is treated as if it were a single variable. The function arguments do not appear in the trend formula. See the Examples.
If data
is a list,
the list entries should have names corresponding to
(some of) the names of covariates in the model formula trend
.
The variable names x
, y
and marks
are reserved for the Cartesian
coordinates and the mark values,
and these should not be used for variables in data
.
Alternatively, data
may be a data frame
giving the values of the covariates at specified locations.
Then pattern
should be a quadrature scheme (object of class
"quad"
) giving the corresponding locations.
See ppm.quad
for details.
Value
- An object of class
"ppm"
describing a fitted point process model. Seeppm.object
for details of the format of this object and methods available for manipulating it.
References
Baddeley, A., Coeurjolly, J.-F., Rubak, E. and Waagepetersen, R. (2013)
A logistic regression estimating function
for spatial Gibbs point processes.
Research Report, Centre for Stochastic Geometry and Bioimaging,
Denmark.
Baddeley, A. and Turner, R. (2000) Practical maximum pseudolikelihood for spatial point patterns. Australian and New Zealand Journal of Statistics 42 283--322. Berman, M. and Turner, T.R. (1992) Approximating point process likelihoods with GLIM. Applied Statistics 41, 31--38. Besag, J. (1975) Statistical analysis of non-lattice data. The Statistician 24, 179-195. Diggle, P.J., Fiksel, T., Grabarnik, P., Ogata, Y., Stoyan, D. and Tanemura, M. (1994) On parameter estimation for pairwise interaction processes. International Statistical Review 62, 99-117.
Huang, F. and Ogata, Y. (1999) Improvements of the maximum pseudo-likelihood estimators in various spatial statistical models. Journal of Computational and Graphical Statistics 8, 510--530. Jensen, J.L. and Moeller, M. (1991) Pseudolikelihood for exponential family models of spatial point processes. Annals of Applied Probability 1, 445--461. Jensen, J.L. and Kuensch, H.R. (1994) On asymptotic normality of pseudo likelihood estimates for pairwise interaction processes, Annals of the Institute of Statistical Mathematics 46, 475--486.
See Also
ppm.ppp
and ppm.quad
for
more details on the fitting technique and edge correction.
ppm.object
for details of how to
print, plot and manipulate a fitted model.
ppp
and quadscheme
for constructing data.
Interactions:
Poisson
,
AreaInter
,
BadGey
,
DiggleGatesStibbard
,
DiggleGratton
,
Geyer
,
Fiksel
,
Hardcore
,
Hybrid
,
LennardJones
,
MultiStrauss
,
MultiStraussHard
,
OrdThresh
,
Ord
,
Pairwise
,
PairPiece
,
Saturated
,
SatPiece
,
Softcore
,
Strauss
and
StraussHard
.
See profilepl
for advice on
fitting nuisance parameters in the interaction,
and ippm
for irregular parameters in the trend.
See valid.ppm
and project.ppm
for
ensuring the fitted model is a valid point process.
Examples
# fit the stationary Poisson process
# to point pattern 'nztrees'
ppm(nztrees ~ 1)
Q <- quadscheme(nztrees)
ppm(Q ~ 1)
# equivalent.
fit1 <- ppm(nztrees ~ x)
# fit the nonstationary Poisson process
# with intensity function lambda(x,y) = exp(a + bx)
# where x,y are the Cartesian coordinates
# and a,b are parameters to be estimated
fit1
coef(fit1)
coef(summary(fit1))
ppm(nztrees ~ polynom(x,2))
<testonly>ppm(nztrees ~ polynom(x,2), nd=16)</testonly>
# fit the nonstationary Poisson process
# with intensity function lambda(x,y) = exp(a + bx + cx^2)
library(splines)
ppm(nztrees ~ bs(x,df=3))
# WARNING: do not use predict.ppm() on this result
# Fits the nonstationary Poisson process
# with intensity function lambda(x,y) = exp(B(x))
# where B is a B-spline with df = 3
ppm(nztrees ~ 1, Strauss(r=10), rbord=10)
<testonly>ppm(nztrees ~ 1, Strauss(r=10), rbord=10, nd=16)</testonly>
# Fit the stationary Strauss process with interaction range r=10
# using the border method with margin rbord=10
ppm(nztrees ~ x, Strauss(13), correction="periodic")
<testonly>ppm(nztrees ~ x, Strauss(13), correction="periodic", nd=16)</testonly>
# Fit the nonstationary Strauss process with interaction range r=13
# and exp(first order potential) = activity = beta(x,y) = exp(a+bx)
# using the periodic correction.
# Compare Maximum Pseudolikelihood and Huang-Ogata fits:
ppm(swedishpines ~ 1, Strauss(9))
ppm(swedishpines ~ 1, Strauss(9), method="ho")
<testonly>ppm(swedishpines ~ 1, Strauss(9), method="ho", nd=16, nsim=8)</testonly>
# COVARIATES
#
X <- rpoispp(42)
weirdfunction <- function(x,y){ 10 * x^2 + 5 * sin(10 * y) }
#
# (a) covariate values as function
ppm(X ~ y + weirdfunction)
#
# (b) covariate values in pixel image
Zimage <- as.im(weirdfunction, unit.square())
ppm(X ~ y + Z, covariates=list(Z=Zimage))
#
# (c) covariate values in data frame
Q <- quadscheme(X)
xQ <- x.quad(Q)
yQ <- y.quad(Q)
Zvalues <- weirdfunction(xQ,yQ)
ppm(Q ~ y + Z, data=data.frame(Z=Zvalues))
# Note Q not X
# COVARIATE FUNCTION WITH EXTRA ARGUMENTS
#
f <- function(x,y,a){ y - a }
ppm(X ~ x + f, covfunargs=list(a=1/2))
# COVARIATE: inside/outside window
b <- owin(c(0.1, 0.6), c(0.1, 0.9))
ppm(X ~ b)
## MULTITYPE POINT PROCESSES ###
# fit stationary marked Poisson process
# with different intensity for each species
ppm(lansing ~ marks, Poisson())
<testonly>a <- ppm(amacrine ~ marks, Poisson(), nd=16)</testonly>
# fit nonstationary marked Poisson process
# with different log-cubic trend for each species
ppm(lansing ~ marks * polynom(x,y,3), Poisson())
<testonly>ppm(amacrine ~ marks * polynom(x,y,2), Poisson(), nd=16)</testonly>