Package 'spatstat.linnet'

Description

The spatstat.linnet package belongs to the spatstat family of packages. It contains the functionality for analysing spatial data on a linear network.

Details

spatstat is a family of R packages for the statistical analysis of spatial data. Its main focus is the analysis of spatial patterns of points in two-dimensional space.

The original spatstat package has now been split into several sub-packages.

This sub-package spatstat.linnet contains the user-level functions from spatstat that are concerned with spatial data on a linear network.

Structure of the spatstat family

The orginal spatstat package grew to be very large. It has now been divided into several sub-packages:

• spatstat.utils containing basic utilities

• spatstat.sparse containing linear algebra utilities

• spatstat.data containing datasets

• spatstat.univar containing functions for estimating probability distributions of random variables

• spatstat.geom containing geometrical objects and geometrical operations

• spatstat.explore containing the main functionality for exploratory and non-parametric analysis of spatial data

• spatstat.model containing the main functionality for statistical modelling and inference for spatial data

• spatstat.linnet containing functions for spatial data on a linear network

• spatstat, which simply loads the other sub-packages listed above, and provides documentation.

When you install spatstat, these sub-packages are also installed. Then if you load the spatstat package by typing library(spatstat), the other sub-packages listed above will automatically be loaded or imported. For an overview of all the functions available in these sub-packages, see the help file for spatstat in the spatstat package,

Additionally there are several extension packages:

• spatstat.gui for interactive graphics

• spatstat.local for local likelihood (including geographically weighted regression)

• spatstat.Knet for additional, computationally efficient code for linear networks

• spatstat.sphere (under development) for spatial data on a sphere, including spatial data on the earth's surface

The extension packages must be installed separately and loaded explicitly if needed. They also have separate documentation.

Overview of spatstat.linnet

A linear network is a subset of the two-dimensional plane composed of straight line segments. It could represent a road network, for example. Our code requires that, if two segments intersect each other, then the intersection is a single point, and the intersection point is treated as a vertex of the network.

The spatstat.linnet package supports spatial data analysis on a linear network. The primary aim is to analyse spatial patterns of points on a network. The points could represent road accidents on a road network, for example.

The spatstat.linnet package provides code for handling

• linear networks

• point patterns on a linear network

• pixel images on a linear network (where the network is divided into small segments and a numerical value is assigned to each segment)

• functions on a linear network (i.e. functions that are defined at every location along the network)

• tessellations of a linear network (where the network is subdivided into disjoint subsets with different labels)

• point process models on a linear network

Here is a list of the main functionality provided in spatstat.linnet.

Linear networks

An object of class "linnet" represents a linear network. Examples of such objects include the dataset simplenet provided in the package.

Linear network objects can be created by the following functions:

Utilities for manipulating networks include:

A network is called a tree if it has no closed loops. The following functions support the creation and manipulation of trees:

Point patterns on a linear network

An object of class "lpp" represents a point pattern on a linear network (for example, road accidents on a road network).

Examples of such objects include the following datasets provided in the spatstat.data package:

There is also a dataset provided in the extension package spatstat.Knet:

Point patterns on a network can be created by the following functions:

Point patterns on a network can be generated randomly using the following functions:

Functions for manipulating a point pattern on a network include the following. An object of class "lpp" also belongs to the class "ppx", for which additional support is available.

Pixel images on a network

An object of class "linim" represents a pixel image on a linear network. Effectively, the network is divided into small segments (lixels) and each small segment is assigned a value, which could be numeric, factor, logical or complex values.

Pixel images on a network can be created using the following functions:

Functions for manipulating a pixel image on a network include:

Functions on a linear network

An object of class "linfun" represents a function defined at any location along the network. Objects of this class are created by the following functions:

The following supporting code is available:

Tessellations of a linear network

An object of class "lintess" represents a tessellation of the network, that is, a subdivision of the network into disjoint subsets called ‘tiles’. Objects of this class are created by the following functions:

The following functions are provided for manipulating a tessellation on a network:

Smoothing a point pattern on a linear network:

Given a point pattern dataset on a linear network, it is often desired to estimate the spatially-varying density or intensity of points along the network. For example if the points represent road accidents, then we may wish to estimate the spatially-varying density of accidents per unit length (over a given period of time).

Related tasks include estimation of relative risk, and smoothing of of values observed at the data points.

Exploration of dependence on a covariate:

Another task is to investigate how the spatially-varying intensity of points depends on an explanatory variable (covariate). The covariate may be given as a pixel image on the network (class "linim") or as a function on the network (class "linfun").

Summary statistics for a point pattern on a linear network:

These are for point patterns on a linear network (class lpp). For unmarked patterns:

For multitype patterns:

Related facilities:

It is also possible to fit point process models to lpp objects.

Point process models on a linear network:

An object of class "lpp" represents a pattern of points on a linear network. Point process models can also be fitted to these objects. Currently only Poisson models can be fitted.

Licence

This library and its documentation are usable under the terms of the "GNU General Public License", a copy of which is distributed with the package.

Acknowledgements

Ottmar Cronie, Tilman Davies, Greg McSwiggan and Suman Rakshit made substantial contributions of code.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

Description

Adds new vertices to a linear network at specified locations outside the network.

Usage

addVertices(L, X, join=NULL, joinmarks=NULL)

Arguments

Details

This function adds new vertices to an existing linear network L, at specified locations X outside the network.

The argument L can be either a linear network (class "linnet") or some other object that includes a linear network.

The new vertex locations are points outside the network, specified as a point pattern X (object of class "ppp").

The argument join specifies how to join the new vertices to the existing network.

• If join=NULL (the default), the new vertices are simply added to the list of network vertices without being joined to the rest of the network.

• If join is a vector of integers, then these are taken to be indices of existing vertices of L in the order given in V = vertices(L). Then each new vertex X[i] will be joined to an existing vertex V[j] where j = join[i]. Each new vertex is joined to exactly one existing vertex.

• If join="vertices" then each new vertex X[i] is joined to the nearest existing vertex V[j]. Each new vertex is joined to exactly one existing vertex.

• If join="nearest" then each new vertex is projected to the nearest location along on the network; these locations are inserted as new vertices of L; and then each vertex X[i] is joined to the corresponding projected point. Each new vertex is joined to exactly one newly-inserted vertex.

• If join is a point pattern on a network (class "lpp"), it must be defined on the same network as L and it must consist of the same number of points as X. The points of join will be inserted as new vertices of L, and then each vertex X[i] is joined to the corresponding point join[i]. Each new vertex is joined to exactly one newly-inserted vertex.

The result is the modified object, with an attribute "id" such that the ith added vertex has become the id[i]th vertex of the new network.

Value

An object of the same class as L representing the result of adding the new vertices. The result also has an attribute "id" as described in Details.

Author(s)

Adrian Baddeley

See Also

insertVertices to insert vertices along an existing network.

as.lpp, linnet, methods.linnet, joinVertices, thinNetwork.

Examples

opa <- par(mfrow=c(1,3))
   L <- simplenet
   X <- runifpoint(20, Window(simplenet))
   plot(L)
   plot(X, add=TRUE, cols="green", pch=16, cex=2)
   plot(addVertices(L, X, "nearest"), col="red")
   plot(L, add=TRUE, col="grey", lwd=3)
   plot(X, add=TRUE, cols="green", pch=16, cex=2)
   plot(addVertices(L, X, "vertices"), col="red")
   plot(L, add=TRUE, col="grey", lwd=3)
   plot(X, add=TRUE, cols="green", pch=16, cex=2)
   par(opa)

Description

Apply geometrical transformations to a linear network.

Usage

## S3 method for class 'linnet'
affine(X, mat=diag(c(1,1)), vec=c(0,0), ...)

  ## S3 method for class 'linnet'
shift(X, vec=c(0,0), ..., origin=NULL)

  ## S3 method for class 'linnet'
rotate(X, angle=pi/2, ..., centre=NULL)

  ## S3 method for class 'linnet'
scalardilate(X, f, ...)

  ## S3 method for class 'linnet'
rescale(X, s, unitname)

Arguments

Details

These functions are methods for the generic functions affine, shift, rotate, rescale and scalardilate applicable to objects of class "linnet".

All of these functions perform geometrical transformations on the object X, except for rescale, which simply rescales the units of length.

Value

Another linear network (of class "linnet") representing the result of applying the geometrical transformation.

Author(s)

Adrian Baddeley [email protected]

and Rolf Turner [email protected]

See Also

linnet and as.linnet.

Generic functions affine, shift, rotate, scalardilate, rescale.

Examples

U <- rotate(simplenet, pi)
  stretch <- diag(c(2,3))
  Y <- affine(simplenet, mat=stretch)
  shear <- matrix(c(1,0,0.6,1),ncol=2, nrow=2)
  Z <- affine(simplenet, mat=shear, vec=c(0, 1))

Description

Apply geometrical transformations to a point pattern on a linear network.

Usage

## S3 method for class 'lpp'
affine(X, mat=diag(c(1,1)), vec=c(0,0), ...)

  ## S3 method for class 'lpp'
shift(X, vec=c(0,0), ..., origin=NULL)

  ## S3 method for class 'lpp'
