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

平凱珂 陳平雁△

·方法介紹·

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

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

平凱珂 陳平雁△

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

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

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

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

1.高級接口

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

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

首先我們在Python中創(chuàng)建R語言的矩陣對象。

表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

可見，使用rpy2.robjects.IntVector將Python中的list列表元素轉(zhuǎn)換成了R語言中常用的數(shù)據(jù)結(jié)構(gòu)vector，即向量，這等同于在R語言中使用data=c(99,75,5,21)創(chuàng)建了一個向量。同樣，rpy2.robjects.r中已經(jīng)集成了R語言里matrix這個創(chuàng)建矩陣的函數(shù)，我們直接在Python中將創(chuàng)建好的向量作為其參數(shù)，定義行數(shù)和列數(shù)后即可在Python中直接使用R語言中的矩陣對象。類似的，可以使用rpy2.robject.DataFrame來定義一個R語言中的數(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語言中的函數(shù)名或變量名，并將其賦值給Python變量，就可以創(chuàng)建對應(yīng)函數(shù)或變量的一個Python對象，這里我們把R語言中進(jìn)行χ2檢驗的chisq.test函數(shù)賦值到變量chisq_test中創(chuàng)建該函數(shù)對象，此后只要將chisq_test當(dāng)做普通的Python函數(shù)使用即可，函數(shù)的參數(shù)設(shè)置與R語言完全一致。創(chuàng)建成功后，我們只要將上一步創(chuàng)建的數(shù)據(jù)放入chisq_test函數(shù)作為參數(shù)，就可以得到對應(yīng)的Pearsonχ2檢驗結(jié)果。

由此可見，使用rpy2在Python中調(diào)用R語言的統(tǒng)計函數(shù)是一個非常簡單的過程，且程序具有較強的可讀性。

2.低級接口

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

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

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

In[4]:import rpy2.rinterface as rinterface

rinterface.initr()

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

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

在這里我們可以看到，低級接口中無法直接將R語言的結(jié)果用print函數(shù)打印出來，因為儲存在Python中的對象實際是指向了內(nèi)存地址，要訪問具體的數(shù)據(jù)則需要訪問該變量的各個下標(biāo)，比如進(jìn)行Pearsonχ2檢驗后使用下標(biāo)表示統(tǒng)計量、自由度、P值等檢驗結(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語言不同，Python中數(shù)組的下標(biāo)是從0開始的。此外，儲存在該對象中的數(shù)據(jù)本質(zhì)上是一個R語言中的list，儲存的順序和R語言中默認(rèn)順序是一樣的，所以可以按照原先R語言函數(shù)輸出的順序方便地在Python中獲取。

實 例

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

1.實例背景

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

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

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

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

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

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

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

5.可視化作圖

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

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

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

小 結(jié)

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

△通信作者