R與SAS軟件在繪制統(tǒng)計(jì)地圖方面的應(yīng)用研究

天津醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)教研室(300070)

鄧曉偉 陳云峰 劉紅偉 段同慶 馬 駿△

統(tǒng)計(jì)地圖是用不同的顏色和花紋表示統(tǒng)計(jì)數(shù)值在地理分布上的變化，適宜描述研究指標(biāo)的地理分布[1]。尤其是在流行病學(xué)研究中，統(tǒng)計(jì)地圖用來(lái)描述疾病在地區(qū)、時(shí)間的分布再合適不過(guò)。

地圖的繪制、處理經(jīng)常出現(xiàn)在地理信息系統(tǒng)(geographical information systems)中。隨著學(xué)科交融和統(tǒng)計(jì)學(xué)的發(fā)展，地圖的繪制，特別是統(tǒng)計(jì)地圖的繪制已經(jīng)不僅限于使用固定的地理信息系統(tǒng)專用軟件，如ArcGIS等。SAS作為統(tǒng)計(jì)分析商用軟件有Graph模塊可以實(shí)現(xiàn)統(tǒng)計(jì)制圖；R作為開(kāi)源統(tǒng)計(jì)編程語(yǔ)言，Base中有plot函數(shù)和maptools、sp、sf等許多包可以實(shí)現(xiàn)相應(yīng)的統(tǒng)計(jì)制圖，其中最廣泛應(yīng)用的當(dāng)屬ggplot2和sf包了。本研究旨在應(yīng)用R與SAS繪制統(tǒng)計(jì)地圖，并比較其優(yōu)缺點(diǎn)。

本研究應(yīng)用地圖數(shù)據(jù)分別是SAS軟件SAS/GRAPH中自帶的中國(guó)地圖，和GADM(global admin-istrativeareas)開(kāi)源數(shù)據(jù)庫(kù)[2]的數(shù)據(jù)，包含中國(guó)(除港澳臺(tái)地區(qū))地圖、中國(guó)香港特別行政區(qū)地圖、中國(guó)澳門特別行政區(qū)地圖和中國(guó)臺(tái)灣地區(qū)地圖。應(yīng)用這兩種軟件，可以實(shí)現(xiàn)除港澳臺(tái)地區(qū)中國(guó)省級(jí)行政區(qū)劃衛(wèi)生機(jī)構(gòu)數(shù)量[3]統(tǒng)計(jì)地圖的繪制。所有制圖均保持基本參數(shù)調(diào)節(jié)一致，不進(jìn)行過(guò)度修飾。

程序介紹

1.SAS數(shù)據(jù)結(jié)構(gòu)和程序

SAS邏輯庫(kù)下MAP路徑下有China1和China2兩個(gè)文件，其中China1中包含繪圖所需要的經(jīng)緯度和省份的ID，China2中包含各省份ID對(duì)應(yīng)的省份的信息(包含現(xiàn)用名、別名、曾用名等)。我們只需要在自己的數(shù)據(jù)中將省份數(shù)據(jù)對(duì)應(yīng)到唯一識(shí)別的ID即可(圖 1)。需要注意的是，沒(méi)有數(shù)據(jù)的省份也需要對(duì)應(yīng)ID，否則繪制的地圖將會(huì)缺少無(wú)ID的省份。代碼介紹如下：

圖1 SAS數(shù)據(jù)結(jié)構(gòu)(左為China2，右為個(gè)人數(shù)據(jù))

/*此處省略導(dǎo)入數(shù)據(jù)的過(guò)程，我們默認(rèn)將數(shù)據(jù)導(dǎo)入到WORK下數(shù)據(jù)集名叫China*/

/*定義圖例標(biāo)題*/

LEGEND1 LABEL=(“衛(wèi)生機(jī)構(gòu)數(shù)”);

/*定義圖例分段和其文字內(nèi)容*/

PROC FORMAT;

VALUE MK LOW-10000=“0～1萬(wàn)”

10000-20000=“1～2萬(wàn)”

20000-30000=“2～3萬(wàn)”

30000-40000=“3～4萬(wàn)”

40000-HIGH=“4萬(wàn)以上”;

RUN;

QUIT;

/*繪制地圖*/

PROC GMAP MAP=MAPS.CHINA DATA=CHINA; /*MAP為地圖數(shù)據(jù)，DATA為統(tǒng)計(jì)數(shù)據(jù)*/

FORMAT V MK.; /*將統(tǒng)計(jì)指標(biāo)V按照上面定義的格式進(jìn)行格式化*/

ID ID;/*標(biāo)識(shí)出ID變量*/

CHORO V/CDEFAULT=RED CEMPTY=GREEN LEGEND=LEGEND1 LEVELS=ALL; /*CHORO為可選的幾種圖形之一，其他還包括BLOCK，PRISM等；CDEFAULT用來(lái)標(biāo)識(shí)無(wú)數(shù)據(jù)區(qū)域的外邊顏色；CEMPTY用來(lái)標(biāo)識(shí)無(wú)數(shù)據(jù)區(qū)域內(nèi)部填充色；LEGEND連接自定義圖例標(biāo)題；LEVELS=ALL顯示全部圖例分級(jí)*/

RUN;

QUIT;

如果希望將圖例按照不同顏色而非漸變色表示，讀者可通過(guò)將需要繪制的變量變成字符型變量，然后進(jìn)行上述操作。

