spatstat (version 1.64-1)

rhohat: Nonparametric Estimate of Intensity as Function of a Covariate

Description

Computes a nonparametric estimate of the intensity of a point process, as a function of a (continuous) spatial covariate.

Usage

rhohat(object, covariate, ...)

# S3 method for ppp rhohat(object, covariate, ..., baseline=NULL, weights=NULL, method=c("ratio", "reweight", "transform"), horvitz=FALSE, smoother=c("kernel", "local", "decreasing", "increasing"), subset=NULL, dimyx=NULL, eps=NULL, n = 512, bw = "nrd0", adjust=1, from = NULL, to = NULL, bwref=bw, covname, confidence=0.95, positiveCI)

# S3 method for quad rhohat(object, covariate, ..., baseline=NULL, weights=NULL, method=c("ratio", "reweight", "transform"), horvitz=FALSE, smoother=c("kernel", "local", "decreasing", "increasing"), subset=NULL, dimyx=NULL, eps=NULL, n = 512, bw = "nrd0", adjust=1, from = NULL, to = NULL, bwref=bw, covname, confidence=0.95, positiveCI)

# S3 method for ppm rhohat(object, covariate, ..., weights=NULL, method=c("ratio", "reweight", "transform"), horvitz=FALSE, smoother=c("kernel", "local", "decreasing", "increasing"), subset=NULL, dimyx=NULL, eps=NULL, n = 512, bw = "nrd0", adjust=1, from = NULL, to = NULL, bwref=bw, covname, confidence=0.95, positiveCI)

# S3 method for lpp rhohat(object, covariate, ..., weights=NULL, method=c("ratio", "reweight", "transform"), horvitz=FALSE, smoother=c("kernel", "local", "decreasing", "increasing"), subset=NULL, nd=1000, eps=NULL, random=TRUE, n = 512, bw = "nrd0", adjust=1, from = NULL, to = NULL, bwref=bw, covname, confidence=0.95, positiveCI)

# S3 method for lppm rhohat(object, covariate, ..., weights=NULL, method=c("ratio", "reweight", "transform"), horvitz=FALSE, smoother=c("kernel", "local", "decreasing", "increasing"), subset=NULL, nd=1000, eps=NULL, random=TRUE, n = 512, bw = "nrd0", adjust=1, from = NULL, to = NULL, bwref=bw, covname, confidence=0.95, positiveCI)

Arguments

object

A point pattern (object of class "ppp" or "lpp"), a quadrature scheme (object of class "quad") or a fitted point process model (object of class "ppm" or "lppm").

covariate

Either a function(x,y) or a pixel image (object of class "im") providing the values of the covariate at any location. Alternatively one of the strings "x" or "y" signifying the Cartesian coordinates.

weights

Optional weights attached to the data points. Either a numeric vector of weights for each data point, or a pixel image (object of class "im") or a function(x,y) providing the weights.

baseline

Optional baseline for intensity function. A function(x,y) or a pixel image (object of class "im") providing the values of the baseline at any location.

method

Character string determining the smoothing method. See Details.

horvitz

Logical value indicating whether to use Horvitz-Thompson weights. See Details.

smoother

Character string determining the smoothing algorithm. See Details.

subset

Optional. A spatial window (object of class "owin") specifying a subset of the data, from which the estimate should be calculated.

dimyx,eps,nd,random

Arguments controlling the pixel resolution at which the covariate will be evaluated. See Details.

bw

Smoothing bandwidth or bandwidth rule (passed to density.default).

adjust

Smoothing bandwidth adjustment factor (passed to density.default).

n, from, to

Arguments passed to density.default to control the number and range of values at which the function will be estimated.

bwref

Optional. An alternative value of bw to use when smoothing the reference density (the density of the covariate values observed at all locations in the window).

Additional arguments passed to density.default or locfit.

covname

Optional. Character string to use as the name of the covariate.

confidence

Confidence level for confidence intervals. A number between 0 and 1.

positiveCI

Logical value. If TRUE, confidence limits are always positive numbers; if FALSE, the lower limit of the confidence interval may sometimes be negative. Default is FALSE if smoother="kernel" and TRUE if smoother="local". See Details.

Value

A function value table (object of class "fv") containing the estimated values of \(\rho\) (and confidence limits) for a sequence of values of \(Z\). Also belongs to the class "rhohat" which has special methods for print, plot and predict.

Smooth estimates

Smooth estimators of \(\rho(z)\) were proposed by Baddeley and Turner (2005) and Baddeley et al (2012). Similar estimators were proposed by Guan (2008) and in the literature on relative distributions (Handcock and Morris, 1999).

The estimated function \(\rho(z)\) will be a smooth function of \(z\).

The smooth estimation procedure involves computing several density estimates and combining them. The algorithm used to compute density estimates is determined by smoother:

  • If smoother="kernel", the smoothing procedure is based on fixed-bandwidth kernel density estimation, performed by density.default.

  • If smoother="local", the smoothing procedure is based on local likelihood density estimation, performed by locfit.

The argument method determines how the density estimates will be combined to obtain an estimate of \(\rho(z)\):

  • If method="ratio", then \(\rho(z)\) is estimated by the ratio of two density estimates, The numerator is a (rescaled) density estimate obtained by smoothing the values \(Z(y_i)\) of the covariate \(Z\) observed at the data points \(y_i\). The denominator is a density estimate of the reference distribution of \(Z\). See Baddeley et al (2012), equation (8). This is similar but not identical to an estimator proposed by Guan (2008).

  • If method="reweight", then \(\rho(z)\) is estimated by applying density estimation to the values \(Z(y_i)\) of the covariate \(Z\) observed at the data points \(y_i\), with weights inversely proportional to the reference density of \(Z\). See Baddeley et al (2012), equation (9).

  • If method="transform", the smoothing method is variable-bandwidth kernel smoothing, implemented by applying the Probability Integral Transform to the covariate values, yielding values in the range 0 to 1, then applying edge-corrected density estimation on the interval \([0,1]\), and back-transforming. See Baddeley et al (2012), equation (10).

If horvitz=TRUE, then the calculations described above are modified by using Horvitz-Thompson weighting. The contribution to the numerator from each data point is weighted by the reciprocal of the baseline value or fitted intensity value at that data point; and a corresponding adjustment is made to the denominator.

Pointwise confidence intervals for the true value of \(\rho(z)\) are also calculated for each \(z\), and will be plotted as grey shading. The confidence intervals are derived using the central limit theorem, based on variance calculations which assume a Poisson point process. If positiveCI=FALSE, the lower limit of the confidence interval may sometimes be negative, because the confidence intervals are based on a normal approximation to the estimate of \(\rho(z)\). If positiveCI=TRUE, the confidence limits are always positive, because the confidence interval is based on a normal approximation to the estimate of \(\log(\rho(z))\). For consistency with earlier versions, the default is positiveCI=FALSE for smoother="kernel" and positiveCI=TRUE for smoother="local".

Monotone estimates

The nonparametric maximum likelihood estimator of a monotone function \(\rho(z)\) was described by Sager (1982). This method assumes that \(\rho(z)\) is either an increasing function of \(z\), or a decreasing function of \(z\). The estimated function will be a step function, increasing or decreasing as a function of \(z\).

This estimator is chosen by specifying smoother="increasing" or smoother="decreasing". The argument method is ignored this case.

To compute the estimate of \(\rho(z)\), the algorithm first computes several primitive step-function estimates, and then takes the maximum of these primitive functions.

If smoother="decreasing", each primitive step function takes the form \(\rho(z) = \lambda\) when \(z \le t\), and \(\rho(z) = 0\) when \(z > t\), where and \(\lambda\) is a primitive estimate of intensity based on the data for \(Z \le t\). The jump location \(t\) will be the value of the covariate \(Z\) at one of the data points. The primitive estimate \(\lambda\) is the average intensity (number of points divided by area) for the region of space where the covariate value is less than or equal to \(t\).

If horvitz=TRUE, then the calculations described above are modified by using Horvitz-Thompson weighting. The contribution to the numerator from each data point is weighted by the reciprocal of the baseline value or fitted intensity value at that data point; and a corresponding adjustment is made to the denominator.

Confidence intervals are not available for the monotone estimators.

Details

This command estimates the relationship between point process intensity and a given spatial covariate. Such a relationship is sometimes called a resource selection function (if the points are organisms and the covariate is a descriptor of habitat) or a prospectivity index (if the points are mineral deposits and the covariate is a geological variable). This command uses nonparametric methods which do not assume a particular form for the relationship.

If object is a point pattern, and baseline is missing or null, this command assumes that object is a realisation of a point process with intensity function \(\lambda(u)\) of the form $$\lambda(u) = \rho(Z(u))$$ where \(Z\) is the spatial covariate function given by covariate, and \(\rho(z)\) is the resource selection function or prospectivity index. A nonparametric estimator of the function \(\rho(z)\) is computed.

If object is a point pattern, and baseline is given, then the intensity function is assumed to be $$\lambda(u) = \rho(Z(u)) B(u)$$ where \(B(u)\) is the baseline intensity at location \(u\). A nonparametric estimator of the relative intensity \(\rho(z)\) is computed.

