Dropout

データ分析が好きです。因果推論とか機械学習の解釈性に興味があります。

バイアス-バリアンスの分解と、アンサンブルの話


\def\calD{\mathcal{D}}
\def\paren#1{\left(#1\right)}
\def\brac#1{\left[#1\right]}
\def\quad#1{\paren{#1}^2}
\def\sumi{\sum_{i=1}^N}
\def\sumk{\sum_{k=1}^K}
\def\E#1{\mathbb{E}_{\calD}\brac{#1}}
\def\var#1{\mathrm{Var}_{\calD}\brac{#1}}
\def\cov#1{\mathrm{Cov}_{\calD}\brac{#1}}
\def\bias#1{\mathrm{Bias}_{\calD}\brac{#1}}
\def\hf{\hat{f}_{\calD}}
\def\hfk{\hf^k}
\def\fx{f(x^*)}
\def\hfx{\hf(x^*)}
\def\hfkx{\hfk(x^*)}
\def\Ehfx{\E{\hfx}}
\def\Ehfkx{\E{\hfkx}}
\def\mse{\mathrm{MSE}_{\calD}}
\def\ss{\sigma^2_*}
\def\cc{\sigma_{**}}

はじめに

この記事では、バイアス(偏り)とバリアンス(分散、ばらつき)の分解について、数式とシミュレーションの両面から確認し、

  • ある種の予測誤差はモデルのバイアスとバリアンスに分解できること
  • バイアスとバリアンスにはトレードオフの関係があること
  • モデルのアンサンブルはバリアンスの減少を通じて予測精度を改善できること
  • アンサンブルはモデル同士の相関が小さいほど有効なこと

を示します。

バイアスとバリアンスの分解

バイアスとバリアンスの分解を考えるため、まずは数式を用いた定式化を行います。

 J個の説明変数を X = (X_1, \dots, X_J)、目的変数を Yとします。
真の関数形は f(x)で、目的変数はこれにノイズ Uが乗って観測されるとします。

\begin{align}
Y_i = f(X_i) + U_i
\end{align}

 N個の観測されたデータの組を \calD = \{(X_i, Y_i)\}_{i=1}^Nとします。
適当なアルゴリズムを使って、データ \calDを学習して予測モデル \hf(x)を構築します。
得られた予測モデルを用いて新しいデータ x^*に対する予測を行います。
真の値 \fxと比べてこの予測値 \hfxがどのくらいずれるのかがこの記事の興味の対象になります。

予測誤差は二乗誤差で評価することにします。
観測データはサンプリングされたものなので、データ \calDのばらつき応じて学習済みモデルも変動し、予測値も変わるはずです。
そこで、二乗誤差の期待値をとることにし、これを平均二乗誤差(MSE)と呼ぶことにします。
このとき、新しいデータ x^*に対するMSEは以下で表現できます。

\begin{align}
\mse(x^*) = \E{\quad{\fx - \hfx}}
\end{align}

ここで、何に対して期待値をとっているかを明示的にするため \E{\cdot}という表記を用いています。

実は、MSEはモデルの偏り(バイアス)とモデルの安定性(バリアンス)に分解できることが知られています。
バイアスとバリアンスの正確な定義は後ほど触れます。
まず、 \Ehfxを足して引くことでMSEを3つのパートに分解できます。

\begin{align}
\mse(x^*)
=& \E{\quad{\fx - \Ehfx - \paren{\hfx - \Ehfx}}}\\
=& \E{\quad{\fx - \Ehfx}} \\
&- \E{2\paren{\fx - \Ehfx}\paren{\hfx - \Ehfx}} \\
&+ \E{\quad{\hfx - \Ehfx}}
\end{align}

ここで、 \fxとか \Ehfxは確率変数ではなく確定的な値なので期待値の外に出せます。

\begin{align}
\mse(x^*)
=& \quad{\fx - \Ehfx} \\
&- 2\paren{\fx - \Ehfx}\paren{\Ehfx - \Ehfx} \\
&+ \E{\quad{\hfx - \Ehfx}}
\end{align}

第二項は0になるので、残りは第一項と第三項です。
これらはそれぞれバイアスとバリアンスと呼ばれています。

\begin{align}
\bias{\hfx} &= \fx - \Ehfx, \\
\var{\hfx} &= \E{\quad{\hfx - \Ehfx}}
\end{align}

バリアンスはモデルの予測値 \Ehfxが訓練データの影響を受けてどのくらいばらつくのかを表す指標です。
僕たちが観測するデータは母集団からのサンプリングなので、仮想的に何度もサンプリングを行うと毎回少し違うデータが観測されるわけですが、そのデータが少し違うことに対してどのくらい予測値が安定的かがバリアンスに反映されます。
ちょっと違うデータが与えられたときにモデルの予測が大きく異なるならバリアンスは大きく、逆にデータが変動しても似たような予測値を返すならバリアンスは小さくなります。

一方で、バイアスはモデルの予測の期待値 \Ehfxが、真の値 \fxからどのくらいずれるのかを表す指標です。
ここで、あくまで期待値の意味でのずれに焦点があたっていることに注意してください。
例えば、バリアンスの大きいモデルは、仮想的に何度もサンプリングを行う状況を考えると毎回違う予測値を出しますが、それを平均すると真の値 \fxと近くなっているのか?ということをバイアスは測定していると解釈できます。


元の式に戻ると、MSEはバイアス(の二乗)とバリアンスに分解できることがわかります。

\begin{align}
\mse(x^*) = \bias{\hfx}^2 + \var{\hfx}
\end{align}

この分解はバイアス-バリアンス分解(bias-variance decomposition)と呼ばれています。

一般に、モデルが単純な場合はバリアンスが小さい代わりにバイアスが大きくなります。
一方で、モデルが複雑な場合はバイアスが小さい代わりにバイアスが大きくなります。
このように、バイアスとバリアンスにはトレードオフが存在し、これはbias-variance tradeoffと呼ばれています。
MSEはバイアスとバリアンスの足し算で決まるので、バイアスとバリアンスがいい塩梅になるようにモデルの複雑さを決めることが予測精度の高いモデルを構築する際には重要です。

なお、上記で計算したMSEはある入力 x^*に対するものなので、もし平均的な予測誤差を知りたい場合は

\begin{align}
\mathbb{E}_{X^*}\brac{{\mse(X^*)}}
\end{align}

のように入力の確率分布を用いて期待値をとる必要があることに注意してください。

シミュレーションでバイアスとバリアンスの関係を確認する

シミュレーションの設定

以下の設定でシミュレーションデータを生成します。

\begin{align}
Y &= \sin(6 \pi X) + U,\\
X &\sim \mathrm{Uniform}(0, 1),\\
U &\sim \mathcal{N}(0, 1)
\end{align}

ここで、説明変数 X区間 [0, 1]の一様分布から生成します。
 Xにsin関数を噛ませて、ノイズ Uを乗せて目的変数としています。
ノイズ Uは平均0、分散1の正規分布から生成します。

この設定で、データを生成してモデルを学習し予測する、ということを何回も繰り返して、モデルのバイアスとバリアンスを確認することにします。
シミュレーションを行いその結果を可視化するために以下のSimulatorクラスを実装しました。

@dataclass
class Simulator:
    """シミュレーションとその結果の可視化を行うクラス

    Args:
        n_grids: Xの理論値をどのくらい細かく作成するか
    """

    n_grids: int = 100

    def __post_init__(self) -> None:
        """理論値を作成"""

        # シミュレーションの関数形を特定
        self.f = lambda X: np.sin(6 * np.pi * X)

        # Xの範囲をグリッドで分割
        self.X_theory = np.linspace(start=0, stop=1, num=self.n_grids)

        # f(X)で理論値を作成
        self.y_theory = self.f(self.X_theory)

    def generate_simulation_data(
        self, n_instances: int
    ) -> tuple[np.ndarray, np.ndarray]:
        """シミュレーション用のデータを生成する

        Args:
            n_instances: シミュレーションで生成するインスタンスの数
        Returns:
            特徴量と目的変数のtuple
        """

        # Xは一様分布から、yは正規分布から生成する
        X = np.random.uniform(low=0, high=1, size=n_instances)
        u = np.random.normal(loc=0, scale=1, size=n_instances)
        # 関数f(X)を噛ませてノイズを足して目的変数とする
        y = self.f(X) + u

        return (X, y)

    def simulate(
        self, estimators: dict[str, Any], n_simulations: int, n_instances: int
    ) -> dict[str, np.ndarray]:
        """シミュレーションデータを生成して、学習と予測を行う

        Args:
            estimators: シミュレーション用のモデル
            n_simulations: シミュレーション回数
            n_instances: シミュレーションで生成するインスタンスの数
        Returns:
            モデルごとの予測値がまとまった辞書
        """

        y_preds = {
            k: np.zeros((self.n_grids, n_simulations)) for k in estimators.keys()
        }
        for s in range(n_simulations):
            X, y = self.generate_simulation_data(n_instances)
            for k, estimator in estimators.items():
                y_preds[k][:, s] = estimator.fit(X.reshape(-1, 1), y).predict(
                    self.X_theory.reshape(-1, 1)
                )

        return y_preds

    def decompose_bias_variance(
        self,
        y_preds: dict[str, np.ndarray],
    ) -> dict[dtr, dict[str, np.ndarray]]:
        """バイアスとバリアンスに分解する

        Args:
            y_preds: モデルごとの予測値
        Returns:
            BiasとVarianceがまとまった辞書
        """

        bias = {
            k: (self.y_theory - y_pred.mean(axis=1)) ** 2
            for k, y_pred in y_preds.items()
        }
        variance = {k: y_pred.var(axis=1) for k, y_pred in y_preds.items()}

        return {"Bias2": bias, "Variance": variance}

    def draw_prediction(
        self,
        y_preds: dict[str, np.ndarray],
        estimator_key: int,
        simulation_ids: list[int],
        ylim: tuple(float, float) | None = None,
    ) -> None:
        """予測結果を可視化

        Args:
            y_preds: 予測結果
            estimator_key: 可視化したいモデルの名前
            simulation_ids: 何回目のシミュレーション結果を可視化するか
            ylim: グラフの上限下限を調整したいときに使う
        """

        if ylim is None:
            y = np.concatenate([y for y in y_preds.values()])
            ylim = (y.min(), y.max())

        fig, axes = plt.subplots(ncols=2, figsize=(8, 4))

        # バイアスの可視化
        axes[0].plot(
            self.X_theory,
            y_preds[estimator_key].mean(axis=1),
            linewidth=2,
            label=f"予測値(シミュレーションの平均)",
        )
        axes[0].plot(self.X_theory, self.y_theory, color=".3", linewidth=2, label="理論値")
        axes[0].legend()
        axes[0].set(xlabel="X", ylabel="Y", ylim=ylim, title="Biasの可視化")

        # バリアンスの可視化
        axes[1].plot(self.X_theory, self.y_theory, color=".3", linewidth=2, label="理論値")

        for s in simulation_ids:
            axes[1].plot(
                self.X_theory,
                y_preds[estimator_key][:, s],
                linewidth=0.1,
                color=sns.color_palette()[0],
                label=f"予測値(シミュレーション{s:02.0f})",
            )
        axes[1].set(xlabel="X", ylabel="Y", ylim=ylim, title="Varianceの可視化")

        fig.suptitle(f"{estimator_key}のBiasとVarianceを可視化")
        fig.savefig(f"output/{estimator_key}のBiasとVarianceを可視化")
        fig.show()

    def draw_bias_variance(
        self, bias_variances: dict[str, dict[str, np.ndarray]], title: str
    ) -> None:
        """バイアスとバリアンスを可視化する

        Args:
            bias_variances: decompose_bias_variance()の結果
        """

        fig, axes = plt.subplots(ncols=2, figsize=(8, 4))
        for ax, (k, bvs) in zip(axes, bias_variances.items()):
            for estimator_key, value in bvs.items():
                # Xは一様分布に従うので単に平均をとればいい
                avg = value.mean()
                ax.plot(
                    self.X_theory,
                    value,
                    linewidth=1,
                    label=f"{estimator_key}(AVG:{avg: .2f})",
                )
                ax.legend()
                ax.set(xlabel="X", title=k)
        fig.suptitle(title)
        fig.savefig(f"output/{title}")
        fig.show()

シミュレーション結果の確認

予測モデルとしては取り扱いの簡単さを重視して決定木を採用することにします。
ハイパーパラメータmax_depthを調整することでモデルの複雑さを簡単に制御することができます。

max_depthは決定木の最大分割回数を制御するハイパーパラメータで、例えば、max_depth=1なら決定木は一度しか分割を行わないので、2種類の予測値しか出しません。
この場合、モデルが非常に単純なので、バイアスが大きくバリアンスが小さくなることが予測されます。
一方で、max_depth=10なら最大で2の10乗つまり1024種類の予測値を出せることになります。
この場合はmax_depth=1よりもモデルが複雑になるので、バイアスが小さくバリアンスが大きくなるでしょう。

それでは、実際にmax_depth=1とmax_depth=10の決定木でシミュレーションを行います。
シミュレーションの回数は100回、一回のシミュレーションで生成するインスタンスの数は1000とします。

# シミュレーターインスタンスの作成
simulator = Simulator()

# 単純なモデルと複雑なモデルを用意
dts = {
    "DecisionTree(01)": DecisionTreeRegressor(max_depth=1),
    "DecisionTree(10)": DecisionTreeRegressor(max_depth=10),
}

# 予測結果
y_preds = simulator.simulate(estimators=dts, n_simulations=100, n_instances=1000)

# バイアスとバリアンスに分解
bias_variances = simulator.decompose_bias_variance(y_preds)

まずは、単純なmax_depth=1の場合を確認します。
モデルはsin関数の波をまともに捉えられておらず、バイアスが大きいことがわかります。
一方で、どのシミュレーションでも左端で分割するか右端で分割するかの2択がばらついているだけでほぼ同じような予測を行っており、バリアンスは小さいことが見て取れます。

# 単純なモデルの予測結果を可視化
# シミュレーション結果をすべて可視化するとうるさいので20個にしぼる
simulator.draw_prediction(
    y_preds=y_preds, simulation_ids=range(20), estimator_key="DecisionTree(01)"
)

f:id:dropout009:20210726192759p:plain

次に、複雑なmax_depth=10の場合を確認します。
モデルは(平均的には)sin関数の波を十分に捉えることができており、バイアスが小さいことがわかります。
一方で、シミュレーションごとに予測値はかなり荒れており、バリアンスは大きいことが見て取れます。

# 複雑なモデルの予測結果を可視化
simulator.draw_prediction(
    y_preds=y_preds, simulation_ids=range(20), estimator_key="DecisionTree(10)"
)

f:id:dropout009:20210726192803p:plain

モデルごとにバイアス(の二乗)とバリアンスを直接比較したグラフが下図になります。
単純なmax_depth=1のモデルはバイアスが大きい反面バリアンスが小さく、複雑なmax_depth=10のモデルはバイアスが小さい反面バリアンスが大きいことが見て取れます。
このように、シミュレーション結果はモデルのバイアスとバリアンスにはトレードオフがあることを示唆しています。

# モデルごとのバイアスとバリアンスを比較
simulator.draw_bias_variance(bias_variances=bias_variances, title="モデルごとのBias-Varianceの比較")

f:id:dropout009:20210726192806p:plain

バイアスとバリアンスのトレードオフを可視化する

単純なモデルと複雑なモデルにはバイアスとバリアンスのトレードオフがあることが示唆されました。
このトレードオフをより詳細に確認するため、max_depthを1から10まで変化させた際にバイアスとバリアンスがどのように変化するかを追跡しましょう。

# max_depthを1から10まで変化させる
max_depth_range = list(range(1, 11))
dts = {f'DT({x:02})': DecisionTreeRegressor(max_depth=x) for x in max_depth_range}

# 予測値を作成
y_preds = simulator.simulate(estimators=dts, n_simulations=100, n_instances=1000)

