はじめに

機械学習を現実の問題に適用する場合、そのモデルに説明性が求められることが少なからず存在すると思います。

その場合、精度を犠牲にして線形回帰を実施するでしょうか？木系モデルの重要度を頑張って説明するでしょうか？それともSHAPやLIMEなど線形近似モデルを利用するでしょうか？（まあ銀の弾丸はないんですが）

今回実験を行う一般化加法モデル（GAM）は、線形モデルの利点（説明性）を保ちつつ精度を高められるモデルであるといわれているもので、実際のところどれくらいの感じになるか確認するための実験を行いました。

参考文献

GAMについて

GAMは式としては以下のようになります。

$y=\alpha+f_{1}\left(x_{1}\right)+f_{2}\left(x_{2}\right)+\cdots+f_{n}\left(x_{n}\right)+\varepsilon$

fが全て線形な場合には、普通の線形回帰と同じ式になります。それぞれの変数に対して関数をおくような形になり、非線形な関係を表せるようになるため、線形回帰より精度が良くなることが多いです。

これら１つ１つの関数をどのように定めているかという問題ですが、これは平滑化スプラインで曲線近似を行っています。この部分は参考１の記事を参照するとわかりやすいです。平滑化スプラインを使う場合、曲げの度合いを強めればいくらでも学習データにフィッティングさせることができてしまうのですが、そこは一般化クロスバリデーション（GCV）という方法で過学習を防いでいるらしいです。

有名どころとしては、時系列解析ツールのProphetでもGAMは用いられています。

実験

今回使用するデータはsklearnに入っているBoston house-pricesです。これはボストン市郊外の地域別の住宅価格とデータセットになります。少し変数を絞り、以下の6つの変数から住宅価格を推定するモデルを作成します。比較のため、線形回帰、lightgbm、GAMの３つの手法を用います。

カラム名	説明
CRIM	1人当たりの犯罪率
ZN	25,000 平方フィート以上の住居区画の占める割合
RM	住居の平均部屋数
AGE	1940 年より前に建てられた物件の割合
DIS	5 つのボストン市の雇用施設からの重み付け距離
LSTAT	低所得者の割合

説明重視のモデル代表として線形回帰を実施します。

from sklearn.datasets import load_boston
boston = load_boston()
df=pd.DataFrame(boston.data, columns=boston.feature_names)

#使用する変数を選択、正規化
X = df[['CRIM', 'ZN','RM', 'AGE', 'DIS', 'LSTAT']]
ss = preprocessing.StandardScaler()
X_ss = pd.DataFrame(ss.fit_transform(X))
X_ss.columns=X.columns
y = boston.target

#trainとtestに分ける
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

#線形回帰の実装（切片あり）
import statsmodels.api as sm
Xols=X.copy()
Xols["coef"]=1
mod = sm.OLS(y, Xols)
res = mod.fit()

実行結果は以下のようになります。

f:id:rmizutaa:20190323091549p:plain

CRIM,AGE,DIS,LSTATはマイナス、ZNとRMはプラスにはたらくようです。 testデータのRMSEは5.174でした。

次に精度重視のモデル代表としてlightgbmを実施します。

#lightgbmの実装
lgbX_train, lgbX_valid, lgby_train, lgby_valid = train_test_split(X_train, y_train, test_size=0.2, random_state=42)
params={
    'random_state' : 1,
    'objective': 'rmse',
    }
dtrain = lgb.Dataset(lgbX_train, label=lgby_train)
dvalid = lgb.Dataset(lgbX_valid, label=lgby_valid)
bst = lgb.train(params, dtrain, num_boost_round=1000,valid_sets=[dtrain,dvalid],early_stopping_rounds=50,verbose_eval=10)

変数重要度は以下のようになりました。

f:id:rmizutaa:20190323091433p:plain

LSTATの影響がもっとも大きくZNは結果与える影響が小さいことがわかります。プラス、マイナスのどちらに寄与するかはわかりません。また、テストデータのRMSEは3.950でした。

最後にGAMを実施します。 pythonだとpyGAMというパッケージがあるのでそれを利用します。 Rだとmgcvというパッケージが使えるみたいです。

from pygam import LinearGAM
gam = LinearGAM().fit(X_train, y_train)

GAMの結果は図示しないとよくわからないのでプロットします。

#それぞれのf(x)を図示
plt.figure(figsize=(12,8))
for i in range(X.shape[1]):
    plt.subplot(2,3,i+1)
    XX = gam.generate_X_grid(term=i)
    plt.plot(XX[:, i],gam.partial_dependence(term=i, X=XX))
    plt.plot(XX[:, i], gam.partial_dependence(term=i, X=XX, width=.95)[1],
                c='r', ls='--')
    plt.title(X.columns[i])
plt.show()

f:id:rmizutaa:20190323091506p:plain

線形モデルとは違ってグネグネした線になります。 DISとLSTATは線形に住宅価格に影響するわけではなく、ある一定の値までは住宅価格と反比例しますが、ある値を越えると価格に与える影響は小さくなるようです。また、RMは0以下の場合には住宅価格に与える影響が小さいですが、 0以上になると値に比例して住宅価格が上がるようです。

テストデータのRMSEは4.116でした。カラムがあまり多くないせいか精度はlightgbmの値(3.950)に近い結果になりました。また各変数が独立しているので予測値の値がなぜこうなったのか？という部分がブラックボックス化することもないかなと思います。

ただ、RM（部屋数）が多すぎると逆に住宅価格が下がるという不自然な結果にもなっている部分もあります。これは平滑化スプラインの外れ値に弱い性質に起因していると考えられます。下記はRMと住宅価格の関係性の散布図ですが、一番右端に部屋数が多いのに価格が低い点が１点あり、その影響を受けている可能性が高い気がします。

f:id:rmizutaa:20190323100343p:plain