Estimation of the Space-Time Inhomogeneous K-function

Description

Compute an estimate of the Space-Time Inhomogeneous K-function.

Usage

1
2
STIKhat(xyt, s.region, t.region, dist, times, lambda, 
correction="isotropic", infectious=FALSE) 

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 K(u,v) is computed. If missing, the maximum of dist is given by min(Sx,Sy)/4, where Sx and Sy represent the maximum width and height of s.region.

times

vector of times v at which K(u,v) is computed.

lambda

vector of values of the space-time intensity function evaluated at the points (x,y,t) in SxT. If lambda is missing, the estimate of the space-time K-function is computed as for the homogeneous case (Diggle et al., 1995), i.e. considering n/|SxT| as an estimate of the space-time intensity.

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

infectious

logical value. If TRUE, only future events are considered and the isotropic edge correction method is used. See Details.

Details

Gabriel (2014) proposes the following unbiased estimator for the STIK-function, based on data giving the locations of events xi: i=1...,n on a spatio-temporal region SxT, where S is an arbitrary polygon and T is a time interval:

K(u,v)=sum_{i=1,...,n} sum_{j != i} 1/wij 1/(lambda(x_i)lambda(x_j)) 1{uij <= u ; tj-ti <= v}

where lambda(xi) is the intensity at xi=(si,ti) and wij is an edge correction factor to deal with spatial-temporal edge effects. The edge correction methods implemented are:

isotropic: wij = |S x T| wij(s) wij(t), where the temporal edge correction factor wij(t)=1 if both ends of the interval of length 2|ti-tj| centred at ti lie within T and wij(t)=1/2 otherwise and wij(s) is the proportion of the circumference of a circle centred at the location si with radius ||si-sj|| lying in S (also called Ripley's edge correction factor).

border: wij=(sum_j 1{d(sj,S)>u ; d(tj,T)>v}/ lambda(xj)) / 1{d(si,S)>u ; d(ti,T)>v}, where d(si,S) denotes the distance between si and the boundary of S and d(ti,T) the distance between ti and the boundary of T.

modified.border: wij = |S(-u) x T(-v)| / 1{d(si,S)>u ; d(ti,T)>v}, where S(-u) and T(-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: |S intersect S(si-sj) x T intersect T(ti-tj)|, where S(si-sj) and T(ti-tj) are the translated spatial and temporal regions.

none: No edge correction is performed and |S x T|.

If parameter infectious = TRUE, ony future events are considered and the estimator is, using an isotropic edge correction factor (Gabriel and Diggle, 2009):

K(u,v) = 1/|SxT| n/nv sum_{i=1,...,nv} sum_{j=1,...,nv; j > i} 1/wij 1/(lambda(x_i)lambda(x_j)) 1{uij <= u} 1{tj - ti <= v}

In this equation, the points xi=(si,ti) are ordered so that ti < t(i+1), with ties due to round-off error broken by randomly unrounding if necessary. To deal with temporal edge-effects, for each v, nv denotes the number of events for which ti <= T1-v, with T=[T0,T1]. To deal with spatial edge-effects, we use Ripley's method.

If lambda is missing in argument, STIKhat computes an estimate of the space-time (homogeneous) K-function:

K(u,v) = |SxT|/(nv(n-1)) sum_{i=1,...,nv} sum_{j=1,...,nv; j>i} 1/wij 1{uij <= u} 1{tj - ti <= v}

Value

A list containing:

Khat

ndist x ntimes matrix containing values of K(u,v).

Ktheo

ndist x ntimes matrix containing theoretical values for a Poisson process; pi u^2 v for K and 2 pi u^2 v) for K^*.

dist, times, infectious

parameters passed in argument.

correction

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

Author(s)

Edith Gabriel <edith.gabriel@univ-avignon.fr>

References

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(1).

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.

Baddeley A, Turner R (2000) Practical maximum pseudolikelihood for spatial point patterns. Aust NZ J Stat 42, 283–322.

Baddeley A., Moller J. and Waagepetersen R. (2000). Non- and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54, 329–350.

Diggle P. , Chedwynd A., Haggkvist R. and Morris S. (1995). Second-order analysis of space-time clustering. Statistical Methods in Medical Research, 4, 124–136.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
## Not run: 
## 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=5000, ny=5000)
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 STIK function
u <- seq(0,10,by=1)
v <- seq(0,15,by=1)
stik1 <- STIKhat(xyt=FMD, s.region=northcumbria/1000,t.region=c(1,200), 
lambda=mhat*mut/dim(FMD)[1], dist=u, times=v, infectious=TRUE)

# plotting the estimation
plotK(stik1)
plotK(stik1,type="persp",theta=-65,phi=35)

## Second example
##################

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

plotK(stik2,type="contour",legend=T,which="translate")

## End(Not run)