rotate(X, angle=pi/2, ..., centre=NULL)

  ## S3 method for class 'lpp'
scalardilate(X, f, ...)

  ## S3 method for class 'lpp'
rescale(X, s, unitname)

Arguments

Details

These functions are methods for the generic functions affine, shift, rotate, rescale and scalardilate applicable to objects of class "lpp".

All of these functions perform geometrical transformations on the object X, except for rescale, which simply rescales the units of length.

Value

Another point pattern on a linear network (object of class "lpp") representing the result of applying the geometrical transformation.

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

See Also

lpp.

Generic functions affine, shift, rotate, scalardilate, rescale.

Examples

X <- rpoislpp(2, simplenet)
  U <- rotate(X, pi)
  V <- shift(X, c(0.1, 0.2))
  stretch <- diag(c(2,3))
  Y <- affine(X, mat=stretch)
  shear <- matrix(c(1,0,0.6,1),ncol=2, nrow=2)
  Z <- affine(X, mat=shear, vec=c(0, 1))

Description

Performs analysis of deviance for two or more fitted point process models on a linear network.

Usage

## S3 method for class 'lppm'
anova(object, ..., test=NULL)

Arguments

Details

This is a method for anova for fitted point process models on a linear network (objects of class "lppm", usually generated by the model-fitting function lppm).

If the fitted models are all Poisson point processes, then this function performs an Analysis of Deviance of the fitted models. The output shows the deviance differences (i.e. 2 times log likelihood ratio), the difference in degrees of freedom, and (if test="Chi") the two-sided p-values for the chi-squared tests. Their interpretation is very similar to that in anova.glm.

If some of the fitted models are not Poisson point processes, then the deviance difference is replaced by the adjusted composite likelihood ratio (Pace et al, 2011; Baddeley et al, 2014).

Value

An object of class "anova", or NULL.

Errors and warnings

models not nested:

There may be an error message that the models are not “nested”. For an Analysis of Deviance the models must be nested, i.e. one model must be a special case of the other. For example the point process model with formula ~x is a special case of the model with formula ~x+y, so these models are nested. However the two point process models with formulae ~x and ~y are not nested.

If you get this error message and you believe that the models should be nested, the problem may be the inability of R to recognise that the two formulae are nested. Try modifying the formulae to make their relationship more obvious.

different sizes of dataset:

There may be an error message from anova.glmlist that “models were not all fitted to the same size of dataset”. This generally occurs when the point process models are fitted on different linear networks.

Author(s)

Adrian Baddeley [email protected]

References

Ang, Q.W. (2010) Statistical methodology for events on a network. Master's thesis, School of Mathematics and Statistics, University of Western Australia.

Ang, Q.W., Baddeley, A. and Nair, G. (2012) Geometrically corrected second-order analysis of events on a linear network, with applications to ecology and criminology. Scandinavian Journal of Statistics 39, 591–617.

Baddeley, A., Turner, R. and Rubak, E. (2015) Adjusted composite likelihood ratio test for Gibbs point processes. Journal of Statistical Computation and Simulation 86 (5) 922–941. DOI: 10.1080/00949655.2015.1044530.

McSwiggan, G., Nair, M.G. and Baddeley, A. (2012) Fitting Poisson point process models to events on a linear network. Manuscript in preparation.

Pace, L., Salvan, A. and Sartori, N. (2011) Adjusting composite likelihood ratio statistics. Statistica Sinica 21, 129–148.

See Also

lppm

Examples

X <- runiflpp(10, simplenet)
 mod0 <- lppm(X ~1)
 modx <- lppm(X ~x)
 anova(mod0, modx, test="Chi")

Description

Converts a tessellation on a linear network into a data frame.

Usage

## S3 method for class 'lintess'
as.data.frame(x, ...)

Arguments

Details

A tessellation on a linear network is a partition of the network into non-overlapping pieces (tiles). Each tile consists of one or more line segments which are subsets of the line segments making up the network. A tile can consist of several disjoint pieces.

This function converts the tessellation x to a data frame. Each row of the data frame specifies one sub-segment of the network, and allocates it to a particular tile. The data frame has the following columns:

• The seg column specifies which line segment of the network contains the sub-segment. Values of seg are integer indices for the network segments in as.psp(as.linnet(x)).

• The t0 and t1 columns specify the start and end points of the sub-segment. They are numeric values between 0 and 1 inclusive, where the values 0 and 1 representing the network vertices that are joined by this network segment.

• The tile column specifies which tile of the tessellation includes this sub-segment. It is a factor whose levels are the names of the tiles.

The tessellation may have marks, which are attached to the tiles of the tessellation. If marks are present, the resulting data frame includes columns containing, for each sub-segment, the mark value of the corresponding tile.

Value

A data frame with columns named seg, t0, t1, tile, and possibly other columns.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

lintess

Examples

X <- lineardirichlet(runiflpp(3, simplenet))
  marks(X) <- letters[1:3]
  as.data.frame(X)

Description

Convert some kind of data to an object of class "linfun" representing a function on a linear network.

Usage

as.linfun(X, ...)

  ## S3 method for class 'linim'
as.linfun(X, ...)

  ## S3 method for class 'lintess'
as.linfun(X, ..., values=marks(X), navalue=NA)

Arguments

Details

An object of class "linfun" represents a function defined on a linear network.

The function as.linfun is generic. The method as.linfun.linim converts objects of class "linim" (pixel images on a linear network) to functions on the network.

The method as.linfun.lintess converts a tessellation on a linear network into a function with a different value on each tile of the tessellation. The function values are specified by the argument values. It should be a vector with one entry for each tile of the tessellation; any point lying in tile number i will return the value v[i]. If values is missing, the marks of the tessellation are taken as the function values. If values is missing and the tessellation has no marks, or if values is given as NULL, then the function returns factor values identifying which tile contains each given point.

Value

Object of class "linfun".

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

linfun

Examples

X <- runiflpp(2, simplenet)
   Y <- runiflpp(5, simplenet)

   # image on network
   D <- density(Y, 0.1)

   f <- as.linfun(D)
   f
   f(X)

   # tessellation on network
   Z <- lineardirichlet(Y)
   g <- as.linfun(Z)
   g(X)
   h <- as.linfun(Z, values = runif(5))
   h(X)

Description

Converts various kinds of data to a pixel image on a linear network.

Usage

as.linim(X, ...)

  ## S3 method for class 'linim'
as.linim(X, ...)

  ## Default S3 method:
as.linim(X, L, ...,
                             eps = NULL, dimyx = NULL, xy = NULL,
                             rule.eps=c("adjust.eps",
                                        "grow.frame", "shrink.frame"),
                             delta=NULL, nd=NULL)

  ## S3 method for class 'linfun'
as.linim(X, L=domain(X), ...,
                            eps = NULL, dimyx = NULL, xy = NULL,
                            rule.eps=c("adjust.eps",
                                        "grow.frame", "shrink.frame"),
                            delta=NULL, nd=NULL)

Arguments

Details

This function converts the data X into a pixel image on a linear network, an object of class "linim" (see linim).

The argument X may be any of the following:

• a function on a linear network, an object of class "linfun".

• a pixel image on a linear network, an object of class "linim".

• a pixel image, an object of class "im".

• any type of data acceptable to as.im, such as a function, numeric value, or window.

First X is converted to a pixel image object Y (object of class "im"). The conversion is performed by as.im. The arguments eps, dimyx, xy and rule.eps determine the pixel resolution.

Next Y is converted to a pixel image on a linear network using linim. The argument L determines the linear network. If L is missing or NULL, then X should be an object of class "linim", and L defaults to the linear network on which X is defined.

In addition to converting the function to a pixel image, the algorithm also generates a fine grid of sample points evenly spaced along each segment of the network (with spacing at most delta coordinate units). The function values at these sample points are stored in the resulting object as a data frame (the argument df of linim). This mechanism allows greater accuracy for some calculations (such as integral.linim).

Value

An image object on a linear network; an object of class "linim".

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected]

See Also

as.im

Examples

f <- function(x,y){ x + y }
  plot(as.linim(f, simplenet))

Description

Given some kind of data on a linear network, the command as.linnet extracts the linear network itself.

Usage

## S3 method for class 'linim'
as.linnet(X, ...)

 ## S3 method for class 'linfun'
as.linnet(X, ...)

 ## S3 method for class 'lintess'
as.linnet(X, ...)

 ## S3 method for class 'lpp'
as.linnet(X, ..., fatal=TRUE, sparse)

Arguments

Details

These are methods for the generic as.linnet for various classes.

The network on which the data are defined is extracted.

Value

A linear network (object of class "linnet").

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

linnet, methods.linnet.

Examples

# make some data
  xcoord <- linfun(function(x,y,seg,tp) { x }, simplenet)
  as.linnet(xcoord)
  X <- as.linim(xcoord)
  as.linnet(X)

Description

Converts a line segment pattern to a linear network.

Usage

## S3 method for class 'psp'
as.linnet(X, ..., eps, sparse=FALSE)

Arguments

Details

This command converts any collection of line segments into a linear network by guessing the connectivity of the network, using the distance threshold eps.

If any segments in X cross over each other, they are first cut into pieces using selfcut.psp.

