Statsmodels formula API と多項式回帰モデル

  • 12
    いいね
  • 0
    コメント
この記事は最終更新日から1年以上が経過しています。

Pythonの統計分析ライブラリであるStatsmodelsには,通常の関数型呼び出しに対応するAPIと,R-styleのモデル式を用いる数式対応API(formula API)の2種類が実装されている.Formula APIは複雑な回帰モデルを使う際に有用であるが,モデル式について理解不足があったので,少し整理してみた.

(記事作成時の環境は,python 2.7.8, patsy 0.3.0, statsmodels 0.61 となります.)

Patsyって何だ

Patsyのドキュメント冒頭には,"It's only a model." と気の利いた一文があるが,統計モデルを表現する数式の扱いをサポートするPythonモジュールである.このPatsyを単独で使用する機会はほとんどなく,StatsmodelsやPyMC3の下請けの形で用いられる.(ドキュメントにはPatsyモジュールのインストール方法が載っていましたが,ほとんどのケースではStatsmodels他のモジュールをインストールする際に一緒に導入しているかと思います.)

使用例 (Statsmodelsドキュメントから引用)

# DataFrame
   Lottery  Literacy  Wealth Region
0       41        37      73      E
1       38        51      22      N
2       66        13      61      C
3       80        46      76      E
4       79        69      83      E

import statsmodels.api as sm
import statsmodels.formula.api as smf

mod = smf.ols(formula='Lottery ~ Literacy + Wealth + Region', data=df)
res = mod.fit()

この'Lottery ~ Literacy + Wealth + Region'の部分がモデル式で,Patsyにより解釈,展開されてStatsmodelsの処理に使われる.容易に想像できる通り,上のモデル式は「被説明変数 Lottery は,独立変数 Literacy,Wealth,Region により説明される」ことを表しており,上のリストは重回帰分析(Multiple Regression)を行っている.

モデル式の文法

Patsyのドキュメントには,さまざまなモデル式について説明されているが,少し分かりにくいな書き方になっている.もともとR-styleの共用を考えての導入なので,Rの解説本が参考になる.以下,代表的なモデル式について示す.

モデル式の内容

