Last chance! 50% off unlimited learning
Sale ends in
Computes a nonparametric estimate of the intensity of a point process, as a function of a (continuous) spatial covariate.
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",
"mountain", "valley", "piecewise"),
subset=NULL,
do.CI=TRUE,
jitter=TRUE, jitterfactor=1, interpolate=TRUE,
dimyx=NULL, eps=NULL,
rule.eps = c("adjust.eps", "grow.frame", "shrink.frame"),
n = 512, bw = "nrd0", adjust=1, from = NULL, to = NULL,
bwref=bw,
covname, confidence=0.95, positiveCI, breaks=NULL)
# S3 method for quad
rhohat(object, covariate, ...,
baseline=NULL, weights=NULL,
method=c("ratio", "reweight", "transform"),
horvitz=FALSE,
smoother=c("kernel", "local", "decreasing", "increasing",
"mountain", "valley", "piecewise"),
subset=NULL,
do.CI=TRUE,
jitter=TRUE, jitterfactor=1, interpolate=TRUE,
dimyx=NULL, eps=NULL,
rule.eps = c("adjust.eps", "grow.frame", "shrink.frame"),
n = 512, bw = "nrd0", adjust=1, from = NULL, to = NULL,
bwref=bw,
covname, confidence=0.95, positiveCI, breaks=NULL)
A function value table (object of class "fv"
)
containing the estimated values of
(and confidence limits) for a sequence of values of "rhohat"
which has special methods for print
, plot
and predict
.
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"
,
"slrm"
or "lppm"
).
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.
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.
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.
Character string determining the estimation method. See Details.
Logical value indicating whether to use Horvitz-Thompson weights. See Details.
Character string determining the smoothing algorithm and the type of curve that will be estimated. See Details.
Optional. A spatial window (object of class "owin"
)
specifying a subset of the data, from which the estimate should
be calculated.
Logical value specifying whether to calculate standard errors and confidence bands.
Logical value. If jitter=TRUE
(the default),
the values of the covariate at the
data points will be jittered (randomly perturbed by adding a small
amount of noise) using the function jitter
.
If jitter=FALSE
, the covariate values at the data points will
not be altered. See the section on Randomisation and discretisation.
Numeric value controlling the scale of noise added to the
covariate values at the data points when jitter=TRUE
.
Passed to the function jitter
as the argument factor
.
Logical value specifying whether to use spatial interpolation
to obtain the values of the covariate at the data points,
when the covariate is a pixel image
(object of class "im"
).
If interpolate=FALSE
, the covariate value for each data point
is simply the value of the covariate image at the pixel centre that
is nearest to the data point. If interpolate=TRUE
, the
covariate value for each data point is obtained by interpolating the
nearest pixel values using interp.im
.
Arguments controlling the pixel resolution at which the covariate will be evaluated. See Details.
Smoothing bandwidth or bandwidth rule
(passed to density.default
).
Smoothing bandwidth adjustment factor
(passed to density.default
).
Arguments passed to density.default
to
control the number and range of values at which the function
will be estimated.
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
.
Optional. Character string to use as the name of the covariate.
Confidence level for confidence intervals. A number between 0 and 1.
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.
Breakpoints for the piecewise-constant function
computed when smoother='piecewise'
.
Either a vector of numeric values specifying the breakpoints,
or a single integer specifying the number of equally-spaced
breakpoints. There is a sensible default.
Smooth estimators of
The estimated function
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
If method="ratio"
, then
If method="reweight"
, then
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
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 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 positiveCI=TRUE
, the confidence limits are always
positive, because the confidence interval is based on a normal
approximation to the estimate of positiveCI=FALSE
for smoother="kernel"
and positiveCI=TRUE
for smoother="local"
.
The nonparametric maximum likelihood estimator
of a monotone function
This estimator is chosen by specifying
smoother="increasing"
or smoother="decreasing"
.
The argument method
is ignored this case.
To compute the estimate of
If smoother="decreasing"
, each primitive step function
takes the form
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.
If smoother="valley"
then we estimate a U-shaped function.
A function
If smoother="mountain"
then we estimate a function which has
an inverted U shape. A function
Confidence intervals are not available for the unimodal estimators.
By default, rhohat
adds a small amount of random noise to the
data. This is designed to suppress the effects of
discretisation in pixel images.
This strategy means that rhohat
does not produce exactly the same result when the computation is
repeated. If you need the results to be exactly reproducible, set
jitter=FALSE
.
By default, the values of the covariate at the data points
will be randomly perturbed by adding a small amount
of noise using the function jitter
. To reduce this
effect, set jitterfactor
to a number smaller than 1. To
suppress this effect entirely, set jitter=FALSE
.
Smoothing algorithm by Adrian Baddeley Adrian.Baddeley@curtin.edu.au, Ya-Mei Chang, Yong Song, and Rolf Turner rolfturner@posteo.net.
Nonparametric maximum likelihood algorithm by Adrian Baddeley Adrian.Baddeley@curtin.edu.au.
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
covariate
, and
If object
is a point pattern, and baseline
is given,
then the intensity function is assumed to be
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
object
. A nonparametric estimator of
the relative intensity
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),
the nonparametric estimator is a kernel smoothing estimator
of do.CI=TRUE
(the default),
confidence bands are also computed, assuming a Poisson point process.
See the section on Smooth estimates.
If smoother="local"
,
the nonparametric estimator is a local regression estimator
of do.CI=TRUE
(the default),
confidence bands are also computed, assuming a Poisson point process.
See the section on Smooth estimates.
If smoother="increasing"
, we assume that
If smoother="decreasing"
, we assume that
If smoother="mountain"
, we assume that
If smoother="valley"
, we assume that
If smoother="piecewise"
, the estimate of
breaks
;
there is a sensible default. The estimate of
do.CI=TRUE
(the default),
confidence bands are computed assuming a Poisson process.
See Baddeley (2018) for a comparison of these estimation techniques
(except for "mountain"
and "valley"
).
If the argument weights
is present, then the contribution
from each data point X[i]
to the estimate of 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,rule.eps
which are passed to as.mask
.
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.
rho2hat
,
methods.rhohat
,
parres
.
See ppm
for a parametric method for the same problem.
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=4, lwd=2, add=TRUE)
rhoB <- rhohat(X, "x", method="reweight")
rhoC <- rhohat(X, "x", method="transform")
rhoI <- rhohat(X, "x", smoother="increasing")
rhoM <- rhohat(X, "x", smoother="mountain")
plot(rhoI, add=TRUE, .y ~ .x, col=6)
legend("top", lty=c(3, 1), col=c(4, 6), lwd=c(2, 1),
legend=c("true", "increasing"))
rh <- rhohat(X, "x", dimyx=32)
Run the code above in your browser using DataLab