# 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後面的資料模型，進而去做預測
每個X值對應出來的是一個分布
右上方表示法的解讀：
常態分佈下(X為某數時，他所有Y的平均值,變異數)
## Geostatistical approach: Random field (隨機場).
地理統計
Z(s):區域化的隨機變數，跟鄰近的數值有關係
s:位置
Z()指定位置的溫度
Raster Data
平均值 + 隨機值 + 雜訊
### Stationary
穩態：每個地方的統計特性都一樣
期望值是0
兩個位置相減的variance是什麼？Varigram(變異元)
差的變異
### Varigram
意義：兩點之間不相似的程度
兩兩差的平方的平均值
### Semi-variogram
除以2：消掉微分後出現的2
### Fitting a Variogram Model
Nugget: 截距的概念，異質特性的差異
Range: 空間自相關的範圍
Sill: 到天花板的高度，相異程度最大可以到哪裡 (Y值: 相異程度)
## Variogram Model
### Parameters
Range:
Sill:
Nugget:
### Models
> Spherical
> Exponential
> Power function
>
#### Spherical model
#### Exponential Model
#### 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')
```
* 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)
```
### 手動計算 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)
```
Sign in with Wallet
Connect another wallet