density.ppp
Kernel Smoothed Intensity of Point Pattern
Compute a kernel smoothed intensity function from a point pattern.
Usage
## S3 method for class 'ppp':
density(x, sigma, \dots,
weights=NULL, edge=TRUE, varcov=NULL,
at="pixels", leaveoneout=TRUE,
adjust=1, diggle=FALSE, se=FALSE)
Arguments
- x
- Point pattern (object of class
"ppp"
). - sigma
- Standard deviation of isotropic Gaussian smoothing kernel.
Either a numerical value, or a function that computes an
appropriate value of
sigma
. - weights
- Optional weights to be attached to the points.
A numeric vector, numeric matrix, or an
expression
. - ...
- Additional arguments passed to
pixellate.ppp
andas.mask
to determine the pixel resolution, or passed tosigma
if it is a function - edge
- Logical value indicating whether to apply edge correction.
- varcov
- Variance-covariance matrix of anisotropic Gaussian kernel.
Incompatible with
sigma
. - at
- String specifying whether to compute the intensity values
at a grid of pixel locations (
at="pixels"
) or only at the points ofx
(at="points"
). - leaveoneout
- Logical value indicating whether to compute a leave-one-out
estimator. Applicable only when
at="points"
. - adjust
- Optional. Adjustment factor for the smoothing parameter.
- diggle
- Logical. If
TRUE
, use Diggle's edge correction, which is more accurate but slower to compute than the correction described under Details. - se
- Logical value indicating whether to compute standard errors as well.
Details
This is a method for the generic function density
.
It computes a fixed-bandwidth kernel estimate
(Diggle, 1985) of the intensity function of the point process
that generated the point pattern x
.
By default it computes the convolution of the
isotropic Gaussian kernel of standard deviation sigma
with point masses at each of the data points in x
.
Anisotropic Gaussian kernels are also supported.
Each point has unit weight, unless the argument weights
is
given.
If edge=TRUE
, the intensity estimate is corrected for
edge effect bias in one of two ways:
- If
diggle=FALSE
(the default) the intensity estimate is correted by dividing it by the convolution of the Gaussian kernel with the window of observation. Thus the intensity value at a point$u$is$$\hat\lambda(u) = e(u) \sum_i k(x_i - u) w_i$$where$k$is the Gaussian smoothing kernel,$e(u)$is an edge correction factor, and$w_i$are the weights. - If
diggle=TRUE
then the method of Diggle (1985) is followed exactly. The intensity value at a point$u$is$$\hat\lambda(u) = \sum_i k(x_i - u) w_i e(x_i)$$where again$k$is the Gaussian smoothing kernel,$e(x_i)$is an edge correction factor, and$w_i$are the weights. This computation is slightly slower but more accurate.
The smoothing kernel is determined by the arguments
sigma
, varcov
and adjust
.
- if
sigma
is a single numerical value, this is taken as the standard deviation of the isotropic Gaussian kernel. - alternatively
sigma
may be a function that computes an appropriate bandwidth for the isotropic Gaussian kernel from the data point pattern by callingsigma(x)
. To perform automatic bandwidth selection using cross-validation, it is recommended to use the functionsbw.diggle
orbw.ppl
. - The smoothing kernel may be chosen to be any Gaussian
kernel, by giving the variance-covariance matrix
varcov
. The argumentssigma
andvarcov
are incompatible. - Alternatively
sigma
may be a vector of length 2 giving the standard deviations of two independent Gaussian coordinates, thus equivalent tovarcov = diag(rep(sigma^2, 2))
. - if neither
sigma
norvarcov
is specified, an isotropic Gaussian kernel will be used, with a default value ofsigma
calculated by a simple rule of thumb that depends only on the size of the window. - The argument
adjust
makes it easy for the user to change the bandwidth specified by any of the rules above. The value ofsigma
will be multiplied by the factoradjust
. The matrixvarcov
will be multiplied byadjust^2
. To double the smoothing bandwidth, setadjust=2
.
...
passed to as.mask
. If at="points"
, the intensity values are computed
to high accuracy at the points of x
only. Computation is
performed by directly evaluating and summing the Gaussian kernel
contributions without discretising the data. The result is a numeric
vector giving the density values.
The intensity value at a point $x_i$ is (if diggle=FALSE
)
$$\hat\lambda(x_i) = e(x_i) \sum_j k(x_j - x_i) w_j$$
or (if diggle=TRUE
)
$$\hat\lambda(x_i) = \sum_j k(x_j - x_i) w_j e(x_j)$$
If leaveoneout=TRUE
(the default), then the sum in the equation
is taken over all $j$ not equal to $i$,
so that the intensity value at a
data point is the sum of kernel contributions from
all other data points.
If leaveoneout=FALSE
then the sum is taken over all $j$,
so that the intensity value at a data point includes a contribution
from the same point.
If weights
is a matrix with more than one column, then the
calculation is effectively repeated for each column of weights. The
result is a list of images (if at="pixels"
) or a matrix of
numerical values (if at="points"
).
The argument weights
can also be an expression
.
It will be evaluated in the data frame as.data.frame(x)
to obtain a vector or matrix of weights. The expression may involve
the symbols x
and y
representing the Cartesian
coordinates, the symbol marks
representing the mark values
if there is only one column of marks, and the names of the columns of
marks if there are several columns.
To select the bandwidth sigma
automatically by
cross-validation, use bw.diggle
or bw.ppl
.
To perform spatial interpolation of values that were observed
at the points of a point pattern, use Smooth.ppp
.
For adaptive nonparametric estimation, see
adaptive.density
.
For data sharpening, see sharpen.ppp
.
To compute a relative risk surface or probability map for
two (or more) types of points, use relrisk
.
Value
- By default, the result is
a pixel image (object of class
"im"
). Pixel values are estimated intensity values, expressed inpoints per unit area .If
at="points"
, the result is a numeric vector of length equal to the number of points inx
. Values are estimated intensity values at the points ofx
.In either case, the return value has attributes
"sigma"
and"varcov"
which report the smoothing bandwidth that was used.If
weights
is a matrix with more than one column, then the result is a list of images (ifat="pixels"
) or a matrix of numerical values (ifat="points"
). Ifse=TRUE
, the result is a list with two elements namedestimate
andSE
, each of the format described above.
Note
This function is often misunderstood.
The result of density.ppp
is not a spatial smoothing
of the marks or weights attached to the point pattern.
To perform spatial interpolation of values that were observed
at the points of a point pattern, use Smooth.ppp
.
The result of density.ppp
is not a probability density.
It is an estimate of the intensity function of the
point process that generated the point pattern data.
Intensity is the expected number of random points
per unit area.
The units of intensity are density.ppp
.
The integral of the intensity function over a spatial region gives the
expected number of points falling in this region.
Inspecting an estimate of the intensity function is usually the first step in exploring a spatial point pattern dataset. For more explanation, see the workshop notes (Baddeley, 2008) or Diggle (2003).
If you have two (or more) types of points, and you want a
probability map or relative risk surface (the spatially-varying
probability of a given type), use relrisk
.
References
Baddeley, A. (2010) Analysing spatial point patterns in R.
Workshop notes. CSIRO online technical publication.
URL: www.uwa.edu.au/resources/pf16h.html
Diggle, P.J. (1985) A kernel method for smoothing point process data. Applied Statistics (Journal of the Royal Statistical Society, Series C) 34 (1985) 138--147.
Diggle, P.J. (2003) Statistical analysis of spatial point patterns, Second edition. Arnold.
See Also
bw.diggle
,
bw.ppl
,
Smooth.ppp
,
sharpen.ppp
,
adaptive.density
,
relrisk
,
ppp.object
,
im.object
Examples
data(cells)
if(interactive()) {
opa <- par(mfrow=c(1,2))
plot(density(cells, 0.05))
plot(density(cells, 0.05, diggle=TRUE))
par(opa)
v <- diag(c(0.05, 0.07)^2)
plot(density(cells, varcov=v))
}
<testonly>Z <- density(cells, 0.05)
Z <- density(cells, 0.05, diggle=TRUE)
Z <- density(cells, varcov=diag(c(0.05^2, 0.07^2)))
Z <- density(cells, 0.05, weights=data.frame(a=1:42,b=42:1))
Z <- density(cells, 0.05, weights=expression(x))</testonly>
# automatic bandwidth selection
plot(density(cells, sigma=bw.diggle(cells)))
# equivalent:
plot(density(cells, bw.diggle))
# evaluate intensity at points
density(cells, 0.05, at="points")