StatsmodelsはPythonというプログラミング言語上で動く統計解析ソフトです。statsmodelsのサンプルを動かすにはPCにPythonがインストールされている必要があります。まだインストールされていない方はJupyter notebookのインストールを参照してください。Jupyter notebookはstatsmodelsを動かすのに大変便利です。また、回帰分析の基礎はstatsmodelsによる線形回帰入門を参考にしてください。

一般化加法モデル Generalized Additive Models (GAM)

Trevor Hastie と Robert Tibshiraniにより提唱された一般化加法モデルは、統計学的には一般化線形モデルの1つであり、予測変数はある予測変数の未知の滑らかな線形関数であるとし、この滑らかな関数の推定に焦点を当てている。一変量の応答変数$Y$の期待値と予測変数$x_i$の関係をリンク関数$g$として

$$ g(E(Y))=\beta_{0}+f_{1}(x_{1})+f_{2}(x_{2})+\cdots +f_{m}(x_{m})$$

という形式で示している。$E$は期待値を表す。関数$f_m$は多項式とかスプライン曲線とかの特定のパラメトリックな形式であったり、ノンパラメトリックであったり、単なる平滑化関数であったりする。$f_1, \cdots, f_m$を説明変数とその係数の積に限定すると

$$ g(E(Y))=\beta_{0}+\alpha_{1}x_{1}+\alpha_{2}x_{2}+\cdots +\alpha_{m}x_{m}$$

という一般化線形モデルの形式になる。

これは1972年にネルダーとウェダーバーンによって提唱された、一般化線形モデル(GLM)に1981年にJerome H. FriedmanとWerner Stuetzleによって提唱された加法モデルを混合したものである。加法モデルは、1次元の滑らかさを用いた制約付きノンパラメトリック回帰モデルであるACE(Alternating conditional expectations)と呼ばれるアルゴリズムを回帰分析の応答変数と予測変数の間の変換に用い、その最適化を達成している。この方法により、たとえば$p$次元の滑らかさを用いて次元の呪いからの影響を削減しているが、一方でその柔軟性ゆえに近似誤差という費用を払っている。

例:自動車の燃費に影響する要因

自動車の燃費についてのトイモデルを構築する。データはUCIの自動車データから取得する。燃費(city_mpg)、燃料(fuel)、駆動方式(drive)、重量(weight)と馬力(hose power(hp))を選択してデータフレームとしてダウンロードする。

カテゴリ変数を線形の項として扱い、2つの説明変数(重量と馬力)の燃費への効果を罰則付きB-スプライン曲線でとらえる。モデルにはガウス過程回帰とポアソン回帰を用いる。

from statsmodels.gam.tests.test_penalized import df_autos
df_autos.head()

city_mpg	fuel	drive	weight	hp
0	21	gas	rwd	2548	111.0
1	21	gas	rwd	2548	111.0
2	19	gas	rwd	2823	154.0
3	24	gas	fwd	2337	102.0
4	18	gas	4wd	2824	115.0

B-スプライン曲線はBasis splineの略でBasisは基底のことである。これは、スプライン、つまり区分多項式を用いて表現する滑らかな曲線の1つである。与えられた複数の制御点(control point)とノットベクトルから定義される。曲線は必ずしも制御点を通らない。また、一部の変更が曲線全体に影響を与えない点が特徴である。一般にスプライン曲線はBスプライン曲線の線形の和として表される。

BSplines

class statsmodels.gam.smooth_basis.BSplines(x, df, degree, include_intercept=False,
　　　　　　　　　　　　　　　　　　　　 constraints=None,variable_names=None, knot_kwds=None)

    additive smooth components using B-Splines

x: array_like - 滑らかな項に用いる1-D または 2-Dの説明変数。もし2次元であれば観測値が行で、説明変数が列である。

df：int - 基底関数の数、または自由度

degree : int - スプラインの次数

include_intercept : bool - もしFalseであれば、基底関数は定数を含まない。

constraints : None, string or array - 制約を満たすために基底関数を変換するために用いられる制約。constraints = ‘center’は定数の削除と基底関数へ集中させるための線形変換。

