PCFhat {stpp}

Estimation of the Space-Time Inhomogeneous Pair Correlation function

Description

Compute an estimate of the space-time pair correlation function.

Usage

PCFhat(xyt, s.region, t.region, dist, times, lambda,
ks="box", hs, kt="box", ht, correction = "isotropic") 

Arguments

xyt

Coordinates and times (x,y,t) of the point pattern.

s.region

Two-column matrix specifying polygonal region containing all data locations. If s.region is missing, the bounding box of xyt[,1:2] is considered.

t.region

Vector containing the minimum and maximum values of the time interval. If t.region is missing, the range of xyt[,3] is considered.

dist

Vector of distances u at which g(u,v) is computed. If missing, the maximum of dist is given by \min(S_x,S_y)/4, where S_x and S_y represent the maximum width and height of the bounding box of s.region.

times

Vector of times v at which g(u,v) is computed. If missing, the maximum of times is given by (T_{\max} - T_{\min})/4, where T_{\min} and T_{\max} are the minimum and maximum of the time interval T.

lambda

Vector of values of the space-time intensity function evaluated at the points (x,y,t) in S\times T. If lambda is missing, the estimate of the space-time pair correlation function is computed as for the homogeneous case, i.e. considering n/|S \times T| as an estimate of the space-time intensity.

ks

Kernel function for the spatial distances. Default is the "box" kernel. Can also be "epanech" for the Epanechnikov kernel or "gaussian" or "biweight".

hs

Bandwidth of the kernel function ks.

kt

Kernel function for the temporal distances. Default is the "box" kernel. Can also be "epanech" for the Epanechnikov kernel or "gaussian" or "biweight".

ht

Bandwidth of the kernel function kt.

correction

A character vector specifying the edge correction(s) to be applied among "isotropic", "border", "modified.border", "translate" and "none" (see Details). The default is "isotropic".

Details

An approximately unbiased estimator for the space-time pair correlation function, based on data giving the locations of events x_i: i=1,...n on a spatio-temporal region S \times T, where S is an arbitrary polygon and T a time interval:

\widehat{g}(u,v)=\frac{1}{4\pi u}\sum_{i=1}^{n}\sum_{j \neq i} \frac{1}{w_{ij}}\frac{k_{s}(u-\|s_i-s_j\|)k_{t}(v-|t_i-t_j|)}{\lambda(x_i) \lambda(x_j)},

where \lambda(x_i) is the intensity at x_i = (s_i,t_i) and w_{ij} is an edge correction factor to deal with spatial-temporal edge effects. The edge correction methods implemented are:

isotropic: w_{ij} = |S \times T| w_{ij}^{(t)} w_{ij}^{(s)}, where the temporal edge correction factor w_{ij}^{(t)} = 1 if both ends of the interval of length 2 |t_i - t_j| centred at t_i lie within T and w_{ij}^{(t)}=1/2 otherwise and w_{ij}^{(s)} is the proportion of the circumference of a circle centred at the location s_i with radius \|s_i -s_j\| lying in S (also called Ripley's edge correction factor).

border: w_{ij}=\frac{\sum_{j=1}^{n}\mathbf{1}\lbrace d(s_j,S)>u \ ; \ d(t_j,T) >v\rbrace/\lambda(x_j)}{\mathbf{1}_{\lbrace d(s_i,S) > u \ ; \ d(t_i,T) >v \rbrace}}, where d(s_i,S) denotes the distance between s_i and the boundary of S and d(t_i,T) the distance between t_i and the boundary of T.

modified.border: w_{ij} = \frac{|S_{\ominus u}|\times|T_{\ominus v}|}{\mathbf{1}_{\lbrace d(s_i,S) > u \ ; \ d(t_i,T) >v \rbrace}}, where S_{\ominus u} and T_{\ominus v} are the eroded spatial and temporal region respectively, obtained by trimming off a margin of width u and v from the border of the original region.

translate: w_{ij} =|S \cap S_{s_i-s_j}| \times |T \cap T_{t_i-t_j}|, where S_{s_i-s_j} and T_{t_i-t_j} are the translated spatial and temporal regions.

none: No edge correction is performed and w_{ij}=|S \times T|.

k_s() and k_t() denotes kernel functions with bandwidth h_s and h_t. Experience with pair correlation function estimation recommends box kernels (the default), see Illian et al. (2008). Epanechnikov, Gaussian and biweight kernels are also implemented. Whatever the kernel function, if the bandwidth is missing, a value is obtain from the function dpik of the package KernSmooth. Note that the bandwidths play an important role and their choice is crucial in the quality of the estimators as they heavily influence their variance.

Value

A list containing:

pcf

ndist x ntimes matrix containing values of \hat{g}(u,v).

pcftheo

ndist x ntimes matrix containing theoretical values for a Poisson process.

dist, times

Parameters passed in argument.

kernel

A vector of names and bandwidths of the spatial and temporal kernels.

correction

The name(s) of the edge correction method(s) passed in argument.

Author(s)

Edith Gabriel <[e-mail]>

References

Baddeley, A., Rubak, E., Turner, R., (2015). Spatial Point Patterns: Methodology and Applications with R. CRC Press, Boca Raton.

Gabriel E., Diggle P. (2009). Second-order analysis of inhomogeneous spatio-temporal point process data. Statistica Neerlandica, 63, 43–51.

Gabriel E., Rowlingson B., Diggle P. (2013). stpp: an R package for plotting, simulating and analyzing Spatio-Temporal Point Patterns. Journal of Statistical Software, 53(2), 1–29.

Gabriel E. (2014). Estimating second-order characteristics of inhomogeneous spatio-temporal point processes: influence of edge correction methods and intensity estimates. Methodology and computing in Applied Probabillity, 16(2), 411–431.

Illian JB, Penttinen A, Stoyan H and Stoyan, D. (2008). Statistical Analysis and Modelling of Spatial Point Patterns. John Wiley and Sons, London.

Examples

# First example

data(fmd)
data(northcumbria)
FMD<-as.3dpoints(fmd[,1]/1000,fmd[,2]/1000,fmd[,3])
Northcumbria=northcumbria/1000

# estimation of the temporal intensity
Mt<-density(FMD[,3],n=1000)
mut<-Mt$y[findInterval(FMD[,3],Mt$x)]*dim(FMD)[1]

# estimation of the spatial intensity
h<-mse2d(as.points(FMD[,1:2]), Northcumbria, nsmse=50, range=4)
h<-h$h[which.min(h$mse)]
Ms<-kernel2d(as.points(FMD[,1:2]), Northcumbria, h, nx=500, ny=500)
atx<-findInterval(x=FMD[,1],vec=Ms$x)
aty<-findInterval(x=FMD[,2],vec=Ms$y)
mhat<-NULL
for(i in 1:length(atx)) mhat<-c(mhat,Ms$z[atx[i],aty[i]])

# estimation of the pair correlation function
g1 <- PCFhat(xyt=FMD, dist=1:15, times=1:15, lambda=mhat*mut/dim(FMD)[1],
 s.region=northcumbria/1000,t.region=c(1,200))

# plotting the estimation 

plotPCF(g1)
plotPCF(g1,type="persp",theta=-65,phi=35) 


# Second example

xyt=rpp(lambda=200)
g2=PCFhat(xyt$xyt,dist=seq(0,0.16,by=0.02),
times=seq(0,0.16,by=0.02),correction=c("border","translate"))

plotPCF(g2,type="contour",which="border")