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, edge=TRUE, varcov=NULL,
        at="pixels", leaveoneout=TRUE,
        adjust=1, diggle=FALSE)

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 (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.

The smoothing kernel is determined by the arguments sigma, varcov and adjust.

• ifsigmais a single numerical value, this is taken as the standard deviation of the isotropic Gaussian kernel.
• alternativelysigmamay 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 functionbw.diggle.
• The smoothing kernel may be chosen to be any Gaussian kernel, by giving the variance-covariance matrixvarcov. The argumentssigmaandvarcovare incompatible.
• Alternativelysigmamay be a vector of length 2 giving the standard deviations of two independent Gaussian coordinates, thus equivalent tovarcov = diag(rep(sigma^2, 2)).
• if neithersigmanorvarcovis specified, an isotropic Gaussian kernel will be used, with a default value ofsigmacalculated by a simple rule of thumb that depends only on the size of the window.
• The argumentadjustmakes it easy for the user to change the bandwidth specified by any of the rules above. The value ofsigmawill be multiplied by the factoradjust. The matrixvarcovwill be multiplied byadjust^2. To double the smoothing bandwidth, setadjust=2.

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 select the bandwidth sigma automatically by cross-validation, use bw.diggle. 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.

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 points per unit area. Intensity is usually a function of spatial location, and it is this function which is estimated by 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.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

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

Aliases

• density.ppp

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)))</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")