Using R package spatstat for modeling point process

spatstat: Spatial Point Pattern Analysis, Model-Fitting, Simulation, Tests The Comprehensive toolbox for analyzing spatial data, mainly Spatial Point Patterns, including multi-type/marked points and spatial covariates, in any two-dimensional spatial region. Also supports three-dimensional point patterns, space-time point patterns in any number of dimensions, and point patterns on a linear network.

setwd("D:/_Lab_Data/R/R_Labs")
library(sp);library(rgdal); library(splancs)
library(spatstat) 
data <- read.table("pts_data/tpedata.csv", header=TRUE, sep=",")
data_bnd <- read.table("pts_data/tpe_sqr_bnd2.csv", header=TRUE, sep=",")
Pts_bnd <- as.points(data_bnd[,2], data_bnd[,3])
asthma.pts <- readOGR(dsn = "Colorado_Geog_Data/asthma", layer = "spasthma")
## OGR data source with driver: ESRI Shapefile 
## Source: "Colorado_Geog_Data/asthma", layer: "spasthma"
## with 1291 features
## It has 9 fields
asthma.bdry <- readOGR(dsn = "Colorado_Geog_Data/asthma", layer = "spbdry")
## OGR data source with driver: ESRI Shapefile 
## Source: "Colorado_Geog_Data/asthma", layer: "spbdry"
## with 1 features
## It has 1 fields

The homogeneous Poisson process counts events that occur at a constant rat. This process is characterized by a rate parameter λ, also known as intensity, such that the number of events in time interval (t, t + τ] follows a Poisson distribution with associated parameter λτ.

Homogeneous Poisson process: * Complete spatial randomness (CSR); * Number of points in disjoint areas independent; * Number of points follows a Poisson distribution;

npts(Pts_bnd) # 4159 points
## [1] 4159
x<-csr(Pts_bnd, 100); plot(x)

Inhomogeneous Poisson process: The inhomogeneous case differs from the homogeneous case by the intensity, which λ is not constant anymore.

### simulating Inhomogeneous Poisson process
intenfun <- function(x, y) {200 * x^2}
plot(rpoispp(intenfun, lmax = 1000), main = "Inhomogeneous Poisson process")

plot(density(rpoispp(intenfun, lmax = 1000)), main = "Density")

Capabilities of ‘spatstat’:

1. Creation and plotting of point patterns

# Importing Point data
x_min<-min(data[,2]);x_max<-max(data[,2])
y_min<-min(data[,3]);y_max<-max(data[,3])

# ppp(x.coordinates, y.coordinates, x.range, y.range)
mypattern<-ppp(data[,2], data[,3], c(x_min,x_max), c(y_min,y_max)) 
## Warning in ppp(data[, 2], data[, 3], c(x_min, x_max), c(y_min, y_max)):
## data contain duplicated points
plot(mypattern)

# Using administration boundary of Taipei City as boundary file
x1<-rev(data_bnd[,2]) # reverse the vector of X coord
y1<-rev(data_bnd[,3]) # reverse the vector of Y coord
xx<-cbind(x1, y1)
# building WINDOW file
PTS_bnd2 <- owin(poly=xx) # owin: The vertices must be listed anticlockwise.
# ppp(x.coordinates, y.coordinates, window)
mypattern2<-ppp(data[,2], data[,3], window=PTS_bnd2) 
## Warning in ppp(data[, 2], data[, 3], window = PTS_bnd2): data contain
## duplicated points
plot(mypattern2)

2. Exploratory data analysis:

standard empirical summaries of the data, such as the average intensity, the K function Ripley (1977) and the kernel-smoothed intensity map, can easily be generated and displayed. Many other empirical statistics are implemented in the package, including the empty space function F, nearest neighbor distance function G, pair correlation function g, inhomogeneous K function, second moment measure, Bartlett spectrum, cross-K function, cross-G function, J-function, and mark correlation function.

# Exploratory data analysis: Summary functions: G/F/K

dismatrix<-pairdist(mypattern2)  #matrix of pairwise distances
nnd<-nndist(mypattern2)  # vector of nearest neighbour distances
emp<-distmap(mypattern2) # distances to the point pattern 
plot(emp)

