Package 'spatstat.explore'

Description

The spatstat.explore package belongs to the spatstat family of packages. It contains the core functionality for statistical analysis and modelling of spatial data.

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.explore contains the user-level functions that perform exploratory data analysis and nonparametric data analysis of spatial data.

(The main exception is that functions for linear networks are in the separate sub-package spatstat.linnet.)

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 functionality for exploratory data analysis and nonparametric analysis of spatial data.

• spatstat.model containing the functionality for statistical modelling, model-fitting, formal statistical inference and informal model diagnostics.

• 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 the sub-packages of spatstat, see the help file for "spatstat-package" 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 Functionality in spatstat.explore

The spatstat family of packages is designed to support a complete statistical analysis of spatial data. It supports

• creation, manipulation and plotting of point patterns;

• exploratory data analysis;

• spatial random sampling;

• simulation of point process models;

• parametric model-fitting;

• non-parametric smoothing and regression;

• formal inference (hypothesis tests, confidence intervals);

• model diagnostics.

For an overview, see the help file for "spatstat-package" in the spatstat package.

Following is a list of the functionality provided in the spatstat.explore package only.

To simulate a random point pattern:

Functions for generating random point patterns are now contained in the spatstat.random package.

To interrogate a point pattern:

Manipulation of pixel images:

An object of class "im" represents a pixel image.

Line segment patterns

An object of class "psp" represents a pattern of straight line segments.

Tessellations

An object of class "tess" represents a tessellation.

Three-dimensional point patterns

An object of class "pp3" represents a three-dimensional point pattern in a rectangular box. The box is represented by an object of class "box3".

Multi-dimensional space-time point patterns

An object of class "ppx" represents a point pattern in multi-dimensional space and/or time.

Classical exploratory tools:

Smoothing:

Modern exploratory tools:

Summary statistics for a point pattern:

Related facilities:

Summary statistics for a multitype point pattern: A multitype point pattern is represented by an object X of class "ppp" such that marks(X) is a factor.

Summary statistics for a marked point pattern: A marked point pattern is represented by an object X of class "ppp" with a component X$marks. The entries in the vector X$marks may be numeric, complex, string or any other atomic type. For numeric marks, there are the following functions:

For marks of any type, there are the following:

Alternatively use cut.ppp to convert a marked point pattern to a multitype point pattern.

Programming tools:

Summary statistics for a three-dimensional point pattern:

These are for 3-dimensional point pattern objects (class pp3).

Related facilities:

Summary statistics for random sets:

These work for point patterns (class ppp), line segment patterns (class psp) or windows (class owin).

Model fitting

Functions for fitting point process models are now contained in the spatstat.model package.

Simulation

There are many ways to generate a random point pattern, line segment pattern, pixel image or tessellation in spatstat.

Random point patterns: Functions for random generation are now contained in the spatstat.random package.

See also varblock for estimating the variance of a summary statistic by block resampling, and lohboot for another bootstrap technique.

Fitted point process models:

If you have fitted a point process model to a point pattern dataset, the fitted model can be simulated.

Methods for simulating a fitted model are now contained in the spatstat.model package.

Other random patterns: Functions for random generation are now contained in the spatstat.random package.

Simulation-based inference

Hypothesis tests:

More recently-developed tests:

Model diagnostics:

Classical measures of model sensitivity such as leverage and influence, and classical model diagnostic tools such as residuals, partial residuals, and effect estimates, have been adapted to point process models. These capabilities are now provided in the spatstat.model package.

Resampling and randomisation procedures

You can build your own tests based on randomisation and resampling using the following capabilities:

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

Kasper Klitgaard Berthelsen, Ottmar Cronie, Tilman Davies, Julian Gilbey, Yongtao Guan, Ute Hahn, Kassel Hingee, Abdollah Jalilian, Marie-Colette van Lieshout, Greg McSwiggan, Tuomas Rajala, Suman Rakshit, Dominic Schuhmacher, Rasmus Waagepetersen and Hangsheng Wang made substantial contributions of code.

For comments, corrections, bug alerts and suggestions, we thank Monsuru Adepeju, Corey Anderson, Ang Qi Wei, Ryan Arellano, Jens Astrom, Robert Aue, Marcel Austenfeld, Sandro Azaele, Malissa Baddeley, Guy Bayegnak, Colin Beale, Melanie Bell, Thomas Bendtsen, Ricardo Bernhardt, Andrew Bevan, Brad Biggerstaff, Anders Bilgrau, Leanne Bischof, Christophe Biscio, Roger Bivand, Jose M. Blanco Moreno, Florent Bonneu, Jordan Brown, Ian Buller, Julian Burgos, Simon Byers, Ya-Mei Chang, Jianbao Chen, Igor Chernayavsky, Y.C. Chin, Bjarke Christensen, Lucia Cobo Sanchez, Jean-Francois Coeurjolly, Kim Colyvas, Hadrien Commenges, Rochelle Constantine, Robin Corria Ainslie, Richard Cotton, Marcelino de la Cruz, Peter Dalgaard, Mario D'Antuono, Sourav Das, Peter Diggle, Patrick Donnelly, Ian Dryden, Stephen Eglen, Ahmed El-Gabbas, Belarmain Fandohan, Olivier Flores, David Ford, Peter Forbes, Shane Frank, Janet Franklin, Funwi-Gabga Neba, Oscar Garcia, Agnes Gault, Jonas Geldmann, Marc Genton, Shaaban Ghalandarayeshi, Jason Goldstick, Pavel Grabarnik, C. Graf, Ute Hahn, Andrew Hardegen, Martin Bogsted Hansen, Martin Hazelton, Juha Heikkinen, Mandy Hering, Markus Herrmann, Maximilian Hesselbarth, Paul Hewson, Hamidreza Heydarian, Kurt Hornik, Philipp Hunziker, Jack Hywood, Ross Ihaka, Cenk Icos, Aruna Jammalamadaka, Robert John-Chandran, Devin Johnson, Mahdieh Khanmohammadi, Bob Klaver, Lily Kozmian-Ledward, Peter Kovesi, Mike Kuhn, Jeff Laake, Robert Lamb, Frederic Lavancier, Tom Lawrence, Tomas Lazauskas, Jonathan Lee, George Leser, Angela Li, Li Haitao, George Limitsios, Andrew Lister, Nestor Luambua, Ben Madin, Martin Maechler, Kiran Marchikanti, Jeff Marcus, Robert Mark, Peter McCullagh, Monia Mahling, Jorge Mateu Mahiques, Ulf Mehlig, Frederico Mestre, Sebastian Wastl Meyer, Mi Xiangcheng, Lore De Middeleer, Robin Milne, Enrique Miranda, Jesper Moller, Annie Mollie, Ines Moncada, Mehdi Moradi, Virginia Morera Pujol, Erika Mudrak, Gopalan Nair, Nader Najari, Nicoletta Nava, Linda Stougaard Nielsen, Felipe Nunes, Jens Randel Nyengaard, Jens Oehlschlaegel, Thierry Onkelinx, Sean O'Riordan, Evgeni Parilov, Jeff Picka, Nicolas Picard, Tim Pollington, Mike Porter, Sergiy Protsiv, Adrian Raftery, Ben Ramage, Pablo Ramon, Xavier Raynaud, Nicholas Read, Matt Reiter, Ian Renner, Tom Richardson, Brian Ripley, Ted Rosenbaum, Barry Rowlingson, Jason Rudokas, Tyler Rudolph, John Rudge, Christopher Ryan, Farzaneh Safavimanesh, Aila Sarkka, Cody Schank, Katja Schladitz, Sebastian Schutte, Bryan Scott, Olivia Semboli, Francois Semecurbe, Vadim Shcherbakov, Shen Guochun, Shi Peijian, Harold-Jeffrey Ship, Tammy L Silva, Ida-Maria Sintorn, Yong Song, Malte Spiess, Mark Stevenson, Kaspar Stucki, Jan Sulavik, Michael Sumner, P. Surovy, Ben Taylor, Thordis Linda Thorarinsdottir, Leigh Torres, Berwin Turlach, Torben Tvedebrink, Kevin Ummer, Medha Uppala, Andrew van Burgel, Tobias Verbeke, Mikko Vihtakari, Alexendre Villers, Fabrice Vinatier, Maximilian Vogtland, Sasha Voss, Sven Wagner, Hao Wang, H. Wendrock, Jan Wild, Carl G. Witthoft, Selene Wong, Maxime Woringer, Luke Yates, Mike Zamboni and Achim Zeileis.

