南方醫(yī)科大學(xué)公共衛(wèi)生學(xué)院 生物統(tǒng)計(jì)學(xué)系(510515)

平凱珂 陳平雁△

·方法介紹·

Python與R語(yǔ)言聯(lián)合應(yīng)用的實(shí)現(xiàn)*

南方醫(yī)科大學(xué)公共衛(wèi)生學(xué)院 生物統(tǒng)計(jì)學(xué)系(510515)

平凱珂 陳平雁△

R語(yǔ)言作為主流趨勢(shì)的統(tǒng)計(jì)軟件正在發(fā)揮越來(lái)越重要的作用。而Python語(yǔ)言作為一種面向?qū)ο蟮母呒?jí)編程語(yǔ)言[1]，易學(xué)易用，特別是近年來(lái)以Anaconda為首的Python發(fā)行版整合了大量的數(shù)據(jù)科學(xué)運(yùn)算工具，使得Python能在數(shù)據(jù)處理領(lǐng)域發(fā)揮更大的作用。在實(shí)際應(yīng)用中，一方面，R語(yǔ)言統(tǒng)計(jì)分析功能強(qiáng)大，但運(yùn)算速度較慢且內(nèi)存管理效率不高[2]，導(dǎo)致其在大型項(xiàng)目的分析與管理方面處于劣勢(shì)；另一方面，雖然Python有完善的內(nèi)存優(yōu)化系統(tǒng)，且可以方便地完成數(shù)據(jù)獲取、軟件接口對(duì)接和數(shù)據(jù)庫(kù)交換等數(shù)據(jù)管理工作，但在數(shù)據(jù)分析方面又需借助R語(yǔ)言得以實(shí)現(xiàn)。鑒于此，本研究將結(jié)合兩者的優(yōu)勢(shì)，介紹Python對(duì)R語(yǔ)言的調(diào)用方法，并通過(guò)實(shí)例展示其可行性。

Python與R語(yǔ)言合并應(yīng)用方法介紹

Python對(duì)R語(yǔ)言的調(diào)用使用的是Python中的一個(gè)程序庫(kù)rpy2[3]，它可以實(shí)現(xiàn)Python對(duì)R語(yǔ)言對(duì)象、函數(shù)、方法的讀取與調(diào)用。當(dāng)操作系統(tǒng)中安裝了Python 2.7和R 3.0以上版本時(shí)，在命令行中運(yùn)行pip install rpy2即可安裝該程序庫(kù)[4]。

rpy2提供了兩種調(diào)用方式，高級(jí)接口(high-level interface)和低級(jí)接口(low-level interface)，高級(jí)接口是對(duì)R語(yǔ)言的一個(gè)高級(jí)封裝，它將R語(yǔ)言中常見(jiàn)的一些函數(shù)和對(duì)象封裝成了原生Python對(duì)象或類(lèi)，使用戶可以像原生語(yǔ)言一樣在Python中使用R語(yǔ)言。

1.高級(jí)接口

在Python中引入rpy2.robjects模塊即可使用該高級(jí)接口，例如：

此后rpy2.robjects.r中封裝好的R語(yǔ)言對(duì)象就可以簡(jiǎn)寫(xiě)成robjects.r了。下面我們以兩樣本率比較的χ2檢驗(yàn)為例，介紹Python調(diào)用R語(yǔ)言的高級(jí)接口使用方法。數(shù)據(jù)如表1所示。

首先我們?cè)赑ython中創(chuàng)建R語(yǔ)言的矩陣對(duì)象。

表1 兩組降低顱內(nèi)壓有效率的比較[5]

In[2]:data=robjects.IntVector([99,75,5,21])

test_data=robjects.r.matrix(data,nrow=2,ncol=2)

print test_data

Out[2]:[,1][,2]

[1,]99 5

[2,]75 21