variable_names : None or list of strings - 説明変数$x$の列名と基底関数のパラメータ名。$x$がpandasのオブジェクトであれば自動取得。

knot_kwds : None or list of dict - ノットオプション。ディフォルトではpatsyと同様の方法でノットは選択される。しかし、ノットの数は定数項のあるなしによらない。ノットの内部選択はデータの量により、patsyとmgcvと同様である。境界点はデータの範囲の極限である。提供されている選択肢は：

knots : None or array interior knots
spacing : ‘quantile’ or ‘equal’
lower_bound : None or float location of lower boundary knots,
all boundary knots are at the same point
upper_bound : None or float location of upper boundary knots,
all boundary knots are at the same point
all_knots : None or array If all knots are provided, then those
will be taken as given and all other options will be ignored.

from statsmodels.gam.api import GLMGam, BSplines
import numpy as np

x_spline = df_autos[['weight', 'hp']]

bs = BSplines(x_spline, df=[12, 10], degree=[3, 3])

GLMGam

class statsmodels.gam.generalized_additive_model.GLMGam(endog, exog=None, 
           smoother=None, alpha=0, family=None, offset=None, exposure=None, 
           missing='none', **kwargs)

一般化加法モデル

パラメータ

endog : array_like
exog : array_like or None
この説明変数は線形として扱われる。このクラスのモデルは部分的線形モデルである( partial linear model).
smoother : instance of additive smoother class
BSplineまたはCyclicCubicSpline。
alpha : list of floats　
滑らかな項の罰則付き重み. リストの長さは滑らかな項の数と同じである必要がある。
family: family class instance
ディフォルトでは正規分布。2項分布を指定するには、family = sm.family.Binomial() それぞれの分布族は引数にlink instanceをもつことができる。
offset: array-like or None
オフセットをモデルに組み込む。その場合、配列の長さは外生変数の行の数になる必要がある。
exposure: array-like or None
Log(exposure) をモデルの線形予測に加えることができる。リンク関数にlogを用いたときにのみExposureは有効になる。その際には配列の長さは内生変数と同じになる必要がある。
missing : ‘none’
欠損値処理は含まれていない。
kwargs
スーパークラスを呼ぶために使われるキーワード

Kolmogorov–Arnoldの表現定理によると、多変量関数は一変量の複数の和を用いた構造

$$f({\vec{x}})=\sum_{q=0}^{2n}\Phi_{q}\left(\sum_{p=1}^{n}\phi_{q,p}(x_{p})\right)$$

により表現できる。しかし、これは非常に複雑な構造なので、

$$f({\vec{x}})=\Phi \left(\sum_{p=1}^{n}\phi_{p}(x_{p})\right)$$

が一般に用いられる。$\Phi$は滑らかな単調な関数である。$g$は$\Phi$の逆関数で

$$g(f({\vec{x}}))=\sum_{i}f_{i}(x_{i})$$

と書くことができる。これは一般的に

$$ g(E(Y))=\beta_{0}+f_{1}(x_{1})+f_{2}(x_{2})+\cdots +f_{m}(x_{m})$$

と表現される。このときこの関数はある観測値の期待値の近似である。

また、$f_i$は

$$ f_{j}(x_{j})=\sum_{k=1}^{K_{j}}\beta_{jk}b_{jk}(x_{j})$$

