Rスタイルのフォーミュラを用いたモデルの適合¶
バージョン0.5.0以降、statsmodelsでは、Rスタイルのフォーミュラを使用して統計モデルを適合させることができます。statsmodelsは内部的にpatsyパッケージを使用して、フォーミュラとデータをモデル適合で使用される行列に変換します。フォーミュラフレームワークは非常に強力ですが、このチュートリアルでは表面をなぞるだけです。フォーミュラ言語の完全な説明はpatsyのドキュメントにあります。
モジュールと関数の読み込み¶
In [1]: import statsmodels.api as sm
In [2]: import statsmodels.formula.api as smf
In [3]: import numpy as np
In [4]: import pandas
通常のstatsmodels.apiに加えてstatsmodels.formula.apiを呼び出したことに注意してください。statsmodels.apiはここではデータセットの読み込みにのみ使用されます。formula.apiはapiにある多くの関数（例：OLS、GLM）と同じ関数の大部分を保持していますが、これらのモデルの大部分の対応する小文字の関数も保持しています。一般的に、小文字のモデルはformulaとdf引数を受け入れ、大文字のモデルはendogとexog計画行列を受け入れます。formulaはpatsyフォーミュラでモデルを記述する文字列を受け入れます。dfはpandasデータフレームを受け入れます。
dir(smf)は使用可能なモデルのリストを出力します。
フォーミュラ対応モデルは、次の一般的な呼び出しシグネチャを持ちます：(formula, data, subset=None, *args, **kwargs)
フォーミュラを用いたOLS回帰¶
まず、はじめにページで説明されている線形モデルを適合させます。データをダウンロードし、列のサブセットを作成し、欠損値を除去するためにリストワイズ削除を行います。
In [5]: df = sm.datasets.get_rdataset("Guerry", "HistData").data
In [6]: df = df[['Lottery', 'Literacy', 'Wealth', 'Region']].dropna()
In [7]: df.head()
Out[7]:
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
モデルの適合
In [8]: mod = smf.ols(formula='Lottery ~ Literacy + Wealth + Region', data=df)
In [9]: res = mod.fit()
In [10]: print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Lottery R-squared: 0.338
Model: OLS Adj. R-squared: 0.287
Method: Least Squares F-statistic: 6.636
Date: Thu, 03 Oct 2024 Prob (F-statistic): 1.07e-05
Time: 16:08:49 Log-Likelihood: -375.30
No. Observations: 85 AIC: 764.6
Df Residuals: 78 BIC: 781.7
Df Model: 6
Covariance Type: nonrobust
===============================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------
Intercept 38.6517 9.456 4.087 0.000 19.826 57.478
Region[T.E] -15.4278 9.727 -1.586 0.117 -34.793 3.938
Region[T.N] -10.0170 9.260 -1.082 0.283 -28.453 8.419
Region[T.S] -4.5483 7.279 -0.625 0.534 -19.039 9.943
Region[T.W] -10.0913 7.196 -1.402 0.165 -24.418 4.235
Literacy -0.1858 0.210 -0.886 0.378 -0.603 0.232
Wealth 0.4515 0.103 4.390 0.000 0.247 0.656
==============================================================================
Omnibus: 3.049 Durbin-Watson: 1.785
Prob(Omnibus): 0.218 Jarque-Bera (JB): 2.694
Skew: -0.340 Prob(JB): 0.260
Kurtosis: 2.454 Cond. No. 371.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
カテゴリ変数¶
上記に出力された要約を見ると、patsyは*Region*の要素がテキスト文字列であると判断し、*Region*をカテゴリ変数として扱ったことに注意してください。patsyのデフォルトでは切片も含まれるため、*Region*のカテゴリの1つを自動的に除外しました。
*Region*が整数変数であり、明示的にカテゴリ変数として扱う必要がある場合は、C()演算子を使用することで行うことができます。
In [11]: res = smf.ols(formula='Lottery ~ Literacy + Wealth + C(Region)', data=df).fit()
In [12]: print(res.params)
Intercept 38.651655
C(Region)[T.E] -15.427785
C(Region)[T.N] -10.016961
C(Region)[T.S] -4.548257
C(Region)[T.W] -10.091276
Literacy -0.185819
Wealth 0.451475
dtype: float64
patsyのカテゴリ変数の機能に関するより高度な機能の例については、こちらをご覧ください：Patsy：カテゴリ変数のためのコントラストコーディングシステム
演算子¶
すでに、「〜」がモデルの左辺と右辺を分離し、「+」が計画行列に新しい列を追加することを見てきました。
変数の削除¶
「-」記号を使用して列/変数を削除できます。例えば、モデルから切片を削除するには、
In [13]: res = smf.ols(formula='Lottery ~ Literacy + Wealth + C(Region) -1 ', data=df).fit()
In [14]: print(res.params)
C(Region)[C] 38.651655
C(Region)[E] 23.223870
C(Region)[N] 28.634694
C(Region)[S] 34.103399
C(Region)[W] 28.560379
Literacy -0.185819
Wealth 0.451475
dtype: float64
乗法項相互作用¶
「:」は、他の2つの列の積を持つ新しい列を計画行列に追加します。「*」は、掛け合わせられた個々の列も含まれます。
In [15]: res1 = smf.ols(formula='Lottery ~ Literacy : Wealth - 1', data=df).fit()
In [16]: res2 = smf.ols(formula='Lottery ~ Literacy * Wealth - 1', data=df).fit()
In [17]: print(res1.params)
Literacy:Wealth 0.018176
dtype: float64
In [18]: print(res2.params)
Literacy 0.427386
Wealth 1.080987
Literacy:Wealth -0.013609
dtype: float64
演算子では他にも多くのことが可能です。詳細については、patsyのドキュメントを参照してください。
関数¶
ベクトル化関数をモデルの変数に適用できます。
In [19]: res = smf.ols(formula='Lottery ~ np.log(Literacy)', data=df).fit()
In [20]: print(res.params)
Intercept 115.609119
np.log(Literacy) -20.393959
dtype: float64
カスタム関数の定義
In [21]: def log_plus_1(x):
....: return np.log(x) + 1.0
....:
In [22]: res = smf.ols(formula='Lottery ~ log_plus_1(Literacy)', data=df).fit()
In [23]: print(res.params)
Intercept 136.003079
log_plus_1(Literacy) -20.393959
dtype: float64
名前空間¶
上記のすべての例では、適用する関数を検索するために呼び出し側の名前空間を使用していることに注意してください。使用される名前空間は、eval_envキーワードを介して制御できます。例えば、patsy:patsy.EvalEnvironmentを使用してカスタム名前空間を指定したり、「クリーン」な名前空間を使用したりすることができます（eval_func=-1を渡すことで実現します）。デフォルトでは、呼び出し側の名前空間が使用されます。これは、例えば、ユーザーの名前空間またはpatsyに渡されたデータ構造内にCという名前の変数があり、カテゴリ変数を処理するためにフォーミュラでCが使用されている場合など、（予期せぬ）結果を招く可能性があります。詳細については、Patsy APIリファレンスを参照してください。
（まだ）フォーミュラをサポートしていないモデルでのフォーミュラの使用¶
特定のstatsmodels関数がフォーミュラをサポートしていなくても、patsyのフォーミュラ言語を使用して計画行列を作成できます。それらの行列は、endogとexog引数として適合関数に供給できます。
numpy配列を生成するには
In [24]: import patsy
In [25]: f = 'Lottery ~ Literacy * Wealth'
In [26]: y, X = patsy.dmatrices(f, df, return_type='matrix')
In [27]: print(y[:5])
[[41.]
[38.]
[66.]
[80.]
[79.]]
In [28]: print(X[:5])
[[ 1. 37. 73. 2701.]
[ 1. 51. 22. 1122.]
[ 1. 13. 61. 793.]
[ 1. 46. 76. 3496.]
[ 1. 69. 83. 5727.]]
yとXはnumpy.ndarrayのサブクラスであるpatsy.DesignMatrixのインスタンスになります。
pandasデータフレームを生成するには
In [29]: f = 'Lottery ~ Literacy * Wealth'
In [30]: y, X = patsy.dmatrices(f, df, return_type='dataframe')
In [31]: print(y[:5])
Lottery
0 41.0
1 38.0
2 66.0
3 80.0
4 79.0
In [32]: print(X[:5])
Intercept Literacy Wealth Literacy:Wealth
0 1.0 37.0 73.0 2701.0
1 1.0 51.0 22.0 1122.0
2 1.0 13.0 61.0 793.0
3 1.0 46.0 76.0 3496.0
4 1.0 69.0 83.0 5727.0
In [33]: print(sm.OLS(y, X).fit().summary())
OLS Regression Results
==============================================================================
Dep. Variable: Lottery R-squared: 0.309
Model: OLS Adj. R-squared: 0.283
Method: Least Squares F-statistic: 12.06
Date: Thu, 03 Oct 2024 Prob (F-statistic): 1.32e-06
Time: 16:08:49 Log-Likelihood: -377.13
No. Observations: 85 AIC: 762.3
Df Residuals: 81 BIC: 772.0
Df Model: 3
Covariance Type: nonrobust
===================================================================================
coef std err t P>|t| [0.025 0.975]
-----------------------------------------------------------------------------------
Intercept 38.6348 15.825 2.441 0.017 7.149 70.121
Literacy -0.3522 0.334 -1.056 0.294 -1.016 0.312
Wealth 0.4364 0.283 1.544 0.126 -0.126 0.999
Literacy:Wealth -0.0005 0.006 -0.085 0.933 -0.013 0.012
==============================================================================
Omnibus: 4.447 Durbin-Watson: 1.953
Prob(Omnibus): 0.108 Jarque-Bera (JB): 3.228
Skew: -0.332 Prob(JB): 0.199
Kurtosis: 2.314 Cond. No. 1.40e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.4e+04. This might indicate that there are
strong multicollinearity or other numerical problems.