莫兰指数计算

加载分析所需安装包

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.4.1

## Loading required package: sp

## rgdal: version: 1.2-8, (SVN revision 663)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-4

library(spdep)

## Loading required package: Matrix

构建邻接矩阵和案例数值

weight <- matrix(data = c(0,1/3,1/3,1/3,1/2,0,1/2,0,1/3,1/3,0,1/3,1/2,0,1/2,0), 
                   ncol = 4, nrow = 4, byrow = TRUE)
print(weight)

##           [,1]      [,2]      [,3]      [,4]
## [1,] 0.0000000 0.3333333 0.3333333 0.3333333
## [2,] 0.5000000 0.0000000 0.5000000 0.0000000
## [3,] 0.3333333 0.3333333 0.0000000 0.3333333
## [4,] 0.5000000 0.0000000 0.5000000 0.0000000

value <- c(1,2,3,4)
print(value)

## [1] 1 2 3 4

计算案例数值的平均值

value.s <- value - mean(value)
print(value.s)

## [1] -1.5 -0.5  0.5  1.5

计算各个数值与平均值之差

value.m <- value.s %*% t(value.s)
print(value.m)

##       [,1]  [,2]  [,3]  [,4]
## [1,]  2.25  0.75 -0.75 -2.25
## [2,]  0.75  0.25 -0.25 -0.75
## [3,] -0.75 -0.25  0.25  0.75
## [4,] -2.25 -0.75  0.75  2.25

计算空间权重与i数值离差的乘积，并求和 （分子）

value.up <- sum(weight * value.m)
print(value.up)

## [1] -1.333333

计算数值方差与权重之和的乘积 （分母）

value.d <- sum((value.s - mean(value.s))^2)/(length(value.s))
weight.sum <- sum(weight)
value.blw <- value.d * weight.sum
print(value.blw)

## [1] 5

得到莫兰指数

moran.s <- value.up/value.blw
print(moran.s)

## [1] -0.2666667

使用moran.test函数计算莫兰指数

sample <- readOGR('/Users/souseki/Documents/CBNweekly/Project_Myselt/algorithm/algorithmProject/空间自相关分析/样例数据/sample2.shp')

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/souseki/Documents/CBNweekly/Project_Myselt/algorithm/algorithmProject/空间自相关分析/样例数据/sample2.shp", layer: "sample2"
## with 4 features
## It has 7 fields

plot(sample)

weight.s <- poly2nb(sample, queen = T)
weight_mat.s <- nb2listw(weight.s, style="W", zero.policy=TRUE)
moran.test(sample$value, listw = weight_mat.s, zero.policy = T, na.action=na.pass)

## 
##  Moran I test under randomisation
## 
## data:  sample$value  
## weights: weight_mat.s  
## 
## Moran I statistic standard deviate = 0.35729, p-value = 0.3604
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.26666667       -0.33333333        0.03481481