---
# System prepended metadata

title: 'Week 12 Measuring spatial dependency: Semi-variogram analysis'
tags: [GlobalAnalysisMethods]

---

# Week 12 Measuring spatial dependency: Semi-variogram analysis
###### tags: `GlobalAnalysisMethods`
## Review:

### Moran's I

利用統計指標的概念量測是否有空間自相關的存在，並利用假設檢定是否顯著
利用自己與鄰居的共變異數，去了解空間自相關的情況

### Moran Correlogram

X軸：定義的空間範圍
Y軸：空間自相關的程度(moran's I)

### Simple Linear Regression Model

簡單線性回歸
了解data後面的資料模型，進而去做預測
![](https://i.imgur.com/aSozJeW.png)

每個X值對應出來的是一個分布
右上方表示法的解讀：
常態分佈下(X為某數時，他所有Y的平均值,變異數)

## Geostatistical approach: Random field (隨機場).

地理統計
Z(s):區域化的隨機變數，跟鄰近的數值有關係
s:位置
Z()指定位置的溫度
Raster Data
![](https://i.imgur.com/lsxdhDR.png)
平均值 + 隨機值 + 雜訊

### Stationary

穩態：每個地方的統計特性都一樣
![](https://i.imgur.com/iSnOqik.png)
期望值是0
兩個位置相減的variance是什麼？Varigram(變異元)
差的變異

### Varigram

意義：兩點之間不相似的程度

兩兩差的平方的平均值
![](https://i.imgur.com/xXbr3Gt.png)

### Semi-variogram

![](https://i.imgur.com/kpil1jZ.png)

除以2：消掉微分後出現的2

### Fitting a Variogram Model

![](https://i.imgur.com/C0yoEAP.png)

Nugget: 截距的概念，異質特性的差異
Range: 空間自相關的範圍
Sill: 到天花板的高度，相異程度最大可以到哪裡 (Y值: 相異程度)

## Variogram Model

![](https://i.imgur.com/PBFMQr9.png)

### Parameters

Range:
Sill:
Nugget:

### Models

> Spherical 
> Exponential
> Power function
> 
#### Spherical model

![](https://i.imgur.com/SKFNnrh.png)

#### Exponential Model

![](https://i.imgur.com/V1g1Rix.png)
![](https://i.imgur.com/pYwUdTs.png)

#### Wave(Hole-Effect) Model

週期性模型
Peaks:
Trough:

### Steps of Variogram Modeling

1. Sampling locationss and measured variable
2. Variogram cloud showing semivariances for all pairs
3. Semivariance aggregated to lags of about 100m
4. The final variogram model fitting

## 應用實例介紹

中國縣域鄉村地域多功能格局及影響因素識別

* 農村地域多功能指數

## R-code

> Package: gstat
> 
### Variogram cloud

* Exploring distance & variance

```r=

x= coordinates(EPA_STN)[,1]
y= coordinates(EPA_STN)[,2]

STNDF = cbind(x,y)
dis_STN= dist(STNDF)

pm= EPA_STN@data[,12]

PMDF= cbind(pm,pm) # 透過 column 合併
dis_PM = dist(PMDF)

plot(dis_PM~sqrt(dis_STN))
abline(lm(dis_PM~sqrt(dis_STN)), lwd=3, col='red')
```
![](https://i.imgur.com/rGYsvEm.png)

* Semivariogram

```r=
pm.vgm = variogram(PM~1, EPA_STN,cutoff=80000, width=2000)
pm.fit = fit.variogram(pm.vgm, model = vgm(500, "Exp",30000,5) )
plot(pm.vgm,pm.fit)
```

![](https://i.imgur.com/j0SxEb0.png)


###  手動計算 moran's I

```r=
# 相鄰矩陣
tp.nb.m = nb2mat(tp.nb)
xx = x - mean(x)
sum(tp.nb.m*(xx%*%t(xx)))/sum(xx^2)
```

### 手動計算 local moran's I

```r=
LISA = localmoran(x, tp.nb.w)
LISA[1]
z = (x-mean(x)) / (sd(x)*sqrt(11/12))
z[1]*sum(tp.nb.m[1,]*z)
```

```r=
LISA = localmoran(x,tp.nb.w, mlvar = F)
LISA[1]
z = (x-mean(x)) / sd(x)
z[1]*sum(tp.nb.m[1,]*z)
```

* 用矩陣方法一次求得所有I_i
```r=
z*(tp.nb.m %*% z)
```

### 手動計算 General G

```r=
G = globalG.test(x, tp.nb.w)
G$estimate[1]

G.num = sum(tp.nb.m*(x %*% t(x)))
G.den = sum(x %*% t(x)) - sum(diag(x %*% t(x)))
G.num / G.den
```

### 手動計算 Gi*

```r=
Gi* = localG(x, tp.nb.w.in, return_internals = T)
attr(Gi*, "internals")[,1]

tp.nb.m.in %*% x/sum(x)
```