Then any pair of segment endpoints lying closer than eps units apart, is treated as a single vertex. The linear network is then constructed using linnet.

It would be wise to check the result by plotting the degree of each vertex, as shown in the Examples.

If X has marks, then these are stored in the resulting linear network Y <- as.linnet(X), and can be extracted as marks(as.psp(Y)) or marks(Y$lines).

Value

A linear network (object of class "linnet").

The result also has an attribute "camefrom" indicating the provenance of each line in the resulting network. For example camefrom[3]=2 means that the third line segment in the result is a piece of the second segment of X.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

linnet, selfcut.psp, methods.linnet.

Examples

# make some data
  A <- psp(0.09, 0.55, 0.79, 0.80, window=owin())
  B <- superimpose(A, as.psp(simplenet))

  # convert to a linear network
  L <- as.linnet(B)

  # check validity
  L
  plot(L)
  text(vertices(L), labels=vertexdegree(L))

  # show the pieces that came from original segment number 1
  S <- as.psp(L)
  (camefrom <- attr(L, "camefrom"))
  parts <- which(camefrom == 1)
  plot(S[parts], add=TRUE, col="green", lwd=2)

Description

Convert various kinds of data to a point pattern on a linear network.

Usage

as.lpp(x=NULL, y=NULL, seg=NULL, tp=NULL, ...,
         marks=NULL, L=NULL, check=FALSE, sparse)

Arguments

Details

This function converts data in various formats into a point pattern on a linear network (object of class "lpp").

The possible formats are:

• x is already a point pattern on a linear network (object of class "lpp"). Then x is returned unchanged.

• x is a planar point pattern (object of class "ppp"). Then x is converted to a point pattern on the linear network L using lpp.

• x,y,seg,tp are vectors of equal length. These specify that the ith point has Cartesian coordinates (x[i],y[i]), and lies on segment number seg[i] of the network L, at a fractional position tp[i] along that segment (with tp=0 representing one endpoint and tp=1 the other endpoint of the segment).

• x,y are missing and seg,tp are vectors of equal length as described above.

• seg,tp are NULL, and x,y are data in a format acceptable to xy.coords specifying the Cartesian coordinates.

• Only the arguments x and L are given, and x is a data frame with one of the following types:

  • two columns labelled seg,tp interpreted as local coordinates on the network.

  • two columns labelled x,y interpreted as Cartesian coordinates.

  • four columns labelled x,y,seg,tp interpreted as Cartesian coordinates and local coordinates.

Value

A point pattern on a linear network (object of class "lpp").

Author(s)

Adrian Baddeley [email protected]

and Rolf Turner [email protected]

See Also

lpp.

Examples

A <- as.psp(simplenet)
   X <- runifpointOnLines(10, A)
   is.ppp(X)
   Y <- as.lpp(X, L=simplenet)

Description

Converts data on a linear network into an object of class "owin".

Usage

## S3 method for class 'lpp'
as.owin(W, ..., fatal=TRUE)

 ## S3 method for class 'lppm'
as.owin(W, ..., fatal=TRUE)

Arguments

Details

The class "owin" is a way of specifying the observation window for a point pattern. See owin.object for an overview.

The function as.owin converts data in any of several formats into an object of class "owin" for use by the spatstat package. The function as.owin is generic, with methods for different classes of objects, and a default method.

A long list of methods for as.owin is documented in the help file for as.owin in the spatstat.geom package.

This help file documents additional methods applicable when W is

• an object of class "lpp" representing a point pattern on a linear network. In this case, as.owin extracts the linear network and returns a window containing this network.

• an object of class "lppm" representing a fitted point process model on a linear network. In this case, as.owin extracts the linear network and returns a window containing this network.

If the argument W cannot be converted to a window, then an error will be generated (if fatal=TRUE) or a value of NULL will be returned (if fatal=FALSE).

Value

An object of class "owin" (see owin.object) specifying an observation window.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

as.owin, owin.object, owin.

Additional methods for as.owin may be provided by other packages outside the spatstat family.

Examples

as.owin(simplenet)

Description

Compute the AUC (area under the Receiver Operating Characteristic curve) for a fitted point process model on a linear network.

Usage

## S3 method for class 'lpp'
auc(X, covariate, ..., high = TRUE)

## S3 method for class 'lppm'
auc(X, ...)

Arguments

Details

This command computes the AUC, the area under the Receiver Operating Characteristic curve. The ROC itself is computed by roc.

The function auc is generic, with methods for "ppp" and "ppm" described in the help file for auc.

This help file describes the methods for classes "lpp" and "lppm".

For a point pattern X and a covariate Z, the AUC is a numerical index that measures the ability of the covariate to separate the spatial domain into areas of high and low density of points. Let $x_i$ be a randomly-chosen data point from X and $U$ a randomly-selected location in the study region. The AUC is the probability that $Z(x_i) > Z(U)$ assuming high=TRUE. That is, AUC is the probability that a randomly-selected data point has a higher value of the covariate Z than does a randomly-selected spatial location. The AUC is a number between 0 and 1. A value of 0.5 indicates a complete lack of discriminatory power.

For a fitted point process model X, the AUC measures the ability of the fitted model intensity to separate the spatial domain into areas of high and low density of points. Suppose $\lambda(u)$ is the intensity function of the model. The AUC is the probability that $\lambda(x_i) > \lambda(U)$. That is, AUC is the probability that a randomly-selected data point has higher predicted intensity than does a randomly-selected spatial location. The AUC is not a measure of the goodness-of-fit of the model (Lobo et al, 2007).

Value

Numeric. For auc.lpp, the result is a single number giving the AUC value. For auc.lppm, the result is a numeric vector of length 2 giving the AUC value and the theoretically expected AUC value for this model.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

References

Lobo, J.M., Jimenez-Valverde, A. and Real, R. (2007) AUC: a misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography 17(2) 145–151.

Nam, B.-H. and D'Agostino, R. (2002) Discrimination index, the area under the ROC curve. Pages 267–279 in Huber-Carol, C., Balakrishnan, N., Nikulin, M.S. and Mesbah, M., Goodness-of-fit tests and model validity, Birkhauser, Basel.

See Also

auc, roc, roc.lpp

Examples

auc(spiders, "x")
  fit <- lppm(spiders ~ x + y)
  auc(fit)

Description

Checks whether a character string begins with a particular prefix.

Usage

begins(x, firstbit)

Arguments

Details

This simple wrapper function checks whether (each entry in) x begins with the string firstbit, and returns a logical value or logical vector with one entry for each entry of x. This function is useful mainly for reducing complexity in model formulae.

Value

Logical vector of the same length as x.

Author(s)

Adrian Baddeley [email protected]

Rolf Turner [email protected]

and Ege Rubak [email protected]

Examples

begins(c("Hello", "Goodbye"), "Hell")
  begins("anything", "")

Description

Tests the goodness-of-fit of a Poisson point process model on a linear network, using the approach of Berman (1986).

Usage

## S3 method for class 'lpp'
berman.test(X, covariate,
                         which = c("Z1", "Z2"),
        alternative = c("two.sided", "less", "greater"), ...)

## S3 method for class 'lppm'
berman.test(model, covariate,
                         which = c("Z1", "Z2"),
               alternative = c("two.sided", "less", "greater"), ...)

Arguments

Details

These functions perform a goodness-of-fit test of a Poisson point process model fitted to point pattern data. The observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same values under the model, are compared using either of two test statistics $Z_1$ and $Z_2$ proposed by Berman (1986). The $Z_1$ test is also known as the Lawson-Waller test.

The function berman.test is generic, with methods for point patterns ("ppp" or "lpp") and point process models ("ppm" or "lppm").

See the help file for berman.test for information on the generic function and the methods for data in two-dimensional space, classes "ppp" and "ppm".

This help file describes the methods for data on a linear network, classes "lpp" and "lppm".

• If X is a point pattern dataset (object of class "ppp" or "lpp"), then berman.test(X, ...) performs a goodness-of-fit test of the uniform Poisson point process (Complete Spatial Randomness, CSR) for this dataset.

• If model is a fitted point process model (object of class "ppm" or "lppm") then berman.test(model, ...) performs a test of goodness-of-fit for this fitted model. In this case, model should be a Poisson point process.

The test is performed by comparing the observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same covariate under the model. Thus, you must nominate a spatial covariate for this test.

