変数を 1 つ足しただけで傾きの符号が反対になる重回帰分析

$y$ を $x_1$ だけで説明すると、傾きは $+3.385$ です。ここへ $x_2$ を加えて同時に当てはめると、$x_1$ の係数は $-2.209$ になります。

同じデータ、同じ最小二乗法。それでも符号が反対になる。

単回帰の傾きと重回帰の係数は、測っているものが違うためです。

モデル

説明変数が $p$ 個ある場合、次の形で書きます^[2]。

$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdots +{\beta }_p{x}_p+\epsilon$$

最小にするのは、単回帰と同じく縦方向のずれの 2 乗和です^[1]。変数が増えても、当てはめる面が直線から平面や超平面に変わるだけ。

行列で書く

観測が $n$ 個、説明変数が $p$ 個なら、次の行列を作ります。

$$X=\left(\begin{array}{cccc}1& x_{11}& \cdots & x_{1p}\\ 1& x_{21}& \cdots & x_{2p}\\ ⋮& ⋮& & ⋮\\ 1& x_{n1}& \cdots & x_{np}\end{array}\right)$$

左端の 1 の列は切片のためのものです。これでモデルが 1 行に収まります。

$$\mathbf{y}=X\mathbf{\beta }+\mathbf{\epsilon }$$

2 乗和を $\mathbf{\beta }$ で微分して $0$ と置くと、正規方程式 $X^{\top }X\mathbf{\beta }=X^{\top }\mathbf{y}$ が出ます。

$$\hat{\mathbf{\beta }}=(X^{\top }X{)}^{-1}{X}^{\top }\mathbf{y}$$

$X^{\top }X$ が逆行列を持つことが条件です。持たないのは、説明変数のあいだに完全な 1 次関係があるとき。

偏回帰係数の意味

${\hat{\beta }}_1$ は「$x_1$ が 1 増えたときの $y$ の増分」ではありません。「ほかの変数を止めたまま $x_1$ が 1 増えたときの増分」です。

単回帰の傾きは、ほかの変数を止めていません。$x_1$ が動くと $x_2$ も一緒に動き、その効果もまとめて拾ってしまう。

止めるかどうかで値が変わるのは当然で、変わらないほうがめずらしい。

例：符号が反転する

12 個の観測で、$x_1$ と $x_2$ の相関が 0.9926 という強い関係にある場合です。

$x_2$ を入れると $x_1$ の係数が $+3.385$ から $-2.209$ へ動きました。$x_2$ の係数は $+5.279$ です。

$x_1$ が増えるとき $x_2$ もほぼ同じだけ増え、$x_2$ の効果が大きいので、合計では $y$ が増えて見えます。$x_2$ を止めれば、$x_1$ 単独の効果は負だった。

どちらが正しい、という問いは立ちません。答えたい問いが「$x_1$ を増やしたら $y$ はどうなるか」なら、$x_2$ も一緒に動くのかどうかを先に決める必要があります。

追加変数プロット

偏回帰係数を絵で見る方法があります。$y$ と $x_1$ の両方から $x_2$ の効果を先にとりのぞき、残ったものどうしを並べる。

この散布図に単回帰を当てた傾きが、重回帰の ${\hat{\beta }}_1$ にちょうど一致します。

左は右上がり、右は右下がりです。$x_2$ をとりのぞくと、$x_1$ と $y$ の関係は反対を向いていました。

面の切り口として見る

2 変数のモデルは、$x_1$-$x_2$-$y$ の空間に平面を当てているのと同じです。

$x_2$ を 1 つの値に止めると、その平面の切り口が $y$-$x_1$ 平面の直線になります。この直線の傾きが $-2.209$ です。

$x_2$ を上げていくと、切り口の直線が上へ平行移動します。データ点は $x_1$ が大きいところほど $x_2$ も大きいので、高い位置の直線の上に乗る。

各点を結んで見ると右上がりに見えますが、それは直線をまたいで移動しているためです。同じ直線の上を動いているわけではありません。

多重共線性

説明変数どうしが強く相関していると、$X^{\top }X$ が逆行列を持ちにくくなります。完全に一致すれば逆行列が存在しません。

一致しない場合でも、近いだけで係数の推定が不安定になります。どちらの変数に効果を割りふるかを、データが決めきれないためです。

目安になるのが分散拡大係数です。$x_j$ をほかの説明変数で回帰したときの決定係数を $R_j^2$ として、次で定義されます^[3]。

$${\mathrm{VIF}}_j=\frac{1}{1-R_j^2}$$

この例では $x_1$ と $x_2$ の相関が 0.9926 なので、$\mathrm{VIF}=67.6$。標準誤差が $\sqrt{67.6}=8.2$ 倍に膨らんでいることを意味します。

10 を超えたら注意する、という目安がよく使われます。ただし標本の大きさや誤差の分散でも値が動くので、これだけで判断は決まりません^[3]。

決定係数は必ず上がる

まったく無関係な変数を足しても、$R^2$ は下がりません。当てはめの自由度が増えるだけなので、2 乗和は減るか同じかのどちらかです。

先の例に無関係な $x_3$ を足してみます。

$R^2$ は上がりましたが、調整済みのほうは下がりました。$x_3$ の $t$ 値は 0.500、p 値は 0.6303 で、効いている証拠はありません。

調整済み決定係数は、変数の数で罰を与えた形です。

$$R_{\text{adj}}^2=1-(1-R^2)\frac{n-1}{n-p-1}$$

変数を足したときの $R^2$ の上がり方が、失う自由度に見合わなければ値が下がります。モデルを比べるときは、こちらを見る。

係数の検定

各係数について、$t={\hat{\beta }}_j/\mathrm{SE}({\hat{\beta }}_j)$ を自由度 $n-p-1$ の t 分布と比べます。

例では $x_1$ の係数が $-2.209$、標準誤差が 0.293 なので $t=-7.54$。自由度 9 で p 値は $3.6\times 10^{-5}$ です。

95% 信頼区間は $[-2.87,-1.55]$ で、全体が負の側にあります。符号が負であること自体は、偶然では説明しにくい。

各係数の検定は、ほかの変数がモデルに入っている状態での話です。入れる変数を変えれば、同じ $x_1$ でも t 値が変わります。

単回帰の傾きは合計の効果、重回帰の係数はほかを止めたときの効果。どちらの問いに答えたいかを決めてから、モデルを組みます。

参考