式: Rスタイルの式を使用したモデルのフィット¶
バージョン0.5.0以降、statsmodelsはユーザーがRスタイルの式を使用して統計モデルをフィットさせることを可能にしました。内部的に、statsmodelsはpatsyパッケージを使用して、式とデータをモデルフィッティングに使用される行列に変換しています。この式フレームワークは非常に強力であり、このチュートリアルはその表面的な部分に触れているに過ぎません。式言語の完全な説明は、patsyのドキュメントに記載されています：
モジュールと関数の読み込み¶
[1]:
import numpy as np # noqa:F401 needed in namespace for patsy
import statsmodels.api as sm
規約をインポートする¶
statsmodels.formula.api から明示的にインポートできます。
[2]:
from statsmodels.formula.api import ols
または、statsmodels.api の主要な formula 名前空間を使用することもできます。
[3]:
sm.formula.ols
[3]:
<bound method Model.from_formula of <class 'statsmodels.regression.linear_model.OLS'>>
または、以下の規約を使用することができます。
[4]:
import statsmodels.formula.api as smf
これらの名前は、各モデルの from_formula クラスメソッドにアクセスするための便利な方法です。例えば、
[5]:
sm.OLS.from_formula
[5]:
<bound method Model.from_formula of <class 'statsmodels.regression.linear_model.OLS'>>
すべての小文字のモデルは「formula」と「data」引数を受け取りますが、大文字のモデルは「endog」と「exog」のデザイン行列を受け取ります。「formula」には、モデルを「patsy」フォーミュラの観点で記述する文字列を渡します。「data」には、pandas データフレームまたは、変数名を定義する「getitem」を持つ他のデータ構造（構造化配列や変数の辞書など）を渡します。
「dir(sm.formula)」を実行すると、利用可能なモデルのリストが表示されます。
フォーミュラ対応のモデルは、次の汎用呼び出しシグネチャを持ちます: (formula, data, subset=None, *args, **kwargs)
式を使用したOLS回帰¶
まず、Getting Started ページで説明された線形モデルをフィットさせます。データをダウンロードし、列をサブセット化し、欠損値のある観測値をリスト単位で削除します。
[6]:
dta = sm.datasets.get_rdataset("Guerry", "HistData", cache=True)
[7]:
df = dta.data[["Lottery", "Literacy", "Wealth", "Region"]].dropna()
df.head()
[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 |
モデルをフィットさせます:
[8]:
mod = ols(formula="Lottery ~ Literacy + Wealth + Region", data=df)
res = mod.fit()
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: Wed, 28 Aug 2024 Prob (F-statistic): 1.07e-05
Time: 10:31:40 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() 演算子を使用することでそれを実現できます。
[9]:
res = ols(formula="Lottery ~ Literacy + Wealth + C(Region)", data=df).fit()
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: カテゴリ変数のコントラスト符号化システム
演算子¶
すでに見たように、「~」はモデルの左辺と右辺を区切り、「+」はデザイン行列に新しい列を追加します。
変数の削除¶
「-」記号は列/変数を削除するために使用できます。例えば、インタセプトをモデルから削除するには次のようにします：
[10]:
res = ols(formula="Lottery ~ Literacy + Wealth + C(Region) -1 ", data=df).fit()
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つの列の相互作用を持つ新しい列をデザイン行列に追加します。「*」は、掛け算された個々の列も含めます。
[11]:
res1 = ols(formula="Lottery ~ Literacy : Wealth - 1", data=df).fit()
res2 = ols(formula="Lottery ~ Literacy * Wealth - 1", data=df).fit()
print(res1.params, "\n")
print(res2.params)
Literacy:Wealth 0.018176
dtype: float64
Literacy 0.427386
Wealth 1.080987
Literacy:Wealth -0.013609
dtype: float64
演算子を使うことで他にも多くのことが可能です。詳細については、patsyのドキュメントをご覧ください。
関数¶
モデル内の変数にベクトル化された関数を適用することができます。
[12]:
res = smf.ols(formula="Lottery ~ np.log(Literacy)", data=df).fit()
print(res.params)
Intercept 115.609119
np.log(Literacy) -20.393959
dtype: float64
カスタム関数を定義する
[13]:
def log_plus_1(x):
return np.log(x) + 1.0
res = smf.ols(formula="Lottery ~ log_plus_1(Literacy)", data=df).fit()
print(res.params)
Intercept 136.003079
log_plus_1(Literacy) -20.393959
dtype: float64
呼び出し元の名前空間にある関数は、すべてその式で使用できます。
式をサポートしていない（またはまだサポートしていない）モデルでの使用¶
特定の statsmodels 関数が式をサポートしていなくても、patsy の数式言語を使用してデザイン行列を生成することができます。その行列は、endog および exog 引数としてフィッティング関数に渡すことができます。
numpy 配列を生成するには、
[14]:
import patsy
f = "Lottery ~ Literacy * Wealth"
y, X = patsy.dmatrices(f, df, return_type="matrix")
print(y[:5])
print(X[:5])
[[41.]
[38.]
[66.]
[80.]
[79.]]
[[1.000e+00 3.700e+01 7.300e+01 2.701e+03]
[1.000e+00 5.100e+01 2.200e+01 1.122e+03]
[1.000e+00 1.300e+01 6.100e+01 7.930e+02]
[1.000e+00 4.600e+01 7.600e+01 3.496e+03]
[1.000e+00 6.900e+01 8.300e+01 5.727e+03]]
pandas データフレームを生成するために:
[15]:
f = "Lottery ~ Literacy * Wealth"
y, X = patsy.dmatrices(f, df, return_type="dataframe")
print(y[:5])
print(X[:5])
Lottery
0 41.0
1 38.0
2 66.0
3 80.0
4 79.0
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
[16]:
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: Wed, 28 Aug 2024 Prob (F-statistic): 1.32e-06
Time: 10:31:41 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.