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

2023-10-30 21:32:16

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

一、本文主要内容

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

二、Moran’s I简介

嗯……其实有很多资料已经写得很很很清楚了，这里简单介绍一下。
莫兰指数一般是用来度量空间相关性的一个重要指标。

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是区域数量。

五、数据预处理

先将数据进行处理，数据处理的方式有很多种，笔者用了归一化和倾斜纠正两种方式。
归一化的公式如下：

将数据归一化后，笔者画了一个图，发现数据偏斜分布比较明显，于是对数据开平方，做了数据变换，目的是将不具有正态分布的数据变换成具有正态分布的数据，

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

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

六、构建自定义相邻矩阵

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

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

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

第一行红色框内，第一个数表示最小邻居，0就是说有的区域没有邻居，第二个数表示一共有多少研究对象，第三个是shapefile的文件名，第四个是研究区代号表头id

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

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

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

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

排序好的txt文件：

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，大功告成~

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

码农公寓