Spatial Point Process 2: Fitting Spatial Models

setwd("D:/_Lab_Data/R/R_Labs")
library(sp);library(rgdal)
library(spatstat) #Objects: PPP
library(splancs)  #Objects: spatial Points

data <- read.table("pts_data/tpedata.csv", header=TRUE, sep=",")
data_bnd <- read.table("pts_data/tpe_sqr_bnd2.csv", header=TRUE, sep=",")

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.

Model-fitting: spatial trend

# ppp(x.coordinates, y.coordinates, window)
mypattern2<-ppp(data[,2], data[,3], window=PTS_bnd2) 
plot(mypattern2)

fit<-ppm(mypattern2, ~x, Poisson())
summary(fit)
## Point process model
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = mypattern2, trend = ~x, interaction = Poisson())
## Edge correction: "border"
##  [border correction distance r = 0 ]
## ---------------------------------------------------------------------------
## Quadrature scheme = data + dummy + weights
## Data pattern:
## Planar point pattern:  63 points
## Average intensity 2.34e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 4149 vertices
## enclosing rectangle: [295264.9, 316374.2] x [2761740.2, 2789392.5] units
## Window area = 269660000 square units
## 
## Dummy quadrature points:
## (32 x 32 grid, plus 4 corner points)
## Planar point pattern:  531 points
## Average intensity 1.97e-06 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 4149 vertices
## enclosing rectangle: [295264.9, 316374.2] x [2761740.2, 2789392.5] units
## Window area = 269660000 square units
## 
## Quadrature weights:
## (counting weights based on 32 x 32 array of rectangular tiles)
## All weights:
##  range: [26400, 570000]  total: 2.69e+08
## Weights on data points:
##  range: [97700, 285000]  total: 15300000
## Weights on dummy points:
##  range: [26400, 570000]  total: 2.54e+08
## ---------------------------------------------------------------------------
## FITTED MODEL:
## 
## Nonstationary Poisson process
## 
## ---- Intensity: ----
## 
## Log intensity: ~x
## 
## Fitted trend coefficients:
##   (Intercept)             x 
##  8.015193e+00 -7.648514e-05 
## 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept)  8.015193e+00 9.929565e+00 -1.144640e+01  2.747678e+01      
## x           -7.648514e-05 3.267593e-05 -1.405288e-04 -1.244151e-05     *
##                   Zval
## (Intercept)  0.8072048
## x           -2.3407185
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
##   (Intercept)             x 
##  8.015193e+00 -7.648514e-05 
## 
## Fitted exp(theta):
##  (Intercept)            x 
## 3026.5931410    0.9999235
par(mar=c(0,0,1,0)) 
diagnose.ppm(fit)

## Model diagnostics (raw residuals)
## Diagnostics available:
##  four-panel plot
##  mark plot 
##  smoothed residual field
##  x cumulative residuals
##  y cumulative residuals
##  sum of all residuals
## sum of raw residuals in entire window = -1.308e-07 
## area of entire window = 269700000 
## quadrature area = 269500000 
## range of smoothed field = [ -2.494e-07,3.384e-07 ]

Mapping the prediction

cif: conditional intensity lambda(u,X)

cif <- predict(fit,type="cif",ngrid=100)
image(cif)

persp(cif, theta=-30,phi=40,d=4,ticktype="detailed",zlab="z")

Simulating Point Process

par(mfrow=c(1,2))
sim1 <- rmh(fit); plot(sim1)
## Extracting model information...Evaluating trend...done.
## Checking arguments..determining simulation windows...Evaluating trend integral...
sim1.dense<-density(sim1); plot(sim1.dense)

Model-fitting: interaction

# str<-rStrauss(beta=10, gamma = 1, R = 0.1, W = PTS_bnd2)
X <- rStrauss(0.05,0.1,0.3,square(60)); plot(X)