可見(jiàn)，使用rpy2.robjects.IntVector將Python中的list列表元素轉(zhuǎn)換成了R語(yǔ)言中常用的數(shù)據(jù)結(jié)構(gòu)vector，即向量，這等同于在R語(yǔ)言中使用data=c(99,75,5,21)創(chuàng)建了一個(gè)向量。同樣，rpy2.robjects.r中已經(jīng)集成了R語(yǔ)言里matrix這個(gè)創(chuàng)建矩陣的函數(shù)，我們直接在Python中將創(chuàng)建好的向量作為其參數(shù)，定義行數(shù)和列數(shù)后即可在Python中直接使用R語(yǔ)言中的矩陣對(duì)象。類(lèi)似的，可以使用rpy2.robject.DataFrame來(lái)定義一個(gè)R語(yǔ)言中的數(shù)據(jù)框。

In[3]:chisq_test=robjects.r[′chisq.test′]

test_result=chisq_test(test_data,correct=False)

print(test_result)

Out[3]:

Pearson′s Chi-squared test

data:structure(c(99L,75L,5L,21L),.Dim=c(2L,2L))

X-squared=12.857,df=1,p-value=0.0003362

在rpy2.robjects.r[′r_code′]這一句中，若將r_code替換成R語(yǔ)言中的函數(shù)名或變量名，并將其賦值給Python變量，就可以創(chuàng)建對(duì)應(yīng)函數(shù)或變量的一個(gè)Python對(duì)象，這里我們把R語(yǔ)言中進(jìn)行χ2檢驗(yàn)的chisq.test函數(shù)賦值到變量chisq_test中創(chuàng)建該函數(shù)對(duì)象，此后只要將chisq_test當(dāng)做普通的Python函數(shù)使用即可，函數(shù)的參數(shù)設(shè)置與R語(yǔ)言完全一致。創(chuàng)建成功后，我們只要將上一步創(chuàng)建的數(shù)據(jù)放入chisq_test函數(shù)作為參數(shù)，就可以得到對(duì)應(yīng)的Pearsonχ2檢驗(yàn)結(jié)果。

由此可見(jiàn)，使用rpy2在Python中調(diào)用R語(yǔ)言的統(tǒng)計(jì)函數(shù)是一個(gè)非常簡(jiǎn)單的過(guò)程，且程序具有較強(qiáng)的可讀性。

2.低級(jí)接口

R語(yǔ)言存在底層基礎(chǔ)運(yùn)算(如加法、減法、乘方、開(kāi)方、循環(huán)等運(yùn)算)速度較慢的劣勢(shì)，這主要是因?yàn)镽語(yǔ)言是一種解釋型語(yǔ)言，每一句代碼都會(huì)被處理器臨時(shí)編譯并下達(dá)指令再計(jì)算[2]。但是其在矩陣運(yùn)算的速度上卻有著優(yōu)勢(shì)，特別是其中的很多統(tǒng)計(jì)運(yùn)算包，直接從C語(yǔ)言層面實(shí)現(xiàn)了很多復(fù)雜統(tǒng)計(jì)方法的運(yùn)算，有時(shí)能達(dá)到Python中numpy、scipy這些數(shù)值運(yùn)算庫(kù)所達(dá)不到的運(yùn)行效率。

rpy2的低級(jí)接口主要應(yīng)用于高級(jí)接口沒(méi)有提供封裝，或?qū)\(yùn)算性能有特殊要求的情況。

在Python中引入rpy2.rinterface模塊即可使用該低級(jí)接口，例如：

In[4]:import rpy2.rinterface as rinterface

rinterface.initr()

這里注意一定要使用initr()函數(shù)來(lái)初始化R語(yǔ)言，才可以正確地使用后續(xù)接口，還要注意這句代碼只用運(yùn)行一次，以免出現(xiàn)錯(cuò)誤。

低級(jí)接口的調(diào)用方法與高級(jí)接口類(lèi)似，不同之處在于低級(jí)接口可以訪問(wèn)當(dāng)前R語(yǔ)言運(yùn)行環(huán)境中的所有R語(yǔ)言對(duì)象，這一過(guò)程主要是通過(guò)rpy2.rinterface.globalenv類(lèi)中的get函數(shù)來(lái)完成，這個(gè)函數(shù)會(huì)在Python中創(chuàng)建一個(gè)指向R語(yǔ)言環(huán)境的sexp指針對(duì)象。我們同樣以兩樣本率比較的χ2檢驗(yàn)為例來(lái)說(shuō)明這一過(guò)程，數(shù)據(jù)參見(jiàn)表1。

