在上一节中使用改进的熵权 TOPSIS 法对长三角各地级市的创新能力进行评价。本节使用 GeoDa 将地理数据与综合评价结果整合，再使用 R 语言探索各地级市创新能力的空间相关性。

空间统计分析

根据 Tobler 的“地理学第一定律”：所有事物都与其他事物相关联，较近的事物比较远的事物关联更密切。

创新能力的空间分布

Code

library(rgdal)
library(spdep)
library(ggplot2)
library(MetBrewer)
library(ggvoronoi)

Code

## 导入shp
changsanjiao <- st_read("./data/changsanjiao.shp")
## Reading layer `changsanjiao' from data source 
##   `F:\PublishMyBlog\Rui-s-Blog\posts\R\基于综合评价结果的空间相关性分析\data\changsanjiao.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 25 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 115.755 ymin: 27.93153 xmax: 122.9424 ymax: 34.47572
## Geodetic CRS:  WGS 84

Code

# topsis的空间分布
ggplot(changsanjiao) +
  geom_sf(aes(fill = topsis)) +
  scale_fill_gradientn(colors = rev(met.brewer("Cassatt1", 25, "continuous"))) + 
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    legend.title = element_text(size = 20), 
    legend.text = element_text(size = 10)
  ) +
  theme_bw()

设置空间权重矩阵

使用空间权重矩阵来表达空间依赖的结构，即使用变量的空间滞后来表达空间的依赖性，而空间依赖性通常可以表达成空间近邻的拓扑结构。

空间权重矩阵分为三类：基于几何的、基于理论的以及基于数据的。本篇使用基于几何的，以 Queen 邻接规则判断空间单元的相邻关系。最终得到形式如下的空间权重矩阵\(\mathbf{W} = (w_{ij})_{n \times n}\)：