# バイアスとバリアンスを分解
bias_variances = simulator.decompose_bias_variance(y_preds)
def draw_bias_variance_tradeoff(
    bias_variances: dict[str, dict[str, np.ndarray]], max_depth_range: list
) -> None:
    """max_depthを変化させた際のバイアスとバリアンスの可視化"""

    bias = np.array([v.mean() for v in bias_variances["Bias2"].values()])
    variance = np.array([v.mean() for v in bias_variances["Variance"].values()])
    mse = bias + variance

    fig, ax = plt.subplots()
    ax.plot(max_depth_range, bias, label="Bias2")
    ax.plot(max_depth_range, variance, label="Variance")
    ax.plot(max_depth_range, mse, label="MSE")
    ax.legend()
    ax.set_xticks(max_depth_range)
    ax.set(xlabel="max_depth", ylabel="Bias/Variance/MSE")
    fig.suptitle("Bias-Varianceのトレードオフ")
    fig.savefig("Bias-Varianceのトレードオフ")
    fig.show()


draw_bias_variance_tradeoff(bias_variances, max_depth_range)

f:id:dropout009:20210726192948p:plain

モデルを複雑にするつれてバイアスが小さくなる反面バリアンスが大きくなっていく様子が見て取れます。
また、MSEはバイアスの二乗とバリアンスの足し算なので、最も理論値との誤差が小さくなるいい塩梅のハイパーパラメータはmax_depth=5であることがわかります。

モデルのアンサンブルでばらつきを抑える

ここまで、バイアスとバリアンスにはトレードオフがあることを見てきました。
モデルを複雑にすればバイアスは小さく出来ますが、代わりにバリアンスが大きくなってしまいます。
実は、バイアスを変化させずにバリアンスを小さくし得る手法として、モデルのアンサンブルが知られています。

アンサンブルの効果を数式で確認する

 K個のベースになる予測モデル \paren{\hf^1(x), \dots, \hf^K(x)}をアンサンブルして予測を行うことを考えます。
一番単純なアンサンブルとして、ベースモデルの単純平均を考えることにします。

\begin{align}
\hfx = \frac{1}{K}\sumk\hfk(x)
\end{align}

バイアス

予測モデルのアンサンブルはバイアスとバリアンスにどのような影響を与えるでしょうか?
まずはモデルのバイアスについて考えます。
バイアスの計算式に上式を代入して変形します。

\begin{align}
\bias{\hfx}
&= \fx - \Ehfx\\
&= \fx - \frac{1}{K}\sumk\hfkx\\
&= \frac{1}{K}\sumk\paren{\fx - \hfkx}\\
&= \frac{1}{K}\sumk\bias{\hfkx}\\
\end{align}

ここから、アンサンブルモデルの予測値のバイアスは各ベースモデルのバイアスの平均になることがわかります。
話を簡単にするために、ベースモデルのバイアスがすべて同じく \deltaだとします。
実際、バギングのように単に復元抽出を繰り返して予測モデルを複数作成しその結果を単純平均するようなアルゴリズムでは、ベースモデルのバイアスは一定だと考えられます。
このとき、アンサンブルモデルのバイアスは

\begin{align}
\bias{\hfx} &= \frac{1}{K}\sumk \delta = \delta
\end{align}

となり、アンサンブルによってバイアスが小さくなるなどの改善は見られないことがわかります。

バリアンス

次に、モデルのバリアンスについて考えます。

 
\begin{align}
    \var{\hfx} 
    &= \var{\frac{1}{K}\sumk\hfkx}\\
    &= \quad{\frac{1}{K}}\paren{\sumk\var{\hfkx} + \sum_{k\neq k'}\cov{\hfkx, \hf^{k'}(x^*)}}
\end{align}