Author(s)

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

Description

Extract a subset of the data for a spatially sampled function.

Usage

## S3 method for class 'ssf'
x[i, j, ..., drop]

Arguments

Details

This is the subset operator for the class "ssf".

Value

Another object of class "ssf".

Author(s)

Adrian Baddeley [email protected].

See Also

ssf, with.ssf

Examples

f <- ssf(cells, data.frame(d=nndist(cells), i=1:42))
  f
  f[1:10,]
  f[ ,1]

Description

Computes an adaptive estimate of the intensity function of a point pattern.

Usage

adaptive.density(X, ..., method=c("voronoi","kernel", "nearest"))

Arguments

Details

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

If method="voronoi" the data are passed to the function densityVoronoi which estimates the intensity using the Voronoi-Dirichlet tessellation.

If method="kernel" the data are passed to the function densityAdaptiveKernel.ppp which estimates the intensity using a variable-bandwidth kernel estimator. (This is not yet supported when X has class "lpp".)

If method="nearest" the data are passed to the function nndensity.ppp which estimates the intensity using the distance to the k-th nearest data point. (This is not yet supported when X has class "lpp".)

Value

A pixel image (object of class "im" or "linim") whose values are estimates of the intensity of X.

Author(s)

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

See Also

density.ppp, densityVoronoi, densityAdaptiveKernel.ppp, nndensity.ppp, im.object.

Examples

plot(adaptive.density(nztrees, 1), main="Voronoi estimate")

Description

Calculates the $F$, $G$, $J$, and $K$ summary functions for an unmarked point pattern. Returns them as a function array (of class "fasp", see fasp.object).

Usage

allstats(pp, ..., dataname=NULL, verb=FALSE)

Arguments

Details

This computes four standard summary statistics for a point pattern: the empty space function $F(r)$, nearest neighbour distance distribution function $G(r)$, van Lieshout-Baddeley function $J(r)$ and Ripley's function $K(r)$. The real work is done by Fest, Gest, Jest and Kest respectively. Consult the help files for these functions for further information about the statistical interpretation of $F$, $G$, $J$ and $K$.

If verb is TRUE, then “progress reports” (just indications of completion) are printed out when the calculations are finished for each of the four function types.

The overall title of the array of four functions (for plotting by plot.fasp) will be formed from the argument dataname. If this is not given, it defaults to the expression for pp given in the call to allstats.

Value

A list of length 4 containing the $F$, $G$, $J$ and $K$ functions respectively.

The list can be plotted directly using plot (which dispatches to plot.anylist).

Each list entry retains the format of the output of the relevant estimating routine Fest, Gest, Jest or Kest. Thus each entry in the list is a function value table (object of class "fv", see fv.object).

The default formulae for plotting these functions are cbind(km,theo) ~ r for F, G, and J, and cbind(trans,theo) ~ r for K.

Author(s)

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

See Also

plot.anylist, plot.fv, fv.object, Fest, Gest, Jest, Kest

Examples

a <- allstats(swedishpines,dataname="Swedish Pines")
        if(interactive()) {
        plot(a)
        plot(a, subset=list("r<=15","r<=15","r<=15","r<=50"))
        }

Description

Given a marked point pattern, this computes the estimates of a selected summary function ($F$,$G$, $J$, $K$ etc) of the pattern, for all possible combinations of marks, and returns these functions in an array.

Usage

alltypes(X, fun="K", ...,
           dataname=NULL,verb=FALSE,envelope=FALSE,reuse=TRUE)

Arguments

Details

This routine is a convenient way to analyse the dependence between types in a multitype point pattern. It computes the estimates of a selected summary function of the pattern, for all possible combinations of marks. It returns these functions in an array (an object of class "fasp") amenable to plotting by plot.fasp().

The argument fun specifies the summary function that will be evaluated for each type of point, or for each pair of types. It may be either an R function or a character string.

Suppose that the points have possible types $1,2,\ldots,m$ and let $X_i$ denote the pattern of points of type $i$ only.

If fun="F" then this routine calculates, for each possible type $i$, an estimate of the Empty Space Function $F_i(r)$ of $X_i$. See Fest for explanation of the empty space function. The estimate is computed by applying Fest to $X_i$ with the optional arguments ....

If fun is "Gcross", "Jcross", "Kcross" or "Lcross", the routine calculates, for each pair of types $(i,j)$, an estimate of the “i-toj” cross-type function $G_{ij}(r)$, $J_{ij}(r)$, $K_{ij}(r)$ or $L_{ij}(r)$ respectively describing the dependence between $X_i$ and $X_j$. See Gcross, Jcross, Kcross or Lcross respectively for explanation of these functions. The estimate is computed by applying the relevant function (Gcross etc) to X using each possible value of the arguments i,j, together with the optional arguments ....

If fun is "pcf" the routine calculates the cross-type pair correlation function pcfcross between each pair of types.

If fun is "Gdot", "Jdot", "Kdot" or "Ldot", the routine calculates, for each type $i$, an estimate of the “i-to-any” dot-type function $G_{i\bullet}(r)$, $J_{i\bullet}(r)$ or $K_{i\bullet}(r)$ or $L_{i\bullet}(r)$ respectively describing the dependence between $X_i$ and $X$. See Gdot, Jdot, Kdot or Ldot respectively for explanation of these functions. The estimate is computed by applying the relevant function (Gdot etc) to X using each possible value of the argument i, together with the optional arguments ....

The letters "G", "J", "K" and "L" are interpreted as abbreviations for Gcross, Jcross, Kcross and Lcross respectively, assuming the point pattern is marked. If the point pattern is unmarked, the appropriate function Fest, Jest, Kest or Lest is invoked instead.

If envelope=TRUE, then as well as computing the value of the summary function for each combination of types, the algorithm also computes simulation envelopes of the summary function for each combination of types. The arguments ... are passed to the function envelope to control the number of simulations, the random process generating the simulations, the construction of envelopes, and so on.

When envelope=TRUE it is possible that errors could occur because the simulated point patterns do not satisfy the requirements of the summary function (for example, because the simulated pattern is empty and fun requires at least one point). If the number of such errors exceeds the maximum permitted number maxnerr, then the envelope algorithm will give up, and will return the empirical summary function for the data point pattern, fun(X), in place of the envelope.

Value

A function array (an object of class "fasp", see fasp.object). This can be plotted using plot.fasp.

If the pattern is not marked, the resulting “array” has dimensions $1 \times 1$. Otherwise the following is true:

If fun="F", the function array has dimensions $m \times 1$ where $m$ is the number of different marks in the point pattern. The entry at position [i,1] in this array is the result of applying Fest to the points of type i only.