fit.tpe<-ppm(mypattern2, ~1, Strauss(r=1000))
summary(fit.tpe);plot(fit.tpe)
## Point process model
## Fitting method: maximum pseudolikelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = mypattern2, trend = ~1, interaction = Strauss(r = 1000))
## Edge correction: "border"
##  [border correction distance r = 1000 ]
## ---------------------------------------------------------------------------
## Quadrature scheme = data + dummy + weights
## Data pattern:
## Planar point pattern:  63 points
## Average intensity 2.34e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 4149 vertices
## enclosing rectangle: [295264.9, 316374.2] x [2761740.2, 2789392.5] units
## Window area = 269660000 square units
## 
## Dummy quadrature points:
## (32 x 32 grid, plus 4 corner points)
## Planar point pattern:  531 points
## Average intensity 1.97e-06 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 4149 vertices
## enclosing rectangle: [295264.9, 316374.2] x [2761740.2, 2789392.5] units
## Window area = 269660000 square units
## 
## Quadrature weights:
## (counting weights based on 32 x 32 array of rectangular tiles)
## All weights:
##  range: [26400, 570000]  total: 2.69e+08
## Weights on data points:
##  range: [97700, 285000]  total: 15300000
## Weights on dummy points:
##  range: [26400, 570000]  total: 2.54e+08
## ---------------------------------------------------------------------------
## FITTED MODEL:
## 
## Stationary Strauss process
## 
## ---- Trend: ----
## 
## 
## First order term:
##         beta 
## 1.513965e-07 
## 
##                Estimate      S.E.     CI95.lo     CI95.hi Ztest       Zval
## (Intercept) -15.7033638 0.2580655 -16.2091629 -15.1975647   *** -60.850304
## Interaction   0.4951264 0.1258418   0.2484811   0.7417718   ***   3.934516
## 
##  ---- Interaction: -----
## 
## Interaction:  Strauss process
## interaction distance:    1000
## Fitted interaction parameter gamma:   1.6407057
## 
## Relevant coefficients:
## Interaction 
##   0.4951264 
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
## (Intercept) Interaction 
## -15.7033638   0.4951264 
## 
## Fitted exp(theta):
##  (Intercept)  Interaction 
## 1.513965e-07 1.640706e+00 
## 
## *** Model is not valid ***
## *** Interaction parameters are outside valid range ***

par(mar=c(0,0,1,0))
diagnose.ppm(fit.tpe); coef(fit.tpe)

## Model diagnostics (raw residuals)
## Diagnostics available:
##  four-panel plot
##  mark plot 
##  smoothed residual field
##  x cumulative residuals
##  y cumulative residuals
##  sum of all residuals
## sum of raw residuals in clipped window = -6.551e-08 
## area of clipped window = 181800000 
## quadrature area = 180300000 
## range of smoothed field = [ -1.497e-07,2.551e-07 ]
## (Intercept) Interaction 
## -15.7033638   0.4951264
cif2 <- predict(fit.tpe,type="cif",ngrid=100)
persp(cif2, theta=-30,phi=40,d=4,ticktype="detailed",zlab="z")

Simulating Point Process

par(mfrow=c(1,2))
sim2 <- rmh(fit.tpe); plot(sim2)
## Model is invalid - projecting it
## Extracting model information...Evaluating trend...done.
## Checking arguments..determining simulation windows...
sim2.dense<-density(sim2); plot(sim2.dense)

Spatial trend + Interaction