The argument covariate should be either a function(x,y) or a pixel image (object of class "im" containing the values of a spatial function. If covariate is an image, it should have numeric values, and its domain should cover the observation window of the model. If covariate is a function, it should expect two arguments x and y which are vectors of coordinates, and it should return a numeric vector of the same length as x and y.

First the original data point pattern is extracted from model. The values of the covariate at these data points are collected.

Next the values of the covariate at all locations in the observation window are evaluated. The point process intensity of the fitted model is also evaluated at all locations in the window.

• If which="Z1", the test statistic $Z_1$ is computed as follows. The sum $S$ of the covariate values at all data points is evaluated. The predicted mean $\mu$ and variance $\sigma^2$ of $S$ are computed from the values of the covariate at all locations in the window. Then we compute $Z_1 = (S-\mu)/\sigma$. Closely-related tests were proposed independently by Waller et al (1993) and Lawson (1993) so this test is often termed the Lawson-Waller test in epidemiological literature.

• If which="Z2", the test statistic $Z_2$ is computed as follows. The values of the covariate at all locations in the observation window, weighted by the point process intensity, are compiled into a cumulative distribution function $F$. The probability integral transformation is then applied: the values of the covariate at the original data points are transformed by the predicted cumulative distribution function $F$ into numbers between 0 and 1. If the model is correct, these numbers are i.i.d. uniform random numbers. The standardised sample mean of these numbers is the statistic $Z_2$.

In both cases the null distribution of the test statistic is the standard normal distribution, approximately.

The return value is an object of class "htest" containing the results of the hypothesis test. The print method for this class gives an informative summary of the test outcome.

Value

An object of class "htest" (hypothesis test) and also of class "bermantest", containing the results of the test. The return value can be plotted (by plot.bermantest) or printed to give an informative summary of the test.

Warning

The meaning of a one-sided test must be carefully scrutinised: see the printed output.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

References

Berman, M. (1986) Testing for spatial association between a point process and another stochastic process. Applied Statistics 35, 54–62.

Lawson, A.B. (1993) On the analysis of mortality events around a prespecified fixed point. Journal of the Royal Statistical Society, Series A 156 (3) 363–377.

Waller, L., Turnbull, B., Clark, L.C. and Nasca, P. (1992) Chronic Disease Surveillance and testing of clustering of disease and exposure: Application to leukaemia incidence and TCE-contaminated dumpsites in upstate New York. Environmetrics 3, 281–300.

See Also

cdf.test, quadrat.test, ppm lppm

Examples

#' test of complete randomness
   berman.test(spiders, "x")
   #' test of fitted model
   fit <- lppm(spiders ~ x)
   berman.test(fit, "y", "Z2")

Description

Creates a function which returns the tree branch membership label for any location on a linear network.

Usage

branchlabelfun(L, root = 1)

Arguments

Details

The linear network L must be an acyclic graph (i.e. must not contain any loops) so that it can be interpreted as a tree.

The result of f <- branchlabelfun(L, root) is a function f which gives, for each location on the linear network L, the tree branch label at that location.

Tree branch labels are explained in treebranchlabels.

The result f also belongs to the class "linfun". It can be called using several different kinds of data, as explained in the help for linfun. The values of the function are character strings.

Value

A function (of class "linfun").

Author(s)

Adrian Baddeley [email protected]

Rolf Turner [email protected]

and Ege Rubak [email protected]

See Also

treebranchlabels, linfun

Examples

# make a simple tree
  m <- simplenet$m
  m[8,10] <- m[10,8] <- FALSE
  L <- linnet(vertices(simplenet), m)
  # make function
  f <- branchlabelfun(L, 1)
  plot(f)
  X <- runiflpp(5, L)
  f(X)

Description

Uses likelihood cross-validation to select a smoothing bandwidth for the kernel estimation of point process intensity on a linear network.

Usage

bw.lppl(X, ..., srange=NULL, ns=16, sigma=NULL, weights=NULL,
           distance=c("euclidean", "path"), shortcut=TRUE, warn=TRUE)

Arguments

Details

This function selects an appropriate bandwidth sigma for the kernel estimator of point process intensity computed by density.lpp.

The argument X should be a point pattern on a linear network (class "lpp").

The bandwidth $\sigma$ is chosen to maximise the point process likelihood cross-validation criterion

$\mbox{LCV}(\sigma) = \sum_i \log\hat\lambda_{-i}(x_i) - \int_L \hat\lambda(u) \, {\rm d}u$

where the sum is taken over all the data points $x_i$, where $\hat\lambda_{-i}(x_i)$ is the leave-one-out kernel-smoothing estimate of the intensity at $x_i$ with smoothing bandwidth $\sigma$, and $\hat\lambda(u)$ is the kernel-smoothing estimate of the intensity at a spatial location $u$ with smoothing bandwidth $\sigma$. See Loader(1999, Section 5.3).

The value of $\mbox{LCV}(\sigma)$ is computed directly, using density.lpp, for ns different values of $\sigma$ between srange[1] and srange[2].

The result is a numerical value giving the selected bandwidth. The result also belongs to the class "bw.optim" which can be plotted to show the (rescaled) mean-square error as a function of sigma.

If shortcut=TRUE, the computation is accelerated by omitting the integral term in the equation above. This is valid because the integral is approximately constant.

Value

A single numerical value giving the selected bandwidth. The result also belongs to the class "bw.optim" (see bw.optim.object) which can be plotted to show the bandwidth selection criterion as a function of sigma.

Author(s)

Greg McSwiggan, Suman Rakshit and Adrian Baddeley [email protected].

References

Loader, C. (1999) Local Regression and Likelihood. Springer, New York.

McSwiggan, G., Baddeley, A. and Nair, G. (2019) Estimation of relative risk for events on a linear network. Statistics and Computing 30 (2) 469–484.

See Also

density.lpp, bw.scott.

bw.optim.object.

For point patterns in two-dimensional space, use bw.ppl.

Examples

if(interactive()) {
    b <- bw.lppl(spiders)
    plot(b, main="Likelihood cross validation for spiders")
    plot(density(spiders, b, distance="e"))
  } else {
    b1 <- bw.lppl(spiders, ns=2)
    b2 <- bw.lppl(spiders, ns=2, shortcut=FALSE)
  }

Description

Uses cross-validation to select a smoothing bandwidth for the estimation of relative risk on a linear network.

Usage

## S3 method for class 'lpp'
bw.relrisk(X, ...,
    method = c("likelihood", "leastsquares", "KelsallDiggle", "McSwiggan"),
    distance=c("path", "euclidean"),
    hmin = NULL, hmax = NULL, nh = NULL,
    fast = TRUE, fastmethod = "onestep",
    floored = TRUE, reference = c("thumb", "uniform", "sigma"),
    allow.infinite = TRUE, epsilon = 1e-20, fudge = 0,
    verbose = FALSE, warn = TRUE)

Arguments

Details

This function is a method for the generic bw.relrisk. It computes an optimal value of smoothing bandwidth for the nonparametric estimation of relative risk on a linear network using relrisk.lpp. The optimal value is found by minimising a cross-validation criterion.

The cross-validation criterion is selected by the argument method:

See McSwiggan et al (2019) for details.

The result is a numerical value giving the selected bandwidth sigma. The result also belongs to the class "bw.optim" allowing it to be printed and plotted. The plot shows the cross-validation criterion as a function of bandwidth. The ‘optimal’ bandwidth is the value of bandwidth which minimises the cross-validation criterion.

The range of values for the smoothing bandwidth sigma is set by the arguments hmin, hmax. There is a sensible default, based on the linear network version of Scott's rule bw.scott.iso.

If the optimal bandwidth is achieved at an endpoint of the interval [hmin, hmax], the algorithm will issue a warning (unless warn=FALSE). If this occurs, then it is probably advisable to expand the interval by changing the arguments hmin, hmax.

The cross-validation procedure is based on kernel estimates of intensity, which are computed by density.lpp. Any arguments ... are passed to density.lpp to control the kernel estimation procedure. This includes the argument distance which specifies the type of kernel. The default is distance="path"; the fastest option is distance="euclidean".

Value

A single numerical value giving the selected bandwidth. The result also belongs to the class "bw.optim" (see bw.optim.object) which can be plotted to show the bandwidth selection criterion as a function of sigma.

Author(s)

Greg McSwiggan and Adrian Baddeley [email protected].

References

Kelsall, J.E. and Diggle, P.J. (1995) Kernel estimation of relative risk. Bernoulli 1, 3–16.

McSwiggan, G., Baddeley, A. and Nair, G. (2019) Estimation of relative risk for events on a linear network. Statistics and Computing 30 (2) 469–484.

See Also

relrisk.lpp, bw.relrisk, bw.optim.object

Examples

set.seed(2020)
   X <- superimpose(A=runiflpp(20, simplenet),
                    B=runifpointOnLines(20, as.psp(simplenet)[1]))
   plot(bw.relrisk(X, hmin=0.13, hmax=0.22, method="McSwiggan"))
   plot(bw.relrisk(X, hmin=0.1, hmax=0.2, nh=8, distance="euclidean"))

Description

Uses cross-validation to select a smoothing bandwidth for the Voronoi estimate of point process intensity on a linear network.

Usage

bw.voronoi(X, ..., probrange = c(0.2, 0.8), nprob = 10,
           prob = NULL, nrep = 100, verbose = TRUE, warn=TRUE)

Arguments

Details

This function uses likelihood cross-validation to choose the optimal value of the thinning fraction f (the retention probability) to be used in the smoothed Voronoi estimator of point process intensity densityVoronoi.lpp.

Value

A single numerical value giving the selected bandwidth. The result also belongs to the class "bw.optim" (see bw.optim.object) which can be plotted to show the bandwidth selection criterion as a function of sigma.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected] and Mehdi Moradi [email protected].

References

Moradi, M., Cronie, 0., Rubak, E., Lachieze-Rey, R., Mateu, J. and Baddeley, A. (2019) Resample-smoothing of Voronoi intensity estimators. Statistics and Computing 29 (5) 995–1010.

See Also

densityVoronoi.lpp, bw.optim.object