If fun is "Gdot", "Jdot", "Kdot" or "Ldot", the function array again has dimensions $m \times 1$. The entry at position [i,1] in this array is the result of Gdot(X, i), Jdot(X, i) Kdot(X, i) or Ldot(X, i) respectively.

If fun is "Gcross", "Jcross", "Kcross" or "Lcross" (or their abbreviations "G", "J", "K" or "L"), the function array has dimensions $m \times m$. The [i,j] entry of the function array (for $i \neq j$) is the result of applying the function Gcross, Jcross, Kcross orLcross to the pair of types (i,j). The diagonal [i,i] entry of the function array is the result of applying the univariate function Gest, Jest, Kest or Lest to the points of type i only.

If envelope=FALSE, then each function entry fns[[i]] retains the format of the output of the relevant estimating routine Fest, Gest, Jest, Kest, Lest, Gcross, Jcross ,Kcross, Lcross, Gdot, Jdot, Kdot or Ldot The default formulae for plotting these functions are cbind(km,theo) ~ r for F, G, and J functions, and cbind(trans,theo) ~ r for K and L functions.

If envelope=TRUE, then each function entry fns[[i]] has the same format as the output of the envelope command.

Note

Sizeable amounts of memory may be needed during the calculation.

Author(s)

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

See Also

plot.fasp, fasp.object, Fest, Gest, Jest, Kest, Lest, Gcross, Jcross, Kcross, Lcross, Gdot, Jdot, Kdot, envelope.

Examples

# bramblecanes (3 marks).
   bram <- bramblecanes
   
   bF <- alltypes(bram,"F",verb=TRUE)
   plot(bF)
   if(interactive()) {
     plot(alltypes(bram,"G"))
     plot(alltypes(bram,"Gdot"))
   }
   
   # Swedishpines (unmarked).
  swed <- swedishpines
   
   plot(alltypes(swed,"K"))

   plot(alltypes(amacrine, "pcf"), ylim=c(0,1.3))

   # envelopes
   bKE <- alltypes(bram,"K",envelope=TRUE,nsim=19)
   # global version:
   
     bFE <- alltypes(bram,"F",envelope=TRUE,nsim=19,global=TRUE)
   

   # extract one entry
   as.fv(bKE[1,1])

Description

Converts an envelope object to a data frame.

Usage

## S3 method for class 'envelope'
as.data.frame(x, ..., simfuns=FALSE)

Arguments

Details

This is a method for the generic function as.data.frame for the class of envelopes (see envelope.

The result is a data frame with columns containing the values of the function argument (usually named r), the function estimate for the original point pattern data (obs), the upper and lower envelope limits (hi and lo), and possibly additional columns.

If simfuns=TRUE, the result also includes columns of values of the simulated functions that were used to compute the envelope. This is possible only when the envelope was computed with the argument savefuns=TRUE in the call to envelope.

Value

A data frame.

Author(s)

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

Examples

E <- envelope(cells, nsim=5, savefuns=TRUE)
  tail(as.data.frame(E))
  tail(as.data.frame(E, simfuns=TRUE))

Description

Converts an object of class "fv" to an R language function.

Usage

## S3 method for class 'fv'
as.function(x, ..., value=".y", extrapolate=FALSE)

Arguments

Details

A function value table (object of class "fv") is a convenient way of storing and plotting several different estimates of the same function. Objects of this class are returned by many commands in spatstat, such as Kest, which returns an estimate of Ripley's $K$-function for a point pattern dataset.

Sometimes it is useful to convert the function value table to a function in the R language. This is done by as.function.fv. It converts an object x of class "fv" to an R function f.

If f <- as.function(x) then f is an R function that accepts a numeric argument and returns a corresponding value for the summary function by linear interpolation between the values in the table x.

Argument values lying outside the range of the table yield an NA value (if extrapolate=FALSE) or the function value at the nearest endpoint of the range (if extrapolate = TRUE). To apply different rules to the left and right extremes, use extrapolate=c(TRUE,FALSE) and so on.

Typically the table x contains several columns of function values corresponding to different edge corrections. Auxiliary information for the table identifies one of these columns as the recommended value. By default, the values of the function f <- as.function(x) are taken from this column of recommended values. This default can be changed using the argument value, which can be a character string or character vector of names of columns of x. Alternatively value can be one of the abbreviations used by fvnames.

If value specifies a single column of the table, then the result is a function f(r) with a single numeric argument r (with the same name as the orginal argument of the function table).

If value specifies several columns of the table, then the result is a function f(r,what) where r is the numeric argument and what is a character string identifying the column of values to be used.

The formal arguments of the resulting function are f(r, what=value), which means that in a call to this function f, the permissible values of what are the entries of the original vector value; the default value of what is the first entry of value.

The command as.function.fv is a method for the generic command as.function.

Value

A function(r) or function(r,what) where r is the name of the original argument of the function table.

Author(s)

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

See Also

as.function.rhohat, fv, fv.object, fvnames, plot.fv, Kest

Examples

K <- Kest(cells)
  f <- as.function(K)
  f
  f(0.1)
  g <- as.function(K, value=c("iso", "trans"))
  g
  g(0.1, "trans")

Description

Converts an object of class "rhohat" to an R language function.

Usage

## S3 method for class 'rhohat'
as.function(x, ..., value=".y", extrapolate=TRUE)

Arguments

Details

An object of class "rhohat" is essentially a data frame of estimated values of the function $rho(x)$ as described in the help file for rhohat.

Sometimes it is useful to convert the function value table to a function in the R language. This is done by as.function.rhohat. It converts an object x of class "rhohat" to an R function f.

The command as.function.rhohat is a method for the generic command as.function for the class "rhohat".

If f <- as.function(x) then f is an R function that accepts a numeric argument and returns a corresponding value for the summary function by linear interpolation between the values in the table x.

Argument values lying outside the range of the table yield an NA value (if extrapolate=FALSE) or the function value at the nearest endpoint of the range (if extrapolate = TRUE). To apply different rules to the left and right extremes, use extrapolate=c(TRUE,FALSE) and so on.

Typically the table x contains several columns of function values corresponding to different edge corrections. Auxiliary information for the table identifies one of these columns as the recommended value. By default, the values of the function f <- as.function(x) are taken from this column of recommended values. This default can be changed using the argument value, which can be a character string or character vector of names of columns of x. Alternatively value can be one of the abbreviations used by fvnames.

If value specifies a single column of the table, then the result is a function f(r) with a single numeric argument r (with the same name as the orginal argument of the function table).

If value specifies several columns of the table, then the result is a function f(r,what) where r is the numeric argument and what is a character string identifying the column of values to be used.

The formal arguments of the resulting function are f(r, what=value), which means that in a call to this function f, the permissible values of what are the entries of the original vector value; the default value of what is the first entry of value.

Value

A function(r) or function(r,what) where r is the name of the original argument of the function table.

Author(s)

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

See Also

rhohat, methods.rhohat, as.function.fv.

Examples

g <- rhohat(cells, "x")
  f <- as.function(g)
  f
  f(0.1)

Description

Converts data into a function table (an object of class "fv").

Usage

as.fv(x)

  ## S3 method for class 'fv'
as.fv(x)

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

  ## S3 method for class 'matrix'
as.fv(x)

  ## S3 method for class 'fasp'
as.fv(x)

  ## S3 method for class 'bw.optim'
as.fv(x)

Arguments

Details

This command converts data x, that could be interpreted as the values of a function, into a function value table (object of the class "fv" as described in fv.object). This object can then be plotted easily using plot.fv.

The dataset x may be any of the following:

• an object of class "fv";

• a matrix or data frame with at least two columns;

• an object of class "fasp", representing an array of "fv" objects.

• an object of class "minconfit", giving the results of a minimum contrast fit by the command mincontrast. The

• an object of class "kppm", representing a fitted Cox or cluster point process model, obtained from the model-fitting command kppm;

• an object of class "dppm", representing a fitted determinantal point process model, obtained from the model-fitting command dppm;

• an object of class "bw.optim", representing an optimal choice of smoothing bandwidth by a cross-validation method, obtained from commands like bw.diggle.

The function as.fv is generic, with methods for each of the classes listed above. The behaviour is as follows:

• If x is an object of class "fv", it is returned unchanged.

• If x is a matrix or data frame, the first column is interpreted as the function argument, and subsequent columns are interpreted as values of the function computed by different methods.

• If x is an object of class "fasp" representing an array of "fv" objects, these are combined into a single "fv" object.

• If x is an object of class "minconfit", or an object of class "kppm" or "dppm", the result is a function table containing the observed summary function and the best fit summary function.

• If x is an object of class "bw.optim", the result is a function table of the optimisation criterion as a function of the smoothing bandwidth.

Value

An object of class "fv" (see fv.object).

Author(s)

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

Examples

r <- seq(0, 1, length=101)
  x <- data.frame(r=r, y=r^2)
  as.fv(x)

Description

Converts data specifying an observation window in any of several formats, into an object of class "owin".

Usage

## S3 method for class 'quadrattest'
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 generic 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.

The argument W may be

• an object of class "owin"

• a structure with entries xrange, yrange specifying the $x$ and $y$ dimensions of a rectangle

• a structure with entries named xmin, xmax, ymin, ymax (in any order) specifying the $x$ and $y$ dimensions of a rectangle. This will accept objects of class bbox in the sf package.

• a numeric vector of length 4 (interpreted as (xmin, xmax, ymin, ymax) in that order) specifying the $x$ and $y$ dimensions of a rectangle

• a structure with entries named xl, xu, yl, yu (in any order) specifying the $x$ and $y$ dimensions of a rectangle as (xmin, xmax) = (xl, xu) and (ymin, ymax) = (yl, yu). This will accept objects of class spp used in the Venables and Ripley spatial package.

• an object of class "ppp" representing a point pattern. In this case, the object's window structure will be extracted.

• an object of class "psp" representing a line segment pattern. In this case, the object's window structure will be extracted.

• an object of class "tess" representing a tessellation. In this case, the object's window structure will be extracted.

• an object of class "quad" representing a quadrature scheme. In this case, the window of the data component will be extracted.

• an object of class "im" representing a pixel image. In this case, a window of type "mask" will be returned, with the same pixel raster coordinates as the image. An image pixel value of NA, signifying that the pixel lies outside the window, is transformed into the logical value FALSE, which is the corresponding convention for window masks.

• an object of class "ppm", "kppm", "slrm" or "dppm" representing a fitted point process model. In this case, if from="data" (the default), as.owin extracts the original point pattern data to which the model was fitted, and returns the observation window of this point pattern. If from="covariates" then as.owin extracts the covariate images to which the model was fitted, and returns a binary mask window that specifies the pixel locations.

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

• A data.frame with exactly three columns. Each row of the data frame corresponds to one pixel. Each row contains the $x$ and $y$ coordinates of a pixel, and a logical value indicating whether the pixel lies inside the window.

• A data.frame with exactly two columns. Each row of the data frame contains the $x$ and $y$ coordinates of a pixel that lies inside the window.

• an object of class "distfun", "nnfun" or "funxy" representing a function of spatial location, defined on a spatial domain. The spatial domain of the function will be extracted.

• an object of class "rmhmodel" representing a point process model that can be simulated using rmh. The window (spatial domain) of the model will be extracted. The window may be NULL in some circumstances (indicating that the simulation window has not yet been determined). This is not treated as an error, because the argument fatal defaults to FALSE for this method.

• an object of class "layered" representing a list of spatial objects. See layered. In this case, as.owin will be applied to each of the objects in the list, and the union of these windows will be returned.

• An object of another suitable class from another package. For full details, see vignette('shapefiles').

If the argument W is not in one of these formats and 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).