と書くことができる。$b_{jk}(x_{j})$は既知の基底関数である。$K_j$は基底の次元である。$f_i$を基底関数の形式で表し、入れ替えると、$g(E(Y)$の式は一般化線形モデルとなる。しかし、このままだと過学習に陥りやすいので

$$ D(\beta )+\sum_{j}\lambda_{j}\int f_{j}^{\prime \prime }(x)^{2}dx$$

を最小化する。ここで、$\beta$はすべてのパラメータを1つにまとめたもの、D(Deviance)は飽和モデルの対数尤度とモデルの対数尤度の差を2階微分したものである。この最適化には反復再重み付け最小二乗法 (IRLS法; iteratively reweighted least squares method)が用いられる。これは再重み付けしながら，重み付2乗和誤差を最小化する。

from_formula

次の例ではfrom_formulaの形を用いてモデルを設定している。formulaは'city_mpg ~ fuel + drive'である。smootherはBスプラインである。

alpha = np.array([21833888.8, 6460.38479])
gam_bs = GLMGam.from_formula('city_mpg ~ fuel + drive', data=df_autos,
                             smoother=bs, alpha=alpha)

res_bs = gam_bs.fit()
print(res_bs.summary())

Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:               city_mpg   No. Observations:                  203
Model:                         GLMGam   Df Residuals:                   189.13
Model Family:                Gaussian   Df Model:                        12.87
Link Function:               identity   Scale:                          4.8825
Method:                         PIRLS   Log-Likelihood:                -441.81
Date:                Tue, 22 Oct 2019   Deviance:                       923.45
Time:                        13:05:20   Pearson chi2:                     923.
No. Iterations:                     3                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept       51.9923      1.997     26.034      0.000      48.078      55.906
fuel[T.gas]     -5.8099      0.727     -7.989      0.000      -7.235      -4.385
drive[T.fwd]     1.3910      0.819      1.699      0.089      -0.213       2.995
drive[T.rwd]     1.0638      0.842      1.263      0.207      -0.587       2.715
weight_s0       -3.5556      0.959     -3.707      0.000      -5.436      -1.676
weight_s1       -9.0876      1.750     -5.193      0.000     -12.518      -5.658
weight_s2      -13.0303      1.827     -7.132      0.000     -16.611      -9.450
weight_s3      -14.2641      1.854     -7.695      0.000     -17.897     -10.631
weight_s4      -15.1805      1.892     -8.024      0.000     -18.889     -11.472
weight_s5      -15.9557      1.963     -8.128      0.000     -19.803     -12.108
weight_s6      -16.6297      2.038     -8.161      0.000     -20.624     -12.636
weight_s7      -16.9928      2.045     -8.308      0.000     -21.002     -12.984
weight_s8      -19.3480      2.367     -8.174      0.000     -23.987     -14.709
weight_s9      -20.7978      2.455     -8.472      0.000     -25.609     -15.986
weight_s10     -20.8062      2.443     -8.517      0.000     -25.594     -16.018
hp_s0           -1.4473      0.558     -2.592      0.010      -2.542      -0.353
hp_s1           -3.4228      1.012     -3.381      0.001      -5.407      -1.438
hp_s2           -5.9026      1.251     -4.717      0.000      -8.355      -3.450
hp_s3           -7.2389      1.352     -5.354      0.000      -9.889      -4.589
hp_s4           -9.1052      1.384     -6.581      0.000     -11.817      -6.393
hp_s5           -9.9865      1.525     -6.547      0.000     -12.976      -6.997
hp_s6          -13.3639      2.228     -5.998      0.000     -17.731      -8.997
hp_s7          -13.8902      3.194     -4.349      0.000     -20.150      -7.630
hp_s8          -11.9752      2.556     -4.685      0.000     -16.985      -6.965
================================================================================

Link Functionのidentityは恒等関数のこと。この場合には何もしていない。

MethodのPIRLSは罰則付き反復再重み付け最小二乗法 (PIRLS法) 。

Scaleは分散で、scaleの平方根は回帰の標準誤差である。

Deviance:$ D = \sum_i freq_weights_i * var_weights * resid_dev_i / scale$

GLMGamResults.plot_partial(smooth_index, plot_se=True, cpr=False, 
                           include_constant=True, ax=None)

線形の予測に滑らかな項の寄与をプロットする。plot the contribution of a smooth term to the linear prediction

パラメータ

smooth_index : int - 滑らか項のリスト内の番号。
plot_se : bool - plot_se が trueであれば信頼区間を表示。
cpr : bool - cpr (構成要素＋残差) が trueであれば残差(the partial working residuals)を表示。
include_constant : bool - もしinclude_constantがtrueであれば,予測における推定切片と標準誤差を付加。
ax : None or matplotlib axis instance - もしaxが Noneでなければ,プロットを付加。

%matplotlib inline
res_bs.plot_partial(0, cpr=True, include_constant=True, ax=None)

res_bs.plot_partial(1, cpr=True)

次の例ではポアソン分布をfamilyとして用いる。リンク関数はlogである。

alpha = np.array([8283989284.5829611, 14628207.58927821])
gam_bs = GLMGam.from_formula('city_mpg ~ fuel + drive', data=df_autos,
                             smoother=bs, alpha=alpha,
                             family=sm.families.Poisson())
res_bs = gam_bs.fit()
print(res_bs.summary())

 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:               city_mpg   No. Observations:                  203
Model:                         GLMGam   Df Residuals:                   194.75
Model Family:                 Poisson   Df Model:                         7.25
Link Function:                    log   Scale:                          1.0000
Method:                         PIRLS   Log-Likelihood:                -530.38
Date:                Tue, 22 Oct 2019   Deviance:                       37.569
Time:                        16:09:46   Pearson chi2:                     37.4
No. Iterations:                     6                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept        3.9960      0.130     30.844      0.000       3.742       4.250
fuel[T.gas]     -0.2398      0.057     -4.222      0.000      -0.351      -0.128
drive[T.fwd]     0.0386      0.075      0.513      0.608      -0.109       0.186
drive[T.rwd]     0.0309      0.078      0.395      0.693      -0.122       0.184
weight_s0       -0.0811      0.030     -2.689      0.007      -0.140      -0.022
weight_s1       -0.1938      0.063     -3.067      0.002      -0.318      -0.070
weight_s2       -0.3160      0.082     -3.864      0.000      -0.476      -0.156
weight_s3       -0.3735      0.090     -4.160      0.000      -0.549      -0.198
weight_s4       -0.4187      0.096     -4.360      0.000      -0.607      -0.230
weight_s5       -0.4645      0.103     -4.495      0.000      -0.667      -0.262
weight_s6       -0.5092      0.112     -4.555      0.000      -0.728      -0.290
weight_s7       -0.5469      0.119     -4.598      0.000      -0.780      -0.314
weight_s8       -0.6211      0.137     -4.528      0.000      -0.890      -0.352
weight_s9       -0.6866      0.153     -4.486      0.000      -0.987      -0.387
weight_s10      -0.7370      0.174     -4.228      0.000      -1.079      -0.395
hp_s0           -0.0247      0.010     -2.378      0.017      -0.045      -0.004
hp_s1           -0.0557      0.022     -2.479      0.013      -0.100      -0.012
hp_s2           -0.1046      0.038     -2.719      0.007      -0.180      -0.029
hp_s3           -0.1438      0.050     -2.857      0.004      -0.242      -0.045
hp_s4           -0.1919      0.063     -3.047      0.002      -0.315      -0.068
hp_s5           -0.2567      0.079     -3.231      0.001      -0.412      -0.101
hp_s6           -0.4152      0.120     -3.455      0.001      -0.651      -0.180
hp_s7           -0.4889      0.152     -3.214      0.001      -0.787      -0.191
hp_s8           -0.5470      0.195     -2.810      0.005      -0.928      -0.166
================================================================================

res_bs.plot_partial(0, cpr=True, include_constant=True, ax=None)

res_bs.plot_partial(1, cpr=True, include_constant=True, ax=None)


GLMGam.select_penweight(criterion='aic', start_params=None, 
                        start_model_params=None, method='basinhopping', **fit_kwds)

基準を最小にするalphaを探索する。

gam_bs.select_penweight()[0]

array([8.28380754e+09, 1.46280830e+07])

GLMGam.select_penweight_kfold(alphas=None, cv_iterator=None, 
                              cost=None, k_folds=5, k_grid=11)

alphasを k-交叉検証を用いて探索する。

gam_bs.select_penweight_kfold()[0]

(10000000.0, 630.957344480193)

関連サイト

statsmodels:一般化線形モデル

参考