Examples

np <- if(interactive()) 10 else 3
   nr <- if(interactive()) 100 else 2
   b <- bw.voronoi(spiders, nprob=np, nrep=nr)
   b
   plot(b)

Description

Performs a test of goodness-of-fit of a point process model on a linear network. The observed and predicted distributions of the values of a spatial covariate are compared using either the Kolmogorov-Smirnov test, Cramer-von Mises test or Anderson-Darling test. For non-Poisson models, a Monte Carlo test is used.

Usage

## S3 method for class 'lpp'
cdf.test(X, covariate,  test=c("ks", "cvm", "ad"), ...,
        interpolate=TRUE, jitter=TRUE)

## S3 method for class 'lppm'
cdf.test(model, covariate,  test=c("ks", "cvm", "ad"), ...,
          interpolate=TRUE, jitter=TRUE, nsim=99, verbose=TRUE)

Arguments

Details

These functions perform a goodness-of-fit test of a Poisson point process model fitted to point pattern data on a linear network. The observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same values under the model, are compared using the Kolmogorov-Smirnov test, the Cramer-von Mises test or the Anderson-Darling test. For Gibbs models, a Monte Carlo test is performed using these test statistics.

The function cdf.test is generic, with methods for point patterns ("ppp" or "lpp"), point process models ("ppm" or "lppm") and spatial logistic regression models ("slrm").

See the help file for cdf.test for information on the generic function and the methods for data in two-dimensional space, classes "ppp", "ppm" and "slrm".

This help file describes the methods for data on a linear network, classes "lpp" and "lppm".

• If X is a point pattern on a linear network (object of class "lpp"), then cdf.test(X, ...) performs a goodness-of-fit test of the uniform Poisson point process (Complete Spatial Randomness, CSR) for this dataset. For a multitype point pattern, the uniform intensity is assumed to depend on the type of point (sometimes called Complete Spatial Randomness and Independence, CSRI).

• If model is a fitted point process model on a network (object of class "lppm") then cdf.test(model, ...) performs a test of goodness-of-fit for this fitted model.

The test is performed by comparing the observed distribution of the values of a spatial covariate at the data points, and the predicted distribution of the same covariate under the model, using a classical goodness-of-fit test. Thus, you must nominate a spatial covariate for this test.

If X is a point pattern that does not have marks, the argument covariate should be either a function(x,y) or a pixel image (object of class "im" or "linim") containing the values of a spatial function, or one of the characters "x" or "y" indicating the Cartesian coordinates. If covariate is an image, it should have numeric values, and its domain should cover the observation window of the model. If covariate is a function, it should expect two arguments x and y which are vectors of coordinates, and it should return a numeric vector of the same length as x and y.

If X is a multitype point pattern, the argument covariate can be either a function(x,y,marks), or a pixel image, or a list of pixel images corresponding to each possible mark value, or one of the characters "x" or "y" indicating the Cartesian coordinates.

First the original data point pattern is extracted from model. The values of the covariate at these data points are collected.

The predicted distribution of the values of the covariate under the fitted model is computed as follows. The values of the covariate at all locations in the observation window are evaluated, weighted according to the point process intensity of the fitted model, and compiled into a cumulative distribution function $F$ using ewcdf.

The probability integral transformation is then applied: the values of the covariate at the original data points are transformed by the predicted cumulative distribution function $F$ into numbers between 0 and 1. If the model is correct, these numbers are i.i.d. uniform random numbers. The A goodness-of-fit test of the uniform distribution is applied to these numbers using stats::ks.test, goftest::cvm.test or goftest::ad.test.

This test was apparently first described (in the context of two-dimensional spatial data, and using Kolmogorov-Smirnov) by Berman (1986). See also Baddeley et al (2005).

If model is not a Poisson process, then a Monte Carlo test is performed, by generating nsim point patterns which are simulated realisations of the model, re-fitting the model to each simulated point pattern, and calculating the test statistic for each fitted model. The Monte Carlo $p$ value is determined by comparing the simulated values of the test statistic with the value for the original data.

The return value is an object of class "htest" containing the results of the hypothesis test. The print method for this class gives an informative summary of the test outcome.

The return value also belongs to the class "cdftest" for which there is a plot method plot.cdftest. The plot method displays the empirical cumulative distribution function of the covariate at the data points, and the predicted cumulative distribution function of the covariate under the model, plotted against the value of the covariate.

The argument jitter controls whether covariate values are randomly perturbed, in order to avoid ties. If the original data contains any ties in the covariate (i.e. points with equal values of the covariate), and if jitter=FALSE, then the Kolmogorov-Smirnov test implemented in ks.test will issue a warning that it cannot calculate the exact $p$-value. To avoid this, if jitter=TRUE each value of the covariate will be perturbed by adding a small random value. The perturbations are normally distributed with standard deviation equal to one hundredth of the range of values of the covariate. This prevents ties, and the $p$-value is still correct. There is a very slight loss of power.

Value

An object of class "htest" containing the results of the test. See ks.test for details. The return value can be printed to give an informative summary of the test.

The value also belongs to the class "cdftest" for which there is a plot method.

Warning

The outcome of the test involves a small amount of random variability, because (by default) the coordinates are randomly perturbed to avoid tied values. Hence, if cdf.test is executed twice, the $p$-values will not be exactly the same. To avoid this behaviour, set jitter=FALSE.

Author(s)

Adrian Baddeley [email protected] and Rolf Turner [email protected]

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.

Berman, M. (1986) Testing for spatial association between a point process and another stochastic process. Applied Statistics 35, 54–62.

See Also

plot.cdftest, quadrat.test, berman.test, ks.test, goftest::cvm.test, goftest::ad.test, lppm

Examples

op <- options(useFancyQuotes=FALSE)

   # test of CSR using x coordinate
   cdf.test(spiders, "x")

   # fit inhomogeneous Poisson model and test
   model <- lppm(spiders ~x)
   cdf.test(model, "y")

   # test of CSR using a function of x and y
   fun <- function(x,y){2* x + y}
   cdf.test(spiders, fun)

   # test of CSR using an image covariate
   fim <- as.linim(fun, domain(spiders))
   cdf.test(spiders, fim)

   options(op)

Description

Given a linear network and a set of infinite lines, divide the network into tiles demarcated by the lines. The result is a tessellation of the network.

Usage

chop.linnet(X, L)

Arguments

Details

The first line of L divides X into two tiles. Subsequent lines divide each of these tiles. The result is a tessellation of X. Tiles are not necessarily connected sets.

Value

Tessellation on a linear network (object of class "lintess").

Author(s)

Adrian Baddeley [email protected].

See Also

crossing.linnet to determine the crossing points between the lines and the network.

divide.linnet to divide a network into a tessellation using arbitrary cut points.

Examples

L <- infline(p=runif(3), theta=runif(3, max=pi/2))
   Y <- chop.linnet(simplenet, L)
   plot(Y, main="")
   plot(L, col="red")

Description

Given a point pattern representing a set of vertices, this command gives a point-and-click interface allowing the user to join pairs of selected vertices by edges.

Usage

clickjoin(X, ..., add = TRUE, m = NULL, join = TRUE)

Arguments

Details

This function makes it easier for the user to create a linear network or a planar graph, given a set of vertices.

The function first displays the point pattern X, then repeatedly prompts the user to click on a pair of points in X. Each selected pair of points will be joined by an edge. The function returns a logical matrix which has entries equal to TRUE for each pair of vertices joined by an edge.

The selection of points is performed using identify.ppp which typically expects the user to click the left mouse button. This point-and-click interaction continues until the user terminates it, by pressing the middle mouse button, or pressing the right mouse button and selecting stop.

The return value can be used in linnet to create a linear network.

Value

Logical matrix m with value m[i,j] = TRUE for every pair of vertices X[i] and X[j] that should be joined by an edge.

Author(s)

Adrian Baddeley [email protected].

See Also

linnet, clickppp

Description

Allows the user to create a point pattern on a linear network by point-and-click in the display.

Usage

clicklpp(L, n=NULL, types=NULL, ...,
           add=FALSE, main=NULL, hook=NULL)

Arguments

Details

This function allows the user to create a point pattern on a linear network by interactively clicking on the screen display.

First the linear network L is plotted on the current screen device. Then the user is prompted to point the mouse at any desired locations and click the left mouse button to add each point. Interactive input stops after n clicks (if n was given) or when the middle mouse button is pressed.

The return value is a point pattern on the network L, containing the locations of all the clicked points, after they have been projected onto the network L. Any points that were clicked outside the bounding window of the network will be ignored.

If the argument types is given, then a multitype point pattern will be created. The user is prompted to input the locations of points of type type[i], for each successive index i. (If the argument n was given, there will be n points of each type.) The return value is a multitype point pattern on a linear network.

This function uses the R command locator to input the mouse clicks. It only works on screen devices such as ‘X11’, ‘windows’ and ‘quartz’. Arguments that can be passed to locator through ... include pch (plotting character), cex (character expansion factor) and col (colour). See locator and par.

Value

A point pattern (object of class "lpp").

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected], based on an idea by Dominic Schuhmacher.

See Also

clickppp, identify.lpp, locator, clickpoly, clickbox, clickdist

Description

Find the topologically-connected components of a linear network.

Usage

## S3 method for class 'linnet'
connected(X, ..., what = c("labels", "components"))

Arguments

Details

The function connected is generic. This is the method for linear networks (objects of class "linnet").

Two vertices of the network are connected if they are joined by a path in the network. This function divides the network into subsets, such that all points in a subset are connected to each other.

If what="labels" the return value is a factor with one entry for each vertex of X, identifying which connected component the vertex belongs to.

If what="components" the return value is a list of linear networks, which are the connected components of X.

Value

If what="labels", a factor. If what="components", a list of linear networks.

Author(s)

Adrian Baddeley [email protected] and Suman Rakshit.

See Also

thinNetwork

Examples

# remove some edges from a network to make it disconnected
   plot(simplenet, col="grey", main="", lty=2)
   A <- thinNetwork(simplenet, retainedges=-c(3,5))
   plot(A, add=TRUE, lwd=2)
   # find the connected components
   connected(A)
   cA <- connected(A, what="components")
   plot(cA[[1]], add=TRUE, col="green", lwd=2)
   plot(cA[[2]], add=TRUE, col="blue", lwd=2)

Description

Finds the topologically-connected components of a point pattern on a linear network, when all pairs of points closer than a threshold distance are joined.

Usage

## S3 method for class 'lpp'
connected(X, R=Inf, ..., dismantle=TRUE)

Arguments

Details

The function connected is generic. This is the method for point patterns on a linear network (objects of class "lpp"). It divides the point pattern X into one or more groups of points.

If R=Inf (the default), then X is divided into groups such that any pair of points in the same group can be joined by a path in the network.

If R is a finite number, then two points of X are declared to be R-close if they lie closer than R units apart, measured by the length of the shortest path in the network. Two points are R-connected if they can be reached by a series of steps between R-close pairs of points of X. Then X is divided into groups such that any pair of points in the same group is R-connected.

If dismantle=TRUE (the default) the algorithm first checks whether the network is connected (i.e. whether any pair of vertices can be joined by a path in the network), and if not, the network is decomposed into its connected components.

Value

A point pattern (of class "lpp") with marks indicating the grouping, or a list of such point patterns.

Author(s)

Adrian Baddeley [email protected].

See Also

thinNetwork

Examples

## behaviour like connected.ppp
   U <- runiflpp(20, simplenet)
   plot(connected(U, 0.15, dismantle=FALSE))

   ## behaviour like connected.owin
   ## remove some edges from a network to make it disconnected
   plot(simplenet, col="grey", main="", lty=2)
   A <- thinNetwork(simplenet, retainedges=-c(3,5))
   plot(A, add=TRUE, lwd=2)
   X <- runiflpp(10, A)
   ## find the connected components
   cX <- connected(X)
   plot(cX[[1]], add=TRUE, col="blue", lwd=2)

Description

Computes the distances between pairs of points taken from two different point patterns on the same linear network.

Usage

## S3 method for class 'lpp'
crossdist(X, Y, ..., method="C", check=TRUE)

Arguments

Details

Given two point patterns on a linear network, this function computes the distance from each point in the first pattern to each point in the second pattern, measuring distance by the shortest path along the network.

This is a method for the generic function crossdist for the class of point patterns on a linear network (objects of class "lpp").

This function expects two point pattern objects X and Y on the same linear network, and returns the matrix whose [i,j] entry is the shortest-path distance from X[i] to Y[j].

If two points cannot be joined by a path, the distance between them is infinite (Inf).

The argument method is not normally used. It is retained only for developers to check the validity of the software.

Value

A matrix whose [i,j] entry is the distance from the i-th point in X to the j-th point in Y. Matrix entries are nonnegative numbers or infinity (Inf).

Algorithms and accuracy

Distances are accurate within the numerical tolerance of the network, summary(X)$toler.

For network data stored in the non-sparse representation described in linnet, then pairwise distances are computed using the matrix of path distances between vertices of the network, using R code if method = "interpreted", or using C code if method="C" (the default).

For networks stored in the sparse representation, the argument method has no effect, and the distances are computed using an efficient C algorithm.

Author(s)

Adrian Baddeley [email protected].

See Also

crossdist, crossdist.ppp, pairdist, nndist

Examples

v <- split(chicago)
   X <- v$cartheft
   Y <- v$burglary
   d <- crossdist(X, Y)
   d[1:3,1:4]

Description

Find all the crossing-points between a linear network and another pattern of lines or line segments.

Usage

crossing.linnet(X, Y)

Arguments

Details

All crossing-points between X and Y are determined. The result is a point pattern on the network X.

Value

Point pattern on a linear network (object of class "lpp").

Author(s)

Adrian Baddeley [email protected].

See Also

crossing.psp

Examples

plot(simplenet, main="")
   L <- infline(p=runif(3), theta=runif(3, max=pi/2))
   plot(L, col="red")
   Y <- crossing.linnet(simplenet, L)
   plot(Y, add=TRUE, cols="blue")

Description

For a point pattern on a linear network, classify the points into distinct types according to the numerical marks in the pattern, or according to another variable.

Usage

## S3 method for class 'lpp'
cut(x, z=marks(x), ...)

Arguments

Details

This function has the effect of classifying each point in the point pattern x into one of several possible types. The classification is based on the dataset z, which may be either

• a factor (of length equal to the number of points in z) determining the classification of each point in x. Levels of the factor determine the classification.

• a numeric vector (of length equal to the number of points in z). The range of values of z will be divided into bands (the number of bands is determined by ...) and z will be converted to a factor using cut.default.

• a pixel image on a network (object of class "linim"). The value of z at each point of x will be used as the classifying variable.

• a function on a network (object of class "linfun", see linfun). The value of z at each point of x will be used as the classifying variable.

• a tessellation on a network (object of class "lintess", see lintess). Each point of x will be classified according to the tile of the tessellation into which it falls.

• a character string, giving the name of one of the columns of marks(x), if this is a data frame.

• a character string identifying one of the coordinates: the spatial coordinates "x", "y" or the segment identifier "seg" or the fractional coordinate along the segment, "tp".

The default is to take z to be the vector of marks in x (or the first column in the data frame of marks of x, if it is a data frame). If the marks are numeric, then the range of values of the numerical marks is divided into several intervals, and each interval is associated with a level of a factor. The result is a marked point pattern, on the same linear network, with the same point locations as x, but with the numeric mark of each point discretised by replacing it by the factor level. This is a convenient way to transform a marked point pattern which has numeric marks into a multitype point pattern, for example to plot it or analyse it. See the examples.

To select some points from x, use the subset operators [.lpp or subset.lpp instead.

Value

A multitype point pattern on the same linear network, that is, a point pattern object (of class "lpp") with a marks vector that is a factor.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

cut, lpp, lintess, linfun, linim

Examples

X <- runiflpp(20, simplenet)
  f <- linfun(function(x,y,seg,tp) { x }, simplenet)
  plot(cut(X, f, breaks=4))
  plot(cut(X, "x", breaks=4))
  plot(cut(X, "seg"))

Description

Given a fitted point process model on a linear network, this function extracts the original point pattern dataset to which the model was fitted.

Usage

data.lppm(object)

Arguments

Details

An object of class "lppm" represents a point process model that has been fitted to a point pattern dataset on a linear network. It is typically produced by the model-fitting algorithm lppm. The object contains complete information about the original data point pattern to which the model was fitted. This function extracts the original data pattern.

Value

A point pattern on a linear network (object of class "lpp").

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

lppm, data.ppm

Examples

fit <- lppm(spiders ~ x)
 X <- data.lppm(fit)
 # 'X' is identical to 'spiders'

Description

Computes the edges of the Delaunay triangulation or Dirichlet tessellation of a point pattern, and returns the result as a linear network object.

Usage

delaunayNetwork(X)

dirichletNetwork(X, ...)

Arguments

Details

For delaunayNetwork, points of X which are neighbours in the Delaunay triangulation (see delaunay) will be joined by a straight line. The result will be returned as a linear network (object of class "linnet").

For dirichletNetwork, the Dirichlet tessellation is computed (see dirichlet) and the edges of the tiles of the tessellation are extracted. This is converted to a linear network using as.linnet.psp.

Value

Linear network (object of class "linnet") or NULL.

Author(s)

Adrian Baddeley [email protected]

Rolf Turner [email protected]

and Ege Rubak [email protected]

See Also

delaunay, dirichlet, delaunayDistance

Examples

LE <- delaunayNetwork(cells)
  LI <- dirichletNetwork(cells)

Description

Deletes or extracts a given branch of a tree.

Usage

deletebranch(X, ...)

## S3 method for class 'linnet'
deletebranch(X, code, labels, ...)

## S3 method for class 'lpp'
deletebranch(X, code, labels, ...)

extractbranch(X, ...)

## S3 method for class 'linnet'
extractbranch(X, code, labels, ..., which=NULL)

## S3 method for class 'lpp'
extractbranch(X, code, labels, ..., which=NULL)

Arguments

Details

The linear network L <- X or L <- as.linnet(X) must be a tree, that is, it has no loops.

The argument labels should be a character vector giving tree branch labels for each vertex of the network. It is usually obtained by calling treebranchlabels.

The branch designated by the string code will be deleted or extracted.

The return value is the result of deleting or extracting this branch from X along with any data associated with this branch (such as points or marks).

Value

Another object of the same type as X obtained by deleting or extracting the specified branch.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected]

See Also

treebranchlabels, branchlabelfun, linnet

Examples

# make a simple tree
  m <- simplenet$m
  m[8,10] <- m[10,8] <- FALSE
  L <- linnet(vertices(simplenet), m)
  plot(L, main="")
  # compute branch labels 
  tb <- treebranchlabels(L, 1)
  tbc <- paste0("[", tb, "]")
  text(vertices(L), labels=tbc, cex=2)

  # delete branch B
  LminusB <- deletebranch(L, "b", tb)
  plot(LminusB, add=TRUE, col="green")

  # extract branch B
  LB <- extractbranch(L, "b", tb)
  plot(LB, add=TRUE, col="red")

Description

Compute a kernel smoothed intensity function for the line segments of a linear network.

Usage

## S3 method for class 'linnet'
density(x, ...)

Arguments

Details

This is the method for the generic function density for the class "linnet" (linear networks).

The network x is first converted to a line segment pattern (object of class "psp"). Then the method density.psp is applied to the segment pattern.

A kernel estimate of the intensity of the line segment pattern is computed. The result is the convolution of the isotropic Gaussian kernel, of standard deviation sigma, with the line segments.

The intensity of a line segment pattern is the (spatially-varying) amount of segment length per unit area, expressed in the same units as the coordinates of x. If the units of x are in metres, then an intensity value of 3 means that there are 3 metres of segment length per square metre of spatial domain.

See density.psp for more details.

Value

A pixel image in two dimensions (object of class "im") or a numeric vector.

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

See Also

density.psp, im.object, density.

Examples

D <- density(simplenet, 0.1)
  plot(D)
  plot(simplenet, add=TRUE, col="white")
  ## compare with average intensity
  volume(simplenet)/area(Window(simplenet))

Description

Estimates the intensity of a point process on a linear network by applying kernel smoothing to the point pattern data.

Usage

## S3 method for class 'lpp'
density(x, sigma=NULL, ...,
        weights=NULL,
        distance=c("path", "euclidean"),
        continuous=TRUE,
        kernel="gaussian")

## S3 method for class 'splitppx'
density(x, sigma=NULL, ...)

Arguments

Details

Kernel smoothing is applied to the points of x using either a kernel based on path distances in the network, or a two-dimensional kernel. The result is a pixel image on the linear network (class "linim") which can be plotted.