When W is a data frame, the argument step can be used to specify the pixel grid spacing; otherwise, the spacing will be guessed from the data.

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, as.owin.rmhmodel, as.owin.lpp.

owin.object, owin.

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

Examples

te <- quadrat.test(redwood, nx=3)
  as.owin(te)

Description

Converts data specifying a tessellation, in any of several formats, into an object of class "tess".

Usage

## S3 method for class 'quadrattest'
as.tess(X)

Arguments

Details

A tessellation is a collection of disjoint spatial regions (called tiles) that fit together to form a larger spatial region. This command creates an object of class "tess" that represents a tessellation.

This function converts data in any of several formats into an object of class "tess" for use by the spatstat package. The argument X may be

• an object of class "tess". The object will be stripped of any extraneous attributes and returned.

• a pixel image (object of class "im") with pixel values that are logical or factor values. Each level of the factor will determine a tile of the tessellation.

• a window (object of class "owin"). The result will be a tessellation consisting of a single tile.

• a set of quadrat counts (object of class "quadratcount") returned by the command quadratcount. The quadrats used to generate the counts will be extracted and returned as a tessellation.

• a quadrat test (object of class "quadrattest") returned by the command quadrat.test. The quadrats used to perform the test will be extracted and returned as a tessellation.

• a list of windows (objects of class "owin") giving the tiles of the tessellation.

The function as.tess is generic, with methods for various classes, as listed above.

Value

An object of class "tess" specifying a tessellation.

Author(s)

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

See Also

tess

Examples

h <- quadrat.test(nztrees, nx=4, ny=3)
 as.tess(h)

Description

Compute the AUC (area under the Receiver Operating Characteristic curve) for an observed point pattern.

Usage

auc(X, ...)

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

Arguments

Details

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

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.

Value

Numeric. For auc.ppp and auc.lpp, the result is a single number giving the AUC value.

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

roc

Examples

auc(swedishpines, "x")

Description

Tests the goodness-of-fit of a Poisson point process model using methods of Berman (1986).

Usage

berman.test(...)