If object is a fitted point process model, suppose X is the original data point pattern to which the model was fitted. Then this command assumes X is a realisation of a Poisson point process with intensity function of the form $$ \lambda(u) = \rho(Z(u)) \kappa(u) $$ where \(\kappa(u)\) is the intensity of the fitted model object. A nonparametric estimator of the relative intensity \(\rho(z)\) is computed.

The nonparametric estimation procedure is controlled by the arguments smoother, method and horvitz.

The argument smoother selects the type of estimation technique.

  • If smoother="kernel" (the default) or smoother="local", the nonparametric estimator is a smoothing estimator of \(\rho(z)\), effectively a kind of density estimator (Baddeley et al, 2012). The estimated function \(\rho(z)\) will be a smooth function of \(z\). Confidence bands are also computed, assuming a Poisson point process. See the section on Smooth estimates.

  • If smoother="increasing" or smoother="decreasing", we use the nonparametric maximum likelihood estimator of \(\rho(z)\) described by Sager (1982). This assumes that \(\rho(z)\) is either an increasing function of \(z\), or a decreasing function of \(z\). The estimated function will be a step function, increasing or decreasing as a function of \(z\). See the section on Monotone estimates.

See Baddeley (2018) for a comparison of these estimation techniques.

If the argument weights is present, then the contribution from each data point X[i] to the estimate of \(\rho\) is multiplied by weights[i].

If the argument subset is present, then the calculations are performed using only the data inside this spatial region.

This technique assumes that covariate has continuous values. It is not applicable to covariates with categorical (factor) values or discrete values such as small integers. For a categorical covariate, use intensity.quadratcount applied to the result of quadratcount(X, tess=covariate).

The argument covariate should be a pixel image, or a function, or one of the strings "x" or "y" signifying the cartesian coordinates. It will be evaluated on a fine grid of locations, with spatial resolution controlled by the arguments dimyx,eps,nd,random. In two dimensions (i.e. if object is of class "ppp", "ppm" or "quad") the arguments dimyx, eps are passed to as.mask to control the pixel resolution. On a linear network (i.e. if object is of class "lpp") the argument nd specifies the total number of test locations on the linear network, eps specifies the linear separation between test locations, and random specifies whether the test locations have a randomised starting position.

References

Baddeley, A., Chang, Y.-M., Song, Y. and Turner, R. (2012) Nonparametric estimation of the dependence of a point process on spatial covariates. Statistics and Its Interface 5 (2), 221--236.

Baddeley, A. and Turner, R. (2005) Modelling spatial point patterns in R. In: A. Baddeley, P. Gregori, J. Mateu, R. Stoica, and D. Stoyan, editors, Case Studies in Spatial Point Pattern Modelling, Lecture Notes in Statistics number 185. Pages 23--74. Springer-Verlag, New York, 2006. ISBN: 0-387-28311-0.

Baddeley, A. (2018) A statistical commentary on mineral prospectivity analysis. Chapter 2, pages 25--65 in Handbook of Mathematical Geosciences: Fifty Years of IAMG, edited by B.S. Daya Sagar, Q. Cheng and F.P. Agterberg. Springer, Berlin.

Guan, Y. (2008) On consistent nonparametric intensity estimation for inhomogeneous spatial point processes. Journal of the American Statistical Association 103, 1238--1247.

Handcock, M.S. and Morris, M. (1999) Relative Distribution Methods in the Social Sciences. Springer, New York.

Sager, T.W. (1982) Nonparametric maximum likelihood estimation of spatial patterns. Annals of Statistics 10, 1125--1136.

See Also

rho2hat, methods.rhohat, parres.

See ppm for a parametric method for the same problem.

Examples

Run this code
# NOT RUN {
  X <-  rpoispp(function(x,y){exp(3+3*x)})
  rho <- rhohat(X, "x")
  rho <- rhohat(X, function(x,y){x})
  plot(rho)
  curve(exp(3+3*x), lty=3, col=2, add=TRUE)

  rhoB <- rhohat(X, "x", method="reweight")
  rhoC <- rhohat(X, "x", method="transform")

  rhoM <- rhohat(X, "x", smoother="increasing")
  plot(rhoM, add=TRUE, col=5)

  
# }
# NOT RUN {
  fit <- ppm(X, ~x)
  rr <- rhohat(fit, "y")

# linear network
  Y <- runiflpp(30, simplenet)
  rhoY <- rhohat(Y, "y")
# }

Run the code above in your browser using DataCamp Workspace