モデル式 式の内容
y ~ x 回帰モデル,yはxにより説明される
y ~ x1 + x2 yは,独立変数 x1 と x2 により説明される(重回帰モデル)
y ~ 1 + x1 + x2 y は,(切片 / intercept と)x1,x2 で説明される
y ~ x1:x2 yは,x1 と x2 の交互作用項 (x1 * x2) で説明される
y ~ x1 * x2 yは,x1とx2の各変数,及びその組み合わせで構成される交互作用項により説明される (x1 + x2 + x1* x2 と同じ
y ~ I(x + 5) I()の中は,数式を示す.Patsy解釈をエスケープする

自分の用意したモデル式によって,データがどのように展開されるかを確認したいケースも出てくるが,Patsyモジュールの dmatrices(), dmatrix() は,与えたデータを展開する機能をサポートしている.

dmatrices() は,被説明変数(y)を含めて式を展開する.(第2引数は,テストデータ.)

>>> from patsy import dmatrices, dmatrix
>>> dmatrices("y ~ x1", data)
>>>
(DesignMatrix with shape (8L, 1L)
          y
    1.49408
   -0.20516
    0.31307
   -0.85410
   -2.55299
    0.65362
    0.86444
   -0.74217
   Terms:
     'y' (column 0),
 DesignMatrix with shape (8L, 2L)
   Intercept        x1
           1   1.76405
           1   0.40016
           1   0.97874
           1   2.24089
           1   1.86756
           1  -0.97728
           1   0.95009
           1  -0.15136
   Terms:
     'Intercept' (column 0)
     'x1' (column 1))

dmatrix()は,("y~"を含めない式で)モデル式を展開する.

>>> dmatrix("x1+x2+x1:x2", data)
>>> 
DesignMatrix with shape (8L, 4L)
  Intercept        x1        x2     x1:x2
          1   1.76405  -0.10322  -0.18208
          1   0.40016   0.41060   0.16430
          1   0.97874   0.14404   0.14098
          1   2.24089   1.45427   3.25887
          1   1.86756   0.76104   1.42128
          1  -0.97728   0.12168  -0.11891
          1   0.95009   0.44386   0.42171
          1  -0.15136   0.33367  -0.05050
  Terms:
    'Intercept' (column 0)
    'x1' (column 1)
    'x2' (column 2)
    'x1:x2' (column 3)

上記のようにカラムごとに分けられたデータにより,"x1:x2" が交互作用項(x1 ∗ x2,x1 と x2の積)となっていることが確認できる.ここで,モデル式,"I()"の機能について確認しておく.

>>> dmatrix("x1 * x2", data)
>>> 
DesignMatrix with shape (8L, 4L)
  Intercept        x1        x2     x1:x2
          1   1.76405  -0.10322  -0.18208
          1   0.40016   0.41060   0.16430
          1   0.97874   0.14404   0.14098
          1   2.24089   1.45427   3.25887
          1   1.86756   0.76104   1.42128
          1  -0.97728   0.12168  -0.11891
          1   0.95009   0.44386   0.42171
          1  -0.15136   0.33367  -0.05050
(後略)
>>> dmatrix("I(x1 * x2)", data)
>>>  
DesignMatrix with shape (8L, 2L)
  Intercept  I(x1 * x2)
          1    -0.18208
          1     0.16430
          1     0.14098
          1     3.25887
          1     1.42128
          1    -0.11891
          1     0.42171
          1    -0.05050
(後略)

前半のエスケープ関数 "I()" を使わないケースでは,アスタリスク "∗" をPatsyの「与えた独立変数とその交互作用項に展開する」といった解釈が行われ,"I()"を使うと「乗算演算子」として解釈されている状況が分かる.

興味がわいたので,独立変数3つの"∗" (I() 無し)を調べてみた.

>>> dmatrix("x1 * x2 * x3", mydata)
>>>
DesignMatrix with shape (8L, 8L)
  Columns:
    ['Intercept', 'x1', 'x2', 'x1:x2', 'x3', 'x1:x3', 'x2:x3', 'x1:x2:x3']
  Terms:
    'Intercept' (column 0)
    'x1' (column 1)
    'x2' (column 2)
    'x1:x2' (column 3)
    'x3' (column 4)
    'x1:x3' (column 5)
    'x2:x3' (column 6)
    'x1:x2:x3' (column 7)
  (to view full data, use np.asarray(this_obj))

数値の出力は省略されてしまったが,"x1 * x2 * x3"が(Interceptを含めて)8つの項に展開されることが分かった.('Intercept', 'x1', 'x2', 'x1:x2', 'x3', 'x1:x3', 'x2:x3', 'x1:x2:x3')

因みに累乗(べき乗)も"I()"でエスケープする必要がある.

# I() 無し.
>>> dmatrix("x1 ** 2", data)
>>>
DesignMatrix with shape (8L, 2L)
  Intercept        x1
          1   1.76405
          1   0.40016
          1   0.97874
          1   2.24089
          1   1.86756
          1  -0.97728
          1   0.95009
          1  -0.15136
  Terms:
    'Intercept' (column 0)
    'x1' (column 1)
# x1 のまま.

# I() を使用.
>>> dmatrix("I(x1 **2)", data)
>>>
DesignMatrix with shape (8L, 2L)
  Intercept  I(x1 ** 2)
          1     3.11188
          1     0.16013
          1     0.95793
          1     5.02160
          1     3.48777
          1     0.95507
          1     0.90267
          1     0.02291
  Terms:
    'Intercept' (column 0)
    'I(x1 ** 2)' (column 1)

I()を使わない"x1 ∗∗2" では,x1 はそのままの形が出力され期待する2乗のデータとなっていない.(Patsyモデル式では(a + b + c)**2 のような使い方で,独立変数の組み合わせ項(交互作用項を含む式)を導出するのに使うそうです.)

多項式回帰モデルを使う

ここまで調べたことを踏まえて,多項式回帰モデルを使った回帰分析を行った.対象データは,sine curveのデータ,Fittingに使うモデルは,1次から5次までの(交互作用項を含まない)5種類の多項式モデルである.一般式で表すと次のようになる.

y = a_0 + a_1 x + a_2 x^2 + a_3 x^3 + \cdots + a_n x^n \ \ (n=1, 2, 3, 4, 5)

polynomial_fit.png

上図の通り,フィッティング対象のデータは,2周期分のSineデータである.1次から4次の多項式(緑色)ではあまりよくフィッティングできておらず,5次(赤色)でようやく近づいたという感じである.決定係数(R-squared)は,次の通りになった.

R-squared numbers:
poly_1 : 0.147086633025  # 緑色の線
poly_2 : 0.147086633025  # 緑色の線
poly_3 : 0.261777522262  # 緑色の線
poly_4 : 0.261777522262  # 緑色の線
poly_5 : 0.897860130176  # 赤色の線
sine_1 : 1.0    # 青色の線

sin関数を使ったフィッティングは,周波数も一致していたので,R-squred = 1.0 で完全に一致している.また,1次(poly_1)と2次(poly_w),3次(poly_3) と4次(poly_4)の結果が同一となったのは,"sin"のテーラー展開を見れば,その理由を理解できる.

\sin x = x - \frac{x^3}{3 !} + \frac{x^5}{5 !} - \frac{x^7}{7 !} + \frac{x^9}{9 !} + \cdots 

上記の通り,べき乗の値が奇数の項のみでsinが構成されているので,偶数べき乗項を追加した2次項までの近似(poly_2)や4次項までの近似(poly_4)では,それぞれpoly_1やpoly_3に対してフィッティング精度が上がらないことになる.(但し,これはsinの位相が (x,y)=(0,0) 始まりのsinデータ列に限っていえることであり,位相が多少でもずれれば,偶数べき乗項の影響が発生する状況となるはずである.)

Patsyの機能は,次のPythonモジュールで使われているとのことである.(Patsyドキュメントより引用)
- Statsmodels
- PyMC3 (tutorial)
- HDDM
- rERPy
- UrbanSim

これから使用する機会があるのは,StatsmodelsとPyMC3の2つと予想される.PyMC3では,まだモデル式についてドキュメント化されておらず,どの程度サポートされているか不明であるが,後日,試してみたい.

参考文献