## S3 method for class 'ppp'
berman.test(X, 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").

• 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

Examples

# Berman's data
   X <- copper$SouthPoints
   L <- copper$SouthLines
   D <- distmap(L, eps=1)
   # test of CSR
   berman.test(X, D)
   berman.test(X, D, "Z2")

Description

Advanced Use Only. Combine objects of class "fv", or glue extra columns of data onto an existing "fv" object.

Usage

## S3 method for class 'fv'
cbind(...)

bind.fv(x, y, labl = NULL, desc = NULL, preferred = NULL, clip=FALSE)

Arguments

Details

This documentation is provided for experienced programmers who want to modify the internal behaviour of spatstat.

The function cbind.fv is a method for the generic R function cbind. It combines any number of objects of class "fv" into a single object of class "fv". The objects must be compatible, in the sense that they have identical values of the function argument.

The function bind.fv is a lower level utility which glues additional columns onto an existing object x of class "fv". It has three modes of use:

• If the additional dataset y is an object of class "fv", then x and y must be compatible as described above. Then the columns of y that contain function values will be appended to the object x.

• Alternatively if y is a data frame, then y must have the same number of rows as x. All columns of y will be appended to x.

• Alternatively if y is a function in the R language, then this function will be evaluated at the argument values stored in the object x, and these function values will be appended as a new column to x.

The arguments labl and desc provide plot labels and description strings (as described in fv) for the new columns. If y is an object of class "fv" then labl and desc are optional, and default to the relevant entries in the object y. If y is a data frame then labl and desc should be provided, but there is a default.

For additional flexibility, cbind.fv also accepts arguments which are data frames or functions.

Value

An object of class "fv".

Author(s)

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

See Also

fv for creating objects of class "fv" from raw data.

collapse.fv for combining several "fv" objects with similar columns.

with.fv for evaluating expressions.

fvnames for extracting and assigning the column names of standard components of "fv" objects.

Undocumented functions for modifying an "fv" object include tweak.fv.entry and rebadge.fv.

Examples

K1 <- Kest(cells, correction="border")
   K2 <- Kest(cells, correction="iso")

   # remove column 'theo' to avoid duplication
   K2 <- K2[, names(K2) != "theo"]

   cbind(K1, K2)

   bind.fv(K1, K2, preferred="iso")

   # constrain border estimate to be monotonically increasing
   bm <- cumsum(c(0, pmax(0, diff(K1$border))))
   bind.fv(K1, data.frame(bmono=bm),
               "%s[bmo](r)",
               "monotone border-corrected estimate of %s",
               "bmono")

   # add a column of values defined by a function
   cbind(K1, upper=function(r) { pi * r^2 + 0.1 })

Description

Computes the global envelopes corresponding to the balanced independent two-stage Monte Carlo test of goodness-of-fit.

Usage

bits.envelope(X, ...,
            nsim = 19, nrank = 1,
            alternative=c("two.sided", "less", "greater"),
            leaveout=1, interpolate = FALSE,
            savefuns=FALSE, savepatterns=FALSE,
            verbose = TRUE)

Arguments

Details

Computes global simulation envelopes corresponding to the balanced independent two-stage Monte Carlo test of goodness-of-fit described by Baddeley et al (2017). The envelopes are described in Baddeley et al (2019).

If X is a point pattern, the null hypothesis is CSR.

If X is a fitted model, the null hypothesis is that model.

This command is similar to dg.envelope which corresponds to the Dao-Genton test of goodness-of-fit. It was shown in Baddeley et al (2017) that the Dao-Genton test is biased when the significance level is very small (small $p$-values are not reliable) and we recommend bits.envelope in this case.

Value

An object of class "fv".

Author(s)

Adrian Baddeley, Andrew Hardegen, Tom Lawrence, Robin Milne, Gopalan Nair and Suman Rakshit. Implemented by Adrian Baddeley [email protected].

References

Dao, N.A. and Genton, M. (2014) A Monte Carlo adjusted goodness-of-fit test for parametric models describing spatial point patterns. Journal of Graphical and Computational Statistics 23, 497–517.

Baddeley, A., Hardegen, A., Lawrence, T., Milne, R.K., Nair, G. and Rakshit, S. (2017) On two-stage Monte Carlo tests of composite hypotheses. Computational Statistics and Data Analysis 114, 75–87.

Baddeley, A., Hardegen, A., Lawrence, L., Milne, R.K., Nair, G.M. and Rakshit, S. (2019) Pushing the envelope: extensions of graphical Monte Carlo tests. In preparation.

See Also

dg.envelope, bits.test, mad.test, envelope

Examples

ns <- if(interactive()) 19 else 4
  E <- bits.envelope(swedishpines, Lest, nsim=ns)
  E
  plot(E)
  Eo <- bits.envelope(swedishpines, Lest, alternative="less", nsim=ns)
  Ei <- bits.envelope(swedishpines, Lest, interpolate=TRUE, nsim=ns)

Description

Performs a Balanced Independent Two-Stage Monte Carlo test of goodness-of-fit for spatial pattern.

Usage

bits.test(X, ...,
        exponent = 2, nsim=19, 
        alternative=c("two.sided", "less", "greater"),
        leaveout=1, interpolate = FALSE,
        savefuns=FALSE, savepatterns=FALSE,
        verbose = TRUE)

Arguments

Details

Performs the Balanced Independent Two-Stage Monte Carlo test proposed by Baddeley et al (2017), an improvement of the Dao-Genton (2014) test.

If X is a point pattern, the null hypothesis is CSR.

If X is a fitted model, the null hypothesis is that model.

The argument use.theory passed to envelope determines whether to compare the summary function for the data to its theoretical value for CSR (use.theory=TRUE) or to the sample mean of simulations from CSR (use.theory=FALSE).

The argument leaveout specifies how to calculate the discrepancy between the summary function for the data and the nominal reference value, when the reference value must be estimated by simulation. The values leaveout=0 and leaveout=1 are both algebraically equivalent (Baddeley et al, 2014, Appendix) to computing the difference observed - reference where the reference is the mean of simulated values. The value leaveout=2 gives the leave-two-out discrepancy proposed by Dao and Genton (2014).

Value

A hypothesis test (object of class "htest" which can be printed to show the outcome of the test.

Author(s)

Adrian Baddeley, Andrew Hardegen, Tom Lawrence, Robin Milne, Gopalan Nair and Suman Rakshit. Implemented by Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

References

Dao, N.A. and Genton, M. (2014) A Monte Carlo adjusted goodness-of-fit test for parametric models describing spatial point patterns. Journal of Graphical and Computational Statistics 23, 497–517.

Baddeley, A., Diggle, P.J., Hardegen, A., Lawrence, T., Milne, R.K. and Nair, G. (2014) On tests of spatial pattern based on simulation envelopes. Ecological Monographs 84 (3) 477–489.

Baddeley, A., Hardegen, A., Lawrence, L., Milne, R.K., Nair, G.M. and Rakshit, S. (2017) On two-stage Monte Carlo tests of composite hypotheses. Computational Statistics and Data Analysis 114, 75–87.

See Also

Simulation envelopes: bits.envelope.

Other tests: dg.test, dclf.test, mad.test.

Examples

ns <- if(interactive()) 19 else 4
 bits.test(cells, nsim=ns)
 bits.test(cells, alternative="less", nsim=ns)
 bits.test(cells, nsim=ns, interpolate=TRUE)

Description

Applies a Gaussian blur to a pixel image.

Usage

blur(x, sigma = NULL, ...,
     kernel="gaussian", normalise=FALSE, bleed = TRUE, varcov=NULL)

## S3 method for class 'im'
Smooth(X, sigma = NULL, ...,
                    kernel="gaussian",
                    normalise=FALSE, bleed = TRUE, varcov=NULL)

Arguments

Details

This command applies a Gaussian blur to the pixel image x.

Smooth.im is a method for the generic Smooth for pixel images. It is currently identical to blur, apart from the name of the first argument.

The blurring kernel is the isotropic Gaussian kernel with standard deviation sigma, or the anisotropic Gaussian kernel with variance-covariance matrix varcov. The arguments sigma and varcov are incompatible. Also sigma may be a vector of length 2 giving the standard deviations of two independent Gaussian coordinates, thus equivalent to varcov = diag(sigma^2).

If the pixel values of x include some NA values (meaning that the image domain does not completely fill the rectangular frame) then these NA values are first reset to zero.

The algorithm then computes the convolution $x \ast G$ of the (zero-padded) pixel image $x$ with the specified Gaussian kernel $G$.

If normalise=FALSE, then this convolution $x\ast G$ is returned. If normalise=TRUE, then the convolution $x \ast G$ is normalised by dividing it by the convolution $w \ast G$ of the image domain w with the same Gaussian kernel. Normalisation ensures that the result can be interpreted as a weighted average of input pixel values, without edge effects due to the shape of the domain.

If bleed=FALSE, then pixel values outside the original image domain are set to NA. Thus the output is a pixel image with the same domain as the input. If bleed=TRUE, then no such alteration is performed, and the result is a pixel image defined everywhere in the rectangular frame containing the input image.

Computation is performed using the Fast Fourier Transform.

Value

A pixel image with the same pixel array as the input image x.

Author(s)

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

See Also

interp.im for interpolating a pixel image to a finer resolution, density.ppp for blurring a point pattern, Smooth.ppp for interpolating marks attached to points.

Examples

Z <- as.im(function(x,y) { 4 * x^2 + 3 * y }, letterR)
   opa <- par(mfrow=c(1,3))
   plot(Z)
   plot(letterR, add=TRUE)
   plot(blur(Z, 0.3, bleed=TRUE))
   plot(letterR, add=TRUE)
   plot(blur(Z, 0.3, bleed=FALSE))
   plot(letterR, add=TRUE)
   par(opa)

Description

Calculate the discrete or continuous Boyce index for a spatial point pattern dataset.

Usage

boyce(X, Z, ..., breaks = NULL, halfwidth = NULL)

Arguments

Details

Given a spatial point pattern X and some kind of explanatory information Z, this function computes either the index originally defined by Boyce et al (2002) or the ‘continuous Boyce index’ defined by Hirzel et al (2006).

Boyce et al (2002) defined an index of habitat suitability in which the study region $W$ is first divided into separate subregions $C_1,\ldots,C_m$ based on appropriate scientific considerations. Then we count the number $n_j$ of data points of X that fall in each subregion $C_j$, measure the area $a_j$ of each subregion $C_j$, and calculate the index

$B_j = \frac{n_j/n}{a_j/a}$

where $a$ is the total area and $n$ is the total number of points in X.

Hirzel et al (2006) defined another version of this index which is based on a continuous spatial covariate. For each possible value $z$ of the covariate $Z$, consider the region $C(z)$ where the value of the covariate lies between $z-h$ and $z+h$, where $h$ is the chosen ‘halfwidth’. The ‘continuous Boyce index’ is

$B(z) = \frac{n(z)/n}{a(z)/a}$

where $n(z)$ is the number of points of X falling in $C(z)$, and $a(z)$ is the area of $C(z)$.

If Z is a tessellation (object of class "tess"), the algorithm calculates the original (‘discrete’) Boyce index (Boyce et al, 2002) for each tile of the tessellation. The result is another tessellation, identical to Z except that the mark values are the values of the discrete Boyce index.

If Z is a pixel image whose values are categorical (i.e. factor values), then Z is treated as a tessellation, with one tile for each level of the factor. The discrete Boyce index is then calculated. The result is a tessellation with marks that are the values of the discrete Boyce index.

Otherwise, if Z is a spatial covariate such as a pixel image, a function(x,y) or one of the characters "x" or "y", then exactly one of the arguments breaks or halfwidth must be given.

• if halfwidth is given, it should be a single positive number. The continuous Boyce index (Hirzel et al, 2006) is computed using the specified halfwidth $h$. The result is an object of class "fv" that can be plotted to show $B(z)$ as a function of $z$.

• if breaks is given, it can be either a numeric vector of possible values of Z defining the breakpoints for the bands of values of Z, or a single integer specifying the number of evenly-spaced breakpoints that should be created. The discrete Boyce index (Boyce et al, 2002) is computed. The result is an object of class "fv" that can be plotted to show the discrete Boyce index as a function of $z$.

When Z is a spatial covariate (not factor-valued), the calculation is performed using rhohat.ppp (since the Boyce index is a special case of rhohat). Arguments ... passed to rhohat.ppp control the accuracy of the spatial discretisation and other parameters of the algorithm.

Value

A tessellation (object of class "tess") or a function value table (object of class "fv") as explained above.

Author(s)

Adrian Baddeley [email protected]

References

Boyce, M.S., Vernier, P.R., Nielsen, S.E. and Schmiegelow, F.K.A. (2002) Evaluating resource selection functions. Ecological modelling 157, 281–300.

Hirzel, A.H., Le Lay, V., Helfer, V., Randin, C. and Guisan, A. (2006) Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling 199, 142–152.

See Also

rhohat

Examples

online <- interactive()
  ## a simple tessellation
  V <- quadrats(Window(bei), 4, 3)
  if(online) plot(V)

  ## discrete Boyce index for a simple tessellation
  A <- boyce(bei, V)

  if(online) {
   plot(A, do.col=TRUE)
   marks(A)
   tilenames(A)
  }

  ## spatial covariate: terrain elevation
  Z <- bei.extra$elev

  ## continuous Boyce index for terrain elevation
  BC <- boyce(bei, Z, halfwidth=10)

  if(online) plot(BC)

  ## discrete Boyce index for terrain elevation steps of height 5 metres
  bk <- c(seq(min(Z), max(Z), by=5), Inf)
  BD <- boyce(bei, Z, breaks=bk)

  if(online) plot(BD)

Description

Computes adaptive smoothing bandwidths for a spatial point pattern, according to the inverse-square-root rule of Abramson (1982).

Usage

## S3 method for class 'ppp'
bw.abram(X, h0, 
         ...,
         at=c("points", "pixels"),
         hp = h0, pilot = NULL, trim=5, smoother=density.ppp)

Arguments

Details

This function computes adaptive smoothing bandwidths using the methods of Abramson (1982) and Hall and Marron (1988).

The function bw.abram is generic. The function bw.abram.ppp documented here is the method for spatial point patterns (objects of class "ppp").

If at="points" (the default) a smoothing bandwidth is computed for each point in the pattern X. Alternatively if at="pixels" a smoothing bandwidth is computed for each spatial location in a pixel grid.

Under the Abramson-Hall-Marron rule, the bandwidth at location $u$ is

$h(u) = \mbox{\texttt{h0}} * \mbox{min}[ \frac{\tilde{f}(u)^{-1/2}}{\gamma}, \mbox{\texttt{trim}} ]$

where $\tilde{f}(u)$ is a pilot estimate of the spatially varying probability density. The variable bandwidths are rescaled by $\gamma$, the geometric mean of the $\tilde{f}(u)^{-1/2}$ terms evaluated at the data; this allows the global bandwidth h0 to be considered on the same scale as a corresponding fixed bandwidth. The trimming value trim has the same interpretation as the required ‘clipping’ of the pilot density at some small nominal value (see Hall and Marron, 1988), to necessarily prevent extreme bandwidths (which can occur at very isolated observations).

The pilot density or intensity is determined as follows:

• If pilot is a pixel image, this is taken as the pilot density or intensity.

• If pilot is NULL, then the pilot intensity is computed as a fixed-bandwidth kernel intensity estimate using density.ppp applied to the data pattern X using the pilot bandwidth hp.

• If pilot is a different point pattern on the same spatial domain as X, then the pilot intensity is computed as a fixed-bandwidth kernel intensity estimate using density.ppp applied to pilot using the pilot bandwidth hp.

In each case the pilot density or intensity is renormalised to become a probability density, and then the Abramson rule is applied.

Instead of calculating the pilot as a fixed-bandwidth density estimate, the user can specify another density estimation procedure using the argument smoother. This should be either a function or the character string name of a function. It will replace density.ppp as the function used to calculate the pilot estimate. The pilot estimate will be computed as smoother(X, sigma=hp, ...) if pilot is NULL, or smoother(pilot, sigma=hp, ...) if pilot is a point pattern. If smoother does not recognise the argument name sigma for the smoothing bandwidth, then hp is effectively ignored, as shown in the Examples.

Value

Either a numeric vector of length npoints(X) giving the Abramson bandwidth for each point (when at = "points", the default), or the entire pixel image of the Abramson bandwidths over the relevant spatial domain (when at = "pixels").

Author(s)

Tilman Davies [email protected]. Adapted by Adrian Baddeley [email protected].

References

Abramson, I. (1982) On bandwidth variation in kernel estimates — a square root law. Annals of Statistics, 10(4), 1217-1223.

Davies, T.M. and Baddeley, A. (2018) Fast computation of spatially adaptive kernel estimates. Statistics and Computing, 28(4), 937-956.

Davies, T.M., Marshall, J.C., and Hazelton, M.L. (2018) Tutorial on kernel estimation of continuous spatial and spatiotemporal relative risk. Statistics in Medicine, 37(7), 1191-1221.

Hall, P. and Marron, J.S. (1988) Variable window width kernel density estimates of probability densities. Probability Theory and Related Fields, 80, 37-49.

Silverman, B.W. (1986) Density Estimation for Statistics and Data Analysis. Chapman and Hall, New York.

See Also

bw.abram

Examples

# 'ch' just 58 laryngeal cancer cases
ch <- split(chorley)[[1]]

h <- bw.abram(ch,h0=1,hp=0.7)
length(h)
summary(h)
if(interactive()) hist(h)

# calculate pilot based on all 1036 observations
h.pool <- bw.abram(ch,h0=1,hp=0.7,pilot=chorley)
length(h.pool)
summary(h.pool)
if(interactive()) hist(h.pool)

# get full image used for 'h' above
him <- bw.abram(ch,h0=1,hp=0.7,at="pixels")
plot(him);points(ch,col="grey")

# use Voronoi-Dirichlet pilot ('hp' is ignored)
hvo <- bw.abram(ch, h0=1, smoother=densityVoronoi)

Description

Uses Cronie and van Lieshout's criterion based on Cambell's formula to select a smoothing bandwidth for the kernel estimation of point process intensity.

Usage

bw.CvL(X, ..., srange = NULL, ns = 16, sigma = NULL, warn=TRUE)

Arguments

Details

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

The bandwidth $\sigma$ is chosen to minimise the discrepancy between the area of the observation window and the sum of reciprocal estimated intensity values at the points of the point process

$\mbox{CvL}(\sigma) = (|W| - \sum_i 1/\hat\lambda(x_i))^2$

where the sum is taken over all the data points $x_i$, and where $\hat\lambda(x_i)$ is the kernel-smoothing estimate of the intensity at $x_i$ with smoothing bandwidth $\sigma$.

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

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)

Ottmar Cronie [email protected] and Marie-Colette van Lieshout [email protected]. Adapted for spatstat by Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

References

Cronie, O and Van Lieshout, M N M (2018) A non-model-based approach to bandwidth selection for kernel estimators of spatial intensity functions, Biometrika, 105, 455-462.

See Also

density.ppp, bw.optim.object.

Alternative methods: bw.diggle, bw.scott, bw.ppl, bw.frac.

For adaptive smoothing bandwidths, use bw.CvL.adaptive.

Examples

if(interactive()) {
    b <- bw.CvL(redwood)
    b
    plot(b, main="Cronie and van Lieshout bandwidth criterion for redwoods")
    plot(density(redwood, b))
    plot(density(redwood, bw.CvL))
  }

Description

Uses the Cronie-Van Lieshout criterion to select the global smoothing bandwidth for adaptive kernel estimation of point process intensity.

Usage

bw.CvL.adaptive(X, ..., 
            hrange = NULL, nh = 16, h=NULL,
            bwPilot = bw.scott.iso(X),
            edge = FALSE, diggle = TRUE)

Arguments

Details

This function selects an appropriate value of global bandwidth h0 for adaptive kernel estimation of the intensity function for the point pattern X.

In adaptive estimation, each point in the point pattern is subjected to a different amount of smoothing, controlled by data-dependent or spatially-varying bandwidths. The global bandwidth h0 is a scale factor which is used to adjust all of the data-dependent bandwidths according to the Abramson (1982) square-root rule.

This function considers each candidate value of bandwidth $h$, performs the smoothing steps described above, extracts the adaptively-estimated intensity values $\hat\lambda(x_i)$ at each data point $x_i$, and calculates the Cronie-Van Lieshout criterion

$\mbox{CvL}(h) = \sum_{i=1}^n \frac 1 {\hat\lambda(x_i)}.$

The value of $h$ which minimises the squared difference

$LP2(h) = (CvL(h) - |W|)^2$

(where |W| is the area of the window of X) is selected as the optimal global bandwidth.

Bandwidths h are physical distance values expressed in the same units as the coordinates of X.

Value

A single numerical value giving the selected global 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)

