# Spatial Statistics

## Intro

• Its not always on a map
• Lots of packages
• Lots of techniques
• Lots of subtleties

## Point Patterns

### When the coordinate locations are the data.

• Q: Are these points just 'random'?
• Q: Is that a cluster?
• Q: Are these points like those points?
• Q: Do we get those points where its like this or that?

## Geostatistics

### When measurements at locations are the data.

• Q: Where's the gold?
• Q: How much oil is left?
• Q: Who has all the rare earth metals?
• Q: How do I draw a pollution map?

## Areal Spatial Regression

### Like linear regression models, on areas

• Q: How are these measures related?
• Q: Have we missed something in our model?

## All the same...

• A Point Pattern is a Poisson process from a Gaussian field.
• Geostatistical data are samples from a Gaussian field.
• Areal measurements are integrals over a Gaussian field...

# Point Pattern Analysis

## The `spatstat` package

### Point-pattern analysis package

• Doesn't use `sp` classes
• Point randomness tests
• Point clustering tests
• Point pattern simulation
• Point model fitting

## `spatstat` examples

```data(humberside)
plot(humberside)
```
```##    case control
##       1       2
```
```qH <- quadratcount(humberside, 8, 8)
plot(qH, add = TRUE, col = "blue", cex = 1.5, lwd = 2)
```

### Point pattern simulation and test

```pPois = rpoispp(200)
plot(pPois)
```
```pSSI <- rSSI(0.05, 200)
plot(pSSI)
```
```plot(Kest(pPois, correction = "best"), . - theo ~ r)
```
```##      lty col  key                       label
## iso    1   1  iso hat(K)[iso](r) - K[pois](r)
## theo   2   2 theo     K[pois](r) - K[pois](r)
##                                           meaning
## iso  Ripley isotropic correction estimate of K(r)
## theo                     theoretical Poisson K(r)
```
```plot(Kest(pSSI, correction = "best"), . - theo ~ r)
```
```##      lty col  key                       label
## iso    1   1  iso hat(K)[iso](r) - K[pois](r)
## theo   2   2 theo     K[pois](r) - K[pois](r)
##                                           meaning
## iso  Ripley isotropic correction estimate of K(r)
## theo                     theoretical Poisson K(r)
```

## The `splancs` package

• Very old. Almost prehistoric
• Mostly superceded by `spatstat` and other packages
• One or two things not implemented elsewhere

## Space-time packages

• spacetime - data structures for space-time data
• stpp - space-time point patterns

# GeoStatistics

## Automatic interpolation

```require(automap)
data(meuse)
coordinates(meuse) = ~x + y
data(meuse.grid)
gridded(meuse.grid) = ~x + y
kriging_result = autoKrige(zinc ~ soil + ffreq + dist, meuse, meuse.grid)
```
```## [using universal kriging]
```
```plot(kriging_result)
```

## Kriging with `geoR`

`geoR` uses its own class of objects for its data

```require(geoR)
plot(s100)
```
```grid = expand.grid(seq(0, 1, l = 101), seq(0, 1, l = 101))
vario100 <- variog(s100, max.dist = 1)
```
```## variog: computing omnidirectional variogram
```
```ini.vals <- expand.grid(seq(0, 1, l = 5), seq(0, 1, l = 5))
vfit <- variofit(vario100, ini = ini.vals, fix.nug = TRUE, wei = "equal")
```
```## variofit: covariance model used is matern
## variofit: weights used: equal
## variofit: minimisation function used: optim
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "1"     "0.25" "0"   "0.5"
## status        "est"   "est"  "fix" "fix"
## loss value: 0.167887026209413
```
```plot(vario100)
lines(vfit)
```
```kc <- krige.conv(s100, loc = grid, krige = krige.control(cov.pars = vfit\$cov.pars))
```
```## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
```
```image(kc)
```

## Kriging with `gstat`

`gstat` was written by one of the `sp` developers so it uses `sp` objects everywhere.

```data(meuse)
coordinates(meuse) = ~x + y
data(meuse.grid)
gridded(meuse.grid) = ~x + y
m <- vgm(0.59, "Sph", 874, 0.04)
# ordinary kriging:
x <- krige(log(zinc) ~ 1, meuse, meuse.grid, model = m)
```
```## [using ordinary kriging]
```
```spplot(x["var1.pred"], main = "ordinary kriging predictions")
```
```spplot(x["var1.var"], main = "ordinary kriging variance")
```

# Spatial Regression

## Packages

• spdep - spatial dependency
• spgwr - GWR methodology

## `spdep` functions

### Global Moran Test

```require(spdep)
require(rgdal)
```
```## OGR data source with driver: ESRI Shapefile
## Source: "./Data/", layer: "NY8_utm18"
## with 281 features and 17 fields
## Feature type: wkbPolygon with 2 dimensions
```
```nb = poly2nb(NY8)
moran.test(NY8\$Cases, listw = nb2listw(nb))
```
```##
## 	Moran's I test under randomisation
##
## data:  NY8\$Cases
## weights: nb2listw(nb)
##
## Moran I statistic standard deviate = 4.393, p-value = 5.594e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance
##          0.156020         -0.003571          0.001320
```
```moran.plot(NY8\$Cases, listw = nb2listw(nb, style = "C"))
```

### Local Moran Statistic

```localM = localmoran(NY8\$Cases, listw = nb2listw(nb))
NY8\$moranZ = localM[, "Z.Ii"]
spplot(NY8, "moranZ")
```
```as.character(NY8[NY8\$moranZ > 8, ]\$AREANAME)
```
```## [1] "Vestal town" "Auburn city" "Auburn city"
```

## `spautolm`

```m <- glm(PROPCAS ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8)
sm <- spautolm(PROPCAS ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8, listw = nb2listw(nb),
family = "SAR", method = "full")
summary(m)\$coef
```
```##               Estimate Std. Error t value  Pr(>|t|)
## (Intercept)  4.129e-04  1.524e-04   2.709 0.0071671
## PEXPOSURE    7.112e-05  3.371e-05   2.110 0.0357572
## PCTAGE65P    2.127e-03  5.821e-04   3.654 0.0003084
## PCTOWNHOME  -4.062e-04  1.637e-04  -2.481 0.0136952
```
```summary(sm)\$Coef
```
```##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  4.036e-04  0.0001371   2.944 3.241e-03
## PEXPOSURE    6.632e-05  0.0000285   2.327 1.997e-02
## PCTAGE65P    2.445e-03  0.0005466   4.473 7.712e-06
## PCTOWNHOME  -4.500e-04  0.0001454  -3.096 1.962e-03
```

# Geographically Weighted Regression

## `spgwr`

```require(spgwr)
```
```## Loading required package: spgwr
```
```## NOTE: This package does not constitute approval of GWR as a method of
## spatial analysis
```