統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

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

上記をPythonstatsmodelsを使って同じように書いてみましょう。irisのデータは事前にRから書き出しておいたものを使用することとし、pandasread_csvで読み込みます。そのデータをstatsmodelsGLMに渡します。

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ではFemaleAdultが推定されているのに対して、PythonではMaleChildが推定されていることがわかります。どうやら参照列(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

これで結果が近くなりましたが、まだ一致していません。回帰係数の絶対値は同じとなりましたが符号が異なっています。目的変数であるYesNoという文字列が、RとPythonで別々に解釈された可能性があります。

今度はPython側でYesNoの定義を修正してみましょう。

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
=================================================================================

これで一致しました。

終わりに

今回見てきたように、同一かつシンプルなモデルであっても、言語の仕様によって出力は容易に変わります。異なる結果が得られたときに大慌てしなくて済むよう、各言語や関数についての仕様を把握しておき、本当に異なる結果となっているのかをチェックするというのが大切です。特にモデルの結果を実システムに組み込む際にはエンジニアによるリファクタリングが行われることも多いでしょうから、エンジニアがハマりそうなポイントを押さえておくと役立ちそうですね。