fit.tpe2<-ppm(mypattern2, ~x, Strauss(r=1000))
summary(fit.tpe2);plot(fit.tpe2)
## Point process model
## Fitting method: maximum pseudolikelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = mypattern2, trend = ~x, interaction = Strauss(r = 1000))
## Edge correction: "border"
##  [border correction distance r = 1000 ]
## ---------------------------------------------------------------------------
## Quadrature scheme = data + dummy + weights
## Data pattern:
## Planar point pattern:  63 points
## Average intensity 2.34e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 4149 vertices
## enclosing rectangle: [295264.9, 316374.2] x [2761740.2, 2789392.5] units
## Window area = 269660000 square units
## 
## Dummy quadrature points:
## (32 x 32 grid, plus 4 corner points)
## Planar point pattern:  531 points
## Average intensity 1.97e-06 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 4149 vertices
## enclosing rectangle: [295264.9, 316374.2] x [2761740.2, 2789392.5] units
## Window area = 269660000 square units
## 
## Quadrature weights:
## (counting weights based on 32 x 32 array of rectangular tiles)
## All weights:
##  range: [26400, 570000]  total: 2.69e+08
## Weights on data points:
##  range: [97700, 285000]  total: 15300000
## Weights on dummy points:
##  range: [26400, 570000]  total: 2.54e+08
## ---------------------------------------------------------------------------
## FITTED MODEL:
## 
## Nonstationary Strauss process
## 
## ---- Trend: ----
## 
## Log trend: ~x
## 
## Fitted trend coefficients:
##   (Intercept)             x 
## -4.931495e+00 -3.531182e-05 
## 
##                  Estimate         S.E.       CI95.lo      CI95.hi Ztest
## (Intercept) -4.931495e+00 1.413597e+01 -3.263748e+01 2.277449e+01      
## x           -3.531182e-05 4.602019e-05 -1.255097e-04 5.488609e-05      
## Interaction  4.788408e-01 1.424441e-01  1.996555e-01 7.580260e-01   ***
##                   Zval
## (Intercept) -0.3488615
## x           -0.7673115
## Interaction  3.3616057
## 
##  ---- Interaction: -----
## 
## Interaction:  Strauss process
## interaction distance:    1000
## Fitted interaction parameter gamma:   1.6142021
## 
## Relevant coefficients:
## Interaction 
##   0.4788408 
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
##   (Intercept)             x   Interaction 
## -4.931495e+00 -3.531182e-05  4.788408e-01 
## 
## Fitted exp(theta):
## (Intercept)           x Interaction 
##  0.00721571  0.99996469  1.61420208 
## 
## *** Model is not valid ***
## *** Interaction parameters are outside valid range ***

par(mar=c(0,0,1,0))
diagnose.ppm(fit.tpe2); coef(fit.tpe2)

## Model diagnostics (raw residuals)
## Diagnostics available:
##  four-panel plot
##  mark plot 
##  smoothed residual field
##  x cumulative residuals
##  y cumulative residuals
##  sum of all residuals
## sum of raw residuals in clipped window = -8.423e-08 
## area of clipped window = 181800000 
## quadrature area = 180300000 
## range of smoothed field = [ -1.516e-07,2.322e-07 ]
##   (Intercept)             x   Interaction 
## -4.931495e+00 -3.531182e-05  4.788408e-01
cif2 <- predict(fit.tpe2,type="cif",ngrid=100)
persp(cif2, theta=-30,phi=40,d=4,ticktype="detailed",zlab="z")

Modeling the Effects of Covariates

data(bei); plot(bei)

class(bei)
## [1] "ppp"
plot(bei.extra)

Slope <- bei.extra$grad; plot(Slope)

fit1.bei <- ppm(bei, ~elev + grad, covariates = bei.extra)
fit2.bei <- ppm(bei, ~factor(elev>140), covariates = bei.extra)
fit3.bei <- ppm(bei, ~elev + grad, Poisson(), covariates = bei.extra)
fit4.bei <- ppm(bei, ~elev + grad, Strauss(r=1), covariates = bei.extra)
fit5.bei <- ppm(bei, ~polynom(x, 2) + factor(elev>140), Poisson(), covariates = bei.extra)

Model Comparisons and Simulations

