Diagnostic Plots for Fitted Point Process Model

Description

Given a point process model fitted to a point pattern, produce diagnostic plots based on residuals.

Usage

  diagnose.ppm(object, ..., type="raw", which="all", sigma=NULL, 
               rbord=reach(object), cumulative=TRUE,
               plot.it=TRUE, rv = NULL, compute.sd=TRUE,
               compute.cts=TRUE, typename, check=TRUE, repair=TRUE,
               oldstyle=FALSE)
  ## S3 method for class 'diagppm':
  plot(x, ..., which, 
               plot.neg="image", plot.smooth="imagecontour", 
               plot.sd=TRUE, spacing=0.1,
               srange=NULL, monochrome=FALSE, main=NULL)

Arguments

Details

This function generates several diagnostic plots for a fitted point process model. The plots display the residuals from the fitted model (Baddeley et al, 2005) or alternatively the `exponential energy marks' (Stoyan and Grabarnik, 1991). These plots can be used to assess goodness-of-fit, to identify outliers in the data, and to reveal departures from the fitted model. See also the companion function qqplot.ppm.

The argument object must be a fitted point process model (object of class "ppm") typically produced by the maximum pseudolikelihood fitting algorithm ppm).

The argument type selects the type of residual or weight that will be computed. Current options are:

"eem":

exponential energy marks (Stoyan and Grabarnik, 1991) computed by eem. These are positive weights attached to the data points (i.e. the points of the point pattern dataset to which the model was fitted). If the fitted model is correct, then the sum of these weights for all data points in a spatial region B has expected value equal to the area of B. See eem for further explanation.

"raw", "inverse" or "pearson":

point process residuals (Baddeley et al, 2005) computed by the function residuals.ppm. These are residuals attached both to the data points and to some other points in the window of observation (namely, to the dummy points of the quadrature scheme used to fit the model). If the fitted model is correct, then the sum of the residuals in a spatial region B has mean zero. The options are

• "raw": the raw residuals;
• "inverse": the `inverse-lambda' residuals, a counterpart of the exponential energy weights;
• "pearson": the Pearson residuals.

See residuals.ppm for further explanation.

The argument which selects the type of plot that is produced. Options are:

"marks":

plot the residual measure. For the exponential energy weights (type="eem") this displays circles centred at the points of the data pattern, with radii proportional to the exponential energy weights. For the residuals (type="raw", type="inverse" or type="pearson") this again displays circles centred at the points of the data pattern with radii proportional to the (positive) residuals, while the plotting of the negative residuals depends on the argument plot.neg. If plot.neg="image" then the negative part of the residual measure, which is a density, is plotted as a colour image. If plot.neg="discrete" then the discretised negative residuals (obtained by approximately integrating the negative density using the quadrature scheme of the fitted model) are plotted as squares centred at the dummy points with side lengths proportional to the (negative) residuals. [To control the size of the circles and squares, use the argument maxsize.]

"smooth":

plot a kernel-smoothed version of the residual measure. Each data or dummy point is taken to have a `mass' equal to its residual or exponential energy weight. (Note that residuals can be negative). This point mass is then replaced by a bivariate isotropic Gaussian density with standard deviation sigma. The value of the smoothed residual field at any point in the window is the sum of these weighted densities. If the fitted model is correct, this smoothed field should be flat, and its height should be close to 0 (for the residuals) or 1 (for the exponential energy weights). The field is plotted either as an image, contour plot or perspective view of a surface, according to the argument plot.smooth. The range of values of the smoothed field is printed if the option which="sum" is also selected.

"x":

produce a `lurking variable' plot for the x coordinate. This is a plot of h(x) against x (solid lines) and of E(h(x)) against x (dashed lines), where h(x) is defined below, and E(h(x)) denotes the expectation of h(x) assuming the fitted model is true.

