[線形回帰]説明変数の個数と(自由度調整済み)決定係数について

# 初めに
筆者は、線形回帰の説明変数について、
「目的変数に無関係な変数は学習に用いない方が良い」(※1)
というイメージを持っている。
しかし、線形回帰分析を勉強していく中で、
決定係数について下記の性質があることを知った。

「説明変数を追加した時、追加前と比べて、決定係数が下回ることは起こり得ない。(なので、説明変数の個数が異なるモデルの精度比較は、"自由度調整済み決定係数"を用いる。)」(※2)

その性質を知った筆者は、下記[概要]記載の検証を行ってみたくなった。

# 概要
筆者の想定(※1)と決定係数の性質(※2)に関して、下記2点(①、②)の検証を (実際にモデルを作って) 行ってみた。
本記事は、その検証方法・結果を記載したものである。
※①は理論的に成り立つものであり、②は筆者の想定である。

①決定係数の比較
→目的変数に無関係な説明変数を追加した時、その決定係数は追加前の決定係数を下回ることがない。

②自由度調整済み決定係数の比較
→自由度調整済み決定係数の場合だと、目的変数に無関係な変数を"追加する前"の方が値が高いのではないか？

# 前提
[線形回帰]
　既存の説明変数と目的変数に対して、目的変数を"説明変数の一次結合(線型結合)"で表現する手法。線形結合で表現するため、表現した結果(推定値)は直線のグラフになる。
![Normdist_regression.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/618640/e30a34b1-f0a3-a8cc-8431-9d1473d3f920.png)

[決定係数]
　作成した線形回帰モデルが、"データに対して当てはまりがどれ程良いか"を示す指標であり、下記の計算式で表される。
　$ 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)の数"。

# 検証用コード
```python:experiment.jpynb
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)に関わらず、追加する"目的変数に関係無し"の説明変数が増えれば増える程、決定係数は単調増加している。
つまり、理論通りであることが確認できた。
![plot_R2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/618640/01a4461a-2dd4-1655-12f2-764d619dc4f6.png)
[ グラフ2 ] 自由度調整済み決定係数
"目的変数に関係あり"の説明変数の個数について、
1〜4個の時は、"目的変数に関係無し"の説明変数を追加した方が良くなっている。
(1〜4個だけでは説明不十分で、少しでも情報を追加した方が良いという感じだろう。)
9〜10個の時は、ほとんど変わらない。
(9〜10個あれば十分に目的変数を説明できるという感じだろう。)
それに対して、5〜8個の場合は、"目的変数に関係無し"の説明変数を追加しない方が決定係数が高くなっている。
筆者の想定通りであることが確認できた。
![plot_adjusted_R2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/618640/618f4107-3347-9dbe-f8c4-9b7f84ec3545.png)

# まとめ
決定係数を用いてモデルを比較する時は、学習に使用する説明変数の個数に注意しなければいけない。
検証結果からも確認できたが、決定係数では、説明変数が"他方の分を包含して個数が多い"場合、(使用するに相応しくない変数が追加されているケースでも)決定係数は高くなる、又は、等しい。
なので、その場合に対して、より適切な比較をするには、"自由度調整済み決定係数"を用いた方が良さそうである。