\[ w_{ij}= \left\{\begin{array}{l} 1, \quad 若i与j相邻, \\ 0, \quad 若i与j不相邻. \end{array}\right. \]

在 GeoDa 软件中初步设置邻接关系，然后导入到 R 中：

Code

# 导入空间邻接关系
nb.city <- read.gal("./data/changsanjiao.gal", override.id = TRUE)
summary(nb.city)
## Neighbour list object:
## Number of regions: 25 
## Number of nonzero links: 96 
## Percentage nonzero weights: 15.36 
## Average number of links: 3.84 
## 1 region with no links:
## 25
## Link number distribution:
## 
## 0 1 2 3 4 5 6 7 8 
## 1 1 2 8 4 6 1 1 1 
## 1 least connected region:
## 4 with 1 link
## 1 most connected region:
## 8 with 8 links

从输出结果中可以看出：第 25 个城市舟山没有邻接城市，即为孤岛。需要手动指定邻接城市，这里指定离舟山市最近的宁波市为其邻接：

Code

code <- attr(nb.city, "region.id") # 获得各个城市代码
which(code == "25") # 根据代码找到舟山市在数据中的的位置
## [1] 25
nb.city[25] # 查看舟山市的邻接情况
## [[1]]
## [1] 0
which(code == "22") # 根据代码找到宁波市在数据中的的位置
## [1] 22
nb.city[[25]] <- 22L # 指定宁波市为舟山市的邻接
nb.city[25] # 查看舟山市的邻接情况，有一个邻接为宁波
## [[1]]
## [1] 22
summary(nb.city) # 查看所有地市的邻接信息，无孤岛
## Neighbour list object:
## Number of regions: 25 
## Number of nonzero links: 97 
## Percentage nonzero weights: 15.52 
## Average number of links: 3.88 
## Link number distribution:
## 
## 1 2 3 4 5 6 7 8 
## 2 2 8 4 6 1 1 1 
## 2 least connected regions:
## 4 25 with 1 link
## 1 most connected region:
## 8 with 8 links

将修正后的邻接关系转化为空间权重矩阵并进行行规范化处理：

Code

lw.city <- nb2listw(nb.city) # nb2listw函数默认对权重进行行标准化处理
names(lw.city)
## [1] "style"      "neighbours" "weights"
names(attributes(lw.city))
## [1] "names"     "class"     "region.id" "call"      "GeoDa"
summary(unlist(lw.city$weights)) # 查看权重信息
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1250  0.2000  0.2000  0.2577  0.3333  1.0000

空间拓扑关系可视化

邻接关系可视化

Code

changsanjiao_ <- readOGR("./data/changsanjiao.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "F:\PublishMyBlog\Rui-s-Blog\posts\R\基于综合评价结果的空间相关性分析\data\changsanjiao.shp", layer: "changsanjiao"
## with 25 features
## It has 14 fields
## Integer64 fields read as strings:  ID_0 ID_1 ID_2
center <- coordinates(changsanjiao_)

par(pin = c(4, 3))
plot(changsanjiao_)
plot(nb.city, coords = center, lwd = 2, col = "red", add = TRUE)

Vonoroi图

Code

x <- center[ , 1]
y <- center[ , 2]
dist <- sqrt((x - 200) ^ 2 + (y - 200) ^ 2)
dd <- data.frame(x, y, dist = dist)

ggplot(dd, aes(x, y, fill = dist)) +
  geom_voronoi() +
  stat_voronoi(geom = "path") +
  geom_point() +
  scale_fill_gradient(low = "#20B2AA", high = "#9370DB") +
  theme_bw()

Code

# theme(legend.position = "none") # 去掉右侧图例标签

全局莫兰指数

用于测度所有空间单元是否存在空间自相关性：

\[ I = \frac{n}{\sum_i\sum_jw_{ij}} \times \frac{\sum_i\sum_jw_{ij}(x_i - \bar{x})(x_j - \bar{x})}{\sum_i(x_i-\bar{x})^2} \]

Code

moran.test(
  changsanjiao$topsis, lw.city, 
  randomisation = FALSE, 
  alternative = "two.sided"
) # 基于正态假设的双侧检验
## 
##  Moran I test under normality
## 
## data:  changsanjiao$topsis  
## weights: lw.city    
## 
## Moran I statistic standard deviate = 3.0026, p-value = 0.002677
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.37964813       -0.04166667        0.01968846

p 值小于 0.05，拒绝无空间自相关的原假设，说明总体来看长三角各地的创新创业能力存在空间自相关性。

局部莫兰指数

反映某一个空间单元与其邻近的空间单元之间的相关关系：

\[ I_i = \frac{n(x_i-\bar{x})\sum_{j \neq i}w_{ij}(x_j-\bar{x})}{\sum_i(x_i-\bar{x})^2} \]

Code

lisa <- localmoran(changsanjiao$topsis, lw.city, alternative = "two.sided")
head(lisa, 25) # 查看lisa
##              Ii          E.Ii       Var.Ii        Z.Ii Pr(z != E(Ii))
## 1   3.111080887 -8.991861e-02 0.9784410785  3.23607268    0.001211865
## 2   2.865405946 -4.908579e-01 1.3582414212  2.87983500    0.003978833
## 3   0.187593674 -2.309551e-03 0.0175320837  1.43421933    0.151509683
## 4   0.167467386 -2.108194e-02 0.5159372009  0.26249834    0.792937268
## 5  -0.056383639 -1.188720e-02 0.0893709204 -0.14884257    0.881677857
## 6   1.439813060 -2.688539e-02 0.1990630215  3.28734692    0.001011362
## 7   0.637898457 -9.439703e-03 0.0508184507  2.87157722    0.004084290
## 8   0.024041518 -1.455714e-02 0.0311852768  0.21857343    0.826982349
## 9   0.050366643 -5.794981e-04 0.0023921922  1.04162948    0.297583494
## 10 -0.026806170 -7.978047e-02 0.2393985327  0.10826910    0.913782236
## 11 -0.031431488 -9.473829e-04 0.0072015193 -0.35922072    0.719429978
## 12  0.319655779 -4.707627e-02 0.1852917072  0.85196328    0.394234473
## 13  0.021310407 -2.127374e-04 0.0005614544  0.90834006    0.363698584
## 14  0.200796655 -3.458275e-02 0.1379019414  0.63384530    0.526181786
## 15 -0.036653301 -1.079703e-02 0.0812642900 -0.09070184    0.927729508
## 16  0.091075104 -3.634547e-03 0.0149576985  0.77439377    0.438697964
## 17  0.168283210 -5.407767e-03 0.0292311049  1.01590861    0.309672906
## 18  0.026475344 -6.651924e-02 0.2564769801  0.18362569    0.854307109
## 19  0.253674917 -8.095271e-03 0.0331663064  1.43738011    0.150610018
## 20  0.213229812 -2.226578e-03 0.0169036318  1.65717879    0.097483332
## 21 -0.505434546 -2.349077e-02 0.1745355251 -1.15359790    0.248665060
## 22  0.459362565 -6.154510e-02 0.6905764549  0.62683735    0.530765856
## 23  0.331050122 -2.500999e-02 0.1325243980  0.97808231    0.328033618
## 24 -0.007773641 -1.475573e-05 0.0001122702 -0.73226297    0.464008079
## 25 -0.412895459 -4.809105e-03 0.1196494283 -1.17976840    0.238092341
sgnf <- which(lisa[ , "Pr(z != E(Ii))"] < 0.05) # 查看局部空间相关性显著的区域
sgnf
## 1 2 6 7 
## 1 2 6 7

Code

changsanjiao$pr <- cut(lisa[ , "Pr(z != E(Ii))"], breaks = c(0, 0.01, 0.05, 0.1, 1))

ggplot(changsanjiao) +
  geom_sf(aes(fill = pr)) +
  scale_fill_discrete(type = c("red", "salmon2", "thistle1","gray99")) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    legend.title = element_text(size = 20), 
    legend.text = element_text(size = 10)
  ) +
  theme_bw()

Code

moran.plot(
  changsanjiao$topsis, lw.city, 
  main = "长三角各地市TOPSIS", 
  xlab = "", ylab = ""
)

创新能力在空间上存在溢出效应，且局部出现聚集效应，主要表现为“高-高聚集”和“低-低聚集”。集聚效应较为显著的地级市有：安庆市、池州市、铜陵市和芜湖市，这些城市的创新能力显著地受到其他邻近城市的影响。

结论

使用了改进的熵权 TOPSIS 法克服了传统 TOPSIS 的缺点，对长三角 25 个地级市的创新能力进行了综合评价，排名前 5 的地级市依次为：南京市、杭州市、宁波市、苏州市和无锡市。指标的相对重要程度排名依次为：新建企业进入、商标授权数目、外观专利公开数目、外来投资笔数、vcpe 投资数目、发明专利授权数目和实用新型专利公开数目。通过空间探索性分析发现：创新能力在空间上存在溢出效应，且局部出现聚集效应，主要表现为“高-高聚集”和“低-低聚集”。