[Python]ステップワイズの実行

[Python]ステップワイズの実行
最近、Pythonのsklearnばかりで作業をしております。その際、たまに発生するのが「あれ?sklearnの回帰分析はP値を出力してくれないの?え?ステップワイズもしてくれないの?Rはしてくれるのに。。。」というもう少し推測統計的観点から結果を考察したい場面。statsmodelsというパッケージもあるらしいのですが、なんだかk使い勝手は悪そう。。。ということで今回は、PythonからRを実行して、Python上でステップワイズをしてしまおう!という内容です。

準備

準備として以下をインスールしておきます。Rは言わずとしれた統計ソフトウェア。pyperはPython上からRを利用するモジュールです。

  • R
  • pyper

Rのインストール方法はWeb上にたくさんあるので、そちらをご参考下さい。pyperはpipからインスールすることが可能です。

pip install pyper

では、まずは普通に重回帰分析を行ってみましょう。

データ

モジュールの読み込みとダミーデータを作成します。

In [1]:
import pyper
import pandas as pd
import numpy as np
In [2]:
df = pd.DataFrame({'seg':['segA','segB','segC','segA','segB','segC'],
                  'y':np.random.normal(size=6),
                  'x1':np.random.normal(size=6),
                  'x2':np.random.normal(size=6),
                  'x3':np.random.normal(size=6),
                  'x4':np.random.normal(size=6)
                  })
df
Out[2]:
seg x1 x2 x3 x4 y
0 segA -1.226449 -0.310213 0.712231 0.131545 0.805389
1 segB 1.139922 -1.138528 0.285189 2.123124 0.912482
2 segC -1.138730 -0.265708 0.947262 0.939928 -1.081869
3 segA -0.135973 0.182349 0.346676 2.254102 -0.245571
4 segB -2.145796 0.382971 -0.386438 0.143505 -0.851627
5 segC -0.529889 -0.234263 -0.288878 1.039283 0.090373

データフレームをRインスタンスに渡す

pyper.R()を使ってRのインスタンスを作ります。作成されたら、assign()を使ってデータをR側に渡します。

In [3]:
# R のインスタンスを作る
r = pyper.R(use_pandas='True')
r.assign("data", df)

print(r("data"))
try({data})
   seg         x1         x2         x3        x4           y
0 segA -1.2264491 -0.3102126  0.7122307 0.1315448  0.80538924
1 segB  1.1399218 -1.1385276  0.2851889 2.1231235  0.91248213
2 segC -1.1387303 -0.2657078  0.9472624 0.9399282 -1.08186906
3 segA -0.1359735  0.1823489  0.3466757 2.2541022 -0.24557137
4 segB -2.1457963  0.3829712 -0.3864384 0.1435050 -0.85162688
5 segC -0.5298892 -0.2342628 -0.2888779 1.0392828  0.09037334

基本的にr(‘〇〇’)の〇〇部分にRのスクリプトを書くことで実行されます。無事にデータを渡せていることが確認できます。

重回帰分析

Rのlm関数を使って、重回帰分析を実行します。

In [4]:
r("result <- lm(y~x1+x2+x3+x4, data=data)")
print(r('summary(result)'))
try({summary(result)})

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = data)

Residuals:
       0        1        2        3        4        5 
 0.06057  0.11901 -0.11303  0.03983  0.14917 -0.25555 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   5.7517     1.6418   3.503    0.177
x1            3.3623     0.9167   3.668    0.169
x2            2.6843     1.1110   2.416    0.250
x3            0.4389     0.3650   1.202    0.442
x4           -2.7601     0.7602  -3.631    0.171

Residual standard error: 0.346 on 1 degrees of freedom
Multiple R-squared:  0.965,	Adjusted R-squared:  0.8251 
F-statistic: 6.895 on 4 and 1 DF,  p-value: 0.2773


おお!ちゃんとRでの計算結果が返されました。ちゃんと各変数の有意水準が計算されています。

ステップワイズ

次にsklearnではできないステップワイズを実行してみます。「全く変数を入れないモデル」と「全ての変数を投入したモデル」を作成します。

In [5]:
r("nullModel <- lm(y~1, data=data)")
r("fullModel <- lm(y~x1+x2+x3+x4, data=data)")
Out[5]:
'try({fullModel <- lm(y~x1+x2+x3+x4, data=data)})\n'

次にstep()を用いて変数選択を行っていきます。今回は増減法(both)で行っていきます。

In [6]:
r("stepModel <- step(nullModel,scope=list(lower=nullModel,upper=fullModel),direction='both')")
print(r("summary(stepModel)"))
try({summary(stepModel)})

Call:
lm(formula = y ~ x2, data = data)

Residuals:
         0          1          2          3          4          5 
 0.7817722  0.0005088 -1.0577554  0.2590769 -0.1318140  0.1482114 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.3091     0.3070  -1.007    0.371
x2           -1.0725     0.5770  -1.859    0.137

Residual standard error: 0.6776 on 4 degrees of freedom
Multiple R-squared:  0.4634,	Adjusted R-squared:  0.3293 
F-statistic: 3.455 on 1 and 4 DF,  p-value: 0.1366


結果をPython側へ抽出

結果をまとめる際に、Rで計算させた結果をPython側に持ってくる必要があるでしょう。Rの計算結果を持ってくるにはr.get()を使います。回帰係数、P値、決定係数を取得してみましょう。

In [7]:
#回帰係数
r.get('stepModel$coefficients')
Out[7]:
array([-0.30908162, -1.07248601])

ちなみに結果はarrayの形で帰ってきます。次にP値と決定係数ですが、少し指定が異なります。summaryの結果を指定します。

In [8]:
#P-values
r.get('summary(stepModel)$coefficients[,"Pr(>|t|)"]')
Out[8]:
array([ 0.37093364,  0.13661399])
In [9]:
#決定係数
r.get('summary(stepModel)$r.squared')
Out[9]:
0.463416817602509

無事にpython側に回帰係数、P値、決定係数を持ってこれました。あとは煮るなり約なり(あえて誤変換w)

まとめ

Pythonのsklearnは実装や予測に重きをおいているためか、従来のP値や調整済み決定係数などの出力はありません。ステップワイズするならlasso回帰しろやコラ!って暗に言われている感が。。。その点、Rは実装や処理速度に難はありますが、主力結果はとても見やすく丁寧なものが多いです。よくマーケティングで用いられる因子分析やコンジョイント分析などのパッケージが豊富なのも魅力ですね。双方を使い分けれれば、分析の幅はさらに広がります。

Rのより詳しい使い方は知りたい方は、以下を書籍などを参考にしてみて下さい♪

Pythonカテゴリの最新記事