rpy2を用いてPython上でRを使用した効果検証手法の簡単な実装

カケハシでデータサイエンティストをしている島吉です。カケハシのデータサイエンティストは、AI在庫管理のエンジニアと連携したり、機械学習を使う業務が多いため、データ分析にはPythonを使用することが多いです。

しかし、統計的な手法のライブラリはRに多く存在しています。たとえば、現在の業務では、効果検証に傾向スコアマッチングを使用しており、さまざまな書籍でRを用いた使用例を多く見かけます。

そこで、PythonとRの両方を使用し、Rが適した部分はRで実装し、ほかの処理は使い慣れたPythonで実装しようと考えました。 PythonとRの連携は、たとえば、下記のような方法で実現できます。

Python上でRを使用できるパッケージ（rpy2）を用いる
普段から使用しているDatabricksではノートブック内のセルごとにPython・Rなどの言語を選択できる

今回、Pythonのループや関数のなかでRのライブラリを使いたいケースがあったので、セルごとに言語を変える方法ではなく、rpy2を用いる方法を選択しました。

この記事は秋の技術特集 2024の 14 記事目です。

Python上でRを使用する実装例

rpy2のインストール

!pip install rpy2

パッケージのインポート

import rpy2.robjects as robjects
from rpy2.robjects.packages import importr
from rpy2.robjects import pandas2ri

import pandas as pd
import numpy as np
from sklearn.linear_model import LogisticRegression

サンプルデータ作成

np.random.seed(1)
size = 1000

# 介入群
df_treated = pd.DataFrame({
    "col_1": np.random.normal(50, 100, size),  # 共変量
    "col_2": np.random.normal(150, 100, size), # 共変量
    "col_3": np.random.normal(650, 100, size), # 比較する指標
    "treated_flg": [1]*size # 介入有無
})

# 対照群
df_control = pd.DataFrame({
    "col_1": np.random.normal(100, 100, size), # 共変量
    "col_2": np.random.normal(200, 100, size), # 共変量
    "col_3": np.random.normal(750, 100, size), # 比較する指標
    "treated_flg": [0]*size # 介入有無
})

df = pd.concat([df_treated, df_control], axis=0)
df = df.reset_index(drop=True)

ロジスティック回帰を用いた傾向スコア算出

df_x = df[["col_1", "col_2"]] # 共変量
df_y = df["treated_flg"] # 介入有無（1/0）

model = LogisticRegression()
model.fit(df_x, df_y)

df["propensity_score"] = model.predict_proba(df_x)[:, 1]

Python上でRのMatchingパッケージのインストール

utils = importr("utils") 
utils.install_packages("Matching", repos="https://cloud.r-project.org")

Python上でRのパッケージのインポート

matching = importr("Matching")

Python上でRを用いた傾向スコアマッチングの関数を定義

def r_ps_matching(df):

    # PandasからRにデータを渡す
    pandas2ri.activate()
    robjects.r.assign("r_df", pandas2ri.py2rpy(df))

    # 傾向スコアマッチング
    robjects.r(
    """
    set.seed(1)

    # 傾向スコアが近いもの同士をマッチング（今回は非復元抽出を選択）
    matched <- Match(Y=r_df$col_3, Tr=r_df$treated_flg, X=r_df$propensity_score, caliper=0.2, replace=FALSE)
    
    # 傾向スコアマッチング前後の共変量のバランスを出力
    MatchBalance(treated_flg ~ col_1 + col_2, data=r_df, match.out=matched)

    # 介入効果の出力
    summary(matched)
    """
    )

    # 傾向スコアマッチング後のデータをRからPandasに渡す
    df_matched = robjects.r("rbind(r_df[matched$index.treated,], r_df[matched$index.control,])")
    df_matched = pandas2ri.rpy2py(df_matched)

    return df_matched

Python上でRを用いた傾向スコアマッチングの関数の呼び出し

df_matched = r_ps_matching(df)

Python上でRを使用した効果

使い慣れたPythonにRを連携することで、使用できるライブラリの幅が広がり、傾向スコアマッチングを簡単に実装できました。また、自力で複雑な処理を作る必要がなくなり、ミスを減らすこともできます。

傾向スコアマッチング後は、RからPandasにデータを渡しているので、Pythonでヒストグラムを作成して比較するなどの追加分析・可視化を行うこともできます。

rpy2を用いたので、Pythonの関数内でRを使用できましたが、robjects.r()を記述する必要があり、やや記述が長くなるデメリットもあります。これに対して、セルごとに言語を変える方法では、セルに%rと書くだけで、Rの処理をそのまま記述できます。そのため、用途に応じて使いやすい方法を選択するのがよさそうです。