Pythonを使って自由度調整済み決定係数の大きい説明変数の組み合わせを探索する①(全探索法)

#作業環境
Jupyter Notebook(6.1.4)を用いて作業を進めた。
各versionはpandas(1.1.3), scikit-learn(0.23.2), statsmodels(0.12.0)である。

#やりたいこと
最小二乗法で最適化された線形回帰モデルについて学習し、そのモデルを評価する指標として決定係数や自由度調整済み決定係数を学んだ。学習時には手作業で係数を追加したり削除したりしながら決定係数を比較していたが、これを自動探索するプログラムを練習として組んでみたいと考えた。
つまりやりたいことは、データ全体を渡したら最も決定係数が高くなる説明変数の組み合わせとその決定係数を返してくれるプログラムを考えて作ってみたい。
なお、その説明変数が説明変数として適切かどうかは特に考慮せず、単純に決定係数を比較して説明変数を選択するのが今回の目的である。

今回はまず全探索の方法を考える。これはすべての組み合わせを比べて決定係数が最大のものを選ぶので、本当に最も係数が大きいものを選ぶことができる。しかし、単純に$2^{説明変数の個数}$通りの計算が起こるのでとてもコストの高い探索方法になる。

#探索の方針を考える

以前の記事で決定係数や自由度調整済み決定係数をPythonで求める自由度調整済み決定係数をPythonで求める方法を取り扱った。今回は自由度調整済み決定係数を指標として探索を行いたい。このときsklearn.linear_modelのLinearRegressionを使って求める方法とstatsmodelsを使って求める方法の２通りを扱ったため、探索も２通りの書き方を考える。

パラメータとしては説明変数のデータセット$X$と目的変数のデータセット$y$を渡すほかに、単純な線形回帰モデル以外でも計算できるようにできるようにしたい。また、返り値は最も大きい自由度調整済み決定係数と、そのときに選択した説明変数であるようにする。

大まかな方針としては、説明変数を選んで回帰モデルを作成、自由度調整済み決定係数を求め、もし今までの自由度調整済み決定係数よりも大きい値ならその値と説明変数を保存する、というようにしていく。説明変数の選び方については、繰り返し回数を説明変数の個数の長さをもつ二進数へと変換し、1をフラグにして選択を行う方法をとる。

#検証に使用するデータ
scikit-learnのdatasetsの中のボストン住宅データ(説明変数13個)を使って検証していく。以下のように説明変数データセットと目的変数データセットを作成した。

import pandas as pd
from sklearn.datasets import load_boston
X = pd.DataFrame(load_boston().data, columns=load_boston().feature_names)
y = pd.DataFrame(load_boston().target, columns=['target'])

#全探索(sklearn.linear_modelを使う方法)

まずはsklearn.linear_modelで自由度調整済み決定係数を求める方法で全探索を実装する。前回作った関数をまずは実装する。

from sklearn.metrics import r2_score
def adj_r2_score(y_true, y_pred, p=X.shape[1]):
    return 1-(1-r2_score(y_true, y_pred)) * (len(y_true)-1) / (len(y_true) - p - 1)

今回は説明変数の個数が変化するので第三引数に注意したい。

探索は次のように実装した。

def R2search1(X, y, model):
    """自動調整済み決定係数が最も大きい説明変数を探索する
    パラメータ
    -----------
    X : 説明変数
    y : 目的変数
    model　:インスタンス化したモデル
    
    戻り値
    -----------
    best_r2, best_col :　最も大きい自由度調整済み決定係数とそれを与える説明変数のリスト
    """
    #説明変数名とその個数を取得
    x_col = X.columns
    n = len(x_col)
    #最大値を保存するオブジェクトの作成
    best_col = []
    best_r2 = 0
    for i in range(1, 2 ** n):
        #iを二進数へ変換し、文字数をnに直す
        #例）5 --> 0000000000101
        flg = str(format(i, 'b').zfill(n))
        #フラグに従って特徴量の選択
        col = []
        for j in range(n):
            if flg[j] == '1':
                col.append(x_col[j])
        #回帰モデルを作成し、自由度調整済み決定係数を求める
        X_train = X[col]
        model.fit(X_train, y)
        y_pred = model.predict(X_train)
        r2 = adj_r2_score(y, y_pred, p=len(col))
        #r2が今までのものより大きいのであれば更新
        if r2 >= best_r2:
            best_r2 = r2
            best_col = col.copy()
    return best_r2, best_col

検証データで動作を確かめてみたところ、次のような結果が得られた。

from sklearn.linear_model import LinearRegression
model = LinearRegression()
print(R2search1(X, y, model))
#---> (0.7348057723274566, ['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT'])

from sklearn.linear_model import Ridge
model = Ridge(alpha=1)
print(R2search1(X, y, model))
#---> (0.7331269601298093, ['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT'])

#全探索(statsmodels.apiを使う方法)

次に、statsmodels.apiで自由度調整済み決定係数を求める方法で全探索を実装する。これはモデルを与えられないので線形回帰モデルのみの扱いである。

def R2search2(X, y):
    """自動調整済み決定係数が最も大きい説明変数を探索する
    パラメータ
    -----------
    X : 説明変数
    y : 目的変数
    
    戻り値
    -----------
    best_r2, best_col :　最も大きい自由度調整済み決定係数とそれを与える説明変数のリスト
    """
    #説明変数名とその個数を取得
    x_col = X.columns
    n = len(x_col)
    #最大値を保存するオブジェクトの作成
    best_col = []
    best_r2 = 0
    for i in range(2 ** n - 1):
        #iを二進数へ変換し、文字数をnに直す
        a = str(format(i+1, 'b').zfill(n))
        #フラグに従って特徴量の選択
        col = []
        for j in range(n):
            if a[j] == '1':
                col.append(x_col[j])
        #自由度調整済み決定係数を求める
        X_train = X[col]
        #定数項の追加
        X_train = sm.add_constant(X_train)
        model = sm.OLS(y, X_train)
        result = model.fit()
        r2 = result.rsquared_adj
        #r2が今までのものより大きいのであれば更新
        if r2 >= best_r2:
            best_r2 = r2
            best_col = col.copy()
    return best_r2, best_col

こちらも検証データを使って動作を確かめてみたところ、次の結果が得られた。

print(R2search2(X, y))
#---> (0.7348057723274566, ['CRIM', 'ZN', 'CHAS', 'NOX', 'RM', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT'])

sklearn.linear_modelを使う全探索と同じ結果を得ることができた。

#実行時間の比較
%%timeitを使用して実行時間を確認すると以下の通りであった。

%%timeit
print(R2search1(X, y, model))
#---> 1min 5s ± 4.66 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

%%timeit
print(R2search2(X, y))
#---> 39.2 s ± 3.62 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

これを見るとsklearn.linear_modelはstatsmodelsの倍近い実行時間がかかっていることがわかる。線形回帰モデルのみで考えるのであればstatsmodelsを使用した方がより早く実行できるのでいいのではないだろうか。

#課題
次のような課題が考えられる。

• やはり実行に時間がかかってしまう(特徴量が増えるととても時間がかかってしまう)
• 自由度調整済み決定係数が最大のものしか取得できない
• statsmodelsの方は別の回帰モデルを選択できない

ひとまず一番上の課題を解決できるように、全探索でない方法を今度考えていこうと思う。
→増加法・減少法を実装してみた。

二番目については更新の仕方を工夫することで解決できそうである。
三番目については変数をいじることで解決する方法があるかもしれない。