Marie-Colette Van Lieshout. Modified by Adrian Baddeley [email protected].

References

Abramson, I. (1982) On bandwidth variation in kernel estimates — a square root law. Annals of Statistics, 10(4), 1217-1223.

Cronie, O and Van Lieshout, M N M (2018) A non-model-based approach to bandwidth selection for kernel estimators of spatial intensity functions, Biometrika, 105, 455-462.

Van Lieshout, M.N.M. (2021) Infill asymptotics for adaptive kernel estimators of spatial intensity. Australian and New Zealand Journal of Statistics 63 (1) 159–181.

See Also

bw.optim.object.

adaptive.density, densityAdaptiveKernel.ppp, bw.abram.ppp, density.ppp.

To select a fixed smoothing bandwidth using the Cronie-Van Lieshout criterion, use bw.CvL.

Examples

online <- interactive()
  if(online) {
    h0 <- bw.CvL.adaptive(redwood3)
  } else {
    ## faster computation for package checker
    h0 <- bw.CvL.adaptive(redwood3, nh=8,
                          hrange=c(1/4, 4) * bw.diggle(redwood3))
  }
  plot(h0)
  plot(as.fv(h0), CvL ~ h)
  if(online) {
    Z <- densityAdaptiveKernel(redwood3, h0)
    plot(Z)
  }

