自由度調整済み決定係数をPythonで求める

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

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

その前準備として、決定係数や自由度調整済み決定係数について、その定義とPythonでの求め方をまとめる。

#決定係数、自由度調整済み決定係数とは
**決定係数$R^2$**は以下の式で定義される。

$$
R^2 = \frac{\sum(特徴量の効果)^2}{\sum(偏差)^2} = 1-\frac{\sum(回帰残差)^2}{\sum(偏差)^2}
$$

ここで、偏差とは目的変数の平均と各データ点の差であり、回帰残差とは回帰モデルの予測値と各データ点との差である。つまり、選んだ特徴量で作った回帰モデルがどれくらいよく目的変数を説明できているか、という解釈ができる。

さらに、**自由度調整済み決定係数${R^*}^2$**は以下の式で定義される。

$$
{R^*}^2 = 1-\frac{\frac{\sum(回帰残差)^2}{n-p-1}}{\frac{\sum(偏差)^2}{n-1}}
$$

ここで、$n$はデータ点の個数、$p$は説明変数の個数である。説明変数を増やすと必然的に決定係数は増加してしまうため、その影響を調整した指標である。簡単な式変形から$R^2$と${R^*}^2$の関係は以下の式となることがわかる。

$$
R^{*2}= 1-(1-R^2)\frac{n-1}{n-p-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'])

#Pythonで自由度調整済み決定係数を求める
今回はsklearn.linear_modelとstatsmodels.apiを使用する方法を使用する方法の２通りで試す。
##sklearn.linear_modelを使用して求める
まずは線形回帰モデルを作成して決定係数を求める

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X, y)
model.score(X, y)
#---> 0.7406426641094094

sklearn.metricsのr2_scoreを使った方が決定係数をより早く求められる。

from sklearn.metrics import r2_score
y_pred = model.predict(X)
r2_score(y, y_pred)
#---> 0.7406426641094094

自由度調整済み決定係数については式を定義して求める。

#自由度調整済み決定係数を計算する関数を定義する
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)

adj_r2_score(y, y_pred)
#---> 0.7337897263724629

##statsmodels.apiを使用して求める
statsmodels.apiを使用すると簡単に決定係数や自由度調整済み決定係数を求めることができる。

import statsmodels.api as sm
#定数項の追加
X1 = sm.add_constant(X)
m = sm.OLS(y, X1)
result = m.fit()
#決定係数の表示
print(result.rsquared)
#---> 0.7406426641094095
#自由度調整済み決定係数の表示
print(result.rsquared_adj)
#---> 0.733789726372463

LinearRegressionではデフォルトとして定数項のある回帰モデルが作成されるが、statsmodelsではデフォルトでは定数項のない回帰モデルが作成される。そのためモデル作成前に定数項の追加を行っている。

2通りの結果を見比べるとほとんど同じ結果を得られることが確認できた。