pcf.ppp: Pair Correlation Function of Point Pattern
In spatstat.core: Core Functionality of the 'spatstat' Family
pcf.ppp R Documentation
Pair Correlation Function of Point Pattern
Description
Estimates the pair correlation function of
a point pattern using kernel methods.
Usage
## S3 method for class 'ppp'
pcf(X, ..., r = NULL, kernel="epanechnikov", bw=NULL,
stoyan=0.15,
correction=c("translate", "Ripley"),
divisor = c("r", "d"),
var.approx = FALSE,
domain=NULL,
ratio=FALSE, close=NULL)
Arguments
X
A point pattern (object of class "ppp").
r
Vector of values for the argument r at which g(r)
should be evaluated. There is a sensible default.
kernel
Choice of smoothing kernel,
passed to density.default.
bw
Bandwidth for smoothing kernel,
passed to density.default.
Either a single numeric value giving the standard deviation of the kernel,
or a character string specifying a bandwidth selection rule
recognised by density.default.
If bw is missing or NULL,
the default value is computed using
Stoyan's rule of thumb: see Details.
...
Other arguments passed to the kernel density estimation
function density.default.
stoyan
Coefficient for Stoyan's bandwidth selection rule; see Details.
correction
Edge correction. A character vector specifying the choice
(or choices) of edge correction. See Details.
divisor
Choice of divisor in the estimation formula:
either "r" (the default) or "d". See Details.
var.approx
Logical value indicating whether to compute an analytic
approximation to the variance of the estimated pair correlation.
domain
Optional. Calculations will be restricted to this subset
of the window. See Details.
ratio
Logical.
If TRUE, the numerator and denominator of
each edge-corrected estimate will also be saved,
for use in analysing replicated point patterns.
close
Advanced use only. Precomputed data. See section on Advanced Use.
Details
The pair correlation function g(r)
is a summary of the dependence between points in a spatial point
process. The best intuitive interpretation is the following: the probability
p(r) of finding two points at locations x and y
separated by a distance r is equal to
p(r) = lambda^2 * g(r) dx dy
where lambda is the intensity of the point process.
For a completely random (uniform Poisson) process,
p(r) = lambda^2 dx dy
so g(r) = 1.
Formally, the pair correlation function of a stationary point process
is defined by
g(r) = K'(r)/ ( 2 * pi * r)
where K'(r) is the derivative of K(r), the
reduced second moment function (aka “Ripley's K function”)
of the point process. See Kest for information
about K(r).
For a stationary Poisson process, the
pair correlation function is identically equal to 1. Values
g(r) < 1 suggest inhibition between points;
values greater than 1 suggest clustering.
This routine computes an estimate of g(r)
by kernel smoothing.
-
If divisor="r" (the default), then the standard
kernel estimator (Stoyan and Stoyan, 1994, pages 284–285)
is used. By default, the recommendations of Stoyan and Stoyan (1994)
are followed exactly.
-
If divisor="d" then a modified estimator is used:
the contribution from
an interpoint distance d[ij] to the
estimate of g(r) is divided by d[ij]
instead of dividing by r. This usually improves the
bias of the estimator when r is close to zero.
There is also a choice of spatial edge corrections
(which are needed to avoid bias due to edge effects
associated with the boundary of the spatial window):
-
If correction="translate" or correction="translation"
then the translation correction
is used. For divisor="r" the translation-corrected estimate
is given in equation (15.15), page 284 of Stoyan and Stoyan (1994).
-
If correction="Ripley" or correction="isotropic"
then Ripley's isotropic edge correction
is used. For divisor="r" the isotropic-corrected estimate
is given in equation (15.18), page 285 of Stoyan and Stoyan (1994).
-
If correction="none" then no edge correction is used,
that is, an uncorrected estimate is computed.
Multiple corrections can be selected. The default is
correction=c("translate", "Ripley").
Alternatively correction="all" selects all options;
correction="best" selects the option which has the best
statistical performance; correction="good" selects the option
which is the best compromise between statistical performance and speed
of computation.
The choice of smoothing kernel is controlled by the
argument kernel which is passed to density.default.
The default is the Epanechnikov kernel, recommended by
Stoyan and Stoyan (1994, page 285).
The bandwidth of the smoothing kernel can be controlled by the
argument bw. Bandwidth is defined as the standard deviation
of the kernel; see the documentation for density.default.
For the Epanechnikov kernel with half-width h,
the argument bw is equivalent to h/sqrt(5).
Stoyan and Stoyan (1994, page 285) recommend using the Epanechnikov
kernel with support [-h,h] chosen by the rule of thumn
h = c/sqrt(lambda),
where lambda is the (estimated) intensity of the
point process, and c is a constant in the range from 0.1 to 0.2.
See equation (15.16).
If bw is missing or NULL,
then this rule of thumb will be applied.
The argument stoyan determines the value of c.
The smoothing bandwidth that was used in the calculation is returned
as an attribute of the final result.
The argument r is the vector of values for the
distance r at which g(r) should be evaluated.
There is a sensible default.
If it is specified, r must be a vector of increasing numbers
starting from r[1] = 0,
and max(r) must not exceed half the diameter of
the window.
If the argument domain is given, estimation will be restricted
to this region. That is, the estimate of
g(r) will be based on pairs of points in which the first point lies
inside domain and the second point is unrestricted.
The argument domain
should be a window (object of class "owin") or something acceptable to
as.owin. It must be a subset of the
window of the point pattern X.
To compute a confidence band for the true value of the
pair correlation function, use lohboot.
If var.approx = TRUE, the variance of the
estimate of the pair correlation will also be calculated using
an analytic approximation (Illian et al, 2008, page 234)
which is valid for stationary point processes which are not
too clustered. This calculation is not yet implemented when
the argument domain is given.
Value
A function value table
(object of class "fv").
Essentially a data frame containing the variables
r
the vector of values of the argument r
at which the pair correlation function g(r) has been estimated
theo
vector of values equal to 1,
the theoretical value of g(r) for the Poisson process
trans
vector of values of g(r)
estimated by translation correction
iso
vector of values of g(r)
estimated by Ripley isotropic correction
v
vector of approximate values of the variance of
the estimate of g(r)
as required.
If ratio=TRUE then the return value also has two
attributes called "numerator" and "denominator"
which are "fv" objects
containing the numerators and denominators of each
estimate of g(r).
The return value also has an attribute "bw" giving the
smoothing bandwidth that was used.
Advanced Use
To perform the same computation using several different bandwidths bw,
it is efficient to use the argument close.
This should be the result of closepairs(X, rmax)
for a suitably large value of rmax, namely
rmax >= max(r) + 3 * bw.
Author(s)
\spatstatAuthors
and Martin Hazelton.
References
Illian, J., Penttinen, A., Stoyan, H. and Stoyan, D. (2008)
Statistical Analysis and Modelling of Spatial Point Patterns.
Wiley.
Stoyan, D. and Stoyan, H. (1994)
Fractals, random shapes and point fields:
methods of geometrical statistics.
John Wiley and Sons.
See Also
Kest,
pcf,
density.default,
bw.stoyan,
bw.pcf,
lohboot.
Examples
X <- simdat
p <- pcf(X)
plot(p, main="pair correlation function for X")
# indicates inhibition at distances r < 0.3
pd <- pcf(X, divisor="d")
# compare estimates
plot(p, cbind(iso, theo) ~ r, col=c("blue", "red"),
ylim.covers=0, main="", lwd=c(2,1), lty=c(1,3), legend=FALSE)
plot(pd, iso ~ r, col="green", lwd=2, add=TRUE)
legend("center", col=c("blue", "green"), lty=1, lwd=2,
legend=c("divisor=r","divisor=d"))
# calculate approximate variance and show POINTWISE confidence bands
pv <- pcf(X, var.approx=TRUE)
plot(pv, cbind(iso, iso+2*sqrt(v), iso-2*sqrt(v)) ~ r)
spatstat.core documentation built on May 18, 2022, 9:05 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| pcf.ppp | R Documentation |
Pair Correlation Function of Point Pattern
Description
Estimates the pair correlation function of a point pattern using kernel methods.
Usage
## S3 method for class 'ppp'
pcf(X, ..., r = NULL, kernel="epanechnikov", bw=NULL,
stoyan=0.15,
correction=c("translate", "Ripley"),
divisor = c("r", "d"),
var.approx = FALSE,
domain=NULL,
ratio=FALSE, close=NULL)
Arguments
X |
A point pattern (object of class |
r |
Vector of values for the argument r at which g(r) should be evaluated. There is a sensible default. |
kernel |
Choice of smoothing kernel,
passed to |
bw |
Bandwidth for smoothing kernel,
passed to |
... |
Other arguments passed to the kernel density estimation
function |
stoyan |
Coefficient for Stoyan's bandwidth selection rule; see Details. |
correction |
Edge correction. A character vector specifying the choice (or choices) of edge correction. See Details. |
divisor |
Choice of divisor in the estimation formula:
either |
var.approx |
Logical value indicating whether to compute an analytic approximation to the variance of the estimated pair correlation. |
domain |
Optional. Calculations will be restricted to this subset of the window. See Details. |
ratio |
Logical.
If |
close |
Advanced use only. Precomputed data. See section on Advanced Use. |
Details
The pair correlation function g(r) is a summary of the dependence between points in a spatial point process. The best intuitive interpretation is the following: the probability p(r) of finding two points at locations x and y separated by a distance r is equal to
p(r) = lambda^2 * g(r) dx dy
where lambda is the intensity of the point process. For a completely random (uniform Poisson) process, p(r) = lambda^2 dx dy so g(r) = 1. Formally, the pair correlation function of a stationary point process is defined by
g(r) = K'(r)/ ( 2 * pi * r)
where K'(r) is the derivative of K(r), the
reduced second moment function (aka “Ripley's K function”)
of the point process. See Kest for information
about K(r).
For a stationary Poisson process, the pair correlation function is identically equal to 1. Values g(r) < 1 suggest inhibition between points; values greater than 1 suggest clustering.
This routine computes an estimate of g(r) by kernel smoothing.
-
If
divisor="r"(the default), then the standard kernel estimator (Stoyan and Stoyan, 1994, pages 284–285) is used. By default, the recommendations of Stoyan and Stoyan (1994) are followed exactly. -
If
divisor="d"then a modified estimator is used: the contribution from an interpoint distance d[ij] to the estimate of g(r) is divided by d[ij] instead of dividing by r. This usually improves the bias of the estimator when r is close to zero.
There is also a choice of spatial edge corrections (which are needed to avoid bias due to edge effects associated with the boundary of the spatial window):
-
If
correction="translate"orcorrection="translation"then the translation correction is used. Fordivisor="r"the translation-corrected estimate is given in equation (15.15), page 284 of Stoyan and Stoyan (1994). -
If
correction="Ripley"orcorrection="isotropic"then Ripley's isotropic edge correction is used. Fordivisor="r"the isotropic-corrected estimate is given in equation (15.18), page 285 of Stoyan and Stoyan (1994). -
If
correction="none"then no edge correction is used, that is, an uncorrected estimate is computed.
Multiple corrections can be selected. The default is
correction=c("translate", "Ripley").
Alternatively correction="all" selects all options;
correction="best" selects the option which has the best
statistical performance; correction="good" selects the option
which is the best compromise between statistical performance and speed
of computation.
The choice of smoothing kernel is controlled by the
argument kernel which is passed to density.default.
The default is the Epanechnikov kernel, recommended by
Stoyan and Stoyan (1994, page 285).
The bandwidth of the smoothing kernel can be controlled by the
argument bw. Bandwidth is defined as the standard deviation
of the kernel; see the documentation for density.default.
For the Epanechnikov kernel with half-width h,
the argument bw is equivalent to h/sqrt(5).
Stoyan and Stoyan (1994, page 285) recommend using the Epanechnikov
kernel with support [-h,h] chosen by the rule of thumn
h = c/sqrt(lambda),
where lambda is the (estimated) intensity of the
point process, and c is a constant in the range from 0.1 to 0.2.
See equation (15.16).
If bw is missing or NULL,
then this rule of thumb will be applied.
The argument stoyan determines the value of c.
The smoothing bandwidth that was used in the calculation is returned
as an attribute of the final result.
The argument r is the vector of values for the
distance r at which g(r) should be evaluated.
There is a sensible default.
If it is specified, r must be a vector of increasing numbers
starting from r[1] = 0,
and max(r) must not exceed half the diameter of
the window.
If the argument domain is given, estimation will be restricted
to this region. That is, the estimate of
g(r) will be based on pairs of points in which the first point lies
inside domain and the second point is unrestricted.
The argument domain
should be a window (object of class "owin") or something acceptable to
as.owin. It must be a subset of the
window of the point pattern X.
To compute a confidence band for the true value of the
pair correlation function, use lohboot.
If var.approx = TRUE, the variance of the
estimate of the pair correlation will also be calculated using
an analytic approximation (Illian et al, 2008, page 234)
which is valid for stationary point processes which are not
too clustered. This calculation is not yet implemented when
the argument domain is given.
Value
A function value table
(object of class "fv").
Essentially a data frame containing the variables
r |
the vector of values of the argument r at which the pair correlation function g(r) has been estimated |
theo |
vector of values equal to 1, the theoretical value of g(r) for the Poisson process |
trans |
vector of values of g(r) estimated by translation correction |
iso |
vector of values of g(r) estimated by Ripley isotropic correction |
v |
vector of approximate values of the variance of the estimate of g(r) |
as required.
If ratio=TRUE then the return value also has two
attributes called "numerator" and "denominator"
which are "fv" objects
containing the numerators and denominators of each
estimate of g(r).
The return value also has an attribute "bw" giving the
smoothing bandwidth that was used.
Advanced Use
To perform the same computation using several different bandwidths bw,
it is efficient to use the argument close.
This should be the result of closepairs(X, rmax)
for a suitably large value of rmax, namely
rmax >= max(r) + 3 * bw.
Author(s)
\spatstatAuthorsand Martin Hazelton.
References
Illian, J., Penttinen, A., Stoyan, H. and Stoyan, D. (2008) Statistical Analysis and Modelling of Spatial Point Patterns. Wiley.
Stoyan, D. and Stoyan, H. (1994) Fractals, random shapes and point fields: methods of geometrical statistics. John Wiley and Sons.
See Also
Kest,
pcf,
density.default,
bw.stoyan,
bw.pcf,
lohboot.
Examples
X <- simdat
p <- pcf(X)
plot(p, main="pair correlation function for X")
# indicates inhibition at distances r < 0.3
pd <- pcf(X, divisor="d")
# compare estimates
plot(p, cbind(iso, theo) ~ r, col=c("blue", "red"),
ylim.covers=0, main="", lwd=c(2,1), lty=c(1,3), legend=FALSE)
plot(pd, iso ~ r, col="green", lwd=2, add=TRUE)
legend("center", col=c("blue", "green"), lty=1, lwd=2,
legend=c("divisor=r","divisor=d"))
# calculate approximate variance and show POINTWISE confidence bands
pv <- pcf(X, var.approx=TRUE)
plot(pv, cbind(iso, iso+2*sqrt(v), iso-2*sqrt(v)) ~ r)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.