In[5]:matrix=rinterface.globalenv.get(“matrix”)

data=rinterface.IntSexpVector([99,75,5,21])

test_data=matrix(data,nrow=2,ncol=2)

print test_data

Out[5]:

8/R:0x00000000071CFC30>

在這里我們可以看到，低級(jí)接口中無(wú)法直接將R語(yǔ)言的結(jié)果用print函數(shù)打印出來(lái)，因?yàn)閮?chǔ)存在Python中的對(duì)象實(shí)際是指向了內(nèi)存地址，要訪問(wèn)具體的數(shù)據(jù)則需要訪問(wèn)該變量的各個(gè)下標(biāo)，比如進(jìn)行Pearsonχ2檢驗(yàn)后使用下標(biāo)表示統(tǒng)計(jì)量、自由度、P值等檢驗(yàn)結(jié)果。

In[6]:chisq_test=rinterface.globalenv.get(′chisq.test′)

test_result=chisq_test(test_data,correct=False)

print(“Method:%s Chi-squared=%s,df=%s,p-value=%s”%

(test_result[3][0],test_result[0][0],test_result[1][0],test_result[2][0]))

Out[6]:Method:Pearson′s Chi-squared test

Chi-squared=12.8570699857,df=1,

p-value=0.000336206596885

需注意，與R語(yǔ)言不同，Python中數(shù)組的下標(biāo)是從0開(kāi)始的。此外，儲(chǔ)存在該對(duì)象中的數(shù)據(jù)本質(zhì)上是一個(gè)R語(yǔ)言中的list，儲(chǔ)存的順序和R語(yǔ)言中默認(rèn)順序是一樣的，所以可以按照原先R語(yǔ)言函數(shù)輸出的順序方便地在Python中獲取。

實(shí) 例

下面以實(shí)例來(lái)說(shuō)明如何利用Python+R來(lái)完成數(shù)據(jù)獲取及其后續(xù)統(tǒng)計(jì)分析與繪圖的任務(wù)。

1.實(shí)例背景

利用weather underground網(wǎng)提供的氣象數(shù)據(jù)API接口，獲取廣州市近五年的每日氣溫?cái)?shù)據(jù)(從2011年11月1日至2016年11月1日)，再根據(jù)此數(shù)據(jù)利用空間狀態(tài)模型預(yù)測(cè)未來(lái)的每日平均氣溫，最后利用R語(yǔ)言的ggplot2包進(jìn)行時(shí)間序列作圖。

2.數(shù)據(jù)獲取

使用Python中的urllib2庫(kù)和bs4庫(kù)中的BeautifulSoup模塊來(lái)抓取網(wǎng)頁(yè)數(shù)據(jù)，其中urllib2是Python的一個(gè)獲取URLs的組件，它通過(guò)向指定的URL發(fā)出請(qǐng)求來(lái)獲取數(shù)據(jù)，此處直接使用GET的方式請(qǐng)求該網(wǎng)站的氣象數(shù)據(jù)接口，就可以獲取到指定城市在指定日期的氣象數(shù)據(jù)。而B(niǎo)eautifulSoup模塊提供了一些專(zhuān)用于處理導(dǎo)航、搜索、修改分析樹(shù)的函數(shù)，通過(guò)解析文檔為用戶提供需要抓取的數(shù)據(jù)，我們利用此模塊可以自由地根據(jù)需要提取出urllib2庫(kù)中已經(jīng)獲取的氣象數(shù)據(jù)，如廣州市每日的氣溫、風(fēng)力、風(fēng)向、濕度、日落日出時(shí)間等。如果我們僅需要?dú)鉁財(cái)?shù)據(jù)，可以再利用BeautifulSoup模塊單獨(dú)從urllib2庫(kù)中提取。

3.數(shù)據(jù)整理