Fc <- Fest(mypattern2); plot(Fc)

##      lty col  key                   label
## km     1   1   km   italic(hat(F)[km](r))
## rs     2   2   rs italic(hat(F)[bord](r))
## cs     3   3   cs   italic(hat(F)[cs](r))
## theo   4   4 theo      italic(F[pois](r))
##                                meaning
## km       Kaplan-Meier estimate of F(r)
## rs   border corrected estimate of F(r)
## cs        Chiu-Stoyan estimate of F(r)
## theo          theoretical Poisson F(r)
Gc <- Gest(mypattern2); plot(Gc)

##      lty col  key                   label
## km     1   1   km   italic(hat(G)[km](r))
## rs     2   2   rs italic(hat(G)[bord](r))
## han    3   3  han  italic(hat(G)[han](r))
## theo   4   4 theo      italic(G[pois](r))
##                                meaning
## km       Kaplan-Meier estimate of G(r)
## rs   border corrected estimate of G(r)
## han           Hanisch estimate of G(r)
## theo          theoretical Poisson G(r)
Kc <- Kest(mypattern2); plot(Kc)

##        lty col    key                    label
## iso      1   1    iso   italic(hat(K)[iso](r))
## trans    2   2  trans italic(hat(K)[trans](r))
## border   3   3 border  italic(hat(K)[bord](r))
## theo     4   4   theo       italic(K[pois](r))
##                                             meaning
## iso    Ripley isotropic correction estimate of K(r)
## trans        translation-corrected estimate of K(r)
## border            border-corrected estimate of K(r)
## theo                       theoretical Poisson K(r)
Kenv<-envelope(mypattern2, Kest, nrank=2, nsim=99) # 95 C.I
## Generating 99 simulations of CSR  ...
## 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, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(Kenv)

##      lty col  key                  label
## obs    1   1  obs italic(hat(K)[obs](r))
## theo   2   2 theo     italic(K[theo](r))
## hi     1   8   hi  italic(hat(K)[hi](r))
## lo     1   8   lo  italic(hat(K)[lo](r))
##                                                meaning
## obs            observed value of K(r) for data pattern
## theo                 theoretical value of K(r) for CSR
## hi   upper pointwise envelope of K(r) from simulations
## lo   lower pointwise envelope of K(r) from simulations
Lc <- Lest(mypattern2); plot(Lc)

##        lty col    key                    label
## iso      1   1    iso   italic(hat(L)[iso](r))
## trans    2   2  trans italic(hat(L)[trans](r))
## border   3   3 border  italic(hat(L)[bord](r))
## theo     4   4   theo       italic(L[pois](r))
##                                             meaning
## iso    Ripley isotropic correction estimate of L(r)
## trans        translation-corrected estimate of L(r)
## border            border-corrected estimate of L(r)
## theo                       theoretical Poisson L(r)

The local K-function

The statistical method was developed by Getis (1984). It is similar to the global K-function in analysis, but differs in that the local K-function only considers those pairs of points having a given point i as one of its members. Like the K-function, the Local K-function is a test of the hypothesis of CSR, and the expected value of Li(d) is d. Again, the confidence envelope is generated by performing a specified number of simulations. If for any distance the observed Li(d) falls outside the confidence envelope, the hypothesis of CSR can be rejected at the appropriate significance level. An observed Li(d) below the envelope indicates that the points are dispersed about point i for distance d. Conversely, an observed Li(d) above the envelope indicates that the points are clustered about point i for distance d.

# Local K function
LK<-localK(mypattern2,correction = "Ripley"); plot(LK)
## 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, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62,  63.

