RとPythonの結果を一致させたい
背景
通常、私はデータを分析する際には主にRを使用します。しかし分析結果を既存のシステムに投入するなど実装を考えた場合には、Rではなく別の言語を求められることもあると思います。 今回、RとPythonのそれぞれでGLMを実行したときに結果が一致するのかを検証する機会があったため、その内容を備忘がてら紹介しておきます。
RからPythonを呼び出す
今回の検証では、RからPythonを呼び出すためのライブラリreticulate
を使用しました。RStudioからR Notebookとして使えば、chunkごとにRとPythonを切り替えることができて大変便利です。
なおMacユーザーでPython環境としてAnacondaをインストールしている場合、.Rprofileに以下を記載しておくことでAnaconda配下のPythonを呼び出すことができます。
Sys.setenv(PATH = paste("/anaconda3/bin", Sys.getenv("PATH"), sep=":"))
では早速reticulate
をインストールし、Pythonを呼び出せるか試してみましょう。
install.packages("reticulate") library(reticulate)
> system("python --version") Python 3.6.4 :: Anaconda custom (64-bit)
このような表示が出れば成功です。
線形回帰
まずはじめに、irisのデータを用いて線形回帰をかけてみます。Rでは以下のように書くことができます。
> glm(Sepal.Length ~ Petal.Length, data = iris, family = gaussian) Call: glm(formula = Sepal.Length ~ Petal.Length, family = gaussian, data = iris) Coefficients: (Intercept) Petal.Length 4.3066 0.4089 Degrees of Freedom: 149 Total (i.e. Null); 148 Residual Null Deviance: 102.2 Residual Deviance: 24.53 AIC: 160
上記をPythonのstatsmodels
を使って同じように書いてみましょう。irisのデータは事前にRから書き出しておいたものを使用することとし、pandas
のread_csvで読み込みます。そのデータをstatsmodels
のGLMに渡します。
import statsmodels.api as sm import pandas as pd iris = pd.read_csv("/YourDirectory/iris.csv") res_iris = sm.GLM(iris['Sepal.Length'], iris['Petal.Length'], family = sm.families.Gaussian()).fit() print(res_iris.summary()) Generalized Linear Model Regression Results ============================================================================== Dep. Variable: Sepal.Length No. Observations: 150 Model: GLM Df Residuals: 149 Model Family: Gaussian Df Model: 0 Link Function: identity Scale: 3.521367262403398 Method: IRLS Log-Likelihood: -306.75 Date: Mon, 04 Jun 2018 Deviance: 524.68 Time: 18:30:59 Pearson chi2: 525. No. Iterations: 2 ================================================================================ coef std err z P>|z| [0.025 0.975] -------------------------------------------------------------------------------- Petal.Length 1.3489 0.037 36.530 0.000 1.277 1.421 ================================================================================
ひとまず分析を実行することはできたようですが、Rの結果とは違いますね。出力されたテーブルの下部にある「coef」を見てみると、切片が出力されていません。statsmodels
はデフォルトでは切片が入らないようです。
statsmodels
で切片を追加するにはadd_constantを使います。以下のように修正しましょう。
import statsmodels.api as sm import pandas as pd iris = pd.read_csv("/YourDirectory/iris.csv") res_iris_2 = sm.GLM(iris['Sepal.Length'], sm.add_constant(iris['Petal.Length']), family = sm.families.Gaussian()).fit() print(res_iris_2.summary()) Generalized Linear Model Regression Results =============================================================================== Dep. Variable: Sepal.Length No. Observations: 150 Model: GLM Df Residuals: 148 Model Family: Gaussian Df Model: 1 Link Function: identity Scale: 0.16570968760697133 Method: IRLS Log-Likelihood: -77.020 Date: Thu, 31 May 2018 Deviance: 24.525 Time: 13:37:26 Pearson chi2: 24.5 No. Iterations: 2 ================================================================================ coef std err z P>|z| [0.025 0.975] -------------------------------------------------------------------------------- const 4.3066 0.078 54.939 0.000 4.153 4.460 Petal.Length 0.4089 0.019 21.646 0.000 0.372 0.446 ================================================================================
無事、Rの結果と一致しました。
ロジスティック回帰
続いてTitanicのデータを用いてロジスティック回帰をかけてみます。Titanicのデータは下記のコードによりdataframeに加工し、iris同様に書き出しておきました。
tmp <- data.frame(Titanic) titanic <- data.frame(Class = rep(tmp$Class, tmp$Freq), Sex = rep(tmp$Sex, tmp$Freq), Age = rep(tmp$Age, tmp$Freq), Survived = rep(tmp$Survived, tmp$Freq))
Rでは以下のように書けます。
> glm(Survived ~ ., data = titanic, family = binomial(link = "logit")) Call: glm(formula = Survived ~ ., family = binomial(link = "logit"), data = titanic) Coefficients: (Intercept) Class2nd Class3rd ClassCrew SexMale AgeChild 2.0438 -1.0181 -1.7778 -0.8577 -2.4201 1.0615 Degrees of Freedom: 2200 Total (i.e. Null); 2195 Residual Null Deviance: 2769 Residual Deviance: 2210 AIC: 2222
これをPythonで実行してみましょう。ここではRのformulaと同じように指定するためにstatsmodels.formula.apiを使用します。
import statsmodels.api as sm import statsmodels.formula.api as smf import pandas as pd import numpy as np dat = pd.read_csv("/YourDirectory/titanic.csv") dat = dat.drop('Unnamed: 0', axis = 1) titanic_formula = "Survived ~ Class + Sex + Age" glm_bin = smf.glm(formula = titanic_formula, data = dat, family = sm.families.Binomial(sm.families.links.logit)).fit() print(glm_bin.summary()) Generalized Linear Model Regression Results ============================================================================================= Dep. Variable: ['Survived[No]', 'Survived[Yes]'] No. Observations: 2201 Model: GLM Df Residuals: 2195 Model Family: Binomial Df Model: 5 Link Function: logit Scale: 1.0 Method: IRLS Log-Likelihood: -1105.0 Date: Mon, 04 Jun 2018 Deviance: 2210.1 Time: 16:58:29 Pearson chi2: 2.25e+03 No. Iterations: 5 ================================================================================= coef std err z P>|z| [0.025 0.975] --------------------------------------------------------------------------------- Intercept -2.0438 0.168 -12.171 0.000 -2.373 -1.715 Class[T.2nd] 1.0181 0.196 5.194 0.000 0.634 1.402 Class[T.3rd] 1.7778 0.172 10.362 0.000 1.441 2.114 Class[T.Crew] 0.8577 0.157 5.451 0.000 0.549 1.166 Sex[T.Male] 2.4201 0.140 17.236 0.000 2.145 2.695 Age[T.Child] -1.0615 0.244 -4.350 0.000 -1.540 -0.583 =================================================================================
やはり、結果が異なります。よく見ると、カテゴリカル変数についてRではFemale
とAdult
が推定されているのに対して、PythonではMale
とChild
が推定されていることがわかります。どうやら参照列(Reference)が異なるようです。
Rの方で、Referenceを変更して確認してみましょう。
titanic$Sex <- relevel(titanic$Sex, ref = "Female") titanic$Age <- relevel(titanic$Age, ref = "Adult") > glm(Survived ~ ., data = titanic, family = binomial(link = "logit")) Call: glm(formula = Survived ~ ., family = binomial(link = "logit"), data = titanic) Coefficients: (Intercept) Class2nd Class3rd ClassCrew SexMale AgeChild 2.0438 -1.0181 -1.7778 -0.8577 -2.4201 1.0615 Degrees of Freedom: 2200 Total (i.e. Null); 2195 Residual Null Deviance: 2769 Residual Deviance: 2210 AIC: 2222
これで結果が近くなりましたが、まだ一致していません。回帰係数の絶対値は同じとなりましたが符号が異なっています。目的変数であるYes
とNo
という文字列が、RとPythonで別々に解釈された可能性があります。
今度はPython側でYes
とNo
の定義を修正してみましょう。
import statsmodels.api as sm import statsmodels.formula.api as smf import pandas as pd import numpy as np dat = pd.read_csv("/YourDirectory/titanic.csv") dat = dat.drop('Unnamed: 0', axis = 1) dat['Survived'] = dat['Survived'].str.replace('No', '0') dat['Survived'] = dat['Survived'].str.replace('Yes', '1') dat['Survived'] = dat['Survived'].astype(np.int64) titanic_formula = "Survived ~ Class + Sex + Age" glm_bin_2 = smf.glm(formula = titanic_formula, data = dat, family=sm.families.Binomial(sm.families.links.logit)).fit() print(glm_bin_2.summary()) Generalized Linear Model Regression Results ============================================================================== Dep. Variable: Survived No. Observations: 2201 Model: GLM Df Residuals: 2195 Model Family: Binomial Df Model: 5 Link Function: logit Scale: 1.0 Method: IRLS Log-Likelihood: -1105.0 Date: Mon, 04 Jun 2018 Deviance: 2210.1 Time: 17:03:20 Pearson chi2: 2.25e+03 No. Iterations: 5 ================================================================================= coef std err z P>|z| [0.025 0.975] --------------------------------------------------------------------------------- Intercept 2.0438 0.168 12.171 0.000 1.715 2.373 Class[T.2nd] -1.0181 0.196 -5.194 0.000 -1.402 -0.634 Class[T.3rd] -1.7778 0.172 -10.362 0.000 -2.114 -1.441 Class[T.Crew] -0.8577 0.157 -5.451 0.000 -1.166 -0.549 Sex[T.Male] -2.4201 0.140 -17.236 0.000 -2.695 -2.145 Age[T.Child] 1.0615 0.244 4.350 0.000 0.583 1.540 =================================================================================
これで一致しました。
終わりに
今回見てきたように、同一かつシンプルなモデルであっても、言語の仕様によって出力は容易に変わります。異なる結果が得られたときに大慌てしなくて済むよう、各言語や関数についての仕様を把握しておき、本当に異なる結果となっているのかをチェックするというのが大切です。特にモデルの結果を実システムに組み込む際にはエンジニアによるリファクタリングが行われることも多いでしょうから、エンジニアがハマりそうなポイントを押さえておくと役立ちそうですね。