由于獲取到的數(shù)據(jù)為華氏度，我們需要將數(shù)據(jù)轉(zhuǎn)換為攝氏度。在Python中需要使用循環(huán)語(yǔ)句來(lái)逐個(gè)轉(zhuǎn)換，而在R語(yǔ)言中我們只需要利用它的矩陣特性，在data.frame數(shù)據(jù)框中只要一句代碼即可批量轉(zhuǎn)換。如前所述，我們使用rpy2.robject.DataFrame來(lái)定義一個(gè)R語(yǔ)言中的數(shù)據(jù)框，然后完成數(shù)據(jù)的轉(zhuǎn)換工作。

4.數(shù)據(jù)分析

我們使用R語(yǔ)言中的forecast包來(lái)進(jìn)行時(shí)間序列數(shù)據(jù)的生成與預(yù)測(cè)。首先利用msts函數(shù)來(lái)創(chuàng)建帶有季節(jié)效應(yīng)的時(shí)間序列數(shù)據(jù)，這里為了簡(jiǎn)便起見(jiàn)，我們直接設(shè)定季節(jié)效應(yīng)期為1年，把seasonal.periods參數(shù)設(shè)置成365.25即可。然后我們將該數(shù)據(jù)使用tbats函數(shù)來(lái)進(jìn)行狀態(tài)空間模型中的TBATS模型建模。最后利用forecast函數(shù)來(lái)生成隨后的預(yù)測(cè)數(shù)據(jù)。

5.可視化作圖

分析完成之后，我們可以利用R語(yǔ)言的ggplot2包來(lái)進(jìn)行數(shù)據(jù)可視化作圖，值得一提的是，rpy2中已經(jīng)封裝了R語(yǔ)言中整個(gè)ggplot2包，只要在Python中直接引入rpy2.robjects.lib.ggplot2模塊，即可在Python中直接使用ggplot2的函數(shù)進(jìn)行作圖。

繪制的圖形如圖1所示，此處我們僅繪制了2016年7月1日至2017年7月1日的數(shù)據(jù)，其中2016年11月1日之后的為預(yù)測(cè)數(shù)據(jù)。

圖1 在Python中使用R語(yǔ)言的ggplot2包進(jìn)行繪圖

小 結(jié)

本文介紹了Python對(duì)R語(yǔ)言的調(diào)用方法，并以實(shí)例展現(xiàn)了“Python+R”的模式用作數(shù)據(jù)獲取和統(tǒng)計(jì)運(yùn)算的可行性。整個(gè)過(guò)程中結(jié)合了Python和R語(yǔ)言各自的優(yōu)勢(shì)，兩種語(yǔ)言各司其職，在數(shù)據(jù)獲取和統(tǒng)計(jì)運(yùn)算方面都有著良好表現(xiàn)。該模式現(xiàn)在也已經(jīng)開(kāi)始在一些實(shí)際研究項(xiàng)目中得到應(yīng)用，例如Vandenbulcke等[6]在一個(gè)飲酒可能是肝癌的危險(xiǎn)因素的前瞻性研究中，利用該模式完成了數(shù)據(jù)管理和統(tǒng)計(jì)分析的工作?！癙ython+R”的工作模式可有效結(jié)合二者的優(yōu)勢(shì)，既能較好地完成項(xiàng)目管理和數(shù)據(jù)管理工作，又能利用R語(yǔ)言完成復(fù)雜的統(tǒng)計(jì)分析，有良好的應(yīng)用前景。

[1]赫特蘭，司維，曾軍崴，等.Python基礎(chǔ)教程.北京：人民郵電出版社，2010.

[2]Aloysius L,William T著,唐李洋譯.R高性能編程.電子工業(yè)出版社,2015.

[3]Laurent G.rpy2-R in Python.http://rpy2.bitbucket.org/.

[4]Python Software Foundation.PyPI-the Python Package Index.https://pypi.python.org/pypi.

[5]孫振球,徐勇勇.醫(yī)學(xué)統(tǒng)計(jì)學(xué).北京：人民衛(wèi)生出版社,2014.

[6]Vandenbulcke H,Moreno C,Colle I,et al.Alcohol intake increases the risk of HCC in hepatitis C virus-related compensated cirrhosis:A prospective study.J Hepatol,2016,65(3):543-551

(責(zé)任編輯：張 悅)

*廣東大學(xué)生科技創(chuàng)新培育專(zhuān)項(xiàng)資金資助(pdjh2016a0091)

△通信作者