Description

Selects an optimal bandwidth for diffusion smoothing using the Cronie-van Lieshout rule.

Usage

bw.CvLHeat(X, ..., srange=NULL, ns=16, sigma=NULL,
         leaveoneout=TRUE, verbose = TRUE)

Arguments

Details

This algorithm selects the optimal global bandwidth for kernel estimation of intensity for the dataset X using diffusion smoothing densityHeat.ppp.

If sigma is a numeric value, the algorithm finds the optimal bandwidth tau <= sigma.

If sigma is a pixel image or function, the algorithm finds the optimal fraction 0 < f <= 1 such that smoothing with f * sigma would be optimal.

Value

A numerical value giving the selected bandwidth (if sigma was a numeric value) or the selected fraction of the maximum bandwidth (if sigma was a pixel image or function). The result also belongs to the class "bw.optim" which can be plotted.

Author(s)

Adrian Baddeley.

See Also

bw.pplHeat for an alternative method.

densityHeat.ppp

Examples

online <- interactive()
  if(!online) op <- spatstat.options(npixel=32)
  f <- function(x,y) { dnorm(x, 2.3, 0.1) * dnorm(y, 2.0, 0.2) }
  X <- rpoint(15, f, win=letterR)
  plot(X)
  b <- bw.CvLHeat(X, sigma=0.25)
  b
  plot(b)
  if(!online) spatstat.options(op)

Description

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

Usage

