Intensity estimation

As seen in the notes, there are three main approaches for explaining spatial variation in risk: density ratio, binary regression and generalised additive model. Here we learn how to produce smooth intensity maps using R.

In the package spatstat, the main function for intensity estimation is density.ppp. When you apply this function to an object of class ppp, you can omit ”.ppp” from the call, i.e., instead of writing density.ppp(patt3,positive=TRUE), you would simply write:

den3 <- density(patt3, positive = TRUE)
den3

real-valued pixel image
128 x 128 pixel array (ny, nx)
enclosing rectangle: [0, 1] x [0, 2] units

Here, the option positive=TRUE forces the function to return positive values only for the density (otherwise, this function can return very small negative numbers where the intensity is close to zero). This retuns an object of class im, for which there exist plotting methods:

plot(den3)

For some (strange?) reason, the values of the intensity funtion are stored in this object in transpose (the element den3$\mathrm{v}$), so if we wanted to use the function image.plot from the package fields, then we’d have to transpose back first:

require(fields)
xvals <- den3$xcol  # the x-values
yvals <- den3$yrow  # the y-values
zvals <- t(den3$v)  # the estimated intensity, transposed

image.plot(xvals, yvals, zvals, asp = 1)

The function density.ppp fits a bivariate, isotropic, Gaussian kernel. The word “isotropic” in this context means that the variance matrix of the kernel functions is diagonal with equal elements i.e., a scalar multiple of the identity matrix, for example .

The function density.ppp computes an estimate, $\widehat{\lambda }$ of the true intensity $\lambda$ via,

$$\widehat{\lambda }(s)=e(s)\sum _{i=1}^{n}K(s-s_i),$$

(1.1)

where $e(s)$ is a bias-correction for edge effects,

$$e(s)=1/{\int }_{W}K(s-x)dx.$$

Equation 1.1 is an unbiased estimate of,

$${\lambda }^{*}(s)=e(s){\int }_{W}K(s-x)\lambda (x)dx,$$

a smoothed version of the true intensity, see Baddeley (2010).

Notice that the call to density.ppp did not require the user to specify a bandwidth. This is because, according to the help file, spatstat, uses “ a simple rule of thumb that depends only on the size of the window”. There are two inbuilt funtions for computing an appropriate bandwidth using cross validation: bw.diggle and bw.ppl.

To see what bandwidth was selected by each of the routines, we extract the sigma attribute of each of these objects:

attr(den3a, "sigma")

[1] 0.06090998

attr(den3b, "sigma")

[1] 0.2264535