話を簡単にするために、すべてのベースモデルに関して分散が \ss、共分散が \ccで一定だとします。
実際、バギングだとベースモデルは学習データが異なる以外に差分はないので分散も共分散も一定になります。
このとき、アンサンブルモデルのバリアンスは

 
\begin{align}
    \var{\hfx} 
    &= \quad{\frac{1}{K}}\paren{\sumk\ss + \sum_{k\neq k'}\cc}\\
    &= \quad{\frac{1}{K}}\paren{K\ss + K(K-1)\cc}\\
    &= \paren{\frac{1}{K}}\ss + \paren{\frac{K-1}{K}}\cc
\end{align}

と単純化できます。

ひとつ極端な状況として、仮にベースモデル同士の相関が全く無ケースを考えます。
このとき \cc = 0であり、アンサンブルモデルのバリアンスは

 
\begin{align}
    \var{\hfx} &= \paren{\frac{1}{K}}\ss 
\end{align}

になります。
よって、ベースモデルの数を増やせば増やすほどアンサンブルモデルのバリアンスは小さくなっていき、結果としてMSEつまり予測誤差も小さくなります。

実際は、ベースモデル同士は大なり小なり相関があるので、ここまで極端な状況は起きませんが、それでも相関が小さいほどアンサンブルモデルのバリアンスは小さくなりやすいことは事実です。
これがアンサンブルを行う際のベースモデルは互いに相関が小さい方がいいと言われる理由です。


まとめると、モデルのアンサンブルはバイアスには影響を与えませんが、バリアンスを小さくすることを通じて予測誤差を改善させる効果があると言えます。

アンサンブルの効果をシミュレーションで確認する

モデルのアンサンブルによるバリアンスの改善が数式から確認できたので、さらにシミュレーションでも改善効果を確認しておきます。
まずは、アンサンブルモデルとして、決定木をベースモデルとしてバギングの結果を返すBaggingTreeRegressorクラスを実装しておきます。

バギングのアルゴリズムは以下になります。

  • 訓練データと同じインスタンス数だけデータを復元抽出する
  • 復元抽出したデータを使ってベースモデルを学習し、保存しておく
  • これを任意の回数繰り返す
  • 予測を行う際は、複数作成した学習済みベースモデルを用いて複数の予測を行い、結果を平均して最終的な予測値とする

復元抽出のたびに訓練データがちょっと変化するので、個々のベースモデルは似ているが少し異なる予測値を出すことが予想されます。
このような、互いに相関はするが完全に相関しているわけではない予測結果を平均することで、バリアンスを抑えることがバギングの狙いです。

@dataclass
class BaggingTreeRegressor:
    """
    決定木をベースモデルとしてBootstrap Aggregationを行う

    Args:
        n_estimators: いくつツリーを作成するか
        tree_params: ツリーに渡すハイパーパラメータ
    """

    n_estimators: int
    tree_params: dict

    def fit(self, X: np.ndarray, y: np.ndarray) -> BaggingRegressor:
        """データを復元抽出して学習"""

        # インスタンスの数
        n_instances = X.shape[0]

        # 復元抽出のためのIDを「インスタンスの数 x ベースモデルの数」だけ作成
        ids = np.random.choice(n_instances, size=(self.n_estimators, n_instances))

        # 復元抽出されたデータで個々のベースモデルを学習して保存
        self.estimators = []
        for e in range(self.n_estimators):
            self.estimators.append(
                DecisionTreeRegressor(**self.tree_params).fit(X[ids[e]], y[ids[e]])
            )

        return self

    def predict(self, X: np.ndarray) -> np.ndarray:
        """予測結果を平均して返す"""

        return np.column_stack([e.predict(X) for e in self.estimators]).mean(axis=1)

それでは、単純な決定木とそのバギングでバイアスとバリアンスがどのように異なるのかをシミュレーションします。
max_depthは先程のシミュレーションで最も予測誤差の小さくなる5を指定します。
また、バギングでは100本の決定木をアンサンブルすることにします。

# 決定木とそのバギングバージョンを用意
dts = {
    'DecisionTree(05)': DecisionTreeRegressor(max_depth=5),
    'BaggingTree': BaggingTreeRegressor(n_estimators=100, tree_params={'max_depth': 5}),
}

# 予測結果
y_preds = simulator.simulate(estimators=dts, n_simulations=100, n_instances=1000)

# バイアスとバリアンスの分解
bias_variances = simulator.decompose_bias_variance(y_preds)

決定木の予測結果をバギングの予測結果を比較すると、バイアスは同じようなものですがバギングはバリアンスがかなり抑えられていることが見て取れます。

# 決定木の予測結果を可視化
simulator.draw_prediction(y_preds=y_preds, simulation_ids=range(20), estimator_key='DecisionTree(05)')

f:id:dropout009:20210726192811p:plain

# バギングの予測結果を可視化
simulator.draw_prediction(y_preds=y_preds, simulation_ids=range(20), estimator_key='BaggingTree')

f:id:dropout009:20210726192817p:plain

実際、単純な決定木とバギングの結果を直接比較すると、バリアンスが三分の一程度に抑えられていることがわかります。
モデルのアンサンブルによるバリアンスの改善をシミュレーションからも確認することができました。

# 決定木とバギングでバイアスとバリアンスを比較
simulator.draw_bias_variance(bias_variances, title="決定木とバギングのBias-Varianceの比較")

f:id:dropout009:20210726192823p:plain

まとめ

この記事では、バイアスとバリアンスの分解について、理論とシミュレーションの両面から確認しました。

  • モデルの予測精度を表すMSEは、バイアスとバリアンスに分解できる
  • バイアスとバリアンスにはトレードオフがある。モデルが単純だとバイアスが大きいがバリアンスは小さく、モデルが複雑だとバイアスが小さいがバリアンスは大きい
  • モデルのアンサンブルはバリアンスを減少させることを通じて予測精度を改善することができる。モデル同士の相関が小さいほどアンサンブルの効果が高い

この記事で利用したシミュレーションコードは以下になります。

https://github.com/dropout009/bias_variance

参考文献

  • Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. The elements of statistical learning. Vol. 1. No. 10. New York: Springer series in statistics, 2001.

多重共線性と回帰係数の信頼性の話。あとリッジ回帰。


\def\paren#1{\left(#1\right)}
\def\brac#1{\left[#1\right]}
\def\sumi{\sum_{i=1}^N}
\def\P#1{\Pr\paren{#1}}
\def\E#1{\mathbb{E}\brac{#1}}
\def\var#1{\mathrm{Var}\brac{#1}}
\def\hvar#1{\widehat{\mathrm{Var}}\brac{#1}}
\def\quad#1{\paren{#1}^2}
\def\inv#1{\paren{#1}^{-1}}
\def\Xij{X_{i, j}}
\def\vX{\mathbf{X}}
\def\vXj{\mathbf{X}_j}
\def\vXmj{\mathbf{X}_{\setminus j}}
\def\vY{\mathbf{Y}}
\def\vU{\mathbf{U}}
\def\hvU{\hat{\vU}}
\def\vVj{\mathbf{V}_j}
\def\hvVj{\hat{\mathbf{V}}_j}
\def\mM{\mathbf{M}}
\def\mMj{\mathbf{M}_j}
\def\mMmj{\mathbf{M}_{\setminus j}}
\def\vP{\mathbf{P}}
\def\mW{\mathbf{W}}
\def\vI{\mathbf{I}}
\def\hbeta{\hat{\beta}}
\def\betamj{\beta_{\setminus j}}
\def\hbetamj{\hbeta_{\setminus j}}
\def\hgamma{\hat{\gamma}}
\def\hXij{\hat{X}_{i, j}}
\def\hVij{\hat{V}_{i, j}}
\def\bXj{\bar{X}_{j}}

はじめに

先日、多重共線性に関する @hizakayuさんや@M123Takahashiさんのコメントを目にしました。

恥ずかしながら、ここで紹介されている多重共線性と推定量の信頼性に関する性質のことは普通に忘れていました。
こちらのコメントを見てそういえばこんなこと習ったなと思い出し、大学院のころの資料*1計量経済学の教科書*2を確認し、改めて学びなおしたことを記事としてまとめました。

この記事では、まずは多重共線性がある場合のOLS推定量の分散について数式を用いて確認します。
次に、数式で確認した結果をシミュレーションを用いて実際にみていきます。
最後に、推定量の信頼性を高める手法としてリッジ推定量を紹介し、リッジ推定量は(バイアスという犠牲を払うことで)信頼性を高められるを確認します。

多重共線性とOLS推定量の信頼度

線形回帰モデルの導入

以下のような線形回帰モデル

\begin{align}
Y_i = X_i'\beta + U_i
\end{align}


を考えます。ここで、 Y_iは目的変数、 X_i = (X_{i, 1}, \dots, X_{i, J})'は説明変数、 U_iは誤差項を表します。
 i = 1, \dots, Nはサンプルを表す添字で、 \{X_i, Y_i\}_{i=1}^Nはランダムサンプリングされているとします。
 \beta = (\beta_1, \dots, \beta_J)'は回帰係数です。

以下、表記を単純にするため行列を利用します。

\begin{align}
\vY = \vX\beta + \vU
\end{align}


それぞれの行列は以下で構成されています。

\begin{align}
\vY = \underbrace{\begin{pmatrix}
Y_1 \\
\vdots \\
Y_N
\end{pmatrix}}_{N \times 1},\;\;
\vX = \begin{pmatrix}
X_{1}'\\
\vdots \\
X_{N}'
\end{pmatrix}
= \underbrace{\begin{pmatrix}
X_{1, 1} & \cdots & X_{1, J}\\
\vdots & \ddots & \vdots \\
X_{N, 1} & \cdots & X_{N, J}
\end{pmatrix}}_{N \times J},\;\;
\vU = \underbrace{\begin{pmatrix}
U_1 \\
\vdots \\
U_N
\end{pmatrix}}_{N \times 1}
\end{align}

標準的な仮定として \E{\vU \mid \vX} = 0をおきます。
また、話を簡単にするために、分散均一性 \var{\vU \mid \vX} = \sigma^2\vI_Nも仮定しておきます。

OLS推定量

このとき、最小二乗法(Ordinary Least Squares)を用いた推定量は以下で与えられます。

\begin{align}
\hbeta &= \arg\min_{\beta} \paren{\vY - \vX\beta}'\paren{\vY - \vX\beta}\\
&=\inv{\vX'\vX}\vX'\vY
\end{align}

これをOLS推定量と呼びます。
いま興味があるのは、このOLS推定量の信頼度と多重共線性の関係です。
この記事では、多重共線性があるとOLSを推定量の信頼性が低くなっていくことを示します。

まず第一に、もし完全な多重共線性が発生していると、 \vX'\vXがランク落ちして逆行列が取れません。
よって、完全な多重共線性がある場合はそもそも回帰係数が推定できないということになります。
これは多重共線性が最も問題になるケースです。
この場合、完全な多重共線性を回避すべく余計な説明変数を削除するなどの調整が必要になります。

以下では、完全な多重共線性は起きていないとして話を進めます。
つまり、 \vX'\vXがフルランクであることを仮定します。
 \vY = \vX\beta + \vUだったことを思い出すと、OLS推定量 \hbetaは以下のように変形できます。

\begin{align}
\hbeta &= \inv{\vX'\vX}\vX'\paren{\vX\beta + \vU}\\
&= \inv{\vX'\vX}\vX'\vX\beta + \inv{\vX'\vX}\vX'\vU\\
&= \beta + \inv{\vX'\vX}\vX'\vU
\end{align}

この変形からOLSを推定量に関する重要な指標である期待値と分散を知ることができます。
まずは期待値をとると、 \E{\vU \mid \vX} = 0の仮定のおかげで

\begin{align}
\E{\hbeta \mid \vX} &= \beta + \inv{\vX'\vX}\vX'\E{\vU \mid \vX}\\
&= \beta
\end{align}

となり、OLS推定量の期待値は真の回帰係数 \betaと一致することがわかります。
この性質は不偏性と呼ばれています。

また、OLS推定量の分散は


\begin{align}
\var{\hbeta_j \mid \vX} 
&= \var{\inv{\vX'\vX}\vX'\vU \mid \vX}\\
&= \inv{\vX'\vX}\vX'\var{\vU \mid \vX}\vX\inv{\vX'\vX}\\
&= \inv{\vX'\vX}\vX'\sigma^2\vI_N\vX\inv{\vX'\vX}\\
&= \inv{\vX'\vX}\sigma^2
\end{align}

となることがわかります。
ここで、分散均一性の仮定 \var{\vU \mid \vX} = \sigma^2\vI_Nを使っていることに注意して下さい。

定量の分散はその推定量がどのくらい安定しているのかを表す指標と言えます。
実際、推定量の分散の平方根をとったものは標準誤差と呼ばれていて、信頼区間の計算や統計的仮説検定に利用します。

ところで、いま知りたいことは、多重共線性と推定量の信頼度(=分散)の関係でした。
残念ながら、 \var{\hbeta_j \mid \vX} = \inv{\vX'\vX}\sigma^2という式からは多重共線性がある場合に分散がどうなるのかいまいち想像ができません。
そこで、 \var{\hbeta_j \mid \vX}を別の形で表現することにします。

Annihilator Matrix

後の行列操作を楽にするために、Annihilator Matrix  \mMを導入します(これがなぜAnnihilatorと呼ばれているのかはよくわかりません)。

\begin{align}
\mM = \vI_N - \vX\inv{\vX'\vX}\vX'
\end{align}

行列 \mMは以下の性質を持ちます。
かなり便利なのでこの記事を読む間だけ覚えておいてください。

  •  \mM' = \mM:転地しても同じ行列
  •  \mM\mM = \mM:かけても同じ行列
  •  \mM\vX = 0:説明変数にかけるとゼロ
  •  \mM\vY = \vY - \vX\hbeta = \hvU:目的変数にかけると残差になる

実際、簡単な行列計算から上記の性質が成り立っていることが見て取れます。

\begin{align}
\mM'
&= \paren{ \vI_N - \vX\inv{\vX'\vX}\vX'}' \\
&= \vI_N - \vX\inv{\vX'\vX}\vX' \\
&= \mM\\
\mM\mM
&= \paren{ \vI_N - \vX\inv{\vX'\vX}\vX'}\paren{\vI_N - \vX\inv{\vX'\vX}\vX'} \\
&= \vI_N - \vX\inv{\vX'\vX}\vX' - \vX\inv{\vX'\vX}\vX' + \vX\inv{\vX'\vX}\vX' \\
&= \mM\\
\mM\vX &= X - \vX\inv{\vX'\vX}\vX'\vX \\
&= 0\\
\mM\vY &= Y - \vX\inv{\vX'\vX}\vX'\vY \\
&= \vY - X\hbeta\\
&= \hvU
\end{align}

 \mM\vY = \vY - \vX\hbeta = \hvUという性質から、 \mMは線形回帰モデルの残差を作る行列として知られています。

OLS推定量の別表現

それでは、残差を作る行列 \mMを用いてOLS推定量を別の形で表現できることを示します。
まず、線形回帰モデルを説明変数 \vXjに関する部分と、それ以外の部分 \vXmjに分解します。

\begin{align}
\vY = \vXj\beta_j + \vXmj\betamj + \vU
\end{align}

ここで、それぞれの行列は以下を表しています。

\begin{align}
\vXj &= (X_{1, j}, \dots, X_{N, j})'\\
\vXmj &= (\vX_1, \dots, \vX_{j-1}, \vX_{j+1}, \dots \vX_{J})\\
\betamj &= (\beta_1, \dots, \beta_{j-1}, \beta_{j+1}, \dots \beta_{J})'
\end{align}

いま、説明変数 \vXjに対するOLS推定量 \hbeta_jを知りたいとします。
通常の最小二乗法はすべての \betaを同時に動かして二乗誤差を最小にしますが、
まず \betamjを動かして二乗誤差を最小にし、残った部分に対して改めて \beta_jを動かして二乗誤差を最小にする \hbeta_jを求める、という二段階のステップを踏んでも同じ結果になります。
この操作を数式で表現すると以下になります。

\begin{align}
\hbeta_j &= \arg\min_{\beta_j}\brac{\min_{\betamj}\paren{\vY - \vXj\beta_j - \vXmj\betamj}'\paren{\vY - \vXj\beta_j - \vXmj\betamj}}
\end{align}

内側の最小化問題をよく見ると、これは \vY - \vXj\beta_jに対して \vXmj\betamjで回帰していると解釈できるので、そのOLS推定量は以下で求めることができます。

\begin{align}
\hbetamj &= \inv{\vXmj'\vXmj}\vXmj'\paren{\vY - \vXj\beta_j}
\end{align}

これを内側の最小化問題に戻すことを考えると

\begin{align}
\vY - \vXj\beta_j - \vXmj\hbetamj
&= \vY - \vXj\beta_j - \vXmj\inv{\vXmj'\vXmj}\vXmj'\paren{\vY - \vXj\beta_j}\\
&= \paren{\vI_N - \vXmj\inv{\vXmj'\vXmj}\vXmj'}\paren{\vY - \vXj\beta_j}\\
& = \mMmj\paren{\vY - \vXj\beta_j}\\
&= \mMmj\vY - \mMmj\vXj\beta_j
\end{align}

なので、元の最小化問題は

\begin{align}
\hbeta_j &= \arg\min_{\beta_j}\paren{\mMmj\vY - \mMmj\vXj\beta_j}\paren{\mMmj\vY - \mMmj\vXj\beta_j}
\end{align}

となることが分かります。

ここから、OLS推定量

\begin{align}
\hbeta_j &= \inv{\vXj'\mMmj'\mMmj\vXj}\vXj'\mMmj'\mMmj\vY\\
&= \inv{\vXj'\mMmj\vXj}\vXj'\mMmj\vY\\
\end{align}

という形式でも表現できることがわかりました。

補助回帰

OLS推定量を別形式で表現できましたが、ここに出てくる \mMmj\vXjは何を計算しているのでしょうか?
残差を作る行列 \mMmjの意味を考えると、 \mMmj\vXjは説明変数 \vXjを他の説明変数 \vXmjで回帰したときの残差を表しています。
つまり、線形回帰モデル

\begin{align}
\vXj = \vXmj\gamma + \vVj
\end{align}

に対する残差

\begin{align}
\mMmj\vXj = \vXj - \vXmj\inv{\vXmj'\vXmj}\vXmj'\vXj = \vXj - \vXmj\hgamma = \hvVj
\end{align}

を表現しています。
説明変数 \vXjをその他の説明変数 \vXmjで回帰することは補助回帰(auxiliary regression)と呼ばれています。
残差が \hvVj小さければ、説明変数 \vXjは他の説明変数 \vXmjでほぼ表現できてしまうということになります。
この場合、説明変数 \vXjは追加的に与える情報量の少ない説明変数ということになります。
一方で、残差が \hvVjが大きい場合は、説明変数 \vXjは他の説明変数では与えることのできない情報をモデルに与えていると言えます。

これを多重共線性の文脈でいうと、残差 \hvVjが小さいことは(説明変数 \vXjに対する)多重共線性が強いことに対応し、残差 \hvVjが大きいことは多重共線性が弱いことに対応します。
実際、完全な多重共線性がある状態として \vXj = \vX_kとなるような別の特徴量が存在するケースを考えると、残差 \hvVjもゼロになります。

OLS推定量の分散を解釈する

それでは、改めてOLS推定量の分散について考えましょう。
 \vY = vXj\beta_j + \vXmj\betamj + \vUであったことを思い出すと、OLS推定量は以下のように変形できます。

\begin{align}
\hbeta_j &= \inv{\vXj'\mMmj \vXj}\vXj'\mMmj\paren{\vXj\beta_j + \vXmj\betamj + \vU}\\
&= \underbrace{\inv{\vXj'\mMmj \vXj}\vXj'\mMmj\vXj}_{=\vI_J}\beta_j
+ \inv{\vXj'\mMmj \vXj}\vXj'\underbrace{\mMmj\vXmj}_{=0}\betamj
+ \inv{\vXj'\mMmj \vXj}\vXj'\mMmj\vU\\
&= \beta_j + \inv{\vXj'\mMmj \vXj}\vXj'\mMmj\vU\\
&= \beta_j + \inv{\vXj'\mMmj\mMmj\vXj}\vXj'\mMmj\mMmj\vU\\
&= \beta_j + \inv{\hvVj'\hvVj}\hvVj'\vU
\end{align}

これの分散をとると


\begin{align}
\var{\hbeta_j \mid \vX} 
&= \var{\inv{\hvVj'\hvVj}\hvVj'\vU \mid \vX}\\
&= \inv{\hvVj'\hvVj}\hvVj'\var{\vU \mid \vX}\hvVj\inv{\hvVj'\hvVj}\\
&= \inv{\hvVj'\hvVj}\hvVj'\sigma^2\vI_N\hvVj\inv{\hvVj'\hvVj}\\
&= \inv{\hvVj'\hvVj}\sigma^2
\end{align}

となります。
さらに、 \hvVj'\hvVj = \sumi \hVij^2 = \sumi\quad{\Xij - \hXij}であることを利用すると、


\begin{align}
\var{\hbeta_j \mid \vX} 
&= \frac{\sigma^2}{\sumi\quad{\Xij - \hXij}}
\end{align}

と表現できます。
説明変数 \vXjの平均を \bXj = \frac{1}{N}\sumi X_{i, j}として、分母に \sumi \quad{\Xij - \bXj}をかけて割ると


\begin{align}
\var{\hbeta_j \mid \vX} &= \frac{\sigma^2}{\sumi \quad{\Xij - \hXij} / \sumi \quad{\Xij - \bXj}} \times \frac{1}{N \times \frac{1}{N}\sumi \quad{\Xij - \bXj}}\\
&= \frac{\sigma^2}{N\paren{1 - R^2_j}\hvar{X_j}}
\end{align}

を得ます。
ここで、
\begin{align}
R^2_j = 1 - \frac{\sumi \quad{\Xij - \hXij}}{\sumi \quad{\Xij - \bXj}}
\end{align}
は説明変数 \vXjをその他の説明変数 \vXmjで回帰したときの決定係数です。 \vXjのばらつきのうち平均では説明できなかった部分を \vXmjでどのくらい説明できるかを表しています。つまり、多重共線性の強さを表す指標です。ちなみに、 \vXmjが一変数の場合は相関係数の二乗に相当します。

また、 \hvar{X_j}は説明変数 \vX_jの標本分散を表しています。

さて、ここまで長い道のりでしたが、この表現から多重共線性がOLS推定量の分散(信頼度)にどう影響するかを解釈することができます。

1. 説明変数 \vXでは説明できない目的変数 \vYの変動 \sigma^2が大きければ大きいほど、分散が大きくなる
2. サンプルサイズ Nが大きくなると、分散が小さくなる
3. 説明変数 \vX_jの標本分散 \hvar{X_j}が大きくなると、分散が小さくなる
4. 決定係数 R^2_jが大きくなると、分散が大きくなる

1.はデータの限界みたいなところがあるので言っても仕方ないかなと思っています。

2.は当たり前といえば当たり前で、データが多いと推定量の信頼性が上がっていくのはよく言われていることですね。
ただ、 \var{\hbeta_j \mid \vX} = \inv{\vX'\vX}\sigma^2という表現からは明示的には見えてこなかった部分ではあります。

3.は説明変数の取りうる値のレンジが大きいと回帰係数が小さくなることからくる性質なので、これもあんまり気にしなくていいと思います。
実際、説明変数をあらかじめ標準化すると分散が1になるので無視できます。
これも \var{\hbeta_j \mid \vX} = \inv{\vX'\vX}\sigma^2という表現からは明示的には見えてこなかった部分ではあります。

やはり重要なのは4.で、説明変数 \vXjをその他の説明変数 \vXmjで説明できてしまう割合が大きいと、OLS推定量 \hbeta_jの信頼性が下がることが見て取れます。
つまり、多重共線性があるケースで推定量の信頼性が下がることがわかります。
ここで重要なことは、 \hbeta_jの信頼性に影響するのはあくまで \hbeta_jに対する決定係数 R_j^2であることです。
言い換えると、 \vXjとは関係ない説明変数同士で強い相関があったとしても、それらの説明変数が \vXjに対して説明力を持たないのであれば、 \hbeta_jの信頼性に特に影響は与えません。
ですので、多重共線性を考える際には、あくまでいま興味のあるターゲットの説明変数 \vXjに対する多重共線性のみに注意すればよく、他の説明変数で多重共線性が起きていても(その説明変数の回帰係数に興味がないなら)無視しても良いと言えます。

これが特に重要になってくるのは因果関係を考えるケースで、多重共線性を気にしてコントロール変数をむやみにモデルから外すと欠落変数バイアスが発生してしまいます。
ターゲットの説明変数と無関係の多重共線性は気にせず、論理的に入れるべきコントロール変数は入れるべきだと思います。
ターゲットの説明変数に対して強い多重共線性が発生しているケースは問題といえば問題だとは思うのですが、そもそもコントロール変数でほとんどの部分が説明できてしまう変数をターゲットにするのはどうなんだという話にもなるかなと思っています。

まとめると、完全な多重共線性はそもそも計算ができないので問題ですが、そうでなければ多重共線性はあまり気にしなくていいケースがほとんどだと思っています。

シミュレーションによる信頼度の確認

シミュレーションの設定

OLS推定量の分散の振る舞いを数式からも確認しましたが、ここではシミュレーションを用いて理解を深めることにします。
以下の設定でシミュレーションデータを生成します。


\begin{align}
    Y &= X_0 + X_1 + X_2 + U\\
    (X_0, X_1, X_2) &\sim \mathcal{N}\paren{
    \begin{pmatrix}
    0 \\ 0 \\ 0
    \end{pmatrix},\;
    \begin{pmatrix}
    1 & \rho & 0\\
    \rho & 1 & 0\\
    0 & 0 & 1 
    \end{pmatrix}
    }\\
    U &\sim \mathcal{N}\paren{0, \sigma^2}
\end{align}

ここで、説明変数 X_0 X_1は相関していて、 X_2は無相関であるような設定を考えています。
話を単純にするために、回帰係数 (\beta_0, \beta_1, \beta_1)はすべて1としています。

多重共線性とOLS推定量の信頼度

まずは多重共線性がOLS推定量の分散に与える影響を見ていきます。
多重共線性の強さによる影響の大きさを知るために、説明変数 X_0 X_1相関係数 \rhoを変化させた場合の分散の大きさを確認します。

f:id:dropout009:20210722194753p:plain

上図を確認すると、説明変数  (X_0, X_1)の相関が強くなるにつれて、  (\beta_0, \beta_1)に関するOLS推定量の分散 \hvar{\beta_0 \mid \vX} \hvar{\beta_1 \mid \vX}が大きくなっていく様子が見て取れます。
よって、多重共線性が強くなると推定量の信頼性は低くなると言えます。

一方で、説明変数  (X_0, X_1)の相関は、  \beta_2に関するOLS推定量の分散 \hvar{\beta_2 \mid \vX}にはなんの影響も与えていないこともわかります。
数式でも確認したように、もし回帰係数  \beta_2にのみ興味がある場合は、説明変数  (X_0, X_1)の多重共線性は無視してモデルに含めてもいいことがわかりました。

サンプルサイズとOLS推定量の信頼度

f:id:dropout009:20210722194757p:plain

次に、相関係数  \rhoを0.9で固定したままサンプルサイズが変化した場合のOLS推定量の分散を可視化したものが上図になります。
強い多重共線性があるぶん  (\beta_0, \beta_1)の信頼度は高めになっていますが、それでもサンプルサイズを大きくれば分散が小さくなっていく様子が見て取れます。
逆に言うと、興味のある説明変数に多重共線性が存在している場合は、できるだけたくさんデータを集める必要があることがわかります。

リッジ推定量で推定量の信頼性を高める

リッジ推定量の導入

ここまで、強い多重共線性が存在すると回帰係数の推定が安定せずOLS推定量 \hbetaの信頼度が下がる(分散が大きくなる)ことを、数式とシミュレーションの両面から見てきました。
定量の信頼度が低いと、回帰分析を用いた意思決定の信頼度が下がることに繋がります。

そこで、推定量の分散を抑える手法として、リッジ回帰(Redge Regression)という手法が知られています。
この手法では、最小二乗法を行う際の損失関数に、 \betaが大きいほど強い罰則を与える罰則項を加えます。
この罰則によって、推定量は極端に大きな値を取らないような制限がかかり、推定が安定します。

\begin{align}
\hbeta(\lambda) &= \arg\min_{\beta} \paren{\vY - \vX\beta}'\paren{\vY - \vX\beta} + \lambda \beta'\beta\\
&= \inv{\vX'\vX + \lambda \vI_J}\vX'\vY
\end{align}

このような罰則を加えた推定量はリッジ推定量と呼ばれています*3
ここで、 \lambdaは罰則の強さを表しています。
 \beta'\beta = \|\beta\|^2_2はL2ノルムによる罰則と解釈できることから、これはL2正則化とも呼ばれています。

OLS推定量と比較して、 \lambda \vI_Jの項目が分母にある分、リッジ推定量のほうが(絶対値でみて)小さい値になります。
 \lambdaが大きいほど強い罰則がかかり、 \betaが0方向に引っ張られます。
実際 \lambdaを無限大に飛ばすと推定量は0になります。
逆に、 \lambda = 0のときは、リッジ推定量はOLS推定量と一致します。

リッジ推定量の性質

それでは、リッジ推定量の性質を確認するため、その期待値と分散を計算します。

変換行列の準備

後の話を簡単にするため、OLS推定量をリッジ推定量に変換する行列 \mWを定義しておきます。

\begin{align}
\mW\hbeta &=\hbeta(\lambda)
\end{align}

行列 \mWの中身を知りたいので、OLS推定量とリッジ推定量を代入します。

\begin{align}
\mW\hbeta &=\hbeta(\lambda) \\
\mW \inv{\vX'\vX}\vX'\vY&= \inv{\vX'\vX + \lambda \vI_J}\vX'\vY
\end{align}

両辺に \inv{\vX'\vY}\paren{\vX'\vX}を右からかけてあげると \mWの中身がわかります。

\begin{align}
\mW \inv{\vX'\vX}\vX'\vY\inv{\vX'\vY}\paren{\vX'\vX} &= \inv{\vX'\vX + \lambda \vI_J}\vX'\vY\inv{\vX'\vY}\paren{\vX'\vX}\\
\mW&= \inv{\vX'\vX + \lambda \vI_J}\paren{\vX'\vX}
\end{align}

リッジ推定量の期待値

さて、この変換行列 \mWを用いてリッジ推定量の期待値を計算します。

\begin{align}
\E{\hbeta(\lambda) \mid \vX}
&= \E{\mW\hbeta \mid \vX}\\
&= \mW\E{\hbeta \mid \vX}\\
&= \inv{\vX'\vX + \lambda \vI_J}\paren{\vX'\vX}\beta\\
&\leqslant \beta
\end{align}

OLS推定量とは異なり、リッジ推定量の期待値は真の回帰係数 \betaと一致しません。
分母の \lambda \vI_Jの分だけ0方向にバイアスがかかることがわかります。
よって、罰則の強さ \lambdaを大きくすればするほどバイアスは大きくなります。

定量にバイアスがかかるのはあまりうれしくない性質です。
ただ、バイアスの方向がゼロに向かうのは、実務的にはそんなに悪くない性質だと思っています。
というのも、ゼロ方向にバイアスがかかるというのは、説明変数が目的変数に与える効果を「過小評価」するという話なので、リッジ回帰で過小評価されていてもなお十分に効果があるなら、保守的に見てもそのくらいの効果がある可能性が高い、という判断が可能になるからです。

リッジ推定量の分散

次に、リッジ推定量の分散について確認します。


\begin{align}
    \var{\hbeta(\lambda) \mid \vX} 
    &= \var{\mW\hbeta \mid \vX}\\
    &= \mW\var{\hbeta \mid \vX}\mW'\\
    &= \inv{\vX'\vX + \lambda \vI_J}\paren{\vX'\vX}\inv{\vX'\vX}\sigma^2\paren{\vX'\vX}\inv{\vX'\vX + \lambda \vI_J}\\
    &= \inv{\vX'\vX + \lambda \vI_J}\paren{\vX'\vX}\inv{\vX'\vX + \lambda \vI_J}\sigma^2
\end{align}

 \mWは分母のほうが大きいので、リッジ推定量の分散はOLS推定量の分散より小さくなることがわかります。また、 \lambdaを大きくすれば分散が小さくなることも見て取れます。

まとめると、リッジ推定量は、バイアスがかかることを代償に推定量の信頼性を高める手法であると解釈できます。
このバイアスと信頼性の関係はBias-Variance tradeoffとして機械学習統計学の分野で知られています。

シミュレーションによる確認

リッジ推定量についても、罰則の強さによって期待値と分散がどのように変化するかシミュレーションで確認しておきます。

まずは期待値について確認します。サンプルサイズを200、相関係数を0.99で固定して、罰則の強さ \lambdaだけを変化させた場合のリッジ推定量の期待値を可視化します。
f:id:dropout009:20210722194802p:plain

回帰係数はすべて1だったことを念頭に置いて上図を確認すると、 \lambdaを大きくするにつれてリッジ推定量は0方向にひっぱられてバイアスが強くなっていくことがわかります(サンプルサイズが小さく多重共線性が強いので、罰則ゼロのときそもそもあんまりうまく推定できていないことも見て取れます)。

同様の設定でリッジ推定量の分散を確認します。
f:id:dropout009:20210722194808p:plain


罰則を強くすると推定量の分散が小さくなっていく様子が見て取れます。
このように、リッジ推定量はバイアスという犠牲を払って信頼性を高めることができます。

まとめ

この記事では、多重共線性がOLS推定量の信頼性に与える影響について、数式とシミュレーションの両面から確認しました。
大雑把に言うと、完全な多重共線性があるケースを除いて、興味のある説明変数以外の多重共線性はあまり気にしなくて良いことがわかりました。

また、推定量の信頼性を高める手法の一つとして、リッジ推定量を紹介しました。
リッジ推定量はゼロ方向にバイアスがかかりますが、推定量の分散を小さくできることを確認しました。
リッジ推定量には、例えば説明変数の数 Jがサンプルサイズ Nよりも大きくても推定が可能になるなど、他にも便利な性質がいくつかあります。リッジ推定量の詳細はvan Wieringen(2015)をご確認ください。


この記事で利用したシミュレーションコードは以下になります。
github.com

参考文献

*1:僕が受講した計量経済学のコースワークは市村英彦先生が担当されていたのですが、講義ノートにはしっかり記述がありました。

*2:Hansen(2011)を特に参考にしています。

*3:リッジ推定量の利点として、完全な多重共線性が存在するケースでも逆行列がとれるので推定が可能になることが挙げられたりしています。それは実際そうなんですが、完全な多重共線性があってしかも回帰係数に興味がある場合はリッジ回帰で対処するのではなく手作業で完全な多重共線性を排除したほうがいいと思います。

Influence Functionでインスタンスの重要度を解釈する


 \def\paren#1{\left(#1\right)}
 \def\sumi{\sum_{i=1}^N}
\def\htheta{\hat{\theta}}
\def\hthetami{\htheta_{- i}}
\def\hthetamie{\hthetami(\epsilon)}
\def\nablat{\nabla_{\theta}}
\def\nablatt{\nabla_{\theta}^2}
\def\evale#1{\left.#1\right|_{\epsilon=0}}

モチベーション

Permutation Feature Importance(PFI)、Partial Dependence(PD)、Accumulated Local Effects(ALE)などの機械学習の解釈手法は、特徴量とモデルの予測値の関係を通じて機械学習モデルを解釈する手法でした*1
たとえば、PFIでは「どの特徴量がモデルにとって重要なのか」という視点でブラックボックスモデルに解釈性を与えることが出来ます。
特徴量重要度の用途の一つはデバッグです。
ドメイン知識から考えると不自然な重要度をもつ特徴量が存在すれば、データやモデルに不備があるのかもしれません。
また、重要度が極端に高い特徴量は、何か情報がリークしている可能性があります。


このように、「どの特徴量がモデルにとって重要なのか」を知ることはモデルの振る舞いを解釈する上で重要です。
ここでは特徴量に注目していますが、視点を変えてインスタンスに注目することもできます。
つまり、「どのインスタンスがモデルにとって重要なのか」という解釈性も、機械学習モデルの振る舞いを理解する上で役に立ちそうです。
このような解釈性を、特徴量重要度に対比してインスタンス重要度とこの記事では呼ぶことにします。
特徴量重要度と同様にモデルのデバッグを考えると、極端に重要度の高いインスタンスが存在すれば、そのインスタンスを外れ値として除外する、そもそもデータとしておかしいのではないか確認するなどの対応が取れます。

インスタンス重要度の定義は?

「どのインスタンスがモデルにとって重要なのか」を知ることには意味がありそうですが、これを知るためにはそもそも何をもって重要とみなすかをきちんと定義する必要があります。
できればどんなモデルに対しても重要度を計算したいので、すべてのモデルが共通に持っている要素である「予測値」に注目します。つまり、「そのインスタンスの有無でモデルの予測がどのくらい大きく変わるか」をもって重要度を定義します。

具体例として、下図を見てください。
右下の外れ値を学習データに入れた場合と入れなかった場合での線形回帰の予測線を可視化しています。
外れ値を入れると回帰直線の傾きが緩やかになることが見て取れます。

f:id:dropout009:20210718020515p:plain

上記画像とシミュレーションデータを作成するためのコードは以下になります。

def generate_simulation_data(
    n_incetances: int,
) -> tuple[np.ndarray, np.ndarray]:
    """シミュレーションデータを生成する"""

    x = np.random.normal(0, 1, n_incetances)
    u = np.random.normal(0, 0.1, n_incetances)

    y = x + u

    # 外れ値で上書きする
    x[0] = 2
    y[0] = -2

    return (x, y)


x, y = generate_simulation_data(n_incetances=20)
def draw_regression_line(x: np.ndarray, y: np.ndarray) -> None:
    """外れ値を含んだ場合と含まない場合の回帰直線を可視化する"""

    fig, ax = plt.subplots()
    sns.regplot(x, y, scatter_kws=dict(alpha=1), ci=None, ax=ax)
    sns.regplot(x[1:], y[1:], scatter_kws=dict(alpha=1), color=".7", ci=None, ax=ax)
    ax.set(xlabel="X", ylabel="Y")
    fig.suptitle("外れ値による回帰直線の変化")
    fig.savefig(f"output/outlier.png")
    fig.show()


draw_regression_line(x, y)

Leave One Outによるインスタンス重要度の計算

このように、そのインスタンスを学習データに含めることでどれだけ予測値が変化するのか、そしてその結果としてどれだけ予測精度が変化するのかを計算し、それをインスタンスの重要度とすることができそうです。
この重要度計算アルゴリズムは、インスタンスをひとつずつ取り除いて(Leave One Out, LOO)繰り返し学習と予測を行うだけなので、簡単に実装することが出来ます。

def calc_leave_one_out_error_diff(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    """LOOでの予測誤差の差分を計算する"""

    n_incetances = x.shape[0]
    X = x.reshape(-1, 1)
    y_pred = LinearRegression().fit(X, y).predict(X)

    error_diff = np.zeros((n_incetances, n_incetances))
    for i in range(n_incetances):
        X_loo = np.delete(X, i, 0)
        y_loo = np.delete(y, i, 0)

        y_pred_loo = LinearRegression().fit(X_loo, y_loo).predict(X)
        error_diff[:, i] = (y - y_pred_loo) ** 2 - (y - y_pred) ** 2

    return error_diff


looed = calc_leave_one_out_error_diff(x, y)
avg_diff = looed.mean(0) # 平均的な誤差を計算

計算結果を可視化してみると、先程の画像で右下に位置していたインスタンス0がモデルの予測に最も大きな影響を与えていたことがわかります。

def draw_importance(imp: np.ndarray) -> None:
    """インスタンスごとの重要度を可視化する"""

    df = (
        pd.DataFrame({"importance": imp})
        .reset_index()
        .rename({"index": "instance"}, axis=1)
        .sort_values("importance")
        .reset_index(drop=True)
        .assign(instance=lambda df: df.instance.astype(str))
    )
    
    fig, ax = plt.subplots()
    ax.barh(df["instance"], df["importance"])
    ax.set(xlabel="重要度", ylabel="インスタンスNo.")
    fig.suptitle("インスタンスごとの重要度")
    fig.savefig(f"output/imp.png")
    fig.show()


draw_importance(avg_diff)

f:id:dropout009:20210718020459p:plain

Influence Functionでインスタンス重要度を近似する

記法の導入

LOOによる予測誤差の変化は直感的でわかりやすいですが、インスタンス重要度の計算に再学習が必要なので計算に時間がかかることが欠点です。
今回のシミュレーションではインスタンスの数が20でモデルは線形回帰だったのですぐに計算は終わりましたが、実務の際にはより大量のデータをより複雑なモデルで分析します。
ですので、より複雑な問題に取り組む際には、再学習を必要としないアルゴリズムが必要になります。

実は、Influence Functionを利用して、再学習を行わずにインスタンス重要度を近似計算するアルゴリズムがKoh and Liang(2017)で提案されています。
ここで、このアルゴリズムを記述するためにいくつかの記法を導入します。

  • 特徴量とアウトカムの組を z = (x, y)とします
  • モデルはパラメトリックなモデルを採用し、そのパラメータを \thetaで表します。パラメータは K個あるとします。つまり、 \theta = (\theta_1, \dots, \theta_K)とします
  • 損失関数を L(z, \theta) = L(y, f(x, \theta))で定義します。モデルが線形回帰で損失関数が二乗誤差なら L(z, \theta)=\paren{y - x^\top\theta}^2です
  • データに対する予測誤差の総和を R(\theta) = \sumi L(z_i, \theta)で表します。モデルが線形回帰で損失関数が二乗誤差なら R(\theta)=\sumi\paren{y_i - x_i^\top\theta}^2です

さて、モデルは予測誤差が最も小さくなるようにパラメータ \thetaを選びます。
つまり、全インスタンスを学習した場合のパラメータ \htheta
\begin{align}
\htheta = \arg\min_{\theta \in \Theta}\sumi L(z_i, \theta)
\end{align}
であり、インスタンス iを取り除いたデータで学習した場合のパラメータ \hthetami
\begin{align}
\hthetami = \arg\min_{\theta \in \Theta} \sum_{j \neq i} L(z_j, \theta)
\end{align}
となります。

この記法を用いると、インスタンス iを使わずにモデルの学習を行った場合に、インスタンス jの予測誤差に与える影響は以下のように表現できます。
\begin{align}
L\paren{z_j, \hthetami} - L\paren{z_i, \htheta}
\end{align}
ここで、

を表します。この差分が大きいほどインスタンスの重要度は高くなります。

Influence Functionの導出

前述のように、このようなLOOによるインスタンス重要度の計算はモデルの再学習が必要なので時間がかかるという問題がありました。
これを克服するためには、モデルの再学習を行わなずに再学習を行った場合の結果を近似できるような手法を考える必要があります。

そこで、Koh and Liang(2017)では以下のような最小化問題を導入しました。
\begin{align}
\hthetamie = \arg\min_{\theta \in \Theta} \sumi L(z_j, \theta) - \epsilon L(z_i, \theta)
\end{align}

ここで、 \epsilon L(z_i, \theta)という項が追加されていることに注意してください。
この項に注目すると、以下の性質があることが見て取れます。

  •  \epsilon = 0のとき、全インスタンスを用いた最小化問題と一致する。つまり、 \hthetami(0) = \htheta
  •  \epsilon = 1のとき、インスタンス iを取り除いた最小化問題と一致する。つまり、 \hthetami(1) = \hthetami

 \epsilon \in (0, 1)の場合は、インスタンス iを完全に取り除いた(重みをゼロにした)状態や、全データを用いた状態ではなく、インスタンス iの重みをちょっとだけ( \epsilonだけ)軽くした状態を考えています。
よって、 \epsilonインスタンス iの重みをコントロールするパラメータと解釈できます。重みが1のときはインスタンス iを取り除いてモデルの学習を行うことに、重みが0のときは全インスタンスを用いてモデルの学習を行うことに相当します。


この最小化問題の解として与えられる \hthetamieを用いた場合、任意の z = (x, y)に対する予測誤差は L(z, \hthetami(\epsilon))です。
ここで、 \epsilonを0から少しだけ変化させたときに、どのくらい誤差が変化するかを考えます。
なんでこんなことを考えるのかを平たく言うと以下になります。

  • ほんとは \epsilon = 0の場合と \epsilon = 1の場合の予測誤差の変化を見たい
  • →それは再学習が必要になって大変
  • →それなら \epsilonが0からちょっとだけ変化したときの予測誤差の変化を使って、 \epsilonを1単位変化させた場合の( \epsilon = 0から \epsilon = 1まで持っていた場合の)予測誤差の変化を近似しよう。 \epsilon = 1のときのパラメータはインスタンス iを取り除いたときの予測誤差なので、この近似によってさも再学習をしたかのようなパラメータを推定できる

上記の操作は予測誤差 L(z, \hthetami(\epsilon)) \epsilon微分して \epsilon=0で評価すれば実行できます。
\begin{align}
\evale{\frac{d L(z, \hthetami(\epsilon))}{d \epsilon}}
&= \evale{\nablat L(z, \hthetamie)^\top \frac{d \hthetamie}{d\epsilon}}\\
&= \evale{\nablat L(z, \htheta)^\top \frac{d \hthetamie}{d\epsilon}}
\end{align}

ここで、 \evale{\hthetami(\epsilon)} = \hthetami(0) = \hthetaであることを利用しています。
上式の \nablat L(z, \htheta)は全データで学習したモデルを用いた予測誤差で求まります。
ですので、残りの \evale{\frac{d \hthetamie}{d\epsilon}}をどうやって求めるかがわかれば、これを計算できます。

パラメータの性質

それでは \hthetamieの性質を調べます。そもそも \hthetamieがどうやって計算されていたかを思い出すと、
\begin{align}
\hthetamie = \arg\min_{\theta \in \Theta} \sumi R(\theta) - \epsilon L(z_i, \theta)
\end{align}
でした。
この後の操作を可能にするために、 R(\theta) \thetaに関してstrongly convexで2階微分可能であることを仮定します*2

 \hthetamieは目的関数を最小化する値なので、1階条件
\begin{align}
\nablat R(\hthetamie) - \nablat L(z_i, \hthetamie)\epsilon = 0
\end{align}
を満たします。

ここで、左辺を \hthetaの周りで1階のテイラー展開すると以下を得ます。

\begin{align}
\paren{\nablat R(\htheta) - \nablat L(z_i, \htheta)\epsilon} + \paren{\nablatt R(\htheta) - \nablatt L(z_i, \htheta)\epsilon}\paren{\hthetamie - \htheta} + \mathrm{residual}
\end{align}

ところで、 \htheta R(\theta)を最小化する値なので、 \nablat R(\htheta) = 0です。ついでに剰余項を無視して近似すると

\begin{align}
0 \approx - \nablat L(z_i, \htheta)\epsilon + \paren{\nablatt R(\htheta) - \nablatt L(z_i, \htheta)\epsilon}\paren{\hthetamie - \htheta}
\end{align}

なので、整理すると

\begin{align}
\paren{\hthetamie - \htheta} \approx \paren{\nablatt R(\htheta) - \nablatt L(z_i, \htheta)\epsilon}^{-1}\nablat L(z_i, \htheta)\epsilon
\end{align}

となります。

 \epsilon 0からちょっとだけ動かしたときに \hthetamie \hthetaからどれだけ離れるかを知りたいので、両辺を \epsilon微分して \epsilon=0で評価します。


\begin{align}
\evale{\frac{d\paren{\hthetamie - \htheta}}{d\epsilon}} \approx& \evale{\paren{\nablatt R(\htheta) - \nablatt L(z_i, \htheta)\epsilon}^{-1}\nablat L(z_i, \htheta)}\\
&+\evale{\paren{\nablatt R(\htheta) - \nablatt L(z_i, \htheta)\epsilon}^{-2}\nablatt L(z_i, \htheta)\nablat L(z_i, \htheta)\epsilon}\\
\evale{\frac{d\hthetamie}{d\epsilon}} &\approx \nablatt R(\htheta)^{-1}\nablat L(z_i, \htheta)
\end{align}


これを元の式に戻すと

\begin{align}
\evale{\frac{d L(z, \hthetami(\epsilon))}{d \epsilon}} \approx \nablat L(z, \htheta)^\top\nablatt R(\htheta)^{-1}\nablat L(z_i, \htheta)
\end{align}

を得ます。上式は \epsilonを1単位動かしたときの予測誤差の変化を1階近似していると解釈できます。 \epsilon = 0 \epsilon = 1の差分は、全データを用いた場合と、インスタンス iを取り除いた場合の予測誤差の差分でした。よって、インスタンス iを取り除いたことによる予測誤差の変化を近似していると解釈できます。

これはInfluence Functionと呼ばれています。インスタンス iを取り除いたときの z=(x, y)に対する予測誤差の変化を

\begin{align}
\mathcal{I}_i(z) = \nablat L(z, \htheta)^\top\nablatt R(\htheta)^{-1}\nablat L(z_i, \htheta)
\end{align}

と表記することにします。

線形回帰モデルの場合

線形回帰モデルのInfluence Function

さて、抽象的な議論が続いてしまったので、具体的に特徴量が1つだけの単純な線形回帰モデル
\begin{align}
y = x\theta + u
\end{align}
を考えます。損失関数は二乗誤差とします。
このとき、Influence Functionに必要なパーツはそれぞれ以下で計算できます。

\begin{align}
L(z, \theta) &= (y - x\theta)^2, \\
\nablat L(z, \theta) &= -2(y - x\theta)x,\\
R(\theta) &= \sumi(y_i - x_i\theta)^2,\\
\nablat R(\theta) &= -2\sumi(y_i - x_i\theta)x_i,\\
\nablatt R(\theta) &= 2\sumi x_i^2,\\
\htheta &= \frac{\sumi x_iy_i}{\sumi x_i^2}
\end{align}

これをInfluence Functionに代入することで、線形回帰モデルに対して再学習無しでインスタンスの重要度を近似することが出来ます*3

\begin{align}
\mathcal{I}_i(z) &= \nablat L(z, \htheta)^\top\nablatt R(\htheta)^{-1}\nablat L(z_i, \htheta)\\
&= \paren{-2\paren{y - x\htheta}x} \paren{2\sumi x_i^2}^{-1} \paren{-2\paren{y_i - x_i\htheta}x_i}\\
&= \frac{2\paren{\paren{y - x\htheta}x}\paren{\paren{y_i - x_i\htheta}x_i}}{\sumi x_i^2}
\end{align}

シミュレーション

Influence Functionが特定できたので、後はLOOで求めた予測誤差の差分をInfluence Functionでうまく近似できるのかをシミュレーションで確認します。

まずはシミュレーションデータを生成します。
できるだけ単純なデータにしたいので、以下のような特徴量が1つの設定を考えます。
\begin{align}
y &= x + u,\\
x &\sim \mathcal{N}(0, 1),\\
u &\sim \mathcal{N}(0, 0.1^2)
\end{align}

def generate_simulation_data(
    n_incetances: int,
) -> tuple[np.ndarray, np.ndarray]:
    """シミュレーションデータを生成する"""
    x = np.random.normal(0, 1, n_incetances)
    u = np.random.normal(0, 0.1, n_incetances)

    y = x + u

    return (x, y)


x, y = generate_simulation_data(n_incetances=100)


シミュレーションデータに対して、LOOによる予測誤差の変化を計算します。

def estimate_theta(x: np.ndarray, y: np.ndarray) -> float:
    """単回帰モデルの係数を推定する"""

    return (x @ y) / (x @ x)


def calc_leave_one_out_error(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    """LOOでの予測誤差の差分を計算する"""
    
    # ベースラインの全データを用いた学習
    n_incetances = x.shape[0]
    theta_hat = estimate_theta(x, y)
    
    # インスタンスのをひとつずつ取り除いて学習
    error = np.zeros((n_incetances, n_incetances))
    for i in range(n_incetances):
        x_loo = np.delete(x, i, 0)
        y_loo = np.delete(y, i, 0)

        theta_hat_loo = estimate_theta(x_loo, y_loo)
        error[i, :] = (y - x * theta_hat_loo) ** 2 - (y - x * theta_hat) ** 2

    return error


looe = calc_leave_one_out_error(x, y)||<

次に、Influence Functionを用いた近似を行います。
>|python|
def calc_influence(
    x: np.ndarray, y: np.ndarray, x_train: np.ndarray, y_train: np.ndarray
) -> np.ndarray:
    """influenceによる近似"""
    
    # 単回帰の係数
    theta_hat = estimate_theta(x_train, y_train)
    
    # パーツごとに計算
    R2 = 2 * (x_train @ x_train)
    L_train = -2 * (y_train - x_train * theta_hat) * x_train
    L = -2 * (y - x * theta_hat) * x

    return L.reshape(-1, 1) @ L_train.reshape(1, -1) / R2


influence = calc_influence(x, y, x, y)

最後に、LOOとInfluenceを比較します。

def draw_scatter_with_line(x: np.ndarray, y: np.ndarray) -> None:
    """LOOとInfluenceの結果を比較する"""
    
    xx = np.linspace(x.min(), x.max())
    
    fig, ax = plt.subplots()
    ax.plot(xx, xx, c='.7', linewidth=2, zorder=1)
    ax.scatter(x, y, zorder=2)
    ax.set(xlabel="Influence", ylabel="LOO")
    fig.suptitle("InfluenceとLOOの比較")
    fig.savefig(f"output/compare.png")
    fig.show()
    
draw_scatter_with_line(influence.flatten(), looe.flatten())


f:id:dropout009:20210719224841p:plain

ひとつひとつのデータ点が、インスタンス iを取り除いた際のインスタンス jに対する予測誤差の変化を表しています。
点はほぼ y = x上の線の上に乗っているので、LOOの結果をInfluence Functionでうまく近似できていることがわかります。

まとめ

この記事では、インスタンス重要度を計算するための手法として、LOOによる正確な計算とInfluence Functionによる近似を紹介しました。

  • インスタンスの重要度を知ることでモデルの振る舞いをより深く理解することができるだけでなく、モデルのデバッグなどに利用することもできる
  • LOOは正確な反面、計算に時間がかる。Influenceは計算コストを抑えられるが、あくまでも近似であることに注意が必要。特に、モデルが複雑になると近似の精度が下がる可能性がある
  • LOOはどんなモデルに対しても適用できるが、Influenceはパラメトリックなモデルに対してのみ適用できる。つまり、線形回帰モデルやNeural Netには適用できるが、Random ForestやGradient Boosting Decision Treeには適用できない

最後に、この記事で使ったコードは以下のgithubにアップロードしています。
github.com

参考文献

*1:特徴量に注目した機械学習の解釈手法については、拙著「機械学習を解釈する技術」で詳しく説明しています(宣伝)。 www.amazon.co.jp

*2:モデルが線形回帰で損失が二乗誤差だとこれを満たしますが、Neural Netなどより複雑なモデルだと条件を満たさない可能性はあります。その場合、パラメータが局所最適解に落ちて近似があんまりうまくいかない可能性があります。

*3:本当のことを言うと、線形回帰モデルは再学習無しで正確なインスタンスの重要度を計算することができます。詳細はHansen(2021)をご確認ください

SHAP(SHapley Additive exPlanations)で機械学習モデルを解釈する


はじめに

ブラックボックスモデルを解釈する手法として、協力ゲーム理論のShapley Valueを応用したSHAP(SHapley Additive exPlanations)が非常に注目されています。SHAPは各インスタンスの予測値の解釈に使えるだけでなく、Partial Dependence Plotのように予測値と変数の関係をみることができ、さらに変数重要度としても解釈が可能であるなど、ミクロ的な解釈からマクロ的な解釈までを一貫して行える点で非常に優れた解釈手法です。

SHAPの論文の作者によって使いやすいPythonパッケージが開発されていることもあり、実際にパッケージを使った実用例はたくさん見かけるので、本記事では協力ゲーム理論の具体例、Shapley Valueのコンセプトと求め方、機械学習モデルを解釈するためのShapley Valueの使われ方を意識してまとめました。


※この記事をベースに2020年1月16日に行われたData Gateway Talk vol.5にて発表した資料は以下になります
speakerdeck.com

この記事で書いていること、書いていないこと

書いていること

  • 協力ゲーム理論の具体例とShapley Valueのコンセプトと求め方
  • 機械学習モデルを解釈するためのShapley Valueの使われ方
  • shapパッケージの簡単な使い方

書いていないこと

アルバイトゲームとShapley Value

ここでは協力ゲーム理論のアルバイトゲームを例にとって、Shapley Valueを直感的に理解することを目指します。

まずはアルバイトゲームについて説明します。*1
アルバイトは単独で行うこともできますし、他の人とチームを組んでやってもいいとします。

  • まずは、1人で働いた場合は、A君が1人でやると6万円、B君が1人でやると4万円、C君が1人でやると2万円がもらえるとしましょう。
  • 次に、2人で働いた場合は、A君とB君が2人でやったときは合計で20万円、A君とC君が2人でやったときは合計で15万円、B君とC君が2人でやったときは合計で10万円がもらえるとしましょう。
  • 最後に、A君B君C君の3人で働いた場合は、合計で24万円がもらえるとしましょう。

これをまとめると以下の表のようになります。

参加者 報酬
A君 6
B君 4
C君 2
A君B君 20
A君C君 15
B君C君 10
A君B君C君 24

今、A君B君C君の3人全員で働いて得た報酬24万円をどうやって分配するのが尤もらしいかを考えます。
直感的には、より貢献度の高い人による多くの報酬を分配するのがフェアな分配のひとつになりそうです。*2
とすると、問題は各人の貢献度をどうやって計るのかということになります。ここで限界貢献度という概念を導入しましょう。これは「各人がアルバイトに参加したときに追加的にどのくらい報酬が増えるか」で計算されます。
たとえば、A君について考えると、限界貢献度は以下のように計算されます。

  • 「誰もいない」→「A君のみ」だと6 - 0 = 6万円
  • 「B君のみ」→「A君とB君」だと20 - 4 = 16万円
  • 「C君のみ」→「A君とC君」だと15 - 2 = 13万円
  • 「B君とC君」→「A君とB君とC君」だと24 - 10 = 14万円

ここでわかるのは、A君の限界報酬はA君が参加する順番に依存するということです。
この影響を打ち消すため、考えられる全ての参加順を用いて平均的な限界貢献度を求めることにしましょう。
発生し得る参加順と、その場合の各人の限界貢献度は以下のようにまとめられます。

参加順 A君の限界貢献度 B君の限界貢献度 C君の限界貢献度
A君→B君→C君 6 14 4
A君→C君→B君 6 9 9
B君→A君→C君 16 4 4
B君→C君→A君 14 4 6
C君→A君→B君 13 9 2
C君→B君→A君 14 8 2

参加順は3!=6通りなので、

  • A君の平均的な限界貢献度は(6 + 6 + 16 + 14 + 13 + 14) / 6 = 11.5万円
  • B君の平均的な限界貢献度は(14 + 9 + 4 + 4 + 9 + 8) / 6 = 8万円
  • C君の平均的な限界貢献度は(4 + 9 + 4 + 6 + 2 + 2) / 6 = 4.5万円

この平均的な限界貢献度のことをShapley Valueと言い、このShapley Valueを用いて報酬を分配しよう、というのがある意味で尤もらしい分配の方法になります。実際、11.5 + 8 + 4.5 = 24万円できれいに分配ができていますし、より貢献度が高い人により多くの報酬が渡るという意味でフェアな分配にもなっています。*3

今回はゲームのプレイヤーがA君B君C君3人のケースを考えました。より一般的なケースとしてN人のプレイヤー\mathcal{N} = \{1, 2, \dots, N\}がゲームに参加するケースを考えると、プレイヤーiのShapley Value\phi_iは以下で計算できます。

\begin{align*}
\phi_i =
\sum_{\mathcal{S} \subset\mathcal{N} \setminus \{i\}}
\frac{S!(N - S - 1)!}{N!}
\bigg(v(\mathcal{S}\cup \{i\}) - v(\mathcal{S})\bigg)
\end{align*}

  • ここで、\mathcal{S} \mathcal{N}からプレイヤーiを除いたプレイヤーの組み合わせです。たとえばA君B君C君3人のケースで、A君のShapley Valueを考える場合は、 \{\emptyset\}, \{B\}, \{C\}, \{B, C\}が該当します。
  • S\mathcal{S}の要素の数、つまりプレイヤーの数になります。先程の例だと、それぞれ0人, 1人, 1人, 2人となります。
  •  v(\cdot)は効用を表す関数です。なので、v(\mathcal{S}\cup \{i\}) - v(\mathcal{S})はプレイヤー iが参加しているときと参加していないときでの効用の差となります。つまり、プレイヤー iが参加したときの限界貢献度となります。

まとめると、プレイヤー iが参加することの限界貢献度を出現する全ての組合わせで求めて平均しています。これはアルバイトゲームの具体例でやったことと全く同じ操作をしています。上式だけみると何を計算しているかよくわからないのですが、具体例をいれて確かめてみると何をやっているかわかりやすいんじゃないかと思います。

機械学習モデルへの応用

最近、協力ゲーム理論のShapley Value機械学習モデルの予測結果を解釈するために利用しよう、という研究が発達してきました。
これらの研究では、モデルに投入したひとつひとつの特徴量をゲームのプレイヤーと見立てて、各特徴量の予測への貢献度をShapley Valueで測ろう、ということをやっています。
具体的なシチュエーションを考えてみましょう。今、特徴量として \mathbf{X} = (X_1, X_2, X_3)の3つがあるとします。
モデルを f(\cdot)とすると、平均的な予測値は E[f( \mathbf{X})]となります。
ここで、ひとつのインスタンスを取り出すと、それぞれ \mathbf{x} = (x_1, x_2, x_3) という値をとっているとします。このとき、予測値としては f(\mathbf{x})が出力されます。

この、平均的な予測値 E[f( \mathbf{X})]と各インスタンスの予測値 f(\mathbf{x})の乖離に対して、各特徴量がどのくらい影響しているのかを調べます。

インスタンスの予測値 f(\mathbf{x})
\begin{align*}
E[f( \mathbf{X} | X_1 = x_1, X_2 = x_2, X_3 = x_3)] & = f( x_1, x_2, x_3 ) = f(\mathbf{x})
\end{align*}
なので、平均的な予測値 E[f( \mathbf{X})]から X_1, X_2, X_3を条件付けていくことで、その特徴量を知ることが各インスタンスの予測に対してどのように影響するかを知ることができます。具体的に、 X_1, X_2, X_3の順で条件づけていくとすると、たとえば下図のような推移が見られます。

f:id:dropout009:20191117175258p:plain

ここで、 \phi_j, j = 1, 2, 3を各特徴量が予測値に与える限界的な効果としています。(なお、 \phi_0は0と平均的な予測値の乖離で、特に重要ではありません。)

まず、特に何の情報もないところから X_1 = x_1という情報を得ると、予測値が \phi_1だけ大きくなります。さらにその状態から X_2 = x_2という情報を得ると、予測値がさらに \phi_2だけ大きくなります。最後に、 X_3 = x_3という情報を得ると、予測値が \phi_3だけ小さくなって、これが最終的なこのインスタンスへの予測結果となります。

ここで、先程のアルバイトゲームとの共通点が現れています。

  • 各特徴量を一つずつ条件付けていくことで予測値に与える限界的な効果を見ていますが、これはアルバイトゲームで言うところの限界貢献度に対応しています。
  • また、今回は X_1, X_2, X_3の順で条件づけていますが、当然別の順番で条件づけていくと予測値に与える限界的な効果は変化します。よって、考えうる全ての順序で限界効果を計算し、それを平均しなければなりません。この平均的な限界効果がShapley Valueに対応します。

詳細は論文をご確認頂きたいのですが、このようにShapley Valueを用いて計算された貢献度は、Shapley Valueが持つ望ましい性質を満たすことが証明されています。
ただ、Shapley Valueを計算するのは非常に計算コストがかかるので、実際の計算ではアルゴリズムごとの近似手法が用いられています。特に、tree系のアルゴリズムに対してはTreeSAHPという高速なアルゴリズムが提案されています。

機械学習モデルに対するShapley Valueを簡単で高速に計算できるPythonパッケージにSHAPがあります。

まずは必要なパッケージを読み込みましょう。

import numpy as np
import pandas as pd

# モデルはRandom Forestを使う
from sklearn.ensemble import RandomForestRegressor

# SHAP(SHapley Additive exPlanations)
import shap
shap.initjs() # いくつかの可視化で必要

データはshapパッケージにあるボストンの不動産価格のデータを用います。
データセットの細かい説明はリンクをご確認下さい。

X, y = shap.datasets.boston()
X.head()

f:id:dropout009:20191117191347p:plain


モデルはRandom Forestを使います。

model = RandomForestRegressor(n_estimators=500, n_jobs=-1)
model.fit(X, y)


ここからがshapの使い方になります。shapにはいくつかのExplainerが用意されていて、まずはExplainerにモデルを渡すします。今回はRandom ForestなのでTreeExplainer()を使います。

explainer = shap.TreeExplainer(model, X)

Explainerにはshap_values()メソッドが用意されています。これにShapley Valueを計算したいインプットの行列を渡すことでShapley Valueを計算できます。shapは高速な計算アルゴリズムが実装されており、特にこのデータセットは500行程度なのですぐに計算は終わりますが、計算コストは依然として高いので大きいデータセットのときは適当にサンプリングして渡す必要があります。

shap_values = explainer.shap_values(X)


まずは一つのインスタンスに対してShapley Valueを確認してみます。これは`force_plot()`を使います。

i = 0
shap.force_plot(explainer.expected_value, shap_values[i,:], X.iloc[i,:])

f:id:dropout009:20191117190407p:plain

basevalueの22.51がが全体平均で、output valueの25.64がこのインスタンスに対する予測値となります。なので、このインスタンスは平均よりも高い値が予測されています。なぜこのような予測になったのかを説明するために、各特徴量がどのくらい大きな因子となっているのかを、Shapley Valueで分解して可視化しています。たとえば、このインスタンスはLSTATが4.98をとっていて、これはこのインスタンスの不動産価格に対して大きなプラスの要因となっています。逆に、RMは6.575となっていて、これはマイナス要因であることも見て取れます。プラスとマイナスを総合するとプラスの方が大きくなっていて、最終的には全体平均より3.13高い予測結果になっていることがわかります。
ちなみにLSTATはその地域に住む低所得層の割合、RMは平均的な部屋の数です。

ひとつひとつのインスタンスでShapley Valueを見ていくことでミクロな分析ができますが、よりマクロな分析として、Shapley Valueを変数ごとに平均して変数重要度のように使うこともできます。
今、データセットインスタンス数が Nとすると、変数 p = 1, \dots, Pの変数重要度 \mathrm{Feature Importance}_pは以下で計算します。

\begin{align*}
\mathrm{Feature Importance}_p = \frac{1}{N}\sum_{i = 1}^{N} \left|\phi_p^i \right|
\end{align*}

ただし、 \phi_p^iインスタンス iでの変数 pのShapley Valueです。また、プラスマイナスの影響を無視するために絶対値をとっています。
この変数重要度を簡単に可視化するための関数として、summry_plot()が用意されています。

shap.summary_plot(shap_values, X, plot_type="bar")

f:id:dropout009:20191117190121p:plain

Shapley Valueの意味で、LSTATとRMが非常に重要な変数であることが見て取れます。


plot_type = "bar"とすると棒グラフが出ますが、指定しないと一つ一つのShapley Valueがそのまま打たれます。*4

shap.summary_plot(shap_values, X)

f:id:dropout009:20191117190134p:plain

上に来るほど先程の棒グラフの意味で重要な変数になります。色が赤いほどその変数の値が高いとき、青いほど低いときのShapley Valueになります。Shapley Valueの分布に加えて、LSTATは低いほうがプラスの要因に、RMは高いほうがプラスの要因になりそうなことが見て取れます。




さらに、特徴量の値とShapley Valueの散布図を書くことで、Partial Dependence Plotのような可視化も可能となります。
これは dependence_plot()を使います。今回は特に重要度の高かったLSTATを見てみましょう。

shap.dependence_plot("LSTAT", shap_values, X)

f:id:dropout009:20191117190116p:plain

LSTATの値が大きくなるほどShapley Valueが小さくなることが見て取れます。
色付けは交互作用が見れるように他の変数の値でされています。デフォルトでは交互作用が一番はっきり現れる変数が自動で選ばれます。今回はDISが選ばれていて、DISの値が大きいほど赤く、小さいほど青い点が打たれます。DISはemployment centreからの距離を表しています。グラフの左半分を見ると、同じLSTATの値でもDISが短いほどShapley Valueが高くなる傾向が見て取れます。


*1:この数値例は岡田章『新板ゲーム理論』の例をそのままお借りしています。

*2:たとえば、山分けで3等分も考えられますが、この設定ではそれはうまくいきません。A君B君C君の3人で働いて3等分すると報酬は8万円になります。この場合、A君B君だけで働いて報酬を2等分するとA君とB君は10万円を稼ぐことができ、C君を外すインセンティブが生まれてしまいます。

*3:Shapley Valueは数学的に証明された望ましい性質をいくつかもっていますが、ここでは具体例と直感的な性質のみ説明しています。詳細を知りたい場合は協力ゲーム理論の教科書か論文をご確認下さい。

*4:僕はこれの正式な名称を知りません。sina plotでいいのかな?

tidymodelsとDALEXによるtidyで解釈可能な機械学習



※この記事をベースにした2020年1月25日に行われた第83回Japan.Rでの発表資料は以下になります。
speakerdeck.com


はじめに

本記事では、tidymodelsを用いて機械学習モデルを作成し、それをDALEXを用いて解釈する方法をまとめています。
DALEX

Collection of tools for Visual Exploration, Explanation and Debugging of Predictive Models

を目的とした複数のパッケージから成り立っています。
具体的にどんなパッケージが所属しているかはリンクを確認して下さい。本記事では、特にメインパッケージであるDALEXingredientsを用います。
ingredientsには、変数重要度、Partial Dependence Plot(PDP)、Individual Conditional Expectation(ICE)、Accumulated Local Effect(ALE)など、特徴量とターゲット変数の関係を見るための関数が用意されています。

なお、本記事ではtidymodelsの使い方には触れないので、過去の記事を参考にして頂ければと思います。

dropout009.hatenablog.com
dropout009.hatenablog.com
dropout009.hatenablog.com

また、本記事では解釈手法そのものには詳しく触れません。Permutationベースの変数重要度やPDP, ICEの詳細な説明は以下の記事も参考にして頂ければと思います。
dropout009.hatenablog.com

パッケージ

本記事で用いるパッケージは以下になります。

library(tidymodels)
library(DALEX) #解釈
library(ingredients) #解釈

library(distributions3) #シミュレーション用
library(colorblindr) #可視化用

set.seed(42)

また、可視化用に以下の関数を用意しておきます

# visualization -----------------------------------------------------------
# ggplot2テーマ
theme_scatter = function() {
  theme_minimal(base_size = 12) %+replace%
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(color = "gray", size = 0.1),
          legend.position = "top",
          axis.title = element_text(size = 15, color = "black"),
          axis.title.x = element_text(margin = margin(10, 0, 0, 0), hjust = 1),
          axis.title.y = element_text(margin = margin(0, 10, 0, 0), angle = 90, hjust = 1),
          axis.text = element_text(size = 12, color = "black"),
          strip.text = element_text(size = 15, color = "black", margin = margin(5, 5, 5, 5)),
          plot.title = element_text(size = 15, color = "black", margin = margin(0, 0, 18, 0)))
}

シミュレーション1

データ

PDPとICEの関係を浮き彫りにするため、今回はシミュレーションデータを用意します。

特徴量として X_1, X_2, X_3の3つが観測されていいて、
 X_2 X_1より強く Yと関係しているが、 X_3は関係してないという状況を考えます。
できるだけ簡単なほうが結果がわかりやすいと思うので、以下の線形で加法な形に特定します。*1

\begin{align}
Y &= X_1 - 5X_2 + \epsilon, \\
X_1, X_2, X_3 &\sim U(-1, 1),\\
\epsilon &\sim \mathcal{N}(0, 0.1^2)
\end{align}

ここで、 X_1, X_2, X_3はそれぞれ独立に区間 [0, 1]に一様分布し、 Yには平均0で標準偏差0.1の正規分布に従うノイズが乗るとしています。

シミュレーションのために乱数を発生させます。
Rには確率分布を扱うための関数がデフォルトで用意されていますが、distributions3を用いると、より直感的に確率分布を扱うことができます。

N = 1000 #サンプルサイズ

U = Uniform(-1, 1) # 分布を指定
Z = Normal(mu = 0, sigma = 0.1)

X1 = random(U, N) # 分布から乱数を生成
X2 = random(U, N)
X3 = random(U, N)
E = random(Z, N)

Y = X1 - 5*X2 + E 

df = tibble(Y, X1, X2, X3) # データフレームに

データを可視化して X Yの関係を確認しておきましょう。

df %>% 
  sample_n(200) %>% 
  pivot_longer(cols = contains("X"), values_to = "X") %>% 
  ggplot(aes(X, Y)) +
  geom_point(size = 3, 
             shape = 21,
             color = "white", fill = palette_OkabeIto[5], 
             alpha = 0.7) +
  facet_wrap(~name) + 
  theme_scatter()

f:id:dropout009:20191116212342p:plain

 Y X_1とは弱い正の関係が、 X_2とは強い負の関係があること、 X_3とは無相関であることが見て取れます。

モデル

データの準備ができたので、モデルを作成します。
モデルはtidymodelsparsnipを使って作成します。

ブラックボックスモデルとして、Random Forestを使うことにします。

fitted = rand_forest(mode = "regression", 
                     trees = 1000,
                     mtry = 3,
                     min_n = 1) %>% 
  set_engine(engine = "ranger", 
             num.threads = parallel::detectCores(), 
             seed = 42) %>% 
  fit(Y ~ ., data = df)
  • rand_forest()でモデルにRandom Forestを利用することと、そのハイパーパラメータを指定
  • set_engin()で利用パッケージとパッケージ固有のパラメータを指定
  • fit()でデータを指定して学習

までを一気に行っています。

DALEXによる解釈

モデルの学習が終わったので、DALEXを用いてモデルの振る舞いを解釈していきましょう。
DALEXによる解釈は、DALEX::explain()でexplainerオブジェクトを作るところから始まります。
嬉しいことに、parsnipで学習したモデルはexplain()にそのまま与えることができます。

> explainer = explain(fitted,# 学習済みモデル
+                     data = df %>% select(-Y), # インプット
+                     y = df %>% pull(Y),# ターゲット
+                     label = "Random Forest")# ラベルをつけておくことができる(なくてもいい)

Preparation of a new explainer is initiated
  -> model label       :  Random Forest 
  -> data              :  1000  rows  3  cols 
  -> data              :  tibbble converted into a data.frame 
  -> target variable   :  1000  values 
  -> predict function  :  yhat.model_fit  will be used (default)
  -> predicted values  :  numerical, min =  -5.886102 , mean =  0.008677733 , max =  5.799964  
  -> residual function :  difference between y and yhat (default)
  -> residuals         :  numerical, min =  -0.1736233 , mean =  -8.650032e-05 , max =  0.1944877  
A new explainer has been created!

あとはexplainerにingredientsで用意された様々な解釈手法を適応するだけです。

変数重要度

まずは変数重要度を見てみましょう。これはexplainerをingredients::feature_importance()に与えるだけで計算できます。

> fi = feature_importance(explainer, 
+                         loss_function = loss_root_mean_square, # 精度の評価関数
+                         type = "raw") # "ratio"にするとフルモデルと比べて何倍悪化するかが出る
> fi
      variable dropout_loss         label
1 _full_model_   0.05063265 Random Forest
2           X3   0.06979976 Random Forest
3           X1   0.80947038 Random Forest
4           X2   4.19405191 Random Forest
5   _baseline_   4.26605636 Random Forest

DALEXは学習アルゴリズムに依存しない手法がベースになっているので、変数重要度の計算もpermutationベースのものが使われます。
つまり、その変数をシャッフルしてどのくらい予測精度が落ちるのかがその変数の重要度として定義とされています。
なお、ingredientsによって作成されたオブジェクトはplot()で簡単に可視化することができます。ggplot2がベースになっているので、ggplot2の設定やレイヤーを+で重ねていくこともできます。

plot(fi)

f:id:dropout009:20191116221221p:plain

変数重要度を見ると、 X_2の重要度が一番高く、 X_1は想定的に重要度が低く、 X_3は全く重要ではないとしていることがわかります。これはそもそものシミュレーションの設定に沿っています。

PDP

次はPDPを見てみましょう。PDPはモデルの平均的な予測結果を可視化したものです。
これもexplainerをingredients::partial_dependency()に与えるだけで計算できます。

pdp = partial_dependency(explainer)

pdp %>% 
  plot() + # ggplot2のレイヤーや設定を重ねていくことができる
  scale_y_continuous(breaks = seq(-6, 6, 2)) + 
  theme_scatter() +
  theme(legend.position = "none")

f:id:dropout009:20191116221810p:plain


モデルが X_1, X_2, X_3 Yの関係をうまく捉えられていることが見て取れます。

シミュレーション2

データの作成

PDPは交互作用がない場合はうまくモデルを解釈する事ができますが、交互作用がある場合は可視化がうまく機能しないことが知られています。
以下のようなシミュレーションを考えてみましょう。

\begin{align}
Y &= X_1 - 5X_2 + 10X_2X_3 + \epsilon, \\
X_1, X_2 &\sim U(-1, 1),\\
X_1, X_2 &\sim Bernoulli(0.5),\\
\epsilon &\sim \mathcal{N}(0, 0.1^2)
\end{align}

今回、 X_3 = 0は0か1をとる変数で、それ単体では効果はありませんが、 X_2との交互作用があり、 X_2 X_3 = 0のときは正の、 X_3 = 1のときは負の効果があるという特定化になっています。

N = 1000

U = Uniform(-1, 1)
B = Bernoulli(p = 0.5)
Z = Normal(mu = 0, sigma = 0.1)

X1 = random(U, N)
X2 = random(U, N)
X3 = random(B, N)
E = random(Z, N)

Y = X1 - 5*X2 + 10*X2*X3 + E

df = tibble(Y, X1, X2, X3)


df %>% 
  sample_n(200) %>% 
  pivot_longer(cols = contains("X"), values_to = "X") %>% 
  ggplot(aes(X, Y)) +
  geom_point(size = 3, 
             shape = 21,
             color = "white", fill = palette_OkabeIto[5], 
             alpha = 0.7) +
  facet_wrap(~name) + 
  theme_scatter()

f:id:dropout009:20191116223552p:plain

DALEXによる解釈

PDP

シミュレーション1と同じくモデルを作成、学習し、explainerオブジェクトを作ります。

fitted = rand_forest(mode = "regression", 
                     trees = 1000,
                     mtry = 3,
                     min_n = 1) %>% 
  set_engine(engine = "ranger", 
             num.threads = parallel::detectCores(), 
             seed = 42) %>% 
  fit(Y ~ ., data = df)


explainer = DALEX::explain(fitted,
                           data = df %>% select(-Y),
                           y = df %>% pull(Y),
                           label = "Random Forest")

PDPも同様に作成します。
このシミュレーションでは Y X_2の関係をうまく可視化できるかにフォーカスすることにします。

pdp = partial_dependency(explainer, variables = "X2") #変数を限定

pdp %>% 
  plot() + 
  geom_abline(slope = -5, color = "gray70", size = 1) + 
  geom_abline(slope = 5, color = "gray70", size = 1) + 
  scale_y_continuous(breaks = seq(-6, 6, 2),
                     limits = c(-6, 6)) + 
  theme_scatter() +
  theme(legend.position = "none")

f:id:dropout009:20191116230316p:plain
青い線がPDP、グレーの線が本来の関係です。
PDPを見ると、 X_2が動いても予測値にはまるで効果がないように見えます。
念の為ですが、モデルの予測自体はうまくいっています。
Random Forestは交互作用をうまく学習できますし、実際、OOBでの R^2は0.96です。

> fitted
parsnip model object

Fit in:  275msRanger result

Call:
 ranger::ranger(formula = formula, data = data, mtry = ~3, num.trees = ~1000,      min.node.size = ~1, num.threads = ~parallel::detectCores(),      seed = ~42, verbose = FALSE) 

Type:                             Regression 
Number of trees:                  1000 
Sample size:                      1000 
Number of independent variables:  3 
Mtry:                             3 
Target node size:                 1 
Variable importance mode:         none 
Splitrule:                        variance 
OOB prediction error (MSE):       0.3297557 
R squared (OOB):                  0.9624996 

ここで起きている現象は、「モデルは交互作用を学習できているが、PDPはそれをうまく可視化できていない」というものです。PDPは交互作用を平均化してしまうため、プラスの交互作用とマイナスの交互作用が相殺して効果がないように見えてしまっています。

ICE Plot

交互作用をうまく捉える手法の一つがICEになります。
ICEは単純にPDPの平均化する前の出力を全てプロットするというものです。
平均化をしていないので交互作用が相殺されることを防ぐことができます。
例によって、ingredients::ceteris_paribus()にexplainerを与えるだけで計算できます。

# 線が多すぎるとわけがわからなくなるので100サンプルだけ抜き出す
# tibbleのまま渡すと警告が出るのでdata.frameにしている。
# 警告の内容を見るとうまく動かなさそうだが、いまのところtibbleのままでもうまく動いているように思う
ice = ceteris_paribus(explainer, 
                      variables = "X2",
                      new_observation = df %>% sample_n(100) %>% as.data.frame()) 

ice %>% 
  plot(alpha = 0.5, size = 0.5, color = colors_discrete_drwhy(1)) + 
  geom_abline(slope = -5, color = "gray70", size = 1) + 
  geom_abline(slope = 5, color = "gray70", size = 1) +
  scale_y_continuous(breaks = seq(-6, 6, 2),
                     limits = c(-6, 6)) + 
  theme_scatter() +
  theme(legend.position = "none")

f:id:dropout009:20191116230534p:plain
青い線の一本一本が各インスタンスの予測結果に対応しています。
PDPとは違って、ICEではモデルが交互作用を学習していることをうまく可視化することができています。
なお、各線のタテ方向のブレはインスタンスごとの X_1の大きさによるものです。仮に X_1の値が全てのインスタンスで共通なら一直線上に並んだきれいなX印になります。

Conditional PDP

PDPは交互作用のある変数についても平均化してしまうのが問題でした。
ICEはひとつの解決策でしたが、これはこれでやたら線が多くなってしまいますし、値の安定性も低くなります。
今回は X_3=0のときと X_3=1のときでX_2の効果が逆になることが問題でした。
ということは、 X_3の値でグループに分けてPDPを計算することで、交互作用の問題を解決できそうです。
これを行うための関数がingredients::aggregate_profiles()として用意されています。

# ICEを与える。グループを指定するとPDPをグループごとに求めることができる。
conditional_pdp = aggregate_profiles(ice, groups = "X3") 

conditional_pdp %>% 
  plot() + 
  geom_abline(slope = -5, color = "gray70", size = 1) + 
  geom_abline(slope = 5, color = "gray70", size = 1) +
  scale_y_continuous(breaks = seq(-6, 6, 2),
                     limits = c(-6, 6)) + 
  theme_scatter() +
  theme(legend.position = "none")

f:id:dropout009:20191116232138p:plain

単純なPDPでは捉えられなかった交互作用を可視化することに成功しました。
現実のシチュエーションでは、たとえば男女でグループ分けすることで、モデルが男女の効果の違いを捉えているのかを確認することができます。

clusterd ICE Plot

今回はシミュレーションなので X_3でグループ分けすればいいことがわかっていましたが、実際にはどの変数でグループ化すればいいかわからない場合も多いかと思います。
ingredients::cluster_profiles()を用いることで、似たようなICEをクラスター化してPDPを計算してくれます。

# ICEを与える。クラスター数を指定する。
clustered_ice = cluster_profiles(ice, k = 2)

clustered_ice %>% 
  plot() + 
  geom_abline(slope = -5, color = "gray70", size = 1) + 
  geom_abline(slope = 5, color = "gray70", size = 1) +
  scale_y_continuous(breaks = seq(-6, 6, 2),
                     limits = c(-6, 6)) + 
  theme_scatter() +
  theme(legend.position = "none")

f:id:dropout009:20191116232848p:plain

自分で X_3でグループ化した場合と同じく、単純なPDPでは捉えられなかった交互作用を可視化することに成功しました。

まとめ

本記事では、tidymodelsを用いて機械学習モデルを作成し、それをDALEXingredientsを用いて解釈する方法をまとめました。
もう一つの重要なパッケージであるiBreakDownは別の記事でまとめたいと思っています。

本記事で使用したコードは以下にまとめてあります。

github.com

*1:実際にこのデータに遭遇したら線形回帰を使ったほうがいいというのはありますが、わかりやすさのためこうしました

tidymodelsによるtidyな機械学習(その3:ハイパーパラメータのチューニング)

はじめに

前回の記事ではハイパーパラメータのチューニングをfor loopを用いたgrid searchでやっっていました。 tidymodels配下のdialstuneを用いることで、より簡単にハイパーパラメータのサーチを行えるので、本記事ではその使い方を紹介したいと思います。 なお、パラメータサーチ以外のtidymodelsの使い方には本記事では言及しないので、以下の記事を参考にして頂ければと思います。

dropout009.hatenablog.com

dropout009.hatenablog.com

前処理

まずは前回の記事と同様、rsampleで訓練/テストデータの分割を行います。なお、例によってデータはdiamondsを用います。

# パッケージ
library(tidyverse)
library(tidymodels)
set.seed(42)

# Train/Testの分割
df_split = initial_split(diamonds,  p = 0.8)

df_train = training(df_split)
df_test  = testing(df_split)

ハイパーパラメータのサーチ

最終的にはtune::tune_grid()でハイパーパラメータを探索しますが、 そのためにTrain/Validationに分割されたデータ、前処理レシピ、学習用モデル、ハイパーパラメータの候補の4つが必要になります。 最初の3つは以前の記事で触れているので、4つ目のハイパーパラメータに関する部分に詳しく触れていきます。

Train/Validationデータ

これはrsamplesの仕事です。普通はCross Validationで評価するので、それに合わせます。

>df_cv = vfold_cv(df_train, v = 5)

> df_cv
#  5-fold cross-validation 
# A tibble: 5 x 2
  splits               id   
  <named list>         <chr>
1 <split [34.5K/8.6K]> Fold1
2 <split [34.5K/8.6K]> Fold2
3 <split [34.5K/8.6K]> Fold3
4 <split [34.5K/8.6K]> Fold4
5 <split [34.5K/8.6K]> Fold5

前処理レシピ

次に、前処理のレシピを作成します。モデルはRandom Forestを使うので、前回の分析同様、前処理は最低限にしておきます。

> rec = recipe(price ~ ., data = df_train) %>% 
+   step_log(price) %>% 
+   step_ordinalscore(all_nominal())
> 
> rec
Data Recipe

Inputs:

      role #variables
   outcome          1
 predictor          9

Operations:

Log transformation on price
Scoring for all_nominal

学習用モデル

parsnipでモデルを設定します。今回もRandom Forestを使いましょう。

> model = rand_forest(mode = "regression",
+                     trees = 50, # 速度重視
+                     min_n = tune(),
+                     mtry = tune()) %>%
+   set_engine("ranger", num.threads = parallel::detectCores(), seed = 42)
> 
> model
Random Forest Model Specification (regression)

Main Arguments:
  mtry = tune()
  trees = 50
  min_n = tune()

Engine-Specific Arguments:
  num.threads = parallel::detectCores()
  seed = 42

Computational engine: ranger 

サーチしたいパラメータはここでは値を決めずtune::tune()を与えます。今回はmin_nmtryをサーチすることにします。

ハイパーパラメータ

ここが本記事での新しい内容になります。 dialesにはparsnipで指定できるハイパーパラメータに関して、探索するレンジを指定するための関数が用意されています。 たとえばmin_n()はRandom Forestのハイパーパラメータmin_nに対応する関数で、デフォルトだと2-40のレンジでハイパーパラメータが探索されます。 なお、これは最終ノードに最低でも必要なインスタンスの数を表していて、これを大きくするとより強い正則化がかかります。

> min_n()
Minimal Node Size  (quantitative)
Range: [2, 40]

自分でレンジを決めたい場合は引数で指定することができます。

> min_n(range = c(1, 10))
Minimal Node Size  (quantitative)
Range: [1, 10]

同様に、もうひとつのハイパーパラメータmtryに関しても見てみましょう。 こちらは各ツリーでの分割の際に用いる特徴量の数で、これを小さくするとより強い正則化がかかります。

> mtry()
# Randomly Selected Predictors  (quantitative)
Range: [1, ?]

min_n()とは違って、レンジの最大値が指定されていません。モデルに投入する特徴量の数よりも大きいmtryを探索しても意味がない*1ので、こちらで直接指定してあげる必要があります。 この際、mtry(range = c(1, 5))のようにレンジを指定することもできますが、実際にモデルに投入するデータフレームを与えてあげることでレンジを指定することもできます。

> # 前処理済み学習用データ
> df_input = rec %>% 
+   prep() %>% 
+   juice() %>% 
+   select(-price)
> 
> finalize(mtry(), df_input)
# Randomly Selected Predictors  (quantitative)
Range: [1, 9]

dials::finalize()の第1引数にレンジを指定したい関数を、第2引数に特徴量のデータフレームを指定することで、適切なレンジが指定されます。 特徴量の数は分析の途中で変動するので、このやり方は柔軟性があって良いんじゃないかと思います。

さて、これで探索範囲の指定ができるようになりました。 次に、探索したいハイパーパラメータのリストを作ってtune::parameters()に与えるとparametersオブジェクトを作ることができます。

> params = list(min_n(),
+               mtry() %>% finalize(rec %>% prep() %>% juice() %>% select(-price))) %>% 
+   parameters()
> 
> params
Collection of 2 parameters for tuning

    id parameter type object class
 min_n          min_n    nparam[+]
  mtry           mtry    nparam[+]

このparametersオブジェクトをdials::grid_*()に渡すことで、実際に探索するハイパーパラメータの値を作ることができます。 たとえばgrid_regular()ならグリッドサーチ、grid_random()ならランダムサーチになります。 今回はランダムサーチにしましょう。サーチの数はsizeで指定できます。

> df_grid = params %>% 
+   grid_random(size = 10) # 実際はもっと多いほうがいい
> 
> df_grid
# A tibble: 10 x 2
   min_n  mtry
   <int> <int>
 1     7     8
 2    19     8
 3    22     3
 4    38     1
 5    12     8
 6    36     1
 7    21     6
 8     7     5
 9    29     2
10    38     3

チューニング

これでハイパーパラメータの探索準備が整いました! これまでに作ったオブジェクトをtune::tune_grid()に渡します。

df_tuned = tune_grid(object = rec,
                     model = model, 
                     resamples = df_cv,
                     grid = df_grid,
                     metrics = metric_set(rmse, mae, rsq),
                     control = control_grid(verbose = T))
  • objectrecipeで作成された前処理レシピを渡します。resamplesに与えたデータがこの前処理を済ませたあとでモデルに投入されます。
  • modelparsnipで定義した学習用のモデルです。
  • resamplesrsamplesで作った学習/評価用のデータを渡します。普通はCross Varidationで評価すると思うのでrsamples::vfold_cv()で作ったデータフレームを渡すのがいいと思います。
  • griddialstuneで作ったハイパーパラメータの候補が格納されたデータフレームを渡します。
  • metrics:精度の評価指標です。yardsticに準備されている関数を指定することができます。
  • control:指定しなくても構いませんが、今回はログが出力されるようにしています。

学習/評価が終わると、以下のようなデータフレームが手に入ります。

> df_tuned
#  5-fold cross-validation 
# A tibble: 5 x 4
  splits               id    .metrics          .notes          
* <list>               <chr> <list>            <list>          
1 <split [34.5K/8.6K]> Fold1 <tibble [20 × 5]> <tibble [0 × 1]>
2 <split [34.5K/8.6K]> Fold2 <tibble [20 × 5]> <tibble [0 × 1]>
3 <split [34.5K/8.6K]> Fold3 <tibble [20 × 5]> <tibble [0 × 1]>
4 <split [34.5K/8.6K]> Fold4 <tibble [20 × 5]> <tibble [0 × 1]>
5 <split [34.5K/8.6K]> Fold5 <tibble [20 × 5]> <tibble [0 × 1]>

.metricsに予測精度が格納されています。unnest()でもとってこれますが、手っ取り早い関数としてtune::collect_metrics()が準備されています。

> df_tuned %>% 
+   collect_metrics()
# A tibble: 30 x 7
    mtry min_n .metric .estimator   mean     n  std_err
   <int> <int> <chr>   <chr>       <dbl> <int>    <dbl>
 1     1    38 mae     standard   0.120      5 0.00131 
 2     1    38 rmse    standard   0.159      5 0.00173 
 3     1    38 rsq     standard   0.979      5 0.000246
 4     2    29 mae     standard   0.0792     5 0.000317
 5     2    29 rmse    standard   0.107      5 0.000518
 6     2    29 rsq     standard   0.989      5 0.000160
 7     4    19 mae     standard   0.0677     5 0.000174
 8     4    19 rmse    standard   0.0932     5 0.000441
 9     4    19 rsq     standard   0.992      5 0.000109
10     5     3 mae     standard   0.0651     5 0.000176
# … with 20 more rows

特に精度の高いハイパーパラメータの候補が知りたい場合は、tune::show_best()で確認することができます。

> df_tuned %>% 
+   show_best(metric = "rmse", n_top = 3, maximize = FALSE)
# A tibble: 3 x 7
   mtry min_n .metric .estimator   mean     n  std_err
  <int> <int> <chr>   <chr>       <dbl> <int>    <dbl>
1     7     8 rmse    standard   0.0910     5 0.000702
2     5     3 rmse    standard   0.0912     5 0.000558
3     7    16 rmse    standard   0.0915     5 0.000679

一番精度の高かったハイパーパラメータを使って全訓練データでモデルを再学習する場合は、tune::select_best()でハイパーパラメータをとってきてupdate()でモデルをアップデートできます。

# 一番精度の良かったハイパーパラメータ
> df_best_param = df_tuned %>% 
+   select_best(metric = "rmse", maximize = FALSE)
> df_best_param
# A tibble: 1 x 2
   mtry min_n
  <int> <int>
1     7     8

# モデルのハイパーパラメータを更新
> model_best = update(model, df_best_param)
> model_best
Random Forest Model Specification (regression)

Main Arguments:
  mtry = 7
  trees = 50
  min_n = 8

Engine-Specific Arguments:
  num.threads = parallel::detectCores()
  seed = 42

Computational engine: ranger 

あとは通常通り学習・評価を行うだけです。

まとめ

本記事ではtidymodels配下のdialstuneを用いたハイパーパラメータのサーチを紹介しました。 前回の記事ではCross Validationデータにどうアクセスするか結構ややこしかったと思うのですが、tuneを用いることでよりスッキリと探索と評価を行うことができます。 dialstuneは公式ドキュメントに個別の使い方はまとまっているのですが、目的がわかりにくい部分や、どの情報が最新かわからないところもあり、典型的なタスクに対しての用途を自分でまとめ直したという経緯になります。

ちなみに、今回はランダムサーチを用いましたが、tuneではbaysian optimizationを用いたパラメータ探索も可能となっています。ぜひ確認してみて下さい*2

本記事で使用したコードは以下にまとめてあります。

github.com

参考文献

*1:実際は特徴量の数よりも大きい値を入れるとエラーを吐きます

*2:Classification Example • tune

Synthetic Difference In Differences(Arkhangelsky et. al., 2019)を読んだ


はじめに

GW中にSynthetic difference in differences(Arkhangelsky, D., Athey, S., Hirshberg, D. A., Imbens, G. W., & Wager, S. (2019)) を読みました。
ランダム化比較試験(RCT)が行えない状況でのパネルデータのから処置効果を推定する際には、バイアスを取り除くためにDifference In DifferencesやSynthetic Controlがよく用いられていますが、
Synthetic Difference In Differencesは名前の通りDifference In DifferencesとSynthetic Controlのいいとこ取りになっています。
パネルデータから因果推論を行う際に非常に強力な武器になると思ったので、本記事でコンセプトをまとめ、論文内の比較実験を再現しました。


コンセプト

セッティング

 i = 1, \dots, Nのユニット(たとえば個人)に対して t = 1, \dots, T時点でのデータがとれているパネルデータを考えます。このとき、アウトカム Y_{it}は以下の N \times T行列で表現できます。*1

\begin{align}
\mathbf{Y}
= \begin{pmatrix}
Y_{11} & \cdots & Y_{1T}\\
\vdots & \ddots & \vdots\\
Y_{N1} & \cdots & Y_{NT}
\end{pmatrix}
\end{align}


簡単のために、一番シンプルな設定として、 i = Nのユニットが  t = T 時点のみ処置を受けるケースを考えます。つまり、 i = 1, \dots, N-1 はコントロール群(c)で i=Nは処置群(t)になりますし、 t = 1, \dots T-1が処置前(pre)、 t = Tが処置後(post)となります。

f:id:dropout009:20190506184615p:plain

ここで、仮想的に i=Nのユニットが t=T時点で処置を受けた場合のアウトカムを Y_{NT}(1)、受けなかった場合のアウトカムを Y_{NT}(0)と書くことにしましょう。このとき、 Y_{NT}(1) - Y_{NT}(0)つまり処置を受けた場合と受けなかった場合の差をみることで、処置がアウトカムに与えるインパクトを測定することができます。ただ、実際に僕たちに観測できるのは処置を受けた Y_{NT}(1)のみで、処置を受けなかったときの Y_{NT}(0)は観測することができません。そこで、なんらかの手法を使って、実際には観測できない Y_{NT}(0)を推定する必要が出てきます。これを反事実と呼びます。以降、反事実 Y_{NT}(0)を推定する手法としてDifference In DIfference(DID)、Synthetic Control(SC)、そしてその融合であるSynthetic Difference In Difference(SDID)を比較していきます。

Difference in Differences (DID)

処置前コントロール群の値にに加えて、(1)コントロール群と処置群のアウトカムにそもそもどのくらい違いがあるのか、(2)処置を受けなくても処置の前後でどのくらいアウトカムが変わるのかを考慮することで、処置後処置群の値 Y_{NT}(0)が推定できるのではというのがDIDの発想です。具体的には、DIDによる推定量 \hat{Y}_{NT}^{\text{DID}}(0) は以下で定義されます。

\begin{align}
\hat{Y}_{NT}^{\text{DID}}(0) &= \bar{Y}^{c, pre} + \left(\bar{Y}^{c, post} - \bar{Y}^{c, pre}\right)+ \left(\bar{Y}^{t, pre} - \bar{Y}^{c, pre}\right)\\
&=\frac{1}{N- 1}\frac{1}{T-1}\sum_{i=1}^{N-1}\sum_{t=1}^{T-1}Y_{it} \\
&\quad+ \left(\frac{1}{N- 1}\sum_{i=1}^{N-1}Y_{iT} - \frac{1}{N- 1}\frac{1}{T-1}\sum_{i=1}^{N-1}\sum_{t=1}^{T-1}Y_{it}\right)\\
&\quad+ \left(\frac{1}{T-1}\sum_{t=1}^{T-1}Y_{Nt} - \frac{1}{N- 1}\frac{1}{T-1}\sum_{i=1}^{N-1}\sum_{t=1}^{T-1}Y_{it}\right)
\end{align}

ここで、 \bar{Y}はそれぞれのグループの平均値を表しています。上式は以下のように解釈できます。まず、単に処置前コントロール群のみを用いて予測を行うと(第一項)、コントロール群と処置群の違いでバイアスがかかり、次に処置前と処置後の違いでもバイアスがかかります。そこで、コントロール群<->処置群の差を第二項で、処置前<->処置後の差を第三項で補正しています。この意味において、DID推定量はバイアスを二重に補正していると言えます。

f:id:dropout009:20190506185121p:plain

最後に、こうしてDIDで推定した \hat{Y}_{NT}^{\text{DID}}(0)と実際の観測値 Y_{NT}(1)の差を見ることで処置のインパクトが推定できます。

\begin{align}
\hat{\tau}^{\text{DID}} = Y_{NT}(1) -\hat{Y}_{NT}^{\text{DID}}(0)
\end{align}

Synthetic Control (SC)

DIDは単純平均でしたが、コントロール群の加重平均で処置群を表現しようというのがSCの発想になります。SCは以下の2ステップを踏みます。まず、加重平均に用いる各コントロール i = 1, \dots, N-1の重み \hat{\omega}_iを、うまく処置群を近似できるように決めます。

\begin{align}
\sum_{i=1}^{N-1} \hat{\omega}_{i} Y_{it} \approx Y_{N t} \text { for all } t=1, \ldots, T-1
\end{align}

f:id:dropout009:20190506185108p:plain

具体的には、以下の二乗誤差を小さくするように \hat{\omega}_iを求めます。

\begin{align}
\hat{\omega}&=\underset{\omega \in \mathbb{W}}{\arg \min }\;\frac{1}{T-1} \sum_{t=1}^{T-1}\left(\sum_{i=1}^{N-1} \omega_{i} Y_{i t}-Y_{N t}\right)^{2} + \frac{1}{2}\zeta\|\omega\|_2\\
\text{where}\quad \mathbb{W}&=\left\{\omega \in \mathbb{R}^{N - 1} \;\bigg|\; \omega_{i} \geq 0, \; \sum_{i=1}^{N-1} \omega_{i}=1\right\}
\end{align}

 \omegaは加重平均の重みなので0以上、足して1の制約がかかっています。また、 L_2正則化を入れることで \omegaの推定を安定させています。

次に、この重みを使って処置後である Y_{NT}(0)を推定します。
\begin{align}
\hat{Y}_{N T}^{\mathrm{SC}}(0)=\frac{1}{T-1} \sum_{i=1}^{N-1} \sum_{t=1}^{T-1} \hat{\omega}_{i} Y_{i t}+\left(\sum_{i=1}^{N-1} \hat{\omega}_{i}Y_{i T}-\frac{1}{T-1} \sum_{i=1}^{N-1} \sum_{t=1}^{T-1} \hat{\omega}_{i}Y_{i t}\right)
\end{align}

DIDの推定量 \hat{Y}_{NT}^{\text{DID}}(0)と見比べると:

  • DIDでは単純平均でしたが、SCは加重平均を用いることで処置群の近似を改善しています。
  • その一方で、SCにはDIDにあった処置前後のバイアス補正(第三項)が存在しません。

Synthetic Difference In Differences (SDID)

DIDとSCを見比べることで、双方の利点と欠点が見えました。そこで、SDIDでは両方のいいとこ取りをします。つまり、

  • 単純平均ではなく加重平均を用いる
  • コントロール/トリートメントのバイアス補正だけでなく、処置前後のバイアス補正を入れる

ことで反事実 Y_{N T}(0)のよりよい推定を目指します。

\begin{align}
\hat{Y}_{N T}^{\mathrm{SDID}}(0)= \sum_{i=1}^{N-1} \sum_{t=1}^{T-1} \hat{\omega}_{i} \hat{\lambda}_{t} Y_{i t}+\left(\sum_{t=1}^{T-1} \hat{\lambda}_{t}Y_{N t}-\sum_{t=1}^{T-1} \sum_{i=1}^{N-1} \hat{\omega}_{i}\hat{\lambda}_{t} Y_{i t}\right)+\left(\sum_{i=1}^{N-1} \hat{\omega}_{i}Y_{i T}-\sum_{i=1}^{N-1} \sum_{t=1}^{T-1} \hat{\omega}_{i}\hat{\lambda}_{t} Y_{i t}\right)
\end{align}


ここで、 \hat{\lambda}_tは時間方向の重みであり、 \hat{\omega}_iと同様に以下の最小化問題を解くことで求めます。

\begin{align}
\hat{\lambda}&=\underset{\lambda_0\in\mathbb{R},\; \lambda \in \mathbb{L}}{\arg \min } \;\frac{1}{N-1}\sum_{i=1}^{N-1}\left(\lambda_0 + \sum_{t=1}^{T-1} \lambda_{t} Y_{i t}-Y_{i T}\right)^{2}+ \frac{1}{2}\xi\|\lambda\|_2\\
\text{where}\quad \mathbb{L}&=\left\{\lambda \in \mathbb{R}^{T - 1} \;\bigg|\; \lambda_{t} \geq 0, \; \sum_{t=1}^{T-1} \lambda_{t}=1\right\}
\end{align}


ただし、トレンドに対応するために \lambda_0が入っています。これは重みとしては用いません。


比較実験

論文ではAbadie, Diamond, and Hainmueller(2010)のデータを用いてDID, SC, SDIDを比較した実験を行っているので、それを再現した結果を紹介します。

ADH(2010)ではカリフォルニアの禁煙法が喫煙に与えたインパクトをSynthetic Controlを用いて推定しています。
使用データはカリフォルニアを含むアメリカ39州の年度別喫煙量のデータで、1970年から2000年にかけて31年分あります。なお、禁煙法は1989年から実施されています。
今回の実験では実際に禁煙法の効果を調べたいのではなく、 Y_{NT}(0)をうまく予測できるかを確かめたいので、処置前(1988年以前)のデータに関して、各州の1期先の喫煙量を他の38州の喫煙量からDID/SC/SDIDで予測し、精度を検証します。*2
よりよい精度で1期先の予測ができるなら、処置の効果 \hat{\tau}_{NT} = Y_{NT}(1) - \hat{Y}_{NT}(0)をより正確に推定できると言えます。

具体的には、各州 i = 1, \dots, 39に対して1980年から1988年までの9年分の予測をDID/SC/SDIDで行い、RMSEを計算します。
\begin{align}
\operatorname{RMSE}_{i}=\sqrt{\frac{1}{9} \sum_{t=1980}^{1988}\left(Y_{it}-\hat{Y}_{it}\right)^{2}}
\end{align}


その結果が以下になります。

f:id:dropout009:20190506172341p:plain


DID/SCと比較して、SDIDはより高い精度で Y_{NT}(0)を予測することに成功しています。中央値で見るとDIDと比較して50%、SCと比較しても15%の改善です。


具体的にどの部分で精度に差が出ているのかプロットを見て確認してみます。

f:id:dropout009:20190506172657p:plain

上図は各州の喫煙量の実測値とDID/SC/SDIDの推定値がプロットされています。横軸が年度で縦軸が喫煙量です。
図が細かいのですが、基本的にSDID(赤)はDID(オレンジ)/SC(緑)よりも実測値Y(青)をうまく予測できています。
特に差が顕著な部分をいくつか抜き出したものが下の図になります。

f:id:dropout009:20190506173610p:plain

DIDは全体的に精度が悪いと思いますが、
Synthetic Controlは非常に特徴的で、New HampshireやUtahなど、他の州の加重平均で表現できない州は予測がうまくいっていません。
これに対して、SDIDは非常にロバストな予測ができています。


まとめ

本記事では、Difference In DifferencesとSynthetic Controlを融合させたSynthetic Difference In Differencesのコンセプトをまとめ、実際に精度が改善することを確認しました。
SDIDはRCTが行えない状況でパネルデータの因果推論を行う際、非常に強力な武器になると思います。


参考文献

Arkhangelsky, D., Athey, S., Hirshberg, D. A., Imbens, G. W., & Wager, S. (2019). Synthetic difference in differences (No. w25532). National Bureau of Economic Research.[1812.09970] Synthetic Difference in Differences

GitHub - swager/sdid-script: Example scripts for synthetic difference in differences

Alberto Abadie, Alexis Diamond, and Jens Hainmueller. Synthetic control methods for comparative case studies: Estimating the effect of California’s tobacco control program. Journal of the American Statistical Association, 105(490):493–505, 2010. https://economics.mit.edu/files/11859

再現コード
GitHub - dropout009/sdid_python: replication python code for synthetic difference in difference simulation. please check https://arxiv.org/abs/1812.09970 and https://github.com/swager/sdid-script

*1:簡単のため、とりあえず共変量 X_{it}は無視することにします。

*2:論文ではやってませんが、カリフォルニアを除いてしまえば2000年まで実験することもできると思います