• If distance="path" (the default) then the smoothing is performed using a kernel based on path distances in the network, as described in described in Okabe and Sugihara (2012) and McSwiggan et al (2016).

  • If continuous=TRUE (the default), smoothing is performed using the “equal-split continuous” rule described in Section 9.2.3 of Okabe and Sugihara (2012). The resulting function is continuous on the linear network.

  • If continuous=FALSE, smoothing is performed using the “equal-split discontinuous” rule described in Section 9.2.2 of Okabe and Sugihara (2012). The resulting function is continuous except at the network vertices.

  • In the default case (where distance="path" and continuous=TRUE and kernel="gaussian", computation is performed rapidly by solving the classical heat equation on the network, as described in McSwiggan et al (2016). The arguments are passed to densityHeat.lpp which performs the computation. Computational time is short, but increases quadratically with sigma.

  • In all other cases, computation is performed by path-tracing as described in Okabe and Sugihara (2012); the arguments are passed to densityEqualSplit which performs the computation. Computation time can be extremely long, and increases exponentially with sigma.

• If distance="euclidean", the smoothing is performed using a two-dimensional kernel. The arguments are passed to densityQuick.lpp to perform the computation. Computation time is very short. See the help for densityQuick.lpp for further details.

There is also a method for split point patterns on a linear network (class "splitppx") which will return a list of pixel images.

The argument sigma specifies the smoothing bandwidth. If sigma is missing or NULL, the default is one-eighth of the length of the shortest side of the bounding box of x. If sigma is a function in the R language, it is assumed to be a bandwidth selection rule, and it will be applied to x to compute the bandwidth value.

Value

A pixel image on the linear network (object of class "linim"), or in some cases, a numeric vector of length equal to npoints(x).

Infinite bandwidth

If sigma=Inf, the resulting density estimate is constant over all locations, and is equal to the average density of points per unit length. (If the network is not connected, then this rule is applied separately to each connected component of the network).

Author(s)

Adrian Baddeley [email protected] and Greg McSwiggan.

References

McSwiggan, G., Baddeley, A. and Nair, G. (2016) Kernel density estimation on a linear network. Scandinavian Journal of Statistics 44, 324–345.

Okabe, A. and Sugihara, K. (2012) Spatial analysis along networks. Wiley.

See Also

lpp, linim, densityQuick.lpp, densityHeat.lpp, densityVoronoi.lpp

Examples

X <- runiflpp(3, simplenet)
  D <- density(X, 0.2, verbose=FALSE)
  plot(D, style="w", main="", adjust=2)
  Dq <- density(X, 0.2, distance="euclidean")
  plot(Dq, style="w", main="", adjust=2)
  Dw <- density(X, 0.2, weights=c(1,2,-1), verbose=FALSE)
  De <- density(X, 0.2, kernel="epanechnikov", verbose=FALSE)
  Ded <- density(X, 0.2, kernel="epanechnikov", continuous=FALSE, verbose=FALSE)

Description

Computes a kernel density estimate on a linear network using the Okabe-Sugihara equal-split algorithms.

Usage

densityEqualSplit(x, sigma = NULL, ...,
                   at = c("pixels", "points"),
                   leaveoneout=TRUE,
                   weights = NULL,
                   kernel = "epanechnikov", continuous = TRUE,
                   epsilon = 1e-06, verbose = TRUE, debug = FALSE, savehistory = TRUE)

Arguments

Details

Kernel smoothing is applied to the points of x using a kernel based on path distances in the network. The result is a pixel image on the linear network (class "linim") which can be plotted.

Smoothing is performed using one of the “equal-split” rules described in Okabe and Sugihara (2012).

• If continuous=TRUE (the default), smoothing is performed using the “equal-split continuous” rule described in Section 9.2.3 of Okabe and Sugihara (2012). The resulting function is continuous on the linear network.

• If continuous=FALSE, smoothing is performed using the “equal-split discontinuous” rule described in Section 9.2.2 of Okabe and Sugihara (2012). The resulting function is not continuous.

Computation is performed by path-tracing as described in Okabe and Sugihara (2012).

It is advisable to choose a kernel with bounded support such as kernel="epanechnikov". With a Gaussian kernel, computation time can be long, and increases exponentially with sigma.

Faster algorithms are available through density.lpp.

The argument sigma specifies the smoothing bandwidth. If sigma is missing or NULL, the default is one-eighth of the length of the shortest side of the bounding box of x. If sigma is a function in the R language, it is assumed to be a bandwidth selection rule, and it will be applied to x to compute the bandwidth value.

Value

If at="pixels" (the default), a pixel image on the linear network (object of class "linim").

If at="points", a numeric vector with one entry for each point of x.

Infinite bandwidth

If sigma=Inf, the resulting density estimate is constant over all locations, and is equal to the average density of points per unit length. (If the network is not connected, then this rule is applied separately to each connected component of the network).

Author(s)

Adrian Baddeley [email protected] and Greg McSwiggan.

References

Okabe, A. and Sugihara, K. (2012) Spatial analysis along networks. Wiley.

See Also

density.lpp

Examples

X <- runiflpp(3, simplenet)
  De <- density(X, 0.2, kernel="epanechnikov", verbose=FALSE)
  Ded <- density(X, 0.2, kernel="epanechnikov", continuous=FALSE, verbose=FALSE)

Description

Computes a kernel estimate of the intensity of a point process on a linear network, and returns the intensity estimate as a function of spatial location.

Usage

## S3 method for class 'lpp'
densityfun(X, sigma, ..., weights=NULL, nsigma=1, verbose=FALSE)

Arguments

Details

Kernel smoothing is applied to the points of X using the diffusion algorithm of McSwiggan et al (2016). The result is a function on the linear network (object of class "linfun") that can be printed, plotted and evaluated at any location.

This is a method for the generic function densityfun for the class "lpp" of point patterns on a linear network.

Value

Function on a linear network (object of class "linfun").

If nsigma=1 (the default), the result is a function giving kernel estimate with bandwidth sigma.

If nsigma > 1, the result is a function with an additional argument k. If k is specified, the function returns the kernel estimate for bandwidth tau = sigma * sqrt(k/nsigma). If k is not specified, results are returned for all k = 1, 2, ..., nsigma.

The result also has attributes

• attr(result, "dt") giving the time step $\Delta t$;

• attr(result, "dx") giving the spacing $\Delta x$ between sample points in the numerical algorithm;

• attr(result, "sigma") giving the smoothing bandwidth $\sigma$ used (or the successive bandwidths used at each sampled time step, if nsigma > 1).

Author(s)

Greg McSwiggan, with tweaks by Adrian Baddeley [email protected].

References

McSwiggan, G., Baddeley, A. and Nair, G. (2016) Kernel Density Estimation on a Linear Network. Scandinavian Journal of Statistics 44, 324–345.

See Also

density.lpp which returns a pixel image on the linear network.

methods.linfun for methods applicable to "linfun" objects.

Examples

X <- unmark(chicago)
  # single bandwidth
  g <- densityfun(X, 30)
  plot(g)
  Y <- X[1:5]
  g(Y)
  # weighted
  gw <- densityfun(X, 30, weights=runif(npoints(X)))
  # sequence of bandwidths 
  g10 <- densityfun(X, 30, nsigma=10)
  g10(Y, k=10)
  g10(Y)
  plot(as.linim(g10, k=5))

Description

Given a point pattern on a linear network, compute a kernel estimate of intensity, by solving the heat equation.

Usage

## S3 method for class 'lpp'
densityHeat(x, sigma=NULL, ...,
              at=c("pixels", "points"), leaveoneout=TRUE,
              weights = NULL,
              dx = NULL, dt = NULL, iterMax = 1e+06,
              finespacing = TRUE, verbose=FALSE)

Arguments

Details

The function densityHeat is generic. This is the method for the class "lpp" of points on a linear network.

Kernel smoothing is applied to the points of x using a kernel based on path distances in the network. If at="pixels" (the default), the result is a pixel image on the linear network (class "linim") which can be plotted. If at="points" the result is a numeric vector giving the density estimates at the data points of x.

The smoothing operation is equivalent to the “equal-split continuous” rule described in Section 9.2.3 of Okabe and Sugihara (2012). However, the actual computation is performed rapidly, by solving the classical time-dependent heat equation on the network, as described in McSwiggan et al (2016). Computational time is short, but increases quadratically with sigma.

If at="points" and leaveoneout=TRUE, a leave-one-out estimate is computed at each data point (that is, the estimate at each data point x[i] is based on all of the points except x[i]) using the truncated series approximation of McSwiggan et al (2019).

The argument sigma specifies the smoothing bandwidth. If sigma is missing or NULL, the default is one-eighth of the length of the shortest side of the bounding box of x. If sigma is a function in the R language, it is assumed to be a bandwidth selection rule, and it will be applied to x to compute the bandwidth value.

Value

If at="pixels" (the default), a pixel image on the linear network (object of class "linim").

If at="points", a numeric vector with one entry for each point of x.

Infinite bandwidth

If sigma=Inf, the resulting density estimate is constant over all locations, and is equal to the average density of points per unit length. (If the network is not connected, then this rule is applied separately to each connected component of the network).

Discretisation and Error Messages

The arguments dx, dt and iterMax determine the discretisation of the network, according to a set of rules. The argument finespacing determines which rule will be applied.

The arguments dx, dt, iterMax are connected by several constraints; specifying one of these arguments will affect the default values of the other two arguments.

The argument finespacing specifies whether a very fine spacing of sample points is required, in order to attain hihg accuracy.

• If finespacing=TRUE (the default), then the sample point spacing dx must not exceed one-third of the length of the shortest segment of the network. This ensures that the discrete approximation is accurate along every segment, no matter how short the segment is. However, this may not be feasible if it implies a very large number of sample points, or a large number of iterations: in such cases, the code may terminate with an error about illegal values of dx, dt or iterMax.

• If finespacing=FALSE, then the sample point spacing dx will be about one-half the width of a pixel in the default pixellation of the window of x. This is usually a much coarser resolution than the one selected by finespacing=TRUE. If it is too coarse, the pixel resolution can be refined using the arguments dimyx, eps or xy passed to as.mask. For example, dimyx=512 would specify a 512 x 512 pixel grid. The default pixel resolution can be changed for the remainder of the R session by spatstat.options('npixel').

Author(s)

Adrian Baddeley [email protected] and Greg McSwiggan.

References

McSwiggan, G., Baddeley, A. and Nair, G. (2016) Kernel density estimation on a linear network. Scandinavian Journal of Statistics 44, 324–345.

McSwiggan, G., Baddeley, A. and Nair, G. (2019) Estimation of relative risk for events on a linear network. Statistics and Computing 30, 469–484.

Okabe, A. and Sugihara, K. (2012) Spatial analysis along networks. Wiley.

See Also

density.lpp

Examples

X <- runiflpp(3, simplenet)
  D <- densityHeat(X, 0.2)
  plot(D, style="w", main="", adjust=2)
  densityHeat.lpp(X, 0.2, at="points")
  Dw <- densityHeat(X, 0.2, weights=c(1,2,-1))

Description

Estimates the intensity of a point process on a linear network using a two-dimensional smoothing kernel.

Usage

densityQuick.lpp(x, sigma=NULL, ...,
         kernel="gaussian",
         at = c("pixels", "points"),
         what = c("estimate", "se", "var"),
         leaveoneout = TRUE,
         diggle = FALSE,
         edge2D = FALSE,
         weights = NULL,
         positive = FALSE)

Arguments

Details

Kernel smoothing is applied to the points of x using a two-dimensional Gaussian kernel, as described in Rakshit et al (2019). The result is a pixel image on the linear network (class "linim") which can be plotted.

Other techniques for kernel smoothing on a network are implemented in density.lpp. The main advantages of using a two-dimensional kernel are very fast computation and insensitivity to changes in the network geometry. The main disadvantage is that it ignores the connectivity of the network. See Rakshit et al (2019) for further explanation.

The argument sigma specifies the smoothing bandwidth. If sigma is missing or NULL, the default is one-eighth of the length of the shortest side of the bounding box of x. If sigma is a function in the R language, it is assumed to be a bandwidth selection rule, and it will be applied to x to compute the bandwidth value.

Value

If at="pixels" (the default), a pixel image on the linear network (object of class "linim").

If at="points", a numeric vector with one entry for each point of x.

Infinite bandwidth

If sigma=Inf, the resulting density estimate is constant over all locations, and is equal to the average density of points per unit length. (If the network is not connected, then this rule is applied separately to each connected component of the network).

Author(s)

Adrian Baddeley, Suman Rakshit and Tilman Davies

References

Rakshit, S., Davies, T., Moradi, M., McSwiggan, G., Nair, G., Mateu, J. and Baddeley, A. (2019) Fast kernel smoothing of point patterns on a large network using 2D convolution. International Statistical Review 87 (3) 531–556. DOI: 10.1111/insr.12327.

See Also

density.lpp, the main function for density estimation on a network.

bw.scott, bw.scott.iso, bw.lpplfor bandwidth selection.

Examples

X <- unmark(chicago)
   plot(densityQuick.lpp(X, 500))
   plot(densityQuick.lpp(X, 500, diggle=TRUE))
   plot(densityQuick.lpp(X, bw.scott.iso))
   plot(densityQuick.lpp(X, 500, what="se"))

Description

Computes an adaptive estimate of the intensity function of a point pattern on a linear network, using the Dirichlet-Voronoi tessellation on the network.

Usage

## S3 method for class 'lpp'
densityVoronoi(X, f = 1, ..., nrep = 1, verbose = TRUE)

Arguments

Details

This function is an alternative to density.lpp. It computes an estimate of the intensity function of a point pattern dataset on a linear network. The result is a pixel image on the network, giving the estimated intensity.

This function is a method for the generic densityVoronoi for the class "lpp" of point patterns on a linear network.

If f=1 (the default), the Voronoi estimate (Barr and Schoenberg, 2010) is computed: the point pattern X is used to construct a Voronoi/Dirichlet tessellation on the network (see lineardirichlet); the lengths of the Dirichlet tiles are computed; the estimated intensity in each tile is the reciprocal of the tile length. The result is a pixel image of intensity estimates which are constant on each tile of the tessellation.

If f=0, the intensity estimate at every location is equal to the average intensity (number of points divided by network length). The result is a pixel image of intensity estimates which are constant.

If f is strictly between 0 and 1, the smoothed Voronoi estimate (Moradi et al, 2019) is computed. The dataset X is randomly thinned by deleting or retaining each point independently, with probability f of retaining a point. The thinned pattern is used to construct a Dirichlet tessellation and form the Voronoi estimate, which is then adjusted by a factor 1/f. This procedure is repeated nrep times and the results are averaged to obtain the smoothed Voronoi estimate.

The value f can be chosen automatically by bandwidth selection using bw.voronoi.

Value

Pixel image on a linear network (object of class "linim").

Author(s)

Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected] and Mehdi Moradi [email protected].

References

Moradi, M., Cronie, 0., Rubak, E., Lachieze-Rey, R., Mateu, J. and Baddeley, A. (2019) Resample-smoothing of Voronoi intensity estimators. Statistics and Computing 29 (5) 995–1010.

See Also

densityVoronoi is the generic, with a method for class "ppp".

lineardirichlet computes the Dirichlet-Voronoi tessellation on a network.

bw.voronoi performs bandwidth selection of the fraction f.

See also density.lpp.

Examples

nr <- if(interactive()) 100 else 3
   plot(densityVoronoi(spiders, 0.1, nrep=nr))

Description

Compute the diameter or bounding radius of a linear network measured using the shortest path distance.

Usage

## S3 method for class 'linnet'
diameter(x)

## S3 method for class 'linnet'
boundingradius(x, ...)