##       lty col   key                        label
## theo    1   1  theo italic(K[loc][","][pois](r))
## iso63   2   2 iso63   italic(K[loc][","][63](r))
## iso62   3   3 iso62   italic(K[loc][","][62](r))
## iso61   4   4 iso61   italic(K[loc][","][61](r))
## iso60   5   5 iso60   italic(K[loc][","][60](r))
## iso59   6   6 iso59   italic(K[loc][","][59](r))
## iso58   7   7 iso58   italic(K[loc][","][58](r))
## iso57   8   8 iso57   italic(K[loc][","][57](r))
## iso56   9   9 iso56   italic(K[loc][","][56](r))
## iso55  10  10 iso55   italic(K[loc][","][55](r))
## iso54  11  11 iso54   italic(K[loc][","][54](r))
## iso53  12  12 iso53   italic(K[loc][","][53](r))
## iso52  13  13 iso52   italic(K[loc][","][52](r))
## iso51  14  14 iso51   italic(K[loc][","][51](r))
## iso50  15  15 iso50   italic(K[loc][","][50](r))
## iso49  16  16 iso49   italic(K[loc][","][49](r))
## iso48  17  17 iso48   italic(K[loc][","][48](r))
## iso47  18  18 iso47   italic(K[loc][","][47](r))
## iso46  19  19 iso46   italic(K[loc][","][46](r))
## iso45  20  20 iso45   italic(K[loc][","][45](r))
## iso44  21  21 iso44   italic(K[loc][","][44](r))
## iso43  22  22 iso43   italic(K[loc][","][43](r))
## iso42  23  23 iso42   italic(K[loc][","][42](r))
## iso41  24  24 iso41   italic(K[loc][","][41](r))
## iso40  25  25 iso40   italic(K[loc][","][40](r))
## iso39  26  26 iso39   italic(K[loc][","][39](r))
## iso38  27  27 iso38   italic(K[loc][","][38](r))
## iso37  28  28 iso37   italic(K[loc][","][37](r))
## iso36  29  29 iso36   italic(K[loc][","][36](r))
## iso35  30  30 iso35   italic(K[loc][","][35](r))
## iso34  31  31 iso34   italic(K[loc][","][34](r))
## iso33  32  32 iso33   italic(K[loc][","][33](r))
## iso32  33  33 iso32   italic(K[loc][","][32](r))
## iso31  34  34 iso31   italic(K[loc][","][31](r))
## iso30  35  35 iso30   italic(K[loc][","][30](r))
## iso29  36  36 iso29   italic(K[loc][","][29](r))
## iso28  37  37 iso28   italic(K[loc][","][28](r))
## iso27  38  38 iso27   italic(K[loc][","][27](r))
## iso26  39  39 iso26   italic(K[loc][","][26](r))
## iso25  40  40 iso25   italic(K[loc][","][25](r))
## iso24  41  41 iso24   italic(K[loc][","][24](r))
## iso23  42  42 iso23   italic(K[loc][","][23](r))
## iso22  43  43 iso22   italic(K[loc][","][22](r))
## iso21  44  44 iso21   italic(K[loc][","][21](r))
## iso20  45  45 iso20   italic(K[loc][","][20](r))
## iso19  46  46 iso19   italic(K[loc][","][19](r))
## iso18  47  47 iso18   italic(K[loc][","][18](r))
## iso17  48  48 iso17   italic(K[loc][","][17](r))
## iso16  49  49 iso16   italic(K[loc][","][16](r))
## iso15  50  50 iso15   italic(K[loc][","][15](r))
## iso14  51  51 iso14   italic(K[loc][","][14](r))
## iso13  52  52 iso13   italic(K[loc][","][13](r))
## iso12  53  53 iso12   italic(K[loc][","][12](r))
## iso11  54  54 iso11   italic(K[loc][","][11](r))
## iso10  55  55 iso10   italic(K[loc][","][10](r))
## iso09  56  56 iso09    italic(K[loc][","][9](r))
## iso08  57  57 iso08    italic(K[loc][","][8](r))
## iso07  58  58 iso07    italic(K[loc][","][7](r))
## iso06  59  59 iso06    italic(K[loc][","][6](r))
## iso05  60  60 iso05    italic(K[loc][","][5](r))
## iso04  61  61 iso04    italic(K[loc][","][4](r))
## iso03  62  62 iso03    italic(K[loc][","][3](r))
## iso02  63  63 iso02    italic(K[loc][","][2](r))
## iso01  64  64 iso01    italic(K[loc][","][1](r))
##                                                              meaning
## theo                                   theoretical Poisson K[loc](r)
## iso63 Ripley isotropic correction estimate of K[loc](r) for point 63
## iso62 Ripley isotropic correction estimate of K[loc](r) for point 62
## iso61 Ripley isotropic correction estimate of K[loc](r) for point 61
## iso60 Ripley isotropic correction estimate of K[loc](r) for point 60
## iso59 Ripley isotropic correction estimate of K[loc](r) for point 59
## iso58 Ripley isotropic correction estimate of K[loc](r) for point 58
## iso57 Ripley isotropic correction estimate of K[loc](r) for point 57
## iso56 Ripley isotropic correction estimate of K[loc](r) for point 56
## iso55 Ripley isotropic correction estimate of K[loc](r) for point 55
## iso54 Ripley isotropic correction estimate of K[loc](r) for point 54
## iso53 Ripley isotropic correction estimate of K[loc](r) for point 53
## iso52 Ripley isotropic correction estimate of K[loc](r) for point 52
## iso51 Ripley isotropic correction estimate of K[loc](r) for point 51
## iso50 Ripley isotropic correction estimate of K[loc](r) for point 50
## iso49 Ripley isotropic correction estimate of K[loc](r) for point 49
## iso48 Ripley isotropic correction estimate of K[loc](r) for point 48
## iso47 Ripley isotropic correction estimate of K[loc](r) for point 47
## iso46 Ripley isotropic correction estimate of K[loc](r) for point 46
## iso45 Ripley isotropic correction estimate of K[loc](r) for point 45
## iso44 Ripley isotropic correction estimate of K[loc](r) for point 44
## iso43 Ripley isotropic correction estimate of K[loc](r) for point 43
## iso42 Ripley isotropic correction estimate of K[loc](r) for point 42
## iso41 Ripley isotropic correction estimate of K[loc](r) for point 41
## iso40 Ripley isotropic correction estimate of K[loc](r) for point 40
## iso39 Ripley isotropic correction estimate of K[loc](r) for point 39
## iso38 Ripley isotropic correction estimate of K[loc](r) for point 38
## iso37 Ripley isotropic correction estimate of K[loc](r) for point 37
## iso36 Ripley isotropic correction estimate of K[loc](r) for point 36
## iso35 Ripley isotropic correction estimate of K[loc](r) for point 35
## iso34 Ripley isotropic correction estimate of K[loc](r) for point 34
## iso33 Ripley isotropic correction estimate of K[loc](r) for point 33
## iso32 Ripley isotropic correction estimate of K[loc](r) for point 32
## iso31 Ripley isotropic correction estimate of K[loc](r) for point 31
## iso30 Ripley isotropic correction estimate of K[loc](r) for point 30
## iso29 Ripley isotropic correction estimate of K[loc](r) for point 29
## iso28 Ripley isotropic correction estimate of K[loc](r) for point 28
## iso27 Ripley isotropic correction estimate of K[loc](r) for point 27
## iso26 Ripley isotropic correction estimate of K[loc](r) for point 26
## iso25 Ripley isotropic correction estimate of K[loc](r) for point 25
## iso24 Ripley isotropic correction estimate of K[loc](r) for point 24
## iso23 Ripley isotropic correction estimate of K[loc](r) for point 23
## iso22 Ripley isotropic correction estimate of K[loc](r) for point 22
## iso21 Ripley isotropic correction estimate of K[loc](r) for point 21
## iso20 Ripley isotropic correction estimate of K[loc](r) for point 20
## iso19 Ripley isotropic correction estimate of K[loc](r) for point 19
## iso18 Ripley isotropic correction estimate of K[loc](r) for point 18
## iso17 Ripley isotropic correction estimate of K[loc](r) for point 17
## iso16 Ripley isotropic correction estimate of K[loc](r) for point 16
## iso15 Ripley isotropic correction estimate of K[loc](r) for point 15
## iso14 Ripley isotropic correction estimate of K[loc](r) for point 14
## iso13 Ripley isotropic correction estimate of K[loc](r) for point 13
## iso12 Ripley isotropic correction estimate of K[loc](r) for point 12
## iso11 Ripley isotropic correction estimate of K[loc](r) for point 11
## iso10 Ripley isotropic correction estimate of K[loc](r) for point 10
## iso09 Ripley isotropic correction estimate of K[loc](r) for point 09
## iso08 Ripley isotropic correction estimate of K[loc](r) for point 08
## iso07 Ripley isotropic correction estimate of K[loc](r) for point 07
## iso06 Ripley isotropic correction estimate of K[loc](r) for point 06
## iso05 Ripley isotropic correction estimate of K[loc](r) for point 05
## iso04 Ripley isotropic correction estimate of K[loc](r) for point 04
## iso03 Ripley isotropic correction estimate of K[loc](r) for point 03
## iso02 Ripley isotropic correction estimate of K[loc](r) for point 02
## iso01 Ripley isotropic correction estimate of K[loc](r) for point 01
LL<-localL(mypattern2,correction = "Ripley"); plot(LL)
## 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, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62,  63.