• if cumulative=TRUE then h(x) is the cumulative sum of the weights or residuals for all points which have X coordinate less than or equal to x. For the residuals E(h(x)) = 0, and for the exponential energy weights E(h(x)) = area of the subset of the window to the left of the line X=x.
• if cumulative=FALSE then h(x) is the marginal integral of the smoothed residual field (see the case which="smooth" described above) on the x axis. This is approximately the derivative of the plot for cumulative=TRUE. The value of h(x) is computed by summing the values of the smoothed residual field over all pixels with the given x coordinate. For the residuals E(h(x)) = 0, and for the exponential energy weights E(h(x)) = length of the intersection between the observation window and the line X=x.

If plot.sd = TRUE, then superimposed on the lurking variable plot are the pointwise two-standard-deviation error limits for h(x) calculated for the inhomogeneous Poisson process. The default is plot.sd = TRUE for Poisson models and plot.sd = FALSE for non-Poisson models.

"y":

produce a similar lurking variable plot for the y coordinate.

"sum":

print the sum of the weights or residuals for all points in the window (clipped by a margin rbord if required) and the area of the same window. If the fitted model is correct the sum of the exponential energy weights should equal the area of the window, while the sum of the residuals should equal zero. Also print the range of values of the smoothed field displayed in the "smooth" case.

"all":

All four of the diagnostic plots listed above are plotted together in a two-by-two display. Top left panel is "marks" plot. Bottom right panel is "smooth" plot. Bottom left panel is "x" plot. Top right panel is "y" plot, rotated 90 degrees.

The argument rbord ensures there are no edge effects in the computation of the residuals. The diagnostic calculations will be confined to those points of the data pattern which are at least rbord units away from the edge of the window. The value of rbord should be greater than or equal to the range of interaction permitted in the model.

By default, the two-standard-deviation limits are calculated from the exact formula for the asymptotic variance of the residuals under the asymptotic normal approximation, equation (37) of Baddeley et al (2006). However, for compatibility with the original paper of Baddeley et al (2005), if oldstyle=TRUE, the two-standard-deviation limits are calculated using the innovation variance, an over-estimate of the true variance of the residuals.

The argument rv would normally be used only by experts. It enables the user to substitute arbitrary values for the residuals or marks, overriding the usual calculations. If rv is present, then instead of calculating the residuals from the fitted model, the algorithm takes the entries of the vector rv to be the values of the residuals or marks, and plots them in the manner appropriate to the type of residual or mark selected by type. The length of rv must equal the number of points in the original data point pattern (if type="eem") or the number of quadrature (data and dummy) points used to fit the model (if type="raw", type="inverse" or type="pearson").

The return value is an object of class "diagppm". There are methods for plot and print for such objects. See the Examples.

See also the companion functions qqplot.ppm, which produces a Q-Q plot of the residuals, and lurking, which produces lurking variable plots for any spatial covariate.

Value

An object of class "diagppm" which contains the coordinates needed to reproduce the selected plots. This object can be plotted using plot.diagppm and printed using print.diagppm.

Author(s)

Adrian Baddeley adrian@maths.uwa.edu.au http://www.maths.uwa.edu.au/~adrian/ and Rolf Turner rolf@math.unb.ca http://www.math.unb.ca/~rolf

References

Baddeley, A., Turner, R., Moller, J. and Hazelton, M. (2005) Residual analysis for spatial point processes. Journal of the Royal Statistical Society, Series B 67, 617–666.

Baddeley, A., Moller, J. and Pakes, A.G. (2006) Properties of residuals for spatial point processes. To appear.

Stoyan, D. and Grabarnik, P. (1991) Second-order characteristics for stochastic structures connected with Gibbs point processes. Mathematische Nachrichten, 151:95–100.

See Also

residuals.ppm, eem, ppm.object, qqplot.ppm, lurking, ppm

Examples

    data(cells)
    fit <- ppm(cells, ~x, Strauss(r=0.15), rbord=0.15)
    diagnose.ppm(fit, rbord=0.15)
    diagnose.ppm(fit, type="pearson", rbord=0.15)
    diagnose.ppm(fit, rbord=0.15, which="marks")
    diagnose.ppm(fit, rbord=0.15, type="raw", plot.neg="discrete")

    # save the diagnostics and plot them later
    u <- diagnose.ppm(fit, rbord=0.15, plot.it=FALSE)
    plot(u)
    plot(u, which="marks")