density.ppp: Kernel Smoothed Intensity of Point Pattern

Description Usage Arguments Details Value Negative Values Note Author(s) References See Also Examples

Description

Compute a kernel smoothed intensity function from a point pattern.

Usage

1
2
3
4
5
6
7
8
  ## S3 method for class 'ppp'
density(x, sigma=NULL, ...,
        weights=NULL, edge=TRUE, varcov=NULL,
        at="pixels", leaveoneout=TRUE,
        adjust=1, diggle=FALSE, se=FALSE,
        kernel="gaussian",
        scalekernel=is.character(kernel), 
        positive=FALSE, verbose=TRUE) 

Arguments

x

Point pattern (object of class "ppp").

sigma

Standard deviation of isotropic 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, an expression, or a pixel image.

...

Additional arguments passed to pixellate.ppp and as.mask to determine the pixel resolution, or passed to sigma if it is a function.

edge

Logical value indicating whether to apply edge correction.

varcov

Variance-covariance matrix of anisotropic smoothing 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 the Jones-Diggle improved edge correction, which is more accurate but slower to compute than the default correction.

kernel

The smoothing kernel. A character string specifying the smoothing kernel (current options are "gaussian", "epanechnikov", "quartic" or "disc"), or a pixel image (object of class "im") containing values of the kernel, or a function(x,y) which yields values of the kernel.

scalekernel

Logical value. If scalekernel=TRUE, then the kernel will be rescaled to the bandwidth determined by sigma and varcov: this is the default behaviour when kernel is a character string. If scalekernel=FALSE, then sigma and varcov will be ignored: this is the default behaviour when kernel is a function or a pixel image.

se

Logical value indicating whether to compute standard errors as well.

positive

Logical value indicating whether to force all density values to be positive numbers. Default is FALSE.

verbose

Logical value indicating whether to issue warnings about numerical problems and conditions.

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:

In both cases, the edge correction term e(u) is the reciprocal of the kernel mass inside the window:

1/e(u) = integral[v in W] k(v-u) dv

where W is the observation window.

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

If at="pixels" (the default), intensity values are computed at every location u in a fine grid, and are returned as a pixel image. The point pattern is first discretised using pixellate.ppp, then the intensity is computed using the Fast Fourier Transform. Accuracy depends on the pixel resolution and the discretisation rule. The pixel resolution is controlled by the arguments ... passed to as.mask (specify the number of pixels by dimyx or the pixel size by eps). The discretisation rule is controlled by the arguments ... passed to pixellate.ppp (the default rule is that each point is allocated to the nearest pixel centre; this can be modified using the arguments fractional and preserve).

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)

λ(x[i]) = e(x[i]) ∑[j] k(x[j] - x[i]) w[j]

or (if diggle=TRUE)

λ(x[i]) = ∑[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.

The argument weights can also be a pixel image (object of class "im"). numerical weights for the data points will be extracted from this image (by looking up the pixel values at the locations of the data points in x).

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

Negative Values

Negative and zero values of the density estimate are possible when at="pixels" because of numerical errors in finite-precision arithmetic.

By default, density.ppp does not try to repair such errors. This would take more computation time and is not always needed. (Also it would not be appropriate if weights include negative values.)

To ensure that the resulting density values are always positive, set positive=TRUE.

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 Baddeley, Rubak and Turner (2015) or Diggle (2003, 2010).

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.

Author(s)

\spatstatAuthors

References

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

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

bw.diggle, bw.ppl, Smooth.ppp, sharpen.ppp, adaptive.density, relrisk, ppp.object, im.object

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
  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))
  }
  
    Z <- density(cells, 0.05)
    Z <- density(cells, 0.05, diggle=TRUE)
    Z <- density(cells, 0.05, se=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))
  
  # 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")

  plot(density(cells, sigma=0.4, kernel="epanechnikov"))

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

Example output

Loading required package: nlme
Loading required package: rpart

spatstat 1.51-0       (nickname: 'Poetic Licence') 
For an introduction to spatstat, type 'beginner' 

 [1]  2.6230115  6.6612002  2.3346718  9.3678806  6.8726363  0.7214982
 [7]  1.4306543  8.4519784  2.2675716  2.6218196  1.1731435  3.9791376
[13]  7.1944496  6.8086060  4.7354784  2.2482630  3.8053247  4.5214383
[19] 12.4922086  4.2571137  7.2928882  1.1176274  2.3007513 22.4488997
[25]  5.0504982 10.9080134 10.0143852  0.8987745  2.1162928 10.1260816
[31]  7.5433476 19.7851133  2.5440620  8.1479922  8.6637017  2.9296855
[37]  2.3009078  2.0798390  1.4667999  3.2263980  3.4377503  1.9701207
attr(,"sigma")
[1] 0.05

spatstat documentation built on Nov. 21, 2017, 9:06 a.m.