初めに
筆者は、線形回帰の説明変数について、
「目的変数に無関係な変数は学習に用いない方が良い」(※1)
というイメージを持っている。
しかし、線形回帰分析を勉強していく中で、
決定係数について下記の性質があることを知った。
「説明変数を追加した時、追加前と比べて、決定係数が下回ることは起こり得ない。(なので、説明変数の個数が異なるモデルの精度比較は、"自由度調整済み決定係数"を用いる。)」(※2)
その性質を知った筆者は、下記[概要]記載の検証を行ってみたくなった。
概要
筆者の想定(※1)と決定係数の性質(※2)に関して、下記2点(①、②)の検証を (実際にモデルを作って) 行ってみた。
本記事は、その検証方法・結果を記載したものである。
※①は理論的に成り立つものであり、②は筆者の想定である。
①決定係数の比較
→目的変数に無関係な説明変数を追加した時、その決定係数は追加前の決定係数を下回ることがない。
②自由度調整済み決定係数の比較
→自由度調整済み決定係数の場合だと、目的変数に無関係な変数を"追加する前"の方が値が高いのではないか？
前提
[線形回帰]
既存の説明変数と目的変数に対して、目的変数を"説明変数の一次結合(線型結合)"で表現する手法。線形結合で表現するため、表現した結果(推定値)は直線のグラフになる。
[決定係数]
作成した線形回帰モデルが、"データに対して当てはまりがどれ程良いか"を示す指標であり、下記の計算式で表される。
$ R^{2}=1-\dfrac {\sum e^{2}}{\sum \left( y-\overline {y}\right) ^{2}} $
$ {\sum e^{2}}:残差変動 (目的変数と推定値の差の二乗和) $
$ {\sum \left( y-\overline {y}\right) ^{2}}:全変動 (目的変数とその平均の差の二乗和) $
[自由度調整済み決定係数]
決定係数に対して、説明変数の個数(k)を加味した指標。決定係数と同様に"データに対して当てはまりがどれ程良いか"を示す。
$ {\widehat {R}^{2}}=1-\dfrac {n-1}{n-k}\left( 1-R^{2}\right) $
$ {n}:データ件数 $
$ {k}:説明変数の個数 $
検証詳細
上記[概要]記載の検証2点(①、②)について、それぞれ下記の方法で検証する。
[用意するデータ]
(1)説明変数(目的変数に関係あり)
→乱数を用いて作成。目的変数は本変数の一次結合で作成する。
(2)説明変数(目的変数に関係無し)
→乱数を用いて作成。
(3)目的変数
→上記(1)の説明変数の一次結合で作成する。
①決定係数の比較
→下記 < モデル1 > と < モデル2 > の決定係数の差(モデル2 - モデル1)を確認する。
< モデル1 > (n 個)
上記変数(1)の変数(1〜n個)で作成したモデル
< モデル2 > (n × m 個)
上記 [モデル1] 1個に対して、
"そのモデルに用いた説明変数"と"上記変数(2)(1〜m個)"で作成したモデル
②自由度調整済み決定係数の比較
→上記検証①に対して、"決定係数"を"自由度調整済み決定係数"にして確認する。
検証手順
※下記[検証用コード]の流れを手順としてまとめたものである。
[1]データ作成
→上記変数(1)〜(3)を持つデータフレーム all_data を作成。
(1)説明変数(目的変数に関係あり)
→ col_rel_[n]（n=0〜10）
(2)説明変数(目的変数に関係無し)
→ col_no_rel_[m]（m=0〜10）
(3)目的変数
→ target
[2]モデル作成
→上記 < モデル1 > と < モデル2 > を作成。
※dict型 dict_models で (valueとして) 保持し、keyは下記の通り。
< モデル1 > (n 個)
→ nothing_ [n] _ [m]
※n、m は学習に用いた説明変数の個数。
n は変数(1)の個数、m は変数(2)の個数
< モデル2 > (n × m 個)
→ contain_ [n] _ [m]
※n、m は上記 < モデル1 > と同じ意味。
[3]決定係数、自由度調整済み決定係数を算出
→両方の値をそれぞれdict型で(valueとして)保持。
各dict名は下記の通り。
決定係数：dict_R2
自由度調整済み決定係数：dict_adjusted_R2
[4]集計結果をデータフレーム1つにまとめる
[5]折れ線グラフ2つ作成
→縦軸は(自由度調整済み)決定係数、
横軸は"モデル学習に使用した上記変数(2)の数"。
検証用コード
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
#【パラメータ】
#シード値
random.seed(10)
#データ件数
data_num = 10 ** 4
#項目数 (rel:目的変数に関連あり、no_rel:目的変数に関連無し)
group_num = 10
dict_group_num = dict()
dict_group_num["rel"] = group_num #目的変数に関連あり
dict_group_num["no_rel"] = group_num #目的変数に関連無し
#【データ作成】
#説明変数
col_format = 'col_{group}_{index}'
all_data = pd.DataFrame({
col_format.format(
group = key
, index = i
):[random.randint(0 , 1000) for _ in range(10 ** 4)]
for key , val in dict_group_num.items()
for i in range(val)
})
#目的変数
w_list = [random.randint(1 , 10) for i in range(dict_group_num["rel"])]
target_series = pd.Series(np.zeros(data_num))
for i in range(dict_group_num["rel"]):
w = w_list[i]
col = col_format.format(
group = "rel"
, index = i
)
#確認
print('-' * 50)
print(target_series.head())
print("w=" , w)
print(all_data[col].head())
print(w * all_data[col].head())
target_series = target_series + w * all_data[col]
#確認
display(target_series.head())
all_data['target'] = target_series
#【モデル作成】
#説明変数の全組み合わせ
dict_features = {}
key_format = "{type_}_{n}_{m}"
#1.「col_rel_x」(n)
type_ = "nothing"
for n in range(dict_group_num["rel"]):
cols = [col_format.format(
group = "rel"
, index = i
) for i in range(n+1)]
key = key_format.format(
type_ = type_
, n = n + 1
, m = 0
)
dict_features[key] = cols
#2.「col_rel_x」(n)と「col_no_rel_x」(m)
type_ = "contain"
for n in range(dict_group_num["rel"]):
cols_rel = [col_format.format(
group = "rel"
, index = i
) for i in range(n+1)]
for m in range(dict_group_num["no_rel"]):
cols_no_rel = [col_format.format(
group = "no_rel"
, index = i
) for i in range(m+1)]
cols = cols_rel + cols_no_rel
key = key_format.format(
type_ = type_
, n = n + 1
, m = m + 1
)
dict_features[key] = cols
#確認
type_ = "nothing"
print("-" * 50 , type_)
_dict = {key:val for key , val in dict_features.items() if key[:len(type_)] == type_}
print(list(_dict.items())[:5])
type_ = "contain"
print("-" * 50 , type_)
_dict = {key:val for key , val in dict_features.items() if key[:len(type_)] == type_}
print(list(_dict.items())[:5])
#モデル作成
dict_models = {}
for key , feature_cols in dict_features.items():
#説明変数と目的変数に分割する
train_X = all_data[feature_cols]
train_y = all_data['target']
#モデル作成
model = LinearRegression()
model.fit(train_X , train_y)
dict_models[key] = model
#確認
print("-" * 50)
print("key={}".format(key))
print(model.intercept_)
print(model.coef_)
#確認
list(dict_models.keys())[:5]
#決定係数 、自由度調整済み決定係数
dict_R2 = {}
dict_adjusted_R2 = {}
for key , feature_cols in dict_models.items():
#モデル取得
model = dict_models[key]
#全データを抽出(numpy)
feature_cols = dict_features[key]
X = np.array(all_data[feature_cols])
y = np.array(all_data["target"])
#決定係数
R2 = model.score(X , y)
dict_R2[key] = R2
#自由度調整済み決定係数
n = data_num
k = int(key.split("_")[1]) + int(key.split("_")[2])
adjusted_R2 = 1 - ((n-1)/(n-k)) * (1 - R2)
dict_adjusted_R2[key] = adjusted_R2
#【検証結果】
R2_df = pd.DataFrame({
"key":list(dict_R2.keys())
, "R2":list(dict_R2.values())
})
adjusted_R2_df = pd.DataFrame({
"key":list(dict_adjusted_R2.keys())
, "adjusted_R2":list(dict_adjusted_R2.values())
})
result_df = pd.merge(R2_df , adjusted_R2_df , on="key" , how='outer')
result_df['rel_num'] = result_df["key"].str.split("_" , expand=True)[1].astype(int)
result_df['no_rel_num'] = result_df["key"].str.split("_" , expand=True)[2].astype(int)
#確認
print(len(R2_df))
print(len(adjusted_R2_df))
print(len(result_df))
display(R2_df.head())
display(adjusted_R2_df.head())
display(result_df.head())
#【グラフ作成】
# 決定係数
value = "R2"
fig = plt.figure(figsize=(10,10))
for i in range(dict_group_num["rel"]):
rel_num = i + 1
df = result_df.query("rel_num == {}".format(rel_num))
base = df.query("no_rel_num == 0")[value]
x = df["no_rel_num"]
y = df[value].apply(lambda x:x - base)
#確認
# print("-" * 50)
# print("base={}".format(base))
# print("x={}".format(x))
# print("y={}".format(y))
plt.plot(x, y, marker = 'o'
, label="rel_num={}".format(rel_num))
plt.title("Diff of {}".format(value))
plt.legend(loc="upper left" , bbox_to_anchor=(1, 1))
plt.xlabel("no_rel_num")
plt.ylabel(value)
plt.grid()
plt.show()
#グラフ保存
fig.savefig("plot_R2.png")
# 自由度調整済み決定係数
value = "adjusted_R2"
fig = plt.figure(figsize=(10,10))
for i in range(dict_group_num["rel"]):
rel_num = i + 1
df = result_df.query("rel_num == {}".format(rel_num))
base = df.query("no_rel_num == 0")[value]
x = df["no_rel_num"]
y = df[value].apply(lambda x:x - base)
#確認
# print("-" * 50)
# print("base={}".format(base))
# print("x={}".format(x))
# print("y={}".format(y))
plt.plot(x, y, marker = 'o'
, label="rel_num={}".format(rel_num))
plt.title("Diff of {}".format(value))
plt.legend(loc="upper left" , bbox_to_anchor=(1, 1))
plt.xlabel("no_rel_num")
plt.ylabel(value)
plt.grid()
plt.show()
#グラフ保存
fig.savefig("plot_adjusted_R2.png")
検証結果
[ グラフ1 ] 決定係数
"目的変数に関係あり"の説明変数の個数(rel_num)に関わらず、追加する"目的変数に関係無し"の説明変数が増えれば増える程、決定係数は単調増加している。
つまり、理論通りであることが確認できた。
[ グラフ2 ] 自由度調整済み決定係数
"目的変数に関係あり"の説明変数の個数について、
1〜4個の時は、"目的変数に関係無し"の説明変数を追加した方が良くなっている。
(1〜4個だけでは説明不十分で、少しでも情報を追加した方が良いという感じだろう。)
9〜10個の時は、ほとんど変わらない。
(9〜10個あれば十分に目的変数を説明できるという感じだろう。)
それに対して、5〜8個の場合は、"目的変数に関係無し"の説明変数を追加しない方が決定係数が高くなっている。
筆者の想定通りであることが確認できた。
まとめ
決定係数を用いてモデルを比較する時は、学習に使用する説明変数の個数に注意しなければいけない。
検証結果からも確認できたが、決定係数では、説明変数が"他方の分を包含して個数が多い"場合、(使用するに相応しくない変数が追加されているケースでも)決定係数は高くなる、又は、等しい。
なので、その場合に対して、より適切な比較をするには、"自由度調整済み決定係数"を用いた方が良さそうである。