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

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

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

# Optuna で最適化

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

```python
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
```

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

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

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

```python
import optuna

study = optuna.create_study()
```

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

```python
n_trials = 100
```

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

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

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

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


```python
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()
```

![スクリーンショット 0005-04-21 12.19.23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/8679f0cd-d6d4-382d-0c7a-f1d4cd10d731.png)

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

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

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

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

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

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

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

```python
n_trials = 1000
```

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


```python
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()
```

以下、実行結果をじっくり長めに眺めてお楽しみください。

## データサイズ 10

![lr10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/3f717810-364f-5fe7-597c-51056392584d.png)


## データサイズ 20

![lr20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/2f49dff1-6df3-8f15-82a2-1b2235fb28af.png)


## データサイズ 30


![lr30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/fcb87add-466e-56e2-caeb-be086a6b8b25.png)

## データサイズ 40


![lr40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/f3b8530f-4002-aaa0-786c-0552973e1001.png)

## データサイズ 50


![lr50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/c1a6ca1a-77d0-f517-a927-48a278135a5b.png)

## データサイズ 60


![lr60.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/1ca2df8a-9b87-2431-b5f9-70fa7e44138a.png)

## データサイズ 70

![lr70.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211162/57479566-5f88-5620-9204-e7f939587e7f.png)

# 良心的なデータ提供者

まあ、そうは言っても、データのプロットを見せてくれる場合などはまだ良心的なんですよ。そのデータの解釈について、ああだこうだと検討できるので。

データを見せずに結論だけ言うタイプの記事などは、検討すらできないので困ったものです。

間違えることはあります。人間だもの。なので、検討できる材料を提供してくれる、良心的なデータ提供者であり続けたいものです。

# 結論

- 相関はありそうか？　ん？
- 相関はそう簡単に分からない