Help us understand the problem. What is going on with this article?

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

More than 3 years have passed since last update.

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では,まだモデル式についてドキュメント化されておらず,どの程度サポートされているか不明であるが,後日,試してみたい.

参考文献

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした