pythonでstepwise regression用の関数を作っておく

はじめに

ほんと、久々の更新になってしまいました。。。
いまだに月間で1000PVほど見られているようでとてもありがたく思いますm(_ _)m

最近も変わらず因果推論の研究を中心に行っておりますが、それ関連の内容はまた機会をみてblogで書いていければと思っています。
また先日、twitterで公開したこちらのスライドもたくさんの方に見ていただけまして、コメントも頂けたりして、とても嬉しく、励みになっています。
speakerdeck.com

また、少しずつではありますが更新いたしますので、たまに覗いていただければ嬉しいです。
では、本題にまいります。

今回の更新

とはいっても、今日の更新は、大した内容ではなく、pythonでstepwise regressionの関数で自分がほしいものがないので、つくりましたという内容です。

Stepwise Regressionについて

特に、回帰モデルでモデルに含める変数をある基準に基づいて選択するという方法です。
有名なのは、AICを用いるものですが、今回はBIC（Bayesian Information Criterion）で実装します。

pythonのstepwise regressionについて

一応、stepwise regression用のモジュールがあるのですが、p-valueをみて変数を選択していて、ちょっと「うーん？」という感じ。
あと、調べると他にもあるのですが、statmodelsに依存した形のものが多かったので、できればそのあたりも依存しないで動かしたくて、結局自作に至りました😣
このあたりは、Rはほんとしっかりしていて、step functionなんかを使うとサクッとできるあたりは嬉しいところですね。

本題 - 関数 -

とりあえずAICとBICを基準に選択できればOKということと、この関数を別の関数の中で使用するために、返り値は回帰係数であってほしいという個人的な理由が存分に入っております！
Rっぽい出力はまた時間があれば作成しますが、やる気が風前の灯なのでこれくらいで妥協いたします...ご勘弁ください！！

# numpy
import numpy as np

#Akaike Information Criterion
def AIC_gaussian(x,y):
    y_pred = np.dot(x,np.linalg.solve(np.dot(x.T,x),np.dot(x.T,y)))
    resid = y - y_pred #residuals
    rss = np.dot(resid.T,resid) #residual sum of square
    AIC = len(y)*np.log(rss/len(y)) + (x.shape[1])*2
    return AIC

#Bayesian Information Criterion
def BIC_gaussian(x,y):
    y_pred = np.dot(x,np.linalg.solve(np.dot(x.T,x),np.dot(x.T,y)))
    resid = y - y_pred #residuals
    rss = np.dot(resid.T,resid) #residual sum of square
    BIC = len(y)*np.log(rss/len(y)) + (x.shape[1])*np.log(len(y))
    return BIC

def forward_stepwise_result(x,y):
    beta = np.linalg.solve(np.dot(x.T,x),np.dot(x.T,y))
    return beta

# stepwise selection
def forward_stepwise(x,y,method='BIC',verbose=True):
    x_tmp = np.concatenate([np.reshape(np.ones(1000),(1000,1)),x],axis=1) #xにinterceptを追加
    initial_list = [0]
    included = list(initial_list)
    while True:
        changed=False
        excluded = list(set(np.arange(x_tmp.shape[1]))-set(included))
        BIC_res = np.zeros(len(excluded))
        base_BIC = BIC_gaussian(x_tmp[:,included],y)
        print('Baseline model with BIC: {:}'.format(base_BIC))
        j = 0
        for new_column in excluded:
            BIC_res[j] = BIC_gaussian(x_tmp[:,included+[new_column]],y)
            print('Add :{:30} with BIC: {:}'.format(excluded[j],BIC_res[j]))
            j = j+1
        if BIC_res.min() < base_BIC:
            best_feature = excluded[BIC_res.argmin()]
            included.append(best_feature)
            changed=True
            if verbose:
                print('Finally Add  {:30} with BIC {:}'.format(best_feature, BIC_res.min()))

        if not changed:
            if verbose:
                print('Any variables does not added, stop forward stepwise regression')
            break
            #end while
    
    #final result の回帰係数を返す
    beta = np.reshape(np.zeros(x_tmp.shape[1]),(x_tmp.shape[1],1))
    beta[included] = forward_stepwise_result(x_tmp[:,included],y)
    
    return beta

#結果
forward_stepwise(x,y,method='BIC',verbose=True)

#Baseline model with BIC: [[7214.43429935]]
#Add :                             1 with BIC: 6564.464745289761
#Add :                             2 with BIC: 7062.695910602746
#Add :                             3 with BIC: 7220.4285978234575
#Add :                             4 with BIC: 6830.525994981331
#Finally Add                               1 with BIC 6564.464745289761
#Baseline model with BIC: [[6564.46474529]]
#Add :                             2 with BIC: 6260.368109310093
#Add :                             3 with BIC: 6571.08187071489
#Add :                             4 with BIC: 5708.315703761152
#Finally Add                               4 with BIC 5708.315703761152
#Baseline model with BIC: [[5708.31570376]]
#Add :                             2 with BIC: 5713.846539819814
#Add :                             3 with BIC: 5710.261885909516
#Any variables does not added, stop forward stepwise regression
#array([[21.86610567],
#       [42.05186786],
#       [ 0.        ],
#       [ 0.        ],
#       [ 0.35454755]])