# Week 7 Nearest Neighbor Analysis

> Outline:
> Analyizing global patterns
> Nearest Neighbor Analysis(NNA)
> 1. R-scale
> 2. 顯著性檢定
> * 跟理論分布比較
> * 蒙地卡羅顯著性檢定
> 3. K-order NNA
> G(d) function
>

## Analyzing global patterns

* First-order property (intensity 強度)

**某一點**event發生的頻率、平均次數

* Second-order property (dependency 相依性)

paired event 在**兩點之間**的關係、特性

## NNA

計算每個點距離他最近點(nearest point)的平均距離(average distance)

與隨機分布下的平均距離比較
大於隨機分布的平均距離 ---> Dispersed (Uniform)
小於隨機分布的平均距離 ---> Clustered

### R-scale

觀察值與理論隨機值的比值: $R=\frac{r_{obs}}{r_{exp}}$

* Generate random points by an<a style="color:red">
indepentdent random process (IRP), also called
complete spatial randomness (CSR)</a>

* 依照IRP/CSR產生的隨機點符合兩個條件：
1. 每個事件點出現在任何一個位置上都有相同的機率
2. 一個點的出現不會影響另一個點出現的機率


### Significance test of NNA

1. 與理論隨機值比較

$$Z=\frac{r_{obs}-r_{exp}}{s.e.},\ r_{exp}=\frac{0.5}{\sqrt{\frac{n}{A}}}; \ s.e.=\frac{0.26136}{\sqrt{\frac{n^2}{A}}}$$

給定信心水準進行假設檢定，計算觀察值的Z-score

以95%的信心水準進行雙尾假設檢定，
H_0: 點分布型態是隨機分布
H_a: 點分布型態不是隨機分布

當$|Z|>1.96$，則統計上顯著，推翻虛無假設

2. 蒙地卡羅顯著性檢定

模擬隨機分布點，做出信賴包絡(confidence envelopes)，取出前後第5%的值比較

### K-order NNA

若只看第一鄰近的點位，就忽視了現實中的實際分布
為了能保有現實中的分布情形，因此使用K-order取出鄰近次序上的點比較

R 實作： 利用 nndist()

nndist(ppp格式的點資料, k=1)--->回傳第一鄰近的平均距離

## G(d) function

學者認為取完平均之後會淡化點與其他點之間的距離分布，因此設計G(d) function試圖解決

核心概念：距離小於d的點有多少%

因此畫出來的圖表會是CDF

### G function的解讀
1. 若一地的點分布是密集的，那麼在短距離內就會上升很快
2. 若點分布是分散的，那圖形在短距離時會緩慢上升，到最後會快速爬升

R 實作 G(d): ecdf(nnd)

### G(d)的顯著性檢定

1. 模擬出隨機點分布並計算這些隨機點的G function
2. 重複模擬很多次，會得到很多條的G function
3. 利用多次模擬結果的信賴包絡去跟觀察值比較

## 實習

以台南市行政區為邊界檢定台南市學校的空間型態分布

library(sp)
library(rgdal)
library(spatstat)
library(maptools)
library(GISTools)
setwd("C:/Users/user/Desktop/SpatialAnalysis")

School <- readOGR(dsn =  "./data", layer = "Schools", encoding="utf8")
TWCounty <- readOGR(dsn = "./data", layer = "Taiwan_county", encoding="utf8")

學校點座標處理：

bnd = bbox(School)
x.coor = School$coords.x1
y.coor = School$coords.x2
x.range = bnd[1,]
y.range = bnd[2,]

取得台南市邊界：

index = TWCounty$COUNTY=="台南市"
TN_BND = TWCounty[index,]

# 加寬500公尺
TN = gBuffer(TN_BND, width = 500) 
area = poly.areas(TN)

格式轉換：ppp

TN_BND2 = as.owin.SpatialPolygons(TN)

school.ppp = ppp(x.coor, y.coor, window = TN_BND2)

### NNA

nnd = nndist(school.ppp2, k=1)
d1 = mean(nnd);d1 # 觀察值

rd = 0.5/sqrt(424/area)
r.scale = d1/rd

模擬1000次

nn2 = c()
for (i in 1:1000) {
  rp = rpoint(424, win = PTS_bnd)
  nn2[i] = mean(nndist(rp, k=1))
}

出圖

hist(nn2, xlim=c(900,1500), main = "NNA", xlab = "平均距離")
abline(v=d1, lwd=4, col='red')


### K-order

# 觀察值
ANN = apply(nndist(school.ppp2,k =1:50),2, FUN=mean) 

## 模擬1000次
df = data.frame(ANN)
for (i in 1:1000) {
  rp = rpoint(424, win=PTS_bnd)
  df[,i+1] = apply(nndist(rp, k=1:50),2,FUN=mean)
}
df = t(df)

for (i in 1:50) {
  df[,i] = sort(df[,i])
}

plot(df[501,], type="l", xlab="K-order nearest", main="K-order NN", col='blue')
lines(df[51,], col='black')
lines(df[951,], col='black')
lines(ANN,col='red')


### G function

G = ecdf(nnd)
plot(G, main = "G function", col='red', cex=0, verticals=T, xaxs="i", yaxs="i")
for(i in 1:1000){
  nnd.sim = nndist(rpoint(424, win=PTS_bnd), k=1)
  G1 = ecdf(nnd.sim)
  lines(G1, col='blue',cex=0, verticals=T, xaxs="i", yaxs="i")
}