# model comparisons
AIC(fit1.bei);AIC(fit2.bei);AIC(fit3.bei);AIC(fit4.bei);AIC(fit5.bei)
## [1] 42296.21
## [1] 42722.64
## [1] 42296.21
## [1] 40719.85
## [1] 42420.78
# model summary
summary(fit4.bei)
## Point process model
## Fitting method: maximum pseudolikelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = bei, trend = ~elev + grad, interaction = Strauss(r = 1), 
##     covariates = bei.extra)
## Edge correction: "border"
##  [border correction distance r = 1 ]
## ---------------------------------------------------------------------------
## Quadrature scheme = data + dummy + weights
## Data pattern:
## Planar point pattern:  3604 points
## Average intensity 0.00721 points per square metre
## Window: rectangle = [0, 1000] x [0, 500] metres
## Window area = 5e+05 square metres
## Unit of length: 1 metre
## 
## Dummy quadrature points:
## (130 x 130 grid, plus 4 corner points)
## Planar point pattern:  16904 points
## Average intensity 0.0338 points per square metre
## Window: rectangle = [0, 1000] x [0, 500] metres
## Window area = 5e+05 square metres
## Unit of length: 1 metre
## 
## Quadrature weights:
## (counting weights based on 130 x 130 array of rectangular tiles)
## All weights:
##  range: [1.64, 29.6] total: 5e+05
## Weights on data points:
##  range: [1.64, 14.8] total: 41000
## Weights on dummy points:
##  range: [1.64, 29.6] total: 459000
## ---------------------------------------------------------------------------
## FITTED MODEL:
## 
## Nonstationary Strauss process
## 
## ---- Trend: ----
## 
## Log trend: ~elev + grad
## Model depends on external covariates 'elev' and 'grad'
## Covariates provided:
##  elev: im
##  grad: im
## 
## Fitted trend coefficients:
## (Intercept)        elev        grad 
## -8.32508801  0.01928238  5.79728448 
## 
##                Estimate        S.E.    CI95.lo     CI95.hi Ztest
## (Intercept) -8.32508801 0.414950865 -9.1383768 -7.51179926   ***
## elev         0.01928238 0.002799172  0.0137961  0.02476866   ***
## grad         5.79728448 0.277383421  5.2536230  6.34094600   ***
## Interaction  1.04627703 0.021631571  1.0038799  1.08867413   ***
##                   Zval
## (Intercept) -20.062828
## elev          6.888601
## grad         20.899895
## Interaction  48.368055
## 
##  ---- Interaction: -----
## 
## Interaction:  Strauss process
## interaction distance:    1
## Fitted interaction parameter gamma:   2.8470319
## 
## Relevant coefficients:
## Interaction 
##    1.046277 
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
## (Intercept)        elev        grad Interaction 
## -8.32508801  0.01928238  5.79728448  1.04627703 
## 
## Fitted exp(theta):
##  (Intercept)         elev         grad  Interaction 
## 2.423596e-04 1.019469e+00 3.294038e+02 2.847032e+00 
## 
## *** Model is not valid ***
## *** Interaction parameters are outside valid range ***
diagnose.ppm(fit4.bei)

## Model diagnostics (raw residuals)
## Diagnostics available:
##  four-panel plot
##  mark plot 
##  smoothed residual field
##  x cumulative residuals
##  y cumulative residuals
##  sum of all residuals
## sum of raw residuals in clipped window = -2.318e-08 
## area of clipped window = 497000 
## quadrature area = 499600 
## range of smoothed field = [ -0.003937,0.006918 ]
coef(fit4.bei)
## (Intercept)        elev        grad Interaction 
## -8.32508801  0.01928238  5.79728448  1.04627703
lam <- predict(fit4.bei, type="cif", ,ngrid=100)
persp(lam, theta=-30,phi=40,d=4,ticktype="detailed",zlab="z")

# Simulation
sim3 <- rmh(fit4.bei); plot(sim3)
## Model is invalid - projecting it
## Extracting model information...Evaluating trend...done.
## Checking arguments..determining simulation windows...Evaluating trend integral...

sim3.dense<-density(sim3); plot(sim3.dense)