在GeoDa中实现自定义相邻矩阵的Moran‘s I计算

末蓝、 2022-11-19 05:27 63阅读 0赞

## 在GeoDa中实现自定义相邻矩阵的Moran’s I计算 ##

## 一、本文主要内容 ##

本文将详细介绍如何在GeoDa中用自己定义的相邻矩阵计算Moran’s I，笔者定义的相邻是按照特征值从大到小进行了排序，排序后与x相邻的前一个和后一个被视为x的邻居（因此第一个和最后一个只有一个邻居）。给出GeoDa下载地址，和数据处理过程中用到的.py文件，希望对各位读者有所帮助，如有不妥之处请批评指正~

## 二、Moran’s I简介 ##

嗯……其实有很多资料已经写得很很很清楚了，这里简单介绍一下。  
莫兰指数一般是用来度量空间相关性的一个重要指标。  
![在这里插入图片描述][20210415091326728.png_pic_center]  
i和j表示区域的索引，如何把两个区域关联起来呢？那就是通过Wij这个空间矩阵，它定义了研究区域的相邻关系，n是总数，Yi和Yj是因变量，Y拔就是均值了。

也就是说，自变量并没有参与计算，只用来生成了相邻矩阵。

Moran’s I >0表示空间正相关性，其值越大，空间相关性越明显，Moran’s I <0表示空间负相关性，其值越小，空间差异越大，否则，Moran’s I = 0，空间呈随机性。

## 三、GeoDa简介及下载地址 ##

Geoda是一款用于空间数据分析的免费软件，最初是由Luc Anselin及其UIUC小组开发的。 Anselin搬到ASU后，由ASU的Geoda中心进行维护。 欲了解更多有关信息和文档，请参见下面的网站：http：//geodacenter.asu.edu/。 下载网址：https：//spatial.uchicago.edu/software。

## 四、数据准备 ##

新冠疫情无人不知无人不晓，笔者利用Moran’s I探索了社会人口特征数据与疫情之间的关系，即疫情和谁相关？相关程度有多少？用到的数据如下：

*  研究地区的shapefile地图
 *  社会人口特征数据csv

其中，每个研究区有唯一代号KEY，csv中对应记录该地区的一系列特征。  
对应到公式里，社会人口特征数据用来生成空间相邻矩阵，Yi，Yj,Y拔带入疫情数据,n是区域数量。

## 五、数据预处理 ##

先将数据进行处理，数据处理的方式有很多种，笔者用了归一化和倾斜纠正两种方式。  
归一化的公式如下：  
![在这里插入图片描述][20210415094658477.png_pic_center]  
将数据归一化后，笔者画了一个图，发现数据偏斜分布比较明显，于是对数据开平方，做了数据变换，目的是将不具有正态分布的数据变换成具有正态分布的数据，

以上处理excel就可以做，熟悉R的读者可以考虑用更高效的R。

读者可以根据需要选用不同的处理方式，没有思路的话可以看看相关主题的论文是如何处理的。

## 六、构建自定义相邻矩阵 ##

好啦，数据处理好可以开始制作自己的相邻矩阵啦。在经典的莫兰指数计算当中，有几种常用的“相邻”，我们首先看一下GeoDa中封装好的相邻矩阵是什么样子：  
打开GeoDa-----Tools----Weights Manager, 需要先打开shapefile,否则此处是灰的。

打开之后点击Create,可以看到有分了两种矩阵，一个是Contiguity，一个是Distance，每个矩阵的意思此处不展开讲啦，毕竟我们是来自定义的；随意点击一个然后点击Create，生成一个.gal文件，有时候可能会提示有的区域没有邻居等，不用理会，然后会提示创建成功。

这个就是空间矩阵，我们用笔记本打开：

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyOTI1MDMx_size_16_color_FFFFFF_t_70_pic_center]  
第一行红色框内，第一个数表示最小邻居，0就是说有的区域没有邻居，第二个数表示一共有多少研究对象，第三个是shapefile的文件名，第四个是研究区代号表头id

下面的就是空间矩阵了，以黄色框为例，112是我的区域代号，4表示这个区域有4个邻居，具体是谁就写在了下面一行113,114,115,141，每个区域都会如此记录一遍，没有邻居的会写0，然后空一行。

好了，接下来要做的就是照猫画虎，前面说过，笔者定义的相邻是按照特征值从小到大排序，排序后与x相邻的前一个和后一个被视为x的邻居，第一个区域和最后一个区域就只有一个邻居。

因此笔者根据特征值大小，把对应的研究区代号在excel中排序，注意，排序是根据想要研究的特征值的大小，而空间矩阵中需要记录的是研究区域的代号，不要搞错了哦。

因为笔者的空间矩阵除了第一个和最后一个元素邻居是1以外，其他的都是2，因此笔者的男朋友帮忙写了一个python文件，方便批量处理，\*\*但是这个py文件写的并不完全，还需要手动加工一下，\*\*这里分享给大家：

排序好的txt文件：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyOTI1MDMx_size_16_color_FFFFFF_t_70_pic_center 1]  
PY代码：

dataFile = open("denf2t6.txt", "r").read()
    dataFile = dataFile.split("\t")
    
    for eachIndex in range(len(dataFile)):
        dataFile[eachIndex] = dataFile[eachIndex]
    
    if dataFile[-1] == [""]: #如果最后一个是空
        dataFile.pop()
    
    filename='geodenf2t6.txt'
    # 新生成的空间矩阵名字，后面需要手动把后缀改成.gal
    len = len(dataFile)
    with open(filename,'w') as f:
        f.write('1 214 dissolve TPU\n') # 矩阵第一行的属性信息 因为笔者的是固定的，所以写死了
        for i,v in enumerate(dataFile):
            # 这里的加了一个限制条件，因为第一个没有前驱，最后一个没有后驱
            if(i > 0 and i < len-1):
                num=2 #两个邻居
                text1 = dataFile[i-1] #前驱
                text2 = dataFile[i+1] #后驱
                # 模板字符串这里改一下，不能直接用+拼接
                text = Template("${v}\t${num}\n${text1}\t${text2}\n") #连起来
                str = text.substitute(v=v,num=num,text1=text1,text2=text2)
                f.write(str)

好了，到这里自定义空间矩阵已经完成！

## 七、在GeoDa中计算Moran’s I ##

回到GeoDa中，TOOLS----WEIGHTS MANAGER-----LOAD，打开创建的.gal文件，如果正确的话下面会显示对应的属性信息。

SPACE----按所需选择计算哪种Moran’s I，大功告成~

后续结果的解读挖个坑吧，祝各位好运！

[20210415091326728.png_pic_center]: /images/20221022/949d8a8390f74a43a8cd0326b729570b.png
[20210415094658477.png_pic_center]: /images/20221022/c90547fb846949d59669ffde764fa139.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyOTI1MDMx_size_16_color_FFFFFF_t_70_pic_center]: /images/20221022/9c60104e7c8a405eb94b997575223e9b.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzMyOTI1MDMx_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221022/40ccb68a47914f83b6f028d4718f4d71.png