spatstat.explore (version 3.7-0)

density.ppp: Kernel Smoothed Intensity of Point Pattern

Description

Compute a kernel smoothed intensity function from a point pattern.

Usage

# S3 method for ppp
density(x, sigma=NULL, ...,
        weights=NULL, edge=TRUE, varcov=NULL,
        at="pixels", leaveoneout=TRUE,
        adjust=1, diggle=FALSE,
        se=FALSE, wtype=c("value", "multiplicity"),
        kernel="gaussian",
        scalekernel=is.character(kernel), 
        positive=FALSE, verbose=TRUE, sameas)

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.

If weights is a matrix with more than one column, then the result is a list of images (if at="pixels") or a matrix of numerical values (if at="points").

If se=TRUE, the result is a list with two elements named

estimate and SE, each of the format described above.

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.

The amount of smoothing is controlled by sigma if it is specified.

By default, smoothing is performed using a Gaussian kernel. The resulting density estimate is the convolution of the isotropic Gaussian kernel, of standard deviation sigma, with point masses at each of the data points in x.

Anisotropic kernels, and non-Gaussian kernels, are also supported. Each point has unit weight, unless the argument weights is given.

If edge=TRUE (the default), the intensity estimate is corrected for edge effect bias.

If at="pixels" (the default), the result is a pixel image giving the estimated intensity at each pixel in a grid. If at="points", the result is a numeric vector giving the estimated intensity at each of the original data points in x.

References

Baddeley, A., Rubak, E. and Turner, R. (2015) Spatial Point Patterns: Methodology and Applications with R. Chapman and Hall/CRC Press.

Bithell, J.F. (1990) An application of density estimation to geographical epidemiology. Statistics in Medicine 9, 691--701.

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.

Diggle, P.J. (2010) Nonparametric methods. Chapter 18, pp. 299--316 in A.E. Gelfand, P.J. Diggle, M. Fuentes and P. Guttorp (eds.) Handbook of Spatial Statistics, CRC Press, Boca Raton, FL.

Jones, M.C. (1993) Simple boundary corrections for kernel density estimation. Statistics and Computing 3, 135--146.

See Also

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

For information about the data structures, see ppp.object, im.object.

Examples

  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))
  }
  # 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")

  # non-Gaussian kernel
  plot(density(cells, sigma=0.4, kernel="epanechnikov"))

  if(interactive()) {
    # see effect of changing pixel resolution
    opa <- par(mfrow=c(1,2))
    plot(density(cells, sigma=0.4))
    plot(density(cells, sigma=0.4, eps=0.05))
    par(opa)
  }

  # relative risk calculation by hand (see relrisk.ppp)
  lung <- split(chorley)$lung
  larynx <- split(chorley)$larynx
  D <- density(lung, sigma=2)
  plot(density(larynx, sigma=2, weights=1/D))