2.R數(shù)據(jù)結(jié)構(gòu)和程序

我們先調(diào)用包，如果沒(méi)有安裝請(qǐng)先安裝。

library(sf)

library(ggplot2)

library(dplyr)

然后導(dǎo)入從GADM下載的地圖(CHN為中國(guó)除港澳臺(tái)地區(qū)；TWN為中國(guó)臺(tái)灣地區(qū)；HKG為中國(guó)香港地區(qū)；MAC為中國(guó)澳門地區(qū))和數(shù)據(jù)。地圖數(shù)據(jù)中包含唯一識(shí)別的ID，變量名為GID，GID后面的數(shù)字編碼為區(qū)劃編碼，這里用省級(jí)GID_1；同時(shí)地圖數(shù)據(jù)中還包含名字、區(qū)劃等級(jí)等信息(圖 2)。

圖2 CHN_1地圖數(shù)據(jù)結(jié)構(gòu)

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”)

同時(shí)把統(tǒng)計(jì)數(shù)據(jù)和中國(guó)(除港澳臺(tái)地區(qū))地圖按照GID_1做左連接。

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

連接完成后，按照分組新建一個(gè)group變量。

total[total$v<=10000，“group”]<-“0～1萬(wàn)”

total[total$v<=20000 &total$v> 10000，“group”] <-“1～2萬(wàn)”

total[total$v<=30000 &total$v> 20000，“group”] <-“2～3萬(wàn)”

total[total$v<=40000 &total$v> 30000，“group”] <-“3～4萬(wàn)”

total[total$v>40000，“group”] <-“4萬(wàn)以上”

最后畫(huà)圖(圖例為漸變色)，其中g(shù)eom_sf中data是數(shù)據(jù)，fill后面是統(tǒng)計(jì)數(shù)據(jù)變量的名字；之后繪制中國(guó)港澳臺(tái)地區(qū)的部分，其中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=“衛(wèi)生機(jī)構(gòu)數(shù)”)

可改變將fill=group，使統(tǒng)計(jì)地圖以不同顏色圖例替換漸變色顯示。

結(jié)果與討論

應(yīng)用以上程序?qū)⒅袊?guó)除港澳臺(tái)地區(qū)的31個(gè)省級(jí)行政區(qū)劃的衛(wèi)生機(jī)構(gòu)總數(shù)繪制成統(tǒng)計(jì)地圖，詳見(jiàn)圖3。

圖3中，①③由R繪制，②④由SAS繪制。港澳臺(tái)地區(qū)由于沒(méi)有數(shù)據(jù)，全部用紅色填充綠色描邊。①②的統(tǒng)計(jì)量為數(shù)值型(FORMAT并不會(huì)影響到變量本身的類型)，R和SAS都默認(rèn)用漸變色填充，其中②是使用FORMAT將漸變色水平定為5，否則SAS默認(rèn)按比例劃分圖例區(qū)間。我們還觀察到ggplot2默認(rèn)數(shù)值由小到大顏色由深到淺，SAS PROC GMAP與此相反。③④的統(tǒng)計(jì)量為字符型變量，SAS和R都采用了不同顏色填充的形式，在默認(rèn)的情況下，ggplot2顏色更鮮艷。SAS的圖例默認(rèn)位于底部，ggplot2的圖例默認(rèn)位于圖外右側(cè)。ggplot2默認(rèn)繪制橫縱參考線，底色灰色且保留橫縱軸及坐標(biāo)，SASPROC GMAP默認(rèn)模板底色白色，并無(wú)參考線橫縱軸及坐標(biāo)。

圖3 SAS與R繪制的統(tǒng)計(jì)地圖

從布局來(lái)看，SAS默認(rèn)參數(shù)下繪制的統(tǒng)計(jì)地圖更符合一般要求，不需要額外去除參考線和坐標(biāo)軸，同時(shí)參數(shù)設(shè)置簡(jiǎn)單，比較適合對(duì)作圖精細(xì)程度要求不高的情況。但是SAS/GRAPH模塊雖然功能強(qiáng)大，可語(yǔ)法并不緊湊，PROCGMAP中的參數(shù)只是可調(diào)整參數(shù)的冰山一角，細(xì)致調(diào)節(jié)需要翻閱PROCGMAP的幫助文件，對(duì)于SAS/GRAPH模塊的新手十分不友好。R的ggplot2包功能強(qiáng)大，語(yǔ)法連續(xù)性和可讀性非常強(qiáng)。但是，對(duì)于完全沒(méi)有基礎(chǔ)的人，去框線作修飾的工作很難入門，SAS PROC GMAP默認(rèn)參數(shù)更符合基本需求。如果有一定基礎(chǔ)，再學(xué)習(xí)的成本較SAS/GRAPH低得多。

小 結(jié)

SAS和R在繪制統(tǒng)計(jì)地圖上的表現(xiàn)難分伯仲。對(duì)于熟悉ggplot2制圖的讀者，R的學(xué)習(xí)成本低；對(duì)于有SAS繪圖基礎(chǔ)的讀者，PROC GMAP是不二之選。針對(duì)全無(wú)基礎(chǔ)的讀者，PROC GMAP方便快捷，再結(jié)合 Adobe Photoshop、Adobe Illustrator或者Corel DRAW，很容易做出令人滿意的統(tǒng)計(jì)地圖。