##       lty col   key                        label
## theo    1   1  theo italic(L[loc][","][pois](r))
## iso63   2   2 iso63   italic(L[loc][","][63](r))
## iso62   3   3 iso62   italic(L[loc][","][62](r))
## iso61   4   4 iso61   italic(L[loc][","][61](r))
## iso60   5   5 iso60   italic(L[loc][","][60](r))
## iso59   6   6 iso59   italic(L[loc][","][59](r))
## iso58   7   7 iso58   italic(L[loc][","][58](r))
## iso57   8   8 iso57   italic(L[loc][","][57](r))
## iso56   9   9 iso56   italic(L[loc][","][56](r))
## iso55  10  10 iso55   italic(L[loc][","][55](r))
## iso54  11  11 iso54   italic(L[loc][","][54](r))
## iso53  12  12 iso53   italic(L[loc][","][53](r))
## iso52  13  13 iso52   italic(L[loc][","][52](r))
## iso51  14  14 iso51   italic(L[loc][","][51](r))
## iso50  15  15 iso50   italic(L[loc][","][50](r))
## iso49  16  16 iso49   italic(L[loc][","][49](r))
## iso48  17  17 iso48   italic(L[loc][","][48](r))
## iso47  18  18 iso47   italic(L[loc][","][47](r))
## iso46  19  19 iso46   italic(L[loc][","][46](r))
## iso45  20  20 iso45   italic(L[loc][","][45](r))
## iso44  21  21 iso44   italic(L[loc][","][44](r))
## iso43  22  22 iso43   italic(L[loc][","][43](r))
## iso42  23  23 iso42   italic(L[loc][","][42](r))
## iso41  24  24 iso41   italic(L[loc][","][41](r))
## iso40  25  25 iso40   italic(L[loc][","][40](r))
## iso39  26  26 iso39   italic(L[loc][","][39](r))
## iso38  27  27 iso38   italic(L[loc][","][38](r))
## iso37  28  28 iso37   italic(L[loc][","][37](r))
## iso36  29  29 iso36   italic(L[loc][","][36](r))
## iso35  30  30 iso35   italic(L[loc][","][35](r))
## iso34  31  31 iso34   italic(L[loc][","][34](r))
## iso33  32  32 iso33   italic(L[loc][","][33](r))
## iso32  33  33 iso32   italic(L[loc][","][32](r))
## iso31  34  34 iso31   italic(L[loc][","][31](r))
## iso30  35  35 iso30   italic(L[loc][","][30](r))
## iso29  36  36 iso29   italic(L[loc][","][29](r))
## iso28  37  37 iso28   italic(L[loc][","][28](r))
## iso27  38  38 iso27   italic(L[loc][","][27](r))
## iso26  39  39 iso26   italic(L[loc][","][26](r))
## iso25  40  40 iso25   italic(L[loc][","][25](r))
## iso24  41  41 iso24   italic(L[loc][","][24](r))
## iso23  42  42 iso23   italic(L[loc][","][23](r))
## iso22  43  43 iso22   italic(L[loc][","][22](r))
## iso21  44  44 iso21   italic(L[loc][","][21](r))
## iso20  45  45 iso20   italic(L[loc][","][20](r))
## iso19  46  46 iso19   italic(L[loc][","][19](r))
## iso18  47  47 iso18   italic(L[loc][","][18](r))
## iso17  48  48 iso17   italic(L[loc][","][17](r))
## iso16  49  49 iso16   italic(L[loc][","][16](r))
## iso15  50  50 iso15   italic(L[loc][","][15](r))
## iso14  51  51 iso14   italic(L[loc][","][14](r))
## iso13  52  52 iso13   italic(L[loc][","][13](r))
## iso12  53  53 iso12   italic(L[loc][","][12](r))
## iso11  54  54 iso11   italic(L[loc][","][11](r))
## iso10  55  55 iso10   italic(L[loc][","][10](r))
## iso09  56  56 iso09    italic(L[loc][","][9](r))
## iso08  57  57 iso08    italic(L[loc][","][8](r))
## iso07  58  58 iso07    italic(L[loc][","][7](r))
## iso06  59  59 iso06    italic(L[loc][","][6](r))
## iso05  60  60 iso05    italic(L[loc][","][5](r))
## iso04  61  61 iso04    italic(L[loc][","][4](r))
## iso03  62  62 iso03    italic(L[loc][","][3](r))
## iso02  63  63 iso02    italic(L[loc][","][2](r))
## iso01  64  64 iso01    italic(L[loc][","][1](r))
##                                                              meaning
## theo                                   theoretical Poisson L[loc](r)
## iso63 Ripley isotropic correction estimate of L[loc](r) for point 63
## iso62 Ripley isotropic correction estimate of L[loc](r) for point 62
## iso61 Ripley isotropic correction estimate of L[loc](r) for point 61
## iso60 Ripley isotropic correction estimate of L[loc](r) for point 60
## iso59 Ripley isotropic correction estimate of L[loc](r) for point 59
## iso58 Ripley isotropic correction estimate of L[loc](r) for point 58
## iso57 Ripley isotropic correction estimate of L[loc](r) for point 57
## iso56 Ripley isotropic correction estimate of L[loc](r) for point 56
## iso55 Ripley isotropic correction estimate of L[loc](r) for point 55
## iso54 Ripley isotropic correction estimate of L[loc](r) for point 54
## iso53 Ripley isotropic correction estimate of L[loc](r) for point 53
## iso52 Ripley isotropic correction estimate of L[loc](r) for point 52
## iso51 Ripley isotropic correction estimate of L[loc](r) for point 51
## iso50 Ripley isotropic correction estimate of L[loc](r) for point 50
## iso49 Ripley isotropic correction estimate of L[loc](r) for point 49
## iso48 Ripley isotropic correction estimate of L[loc](r) for point 48
## iso47 Ripley isotropic correction estimate of L[loc](r) for point 47
## iso46 Ripley isotropic correction estimate of L[loc](r) for point 46
## iso45 Ripley isotropic correction estimate of L[loc](r) for point 45
## iso44 Ripley isotropic correction estimate of L[loc](r) for point 44
## iso43 Ripley isotropic correction estimate of L[loc](r) for point 43
## iso42 Ripley isotropic correction estimate of L[loc](r) for point 42
## iso41 Ripley isotropic correction estimate of L[loc](r) for point 41
## iso40 Ripley isotropic correction estimate of L[loc](r) for point 40
## iso39 Ripley isotropic correction estimate of L[loc](r) for point 39
## iso38 Ripley isotropic correction estimate of L[loc](r) for point 38
## iso37 Ripley isotropic correction estimate of L[loc](r) for point 37
## iso36 Ripley isotropic correction estimate of L[loc](r) for point 36
## iso35 Ripley isotropic correction estimate of L[loc](r) for point 35
## iso34 Ripley isotropic correction estimate of L[loc](r) for point 34
## iso33 Ripley isotropic correction estimate of L[loc](r) for point 33
## iso32 Ripley isotropic correction estimate of L[loc](r) for point 32
## iso31 Ripley isotropic correction estimate of L[loc](r) for point 31
## iso30 Ripley isotropic correction estimate of L[loc](r) for point 30
## iso29 Ripley isotropic correction estimate of L[loc](r) for point 29
## iso28 Ripley isotropic correction estimate of L[loc](r) for point 28
## iso27 Ripley isotropic correction estimate of L[loc](r) for point 27
## iso26 Ripley isotropic correction estimate of L[loc](r) for point 26
## iso25 Ripley isotropic correction estimate of L[loc](r) for point 25
## iso24 Ripley isotropic correction estimate of L[loc](r) for point 24
## iso23 Ripley isotropic correction estimate of L[loc](r) for point 23
## iso22 Ripley isotropic correction estimate of L[loc](r) for point 22
## iso21 Ripley isotropic correction estimate of L[loc](r) for point 21
## iso20 Ripley isotropic correction estimate of L[loc](r) for point 20
## iso19 Ripley isotropic correction estimate of L[loc](r) for point 19
## iso18 Ripley isotropic correction estimate of L[loc](r) for point 18
## iso17 Ripley isotropic correction estimate of L[loc](r) for point 17
## iso16 Ripley isotropic correction estimate of L[loc](r) for point 16
## iso15 Ripley isotropic correction estimate of L[loc](r) for point 15
## iso14 Ripley isotropic correction estimate of L[loc](r) for point 14
## iso13 Ripley isotropic correction estimate of L[loc](r) for point 13
## iso12 Ripley isotropic correction estimate of L[loc](r) for point 12
## iso11 Ripley isotropic correction estimate of L[loc](r) for point 11
## iso10 Ripley isotropic correction estimate of L[loc](r) for point 10
## iso09 Ripley isotropic correction estimate of L[loc](r) for point 09
## iso08 Ripley isotropic correction estimate of L[loc](r) for point 08
## iso07 Ripley isotropic correction estimate of L[loc](r) for point 07
## iso06 Ripley isotropic correction estimate of L[loc](r) for point 06
## iso05 Ripley isotropic correction estimate of L[loc](r) for point 05
## iso04 Ripley isotropic correction estimate of L[loc](r) for point 04
## iso03 Ripley isotropic correction estimate of L[loc](r) for point 03
## iso02 Ripley isotropic correction estimate of L[loc](r) for point 02
## iso01 Ripley isotropic correction estimate of L[loc](r) for point 01

