空间分析与R语言系列之空间扫描统计量
R语言是一门用于统计分析和绘图的语言，在空间数据分析方面有着强大的功能。为了让大家更好地掌握R语言在空间分析方面的应用，在接下来的一段时间内，小编会介绍下几个空间分析经典模型的R语言工具包，供大家选择学习。这次介绍的是空间聚集性探测方法中应用较广的空间扫描统计量。
空间扫描统计量是1997年由哈佛大学的Kulldorff提出的一种聚集性探测检验方法，因此也叫Kulldorffand-Nagarwalla方法，目的是运用一系列扫描圆，在研究区域探测出疾病的空间聚集性。
以圆形空间扫描统计量为例，该方法在开始进行探测时，随机选取研究区域内某一点或小范围中心点，以其为圆心生成一系列扫描圆。
这些扫描圆的半径由0到规定的上限按照一定的步长逐步变化。当扫描圆半径达到规定的上限后，该方法便又以区域内另外一个点为圆心，开始新一轮的圆形扫描。对于每一个扫描圆，利用圆内外病例的实际值和期望值计算一个似然比值。病例概率分布情况不同，所用的似然比求解公式也不同。以泊松分布为例：
式中，c表示扫描圆中实际病例数；μ表示扫描圆中期望病例数；C表示整个研究区的病例数。
整个扫描过程直到遍历完研究区域内所有的点后结束。扫描过程结束后，将所有扫描圆的似然比由大到小排序，选择排在前面的若干个作为聚类的备选区域进入MonteCarlo检验。通过检验的扫描圆便是探测到的疾病聚集高发区。
下面通过一个实例来说明在R语言中如何利用空间扫描统计量实现聚集性探测。
数据采用maptools工具包自带的sids数据，包括美国婴儿猝死综合症病例数据和出生人口数据（Cressie and Chan 1989, and Cressie 1993)。所用R包为RCluster，要求的数据格式包括婴儿猝死综合症病例、预期病例数（各地区风险人口数*卡罗莱纳州病例总数/卡罗莱纳州总风险人口数）、风险人口数和坐标信息。
数据获取代码：
nc_file<- system.file("shapes/sids.shp",package ="maptools")[1] #read ShapeFile from local
nc_file_CRS<- CRS("+proj=longlat +datum=NAD27") #为文件设置投影
#读取多边形文件并设置投影
nc<- readShapePoly(nc_file,ID="FIPSNO",proj4string=nc_file_CRS)
###将实验所需数据存储在dataframe中，
#Observed(number of cases)病例数
sids<- data.frame(Observed=nc$SID74)
#Expected(number of expected cases)预期病例数
sids<- cbind(sids, Expected=nc$BIR74*sum(nc$SID74)/sum(nc$BIR74))
#经纬度坐标（每个地区中心点坐标在nc中不好获取，借助于spdep包中的nc.sids数据，nc和nc.sids通过ID连接）
data("nc.sids")
for(iin 1:100){
sids$x[i] <-nc.sids$lon[which(nc.sids$CNTY.ID==nc@data$CNTY_ID[i])]
sids$y[i] <-nc.sids$lat[which(nc.sids$CNTY.ID==nc@data$CNTY_ID[i])]
}
然后来看一下卡罗莱纳州婴儿猝死综合症病例分布，这样对卡罗莱纳州婴儿猝死综合症高发区域分布有一个整体的认识。
病例分布展示代码：
#设置颜色分级和相应图例属性
brks<- c(0,2,5,10,20,50) # 0 >=sids$Observed <= 44
grps<- as.ordered(cut(sids$Observed,brks,include.lowest = TRUE))
cols<- brewer.pal(5,"Reds")
#展示
plot(nc,col = cols[unclass(grps)],main = "1974年美国北卡罗来纳州婴儿猝死综合症病例分布", axes = TRUE)
#添加图例
legend(title= "病例数目（人）","bottomleft",legend=levels(grps),fill=cols,bty="n",cex=0.7,
y.intersp=0.6)
从图中可以看到卡罗莱纳州婴儿猝死综合症病例多分布于州东北部、西南部、中部偏北和中部偏东地区。
再来看一下卡罗莱纳州婴儿猝死综合症标准化死亡比（SMR）的分布。SMR是实际死亡人数与该地区（或人群）预期死亡人数之比，可以合理比较不同人群之间同一种疾病的死亡率水平，消除因不同人群内部年龄和性别构成差异所带来的影响。
聚集性探测代码：
#计算抽样过程中需要用到的的参数
#opgam是聚集探测的主函数，模型是poison分布，显著性水平为0.05（即只计算P值小于0.05的）
knresults<- opgam(data=sids, thegrid=sids[,c("x","y")],alpha=0.05,
iscluster=kn.iscluster,fractpop=0.15, R=100, model="poisson", mle=mle)
得到结果后，将似然比（likelihood ratio，LR）结果成图，LR是反映真实性的一个指标。
从LR值分布来看，明显有两个聚集区，分别位于卡罗莱纳州东北部和西南部。
##根据聚集性探测的结果取LR值最大的聚集区展示，
#find cluster with the highest likelihood ratio test找到LR最高的聚集区
clusters <- get.knclusters(sids, knresults)
idx <- which.max(knresults$statistic) #聚集区代码是32
#clusters[[idx]] 对应的dataframe中78、9、24、47四个ID所代表的地区组成第一聚集区
#展示聚集区
从图中可以看出，婴儿猝死综合症的热点区域位于卡罗莱纳州西南部，这个区域的LR值为970062.5（P值为0.0099）。