bw.diggle(X, ..., correction="good", hmax=NULL, nr=512, warn=TRUE)

Arguments

Details

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

The bandwidth $\sigma$ is chosen to minimise the mean-square error criterion defined by Diggle (1985). The algorithm uses the method of Berman and Diggle (1989) to compute the quantity

$M(\sigma) = \frac{\mbox{MSE}(\sigma)}{\lambda^2} - g(0)$

as a function of bandwidth $\sigma$, where $\mbox{MSE}(\sigma)$ is the mean squared error at bandwidth $\sigma$, while $\lambda$ is the mean intensity, and $g$ is the pair correlation function. See Diggle (2003, pages 115-118) for a summary of this method.

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.

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.

Definition of bandwidth

The smoothing parameter sigma returned by bw.diggle (and displayed on the horizontal axis of the plot) corresponds to h/2, where h is the smoothing parameter described in Diggle (2003, pages 116-118) and Berman and Diggle (1989). In those references, the smoothing kernel is the uniform density on the disc of radius h. In density.ppp, the smoothing kernel is the isotropic Gaussian density with standard deviation sigma. When replacing one kernel by another, the usual practice is to adjust the bandwidths so that the kernels have equal variance (cf. Diggle 2003, page 118). This implies that sigma = h/2.

Author(s)

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

References

Berman, M. and Diggle, P. (1989) Estimating weighted integrals of the second-order intensity of a spatial point process. Journal of the Royal Statistical Society, series B 51, 81–92.

Diggle, P.J. (1985) A kernel method for smoothing point process data. Applied Statistics (Journal of the Royal Statistical Society, Series C) 34 (1985) 138–147.

Diggle, P.J. (2003) Statistical analysis of spatial point patterns, Second edition. Arnold.

See Also

density.ppp, bw.optim.object.

Alternative methods: bw.ppl, bw.scott, bw.CvL, bw.frac.

Examples

attach(split(lansing))
  b <- bw.diggle(hickory)
  plot(b, ylim=c(-2, 0), main="Cross validation for hickories")
  if(interactive()) {
   plot(density(hickory, b))
  }

Description

Select a smoothing bandwidth for smoothing a point pattern, based only on the geometry of the spatial window. The bandwidth is a specified quantile of the distance between two independent random points in the window.

Usage

bw.frac(X, ..., f=1/4)

Arguments

Details

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

The bandwidth $\sigma$ is computed as a quantile of the distance between two independent random points in the window. The default is the lower quartile of this distribution.

If $F(r)$ is the cumulative distribution function of the distance between two independent random points uniformly distributed in the window, then the value returned is the quantile with probability $f$. That is, the bandwidth is the value $r$ such that $F(r) = f$.

The cumulative distribution function $F(r)$ is computed using distcdf. We then we compute the smallest number $r$ such that $F(r) \ge f$.

Value

A numerical value giving the selected bandwidth. The result also belongs to the class "bw.frac" which can be plotted to show the cumulative distribution function and the selected quantile.

Author(s)

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

See Also

For estimating point process intensity, see density.ppp, bw.diggle, bw.ppl, bw.scott, bw.CvL.

For other smoothing purposes, see bw.stoyan, bw.smoothppp, bw.relrisk.

Examples

h <- bw.frac(letterR)
  h
  plot(h, main="bw.frac(letterR)")

Description

An object of the class "bw.optim" represents a tuning parameter (usually a smoothing bandwidth) that has been selected automatically. The object can be used as if it were a numerical value, but it can also be plotted to show the optimality criterion.

Details

An object of the class "bw.optim" represents the numerical value of a smoothing bandwidth, a threshold, or a similar tuning parameter, that has been selected by optimising a criterion such as cross-validation.

The object is a numerical value, with some attributes that retain information about how the value was selected.

Attributes include the vector of candidate values that were examined, the corresponding values of the optimality criterion, the name of the parameter, the name of the optimality criterion, and the units in which the parameter is measured.

There are methods for print, plot, summary, as.data.frame and as.fv for the class "bw.optim".

The print method simply prints the numerical value of the parameter. The summary method prints this value, and states how this value was selected.

The plot method produces a plot of the optimisation criterion against the candidate value of the parameter. The as.data.frame and as.fv methods extract this graphical information as a data frame or function table, respectively.

Author(s)

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

See Also

Functions which produce objects of class bw.optim include bw.CvL, bw.CvL.adaptive, bw.diggle, bw.lppl, bw.pcf, bw.ppl, bw.relrisk, bw.relrisk.lpp, bw.smoothppp and bw.voronoi

Examples

Ns <- if(interactive()) 32 else 3
  b <- bw.ppl(redwood, srange=c(0.02, 0.07), ns=Ns)
  b
  summary(b)
  plot(b)

Description

Uses composite likelihood or generalized least squares cross-validation to select a smoothing bandwidth for the kernel estimation of pair correlation function.

Usage

bw.pcf(X, rmax=NULL, lambda=NULL, divisor="r", 
         kernel="epanechnikov", nr=10000, bias.correct=TRUE, 
         cv.method=c("compLik", "leastSQ"), simple=TRUE, srange=NULL,
	 ..., verbose=FALSE, warn=TRUE)

Arguments

Details

This function selects an appropriate bandwidth bw for the kernel estimator of the pair correlation function of a point process intensity computed by pcf.ppp (homogeneous case) or pcfinhom (inhomogeneous case).

With cv.method="leastSQ", the bandwidth $h$ is chosen to minimise an unbiased estimate of the integrated mean-square error criterion $M(h)$ defined in equation (4) in Guan (2007a). The code implements the fast algorithm of Jalilian and Waagepetersen (2018).

With cv.method="compLik", the bandwidth $h$ is chosen to maximise a likelihood cross-validation criterion $CV(h)$ defined in equation (6) of Guan (2007b).

$M(b) = \frac{\mbox{MSE}(\sigma)}{\lambda^2} - g(0)$

The result is a numerical value giving the selected bandwidth.

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.

Definition of bandwidth

The bandwidth bw returned by bw.pcf is the standard deviation of the smoothing kernel, following the standard convention in R. As mentioned in the documentation for density.default and pcf.ppp, this differs from other definitions of bandwidth that can be found in the literature. The scale parameter h, which is called the bandwidth in some literature, is defined differently. For example for the Epanechnikov kernel, h is the half-width of the kernel, and bw=h/sqrt(5).

Author(s)

Rasmus Waagepetersen and Abdollah Jalilian. Adapted for spatstat by Adrian Baddeley [email protected], Rolf Turner [email protected] and Ege Rubak [email protected].

References

Guan, Y. (2007a). A composite likelihood cross-validation approach in selecting bandwidth for the estimation of the pair correlation function. Scandinavian Journal of Statistics, 34(2), 336–346.

Guan, Y. (2007b). A least-squares cross-validation bandwidth selection approach in pair correlation function estimations. Statistics & Probability Letters, 77(18), 1722–1729.

Jalilian, A. and Waagepetersen, R. (2018) Fast bandwidth selection for estimation of the pair correlation function. Journal of Statistical Computation and Simulation, 88(10), 2001–2011. https://www.tandfonline.com/doi/full/10.1080/00949655.2018.1428606

See Also

pcf.ppp, pcfinhom, bw.optim.object

Examples

b <- bw.pcf(redwood)
  plot(pcf(redwood, bw=b))

Description

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

Usage

bw.ppl(X, ..., srange=NULL, ns=16, sigma=NULL, varcov1=NULL,
          weights=NULL, shortcut=TRUE, warn=TRUE)

Arguments