The Cross K-function:

Suppose now that we have two (or more) types of events, and the goal is not only to characterize the individual point patterns of the two types, but also to discern the level to which the two point patterns are related to one another. Specifically, it is of interest to know whether the two point patterns are: - independent of one another, - exhibit attraction (the two types of events tend to occur close together), - exhibit repulsion (the two types of events tend to not occur close together).

In the same sense that CSR provided the base case against which clustering or regularity was compared for a univariate point process, independence provides the base case against which attraction or repulsion is compared for a bivariate (or multivariate) point process.

# Cross K function
library("maptools")
asthma.ppp<-as(asthma.pts,"ppp")
asthma.owin<-as(asthma.bdry,"owin")
head(asthma.pts@data)
##   Nsmokers  Asthma HayFever Age Gender roaddist2 d2source1 d2source2
## 1        2 control        0 8.9      1  886552.2  0.000145  0.004850
## 2        2 control        0 8.5      2  889337.6  0.000049  0.006290
## 3        2 control        1 8.5      2  827783.4  0.000050  0.006445
## 4        1 control        0 8.9      2  889337.6  0.000049  0.006290
## 5        1 control        0 6.9      1  889337.6  0.000049  0.006290
## 6        1    case        0 7.1      2 1213531.4  0.000082  0.006373
##   d2source3
## 1  0.021265
## 2  0.023725
## 3  0.024008
## 4  0.023725
## 5  0.023725
## 6  0.023764
m<-asthma.pts$Asthma
asthppp<-ppp(asthma.ppp$x,asthma.ppp$y,marks=m,window=asthma.owin)
## Warning in ppp(asthma.ppp$x, asthma.ppp$y, marks = m, window =
## asthma.owin): data contain duplicated points
asthkcross<-Kcross(asthppp, "case", "control",r=NULL)
asthlcross<-Lcross(asthppp, "case", "control",r=NULL)
plot(asthkcross); plot(asthlcross)

