More than 1 year has passed since last update.

外れ値の影響で回帰直線の結論が真逆になる話

Last updated at 2023-04-21Posted at 2023-04-21

ニュースとかでよくありますよね。散布図が描いてあって、回帰直線が描いてあって、「相関がある」とか「○○なほど〇〇であることが分かった」とか解説してるやつ。

ここでは、データの一部（ここでは「外れ値」と呼ぶことにします）を取り除くだけで「正の相関」「負の相関」が逆転するケースを、人為的に作ってみたいと思います。そしてそれが、データサイズを変えるとどう影響するのかをみてみたいと思います。

Optuna で最適化

はい、いつものように（？）Optunaで最適化したいと思います。

import numpy as np
from sklearn.linear_model import LinearRegression

# Optuna の目的関数
class Objective:
    def __init__(self, n_data=10, n_outlier=1):
        self.X = []
        self.Y = []
        self.n_data = n_data # データサイズ
        self.n_outlier = n_outlier # そのうち外れ値の数
        self.best_score = None
        self.best_X = None
        self.best_Y = None

    def __call__(self, trial):
        self.X = []
        self.Y = []

        # データ点の座標を決める
        for i in range(self.n_data):
            self.X.append(trial.suggest_float("x{}".format(i), 0, 1))
            self.Y.append(trial.suggest_float("y{}".format(i), 0, 1))

        # 全てのデータを用いて線形回帰
        model1 = LinearRegression()
        model1.fit(np.array(self.X).reshape(-1, 1), self.Y)

        # 外れ値を除いたデータで線形回帰
        model2 = LinearRegression()
        model2.fit(
            np.array(self.X[: -self.n_outlier]).reshape(-1, 1),
            self.Y[: -self.n_outlier],
        )

        # ２つの線形回帰の傾きが正負逆で、かつ差が大きい時にスコアが小さくなる
        score = (
            abs(model1.coef_[0] - model2.coef_[0]) * model1.coef_[0] * model2.coef_[0]
        )

        # 歴代ベストスコアを下回った時に、新しいベストスコアとする。
        if self.best_score is None or self.best_score > score:
            self.best_score = score
            self.best_X = self.X
            self.best_Y = self.Y

            self.corrcoef1 = np.corrcoef(self.X, self.Y)[0, 1]
            self.corrcoef2 = np.corrcoef(
                self.X[: -self.n_outlier], self.Y[: -self.n_outlier]
            )[0, 1]

        # このスコアを Optuna の評価値とする。
        return score

上記のように目的関数を設計すれば、次のように、データサイズと外れ値の数をセットできます。

objective = Objective(n_data=10, n_outlier=1)

次のようにしてスタディを立ち上げてみましょう。

import optuna

study = optuna.create_study()

試行回数はとりあえず 100 回とします。

n_trials = 100

次のようにして最適化計算を実行します。

optuna.logging.set_verbosity(optuna.logging.WARN)  # 途中経過メッセージを省略する
study.optimize(objective, n_trials=n_trials, show_progress_bar=True)

外れ値の影響で回帰直線の結論が真逆になる例

計算の結果得られた最適解を図示してみましょう。

import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

plt.scatter(
    objective.best_X[: -objective.n_outlier],
    objective.best_Y[: -objective.n_outlier],
    alpha=0.8,
)
plt.scatter(
    objective.best_X[-objective.n_outlier :],
    objective.best_Y[-objective.n_outlier :],
    alpha=0.8,
)

model2 = LinearRegression()
model2.fit(
    np.array(objective.best_X[: -objective.n_outlier]).reshape(-1, 1),
    objective.best_Y[: -objective.n_outlier],
)
Y2 = model2.predict(np.array([0, 1]).reshape(-1, 1))
plt.plot(np.array([0, 1]), Y2, label="R = {:.3f}".format(objective.corrcoef2))

model1 = LinearRegression()
model1.fit(np.array(objective.best_X).reshape(-1, 1), objective.best_Y)
Y1 = model1.predict(np.array([0, 1]).reshape(-1, 1))
plt.plot(np.array([0, 1]), Y1, label="R = {:.3f}".format(objective.corrcoef1))
plt.legend()
plt.show()

上記の図において、「青い直線」は、「青い点だけを用いて計算した回帰直線」を表しており、凡例にある「R = 0.464」は「青い点だけを用いて計算したピアソンの相関係数の値」を表しています。

一方、「オレンジの直線」は、「全ての点（青い点＋オレンジの点）を用いて計算した回帰直線」を表しており、凡例にある「R = -0.057」は「全ての点（青い点＋オレンジの点）を用いて計算したピアソンの相関係数の値」を表しています。

青い点だけを考えれば、相関係数が正であり、回帰直線の傾きも正なので「Xが大きくなるほどYが大きくなることが分かった」と言えそうな気がするのですが、オレンジの点を入れると結論が逆転してしまう、つまり「Xが大きくなるほどYが小さくなることが分かった」になってしまう例ができあがりました、というわけです。

このような時には、どんな結論を出せば良いでしょうか？少なくとも、「Xが大きくなるほどYが大きくなることが分かった」、「Xが大きくなるほどYが小さくなることが分かった」、どちらの結論も怪しいと言わざるを得ないでしょう。

こういうデータで結論を言おうとしているニュース等、案外多いので気をつけてください。

データサイズと外れ値の数を変えて影響を調べる

それでは、データサイズと外れ値の数を変えて影響を調べてみましょう。今度は

n_trials = 1000

のように試行回数を増やして徹底的にやってみます。

import optuna

for n_data in [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]:
    fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(20, 6))    for n_outlier in [1, 2, 3, 4, 5]:
        objective = Objective(n_data=n_data, n_outlier=n_outlier)
        study = optuna.create_study()
        study.optimize(objective, n_trials=n_trials, show_progress_bar=True)

        axes[n_outlier - 1].set_title(
            "n_data = {}, n_outlier = {}".format(n_data, n_outlier)
        )
        axes[n_outlier - 1].scatter(
            objective.best_X[: -objective.n_outlier],
            objective.best_Y[: -objective.n_outlier],
            alpha=0.8,
        )
        axes[n_outlier - 1].scatter(
            objective.best_X[-objective.n_outlier :],
            objective.best_Y[-objective.n_outlier :],
            alpha=0.8,
        )

        model2 = LinearRegression()
        model2.fit(
            np.array(objective.best_X[: -objective.n_outlier]).reshape(-1, 1),
            objective.best_Y[: -objective.n_outlier],
        )
        Y2 = model2.predict(np.array([0, 1]).reshape(-1, 1))
        axes[n_outlier - 1].plot(
            np.array([0, 1]), Y2, label="R = {:.3f}".format(objective.corrcoef2)
        )

        model1 = LinearRegression()
        model1.fit(np.array(objective.best_X).reshape(-1, 1), objective.best_Y)
        Y1 = model1.predict(np.array([0, 1]).reshape(-1, 1))
        axes[n_outlier - 1].plot(
            np.array([0, 1]), Y1, label="R = {:.3f}".format(objective.corrcoef1)
        )
        axes[n_outlier - 1].legend(
            bbox_to_anchor=(0, -0.1), loc="upper left", borderaxespad=0, fontsize=18
        )
        axes[n_outlier - 1].set_aspect("equal", "box")
        axes[n_outlier - 1].set_xlim([-0.05, 1.05])
        axes[n_outlier - 1].set_ylim([-0.05, 1.05])
    plt.savefig("lr{}.png".format(n_data))
    plt.show()