PythonからRを呼べるrpy2ですが、ちょっとはまったのでメモっておきます。生存時間解析、とくにCOX比例ハザードモデルのためのスクリプトです。

データは、ヘッダ一行で、time[tab]event[tab]変数1[tab]変数2...の順。timeは生存時間。eventは死亡、再発など。以下、Pandas 0.8.1, ry2 2.2.2を利用

import pandas
import pandas.rpy.common as com
_d = pandas.read_table('test_data.txt')
demo = com.convert_to_r_dataframe(_d)
import rpy2.robjects as robjects
# チェック
l = robjects.r['length']
l(demo)
from rpy2.robjects.packages import importr
survival = importr('survival')
robjects.globalenv["demo"] = demo
cox_res = robjects.r('coxph(Surv(time,event)~.,data=demo,method="exact")')
s = robjects.r['summary']
s(cox_res)[8][2]

これで、Likelihood ratio testのP-valueがとれます。その他の値も、適宜Rの出力と見比べれば取得できるはずです。

テーマ:プログラミング - ジャンル:コンピュータ