##        lty col    key                                           label
## iso      1   1    iso   italic({hat(K)[list(case,control)]^{iso}}(r))
## trans    2   2  trans italic({hat(K)[list(case,control)]^{trans}}(r))
## border   3   3 border  italic({hat(K)[list(case,control)]^{bord}}(r))
## theo     4   4   theo       italic({K[list(case,control)]^{pois}}(r))
##                                                                    meaning
## iso    Ripley isotropic correction estimate of Kcross["case","control"](r)
## trans        translation-corrected estimate of Kcross["case","control"](r)
## border            border-corrected estimate of Kcross["case","control"](r)
## theo                       theoretical Poisson Kcross["case","control"](r)

##        lty col    key                                           label
## iso      1   1    iso   italic({hat(L)[list(case,control)]^{iso}}(r))
## trans    2   2  trans italic({hat(L)[list(case,control)]^{trans}}(r))
## border   3   3 border  italic({hat(L)[list(case,control)]^{bord}}(r))
## theo     4   4   theo       italic({L[list(case,control)]^{pois}}(r))
##                                                               meaning
## iso    Ripley isotropic correction estimate of L["case","control"](r)
## trans        translation-corrected estimate of L["case","control"](r)
## border            border-corrected estimate of L["case","control"](r)
## theo                       theoretical Poisson L["case","control"](r)
# Envelope for Cross K function
r=seq(0,0.4,by=0.01)
akmult<-envelope(asthppp, Kcross, i="case", j="control", r=r, simulate=expression(rlabel(asthppp)))
## Generating 99 simulations by evaluating expression  ...
## 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, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(akmult, main="Cross-K Function")

##       lty col   key                                         label
## obs     1   1   obs italic({hat(K)[list(case,control)]^{obs}}(r))
## mmean   2   2 mmean       italic({bar(K)[list(case,control)]}(r))
## hi      1   8    hi  italic({hat(K)[list(case,control)]^{hi}}(r))
## lo      1   8    lo  italic({hat(K)[list(case,control)]^{lo}}(r))
##                                                                        meaning
## obs             observed value of Kcross["case","control"](r) for data pattern
## mmean              sample mean of Kcross["case","control"](r) from simulations
## hi    upper pointwise envelope of Kcross["case","control"](r) from simulations
## lo    lower pointwise envelope of Kcross["case","control"](r) from simulations

Inhomogeneous K function:

In practice the intensity λ(x) needs to be estimated either parametrically or non-parametrically. The K(s) values higher than πs2 will mean that the point pattern shows more aggregation than that shown by λ(x) and values lower than πs2 reflect more relative homogeneity.

#Inhomogeneous K function
mydense<-density(mypattern2); plot(mydense)

Ki <- Kinhom(mypattern2, mydense)
plot(Ki, main = "Inhomogeneous K function")

##            lty col        key                              label
## iso          1   1        iso   italic({hat(K)[inhom]^{iso}}(r))
## trans        2   2      trans italic({hat(K)[inhom]^{trans}}(r))
## bord.modif   3   3 bord.modif italic({hat(K)[inhom]^{bordm}}(r))
## border       4   4     border  italic({hat(K)[inhom]^{bord}}(r))
## theo         5   5       theo       italic({K[inhom]^{pois}}(r))
##                                                        meaning
## iso        Ripley isotropic correction estimate of K[inhom](r)
## trans           translation-correction estimate of K[inhom](r)
## bord.modif   modified border-corrected estimate of K[inhom](r)
## border                border-corrected estimate of K[inhom](r)
## theo                           theoretical Poisson K[inhom](r)
Kenv2<-envelope(mypattern2, Kinhom, simulate = expression(rpoispp(mydense)),nrank=2, nsim=99) # 95 C.I
## Generating 99 simulations by evaluating expression  ...
## 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, 43, 44, 45,
## 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
## 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,
## 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(Kenv2)

##       lty col   key                            label
## obs     1   1   obs italic({hat(K)[inhom]^{obs}}(r))
## mmean   2   2 mmean       italic({bar(K)[inhom]}(r))
## hi      1   8    hi  italic({hat(K)[inhom]^{hi}}(r))
## lo      1   8    lo  italic({hat(K)[inhom]^{lo}}(r))
##                                                        meaning
## obs             observed value of K[inhom](r) for data pattern
## mmean              sample mean of K[inhom](r) from simulations
## hi    upper pointwise envelope of K[inhom](r) from simulations
## lo    lower pointwise envelope of K[inhom](r) from simulations