R与SAS软件在绘制统计地图方面的应用研究

天津医科大学公共卫生学院流行病与卫生统计学教研室(300070)

邓晓伟 陈云峰 刘红伟 段同庆 马 骏△

统计地图是用不同的颜色和花纹表示统计数值在地理分布上的变化，适宜描述研究指标的地理分布[1]。尤其是在流行病学研究中，统计地图用来描述疾病在地区、时间的分布再合适不过。

地图的绘制、处理经常出现在地理信息系统(geographical information systems)中。随着学科交融和统计学的发展，地图的绘制，特别是统计地图的绘制已经不仅限于使用固定的地理信息系统专用软件，如ArcGIS等。SAS作为统计分析商用软件有Graph模块可以实现统计制图；R作为开源统计编程语言，Base中有plot函数和maptools、sp、sf等许多包可以实现相应的统计制图，其中最广泛应用的当属ggplot2和sf包了。本研究旨在应用R与SAS绘制统计地图，并比较其优缺点。

本研究应用地图数据分别是SAS软件SAS/GRAPH中自带的中国地图，和GADM(global admin-istrativeareas)开源数据库[2]的数据，包含中国(除港澳台地区)地图、中国香港特别行政区地图、中国澳门特别行政区地图和中国台湾地区地图。应用这两种软件，可以实现除港澳台地区中国省级行政区划卫生机构数量[3]统计地图的绘制。所有制图均保持基本参数调节一致，不进行过度修饰。

程序介绍

1.SAS数据结构和程序

SAS逻辑库下MAP路径下有China1和China2两个文件，其中China1中包含绘图所需要的经纬度和省份的ID，China2中包含各省份ID对应的省份的信息(包含现用名、别名、曾用名等)。我们只需要在自己的数据中将省份数据对应到唯一识别的ID即可(图 1)。需要注意的是，没有数据的省份也需要对应ID，否则绘制的地图将会缺少无ID的省份。代码介绍如下：

图1 SAS数据结构(左为China2，右为个人数据)

/*此处省略导入数据的过程，我们默认将数据导入到WORK下数据集名叫China*/

/*定义图例标题*/

LEGEND1 LABEL=(“卫生机构数”);

/*定义图例分段和其文字内容*/

PROC FORMAT;

VALUE MK LOW-10000=“0～1万”

10000-20000=“1～2万”

20000-30000=“2～3万”

30000-40000=“3～4万”

40000-HIGH=“4万以上”;

RUN;

QUIT;

/*绘制地图*/

PROC GMAP MAP=MAPS.CHINA DATA=CHINA; /*MAP为地图数据，DATA为统计数据*/

FORMAT V MK.; /*将统计指标V按照上面定义的格式进行格式化*/

ID ID;/*标识出ID变量*/

CHORO V/CDEFAULT=RED CEMPTY=GREEN LEGEND=LEGEND1 LEVELS=ALL; /*CHORO为可选的几种图形之一，其他还包括BLOCK，PRISM等；CDEFAULT用来标识无数据区域的外边颜色；CEMPTY用来标识无数据区域内部填充色；LEGEND连接自定义图例标题；LEVELS=ALL显示全部图例分级*/

RUN;

QUIT;

如果希望将图例按照不同颜色而非渐变色表示，读者可通过将需要绘制的变量变成字符型变量，然后进行上述操作。

2.R数据结构和程序

我们先调用包，如果没有安装请先安装。

library(sf)

library(ggplot2)

library(dplyr)

然后导入从GADM下载的地图(CHN为中国除港澳台地区；TWN为中国台湾地区；HKG为中国香港地区；MAC为中国澳门地区)和数据。地图数据中包含唯一识别的ID，变量名为GID，GID后面的数字编码为区划编码，这里用省级GID_1；同时地图数据中还包含名字、区划等级等信息(图 2)。

图2 CHN_1地图数据结构

x <-st_read(“e：/gadm36_CHN_shp/gadm36_CHN_1.shp”)

y1 <-st_read(“e：/gadm36_TWN_shp/gadm36_TWN_0.shp”)

y2 <-st_read(“e：/gadm36_HKG_shp/gadm36_HKG_0.shp”)

y3 <-st_read(“e：/gadm36_MAC_shp/gadm36_MAC_0.shp”)

y <-read.csv(“e：china.csv”)

同时把统计数据和中国(除港澳台地区)地图按照GID_1做左连接。

total <-left_join(x，y，by=“GID_1”)

连接完成后，按照分组新建一个group变量。

total[total$v<=10000，“group”]<-“0～1万”

total[total$v<=20000 &total$v> 10000，“group”] <-“1～2万”

total[total$v<=30000 &total$v> 20000，“group”] <-“2～3万”

total[total$v<=40000 &total$v> 30000，“group”] <-“3～4万”

total[total$v>40000，“group”] <-“4万以上”

最后画图(图例为渐变色)，其中geom_sf中data是数据，fill后面是统计数据变量的名字；之后绘制中国港澳台地区的部分，其中colour代表轮廓线颜色，fill代表填充颜色。

ggplot()+ geom_sf(data=total，aes(fill=v))+ geom_sf(data=y1，colour=“green”，fill=“red”)+ geom_sf(data=y2，colour=“green”，fill=“red”)+ geom_sf(data=y3，colour=“green”，fill=“red”)+ labs(fill=“卫生机构数”)

可改变将fill=group，使统计地图以不同颜色图例替换渐变色显示。

结果与讨论

应用以上程序将中国除港澳台地区的31个省级行政区划的卫生机构总数绘制成统计地图，详见图3。

图3中，①③由R绘制，②④由SAS绘制。港澳台地区由于没有数据，全部用红色填充绿色描边。①②的统计量为数值型(FORMAT并不会影响到变量本身的类型)，R和SAS都默认用渐变色填充，其中②是使用FORMAT将渐变色水平定为5，否则SAS默认按比例划分图例区间。我们还观察到ggplot2默认数值由小到大颜色由深到浅，SAS PROC GMAP与此相反。③④的统计量为字符型变量，SAS和R都采用了不同颜色填充的形式，在默认的情况下，ggplot2颜色更鲜艳。SAS的图例默认位于底部，ggplot2的图例默认位于图外右侧。ggplot2默认绘制横纵参考线，底色灰色且保留横纵轴及坐标，SASPROC GMAP默认模板底色白色，并无参考线横纵轴及坐标。

图3 SAS与R绘制的统计地图

从布局来看，SAS默认参数下绘制的统计地图更符合一般要求，不需要额外去除参考线和坐标轴，同时参数设置简单，比较适合对作图精细程度要求不高的情况。但是SAS/GRAPH模块虽然功能强大，可语法并不紧凑，PROCGMAP中的参数只是可调整参数的冰山一角，细致调节需要翻阅PROCGMAP的帮助文件，对于SAS/GRAPH模块的新手十分不友好。R的ggplot2包功能强大，语法连续性和可读性非常强。但是，对于完全没有基础的人，去框线作修饰的工作很难入门，SAS PROC GMAP默认参数更符合基本需求。如果有一定基础，再学习的成本较SAS/GRAPH低得多。

小 结

SAS和R在绘制统计地图上的表现难分伯仲。对于熟悉ggplot2制图的读者，R的学习成本低；对于有SAS绘图基础的读者，PROC GMAP是不二之选。针对全无基础的读者，PROC GMAP方便快捷，再结合 Adobe Photoshop、Adobe Illustrator或者Corel DRAW，很容易做出令人满意的统计地图。