Learn R Programming

spatstat (version 1.23-1)

density.ppp: Kernel Smoothed Intensity of Point Pattern

Description

Compute a kernel smoothed intensity function from a point pattern.

Usage

## S3 method for class 'ppp':
density(x, sigma, \dots,
        weights, edge=TRUE, varcov=NULL,
        at="pixels", leaveoneout=TRUE,
        adjust=1, diggle=FALSE)

Arguments

x
Point pattern (object of class "ppp").
sigma
Standard deviation of isotropic Gaussian smoothing kernel.
weights
Optional vector of weights to be attached to the points. May include negative values.
...
Arguments passed to as.mask to determine the pixel resolution.
edge
Logical flag: if TRUE, 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 of x (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.

Value

  • By default, the result is a pixel image (object of class "im"). Pixel values are estimated intensity values, expressed in points per unit area.

    If at="points", the result is a numeric vector of length equal to the number of points in x. Values are estimated intensity values at the points of x.

    In either case, the return value has attributes "sigma" and "varcov" which report the smoothing bandwidth that was used.

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. Each point has unit weight, unless the argument weights is given (it should be a numeric vector; weights can be negative or zero).

If edge=TRUE, the intensity estimate is corrected for edge effect bias in one of two ways:

  • Ifdiggle=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.
  • Ifdiggle=TRUEthen 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.
In both cases, the edge correction term $e(u)$ is the reciprocal of the kernel mass inside the window: $$\frac{1}{e(u)} = \int_W k(v-u) \, {\rm d}v$$ where $W$ is the observation window.

Instead of the isotropic Gaussian kernel with standard deviation sigma, the smoothing kernel may be chosen to be any Gaussian kernel, by giving the variance-covariance matrix varcov. The arguments sigma and varcov are incompatible. Also sigma may be a vector of length 2 giving the standard deviations of two independent Gaussian coordinates, thus equivalent to varcov = diag(rep(sigma^2, 2)). The smoothing parameter sigma has a default value 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 this default. The value of sigma will be multiplied by the factor adjust. To double the smoothing parameter, set adjust=2.

By default the intensity values are computed at every location $u$ in a fine grid, and are returned as a pixel image. Computation is performed using the Fast Fourier Transform. Accuracy depends on the pixel resolution, controlled by the arguments ... 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.

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.

References

Baddeley, A. (2008) Analysing spatial point patterns in R. Workshop notes. CSIRO online technical publication. URL: www.csiro.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

smooth.ppp, sharpen.ppp, adaptive.density, relrisk, ppp.object, im.object

Examples

Run this code
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)))</testonly>
  density(cells, 0.05, at="points")

Run the code above in your browser using DataLab