はじめに
前回の記事の続きである.
https://zenn.dev/miloc/articles/4c1e36b7c8050e
線形重回帰における回帰係数ベクトル
β^=(X⊤X)−1X⊤y=detX⊤X1X⊤X X⊤y
について,前回はdetX⊤Xの変換性について以下の二つが成り立つことを確認した:
- 説明変数データの定数c倍については,detX⊤Xがc2倍となる.
- 説明変数データに定数を加算しても,detX⊤Xは不変である.
(ここでX=(1x1⋯xp)はp個の説明変数データベクトルx1,…xp∈Rnと全ての成分が1である定数ベクトル1∈Rnからなる計画行列である.この記事でも簡単のためp=2として考察する.)
今回の記事ではX⊤X X⊤yの変換性について考察し,説明変数データの(線形)変換に応じた回帰係数の変換の様子を確認する.加えて,予測値の変換の様子も見ることとする.また,Pythonを使った確認も行う.
最初に結論
計画行列X=(1x1x2)から作成した線形回帰モデル
fX(t1,t2)=β^0+β^1t1+β^2t2
があるとする.
ここで,x1→x1new=c1x1+d11, x2→x2new=c2x2+d21という変換を行ったデータからなる計画行列Xnew=(1x1newx2new)で作られる線形回帰モデルは,上式のt1をt1→c1t1new−d1で,t2をt2→c2t2new−d2でそれぞれ置き換えた
fXnew(t1new,t2new)=β^0+β^1(c1t1new−d1)+β^2(c2t1new−d2)
である.
このうえで,変換後のモデルfXnewに対して"学習データ",つまり計画行列Xnewの第i行にあるxi1new,xi2newを入力すると,
fXnew(xi1new,xi2new)=β^0+β^1(c1xi1new−d1)+β^2(c2xi2new−d2)=β^0+β^1(c1(c1xi1+d1)−d1)+β^2(c2(c2xi2+d2)−d2)=β^0+β^1xi1+β^2xi2=fX(xi1,xi2)
となる.つまり入力が学習データである限りにおいては,予測値は学習データの(線形)変換をしても不変である.よって決定係数R2も不変である.
以下,これらのことを確かめていく.
行列式の計算による確認
余因子行列X⊤Xの成分を∗ij(添え字は0から始まるとした)で表すと,
β^=β^0β^1β^2=detX⊤X1X⊤X X⊤y=detX⊤X1∗00∗10∗20∗01∗11∗21∗02∗12∗221⊤yx1⊤yx2⊤y
という状況である.変換性のわかっているdetX⊤X以外の部分について考察すればよい.今
X⊤X=1⊤x1⊤x2⊤(1x1x2)=1⊤1x1⊤1x2⊤11⊤x1x1⊤x1x2⊤x11⊤x2x1⊤x2x2⊤x2
であるので,余因子行列の各成分は以下の通りとなる(余因子行列の定義は例えば前回の記事を参照):
∗00=+x1⊤x1x2⊤x1x1⊤x2x2⊤x2,∗01=−1⊤x1x2⊤x11⊤x2x2⊤x2,∗02=+1⊤x1x1⊤x11⊤x2x1⊤x2,∗10=−x1⊤1x2⊤1x1⊤x2x2⊤x2,∗11=+1⊤1x2⊤11⊤x2x2⊤x2,∗12=−1⊤1x1⊤11⊤x2x1⊤x2,∗20=+x1⊤1x2⊤1x1⊤x1x2⊤x1,∗21=−1⊤1x2⊤11⊤x1x2⊤x1,∗22=+1⊤1x1⊤11⊤x1x1⊤x1.
よって,回帰係数ベクトルβ^=(β^0, β^1, β^2)⊤のそれぞれの成分は
(detX⊤X)β^0(detX⊤X)β^1(detX⊤X)β^2=∗00×1⊤y+∗01×x1⊤y+∗02×x2⊤y=x1⊤x1x2⊤x1x1⊤x2x2⊤x21⊤y−1⊤x1x2⊤x11⊤x2x2⊤x2x1⊤y+1⊤x1x1⊤x11⊤x2x1⊤x2x2⊤y,=∗10×1⊤y+∗11×x1⊤y+∗12×x2⊤y=−x1⊤1x2⊤1x1⊤x2x2⊤x21⊤y+1⊤1x2⊤11⊤x2x2⊤x2x1⊤y−1⊤1x1⊤11⊤x2x1⊤x2x2⊤y,=∗20×1⊤y+∗21×x1⊤y+∗22×x2⊤y=x1⊤1x2⊤1x1⊤x1x2⊤x11⊤y−1⊤1x2⊤11⊤x1x2⊤x1x1⊤y+1⊤1x1⊤11⊤x1x1⊤x1x2⊤y
と計算できる.
定数倍
この表示を見るに,対応する変数の定数c倍は,係数に1/c倍の影響がある.対応しない変数の定数倍は影響なしである.
例えば,β1については
(detX⊤X)β^1=−x1⊤1x2⊤1x1⊤x2x2⊤x21⊤y+1⊤1x2⊤11⊤x2x2⊤x2x1⊤y−1⊤1x1⊤11⊤x2x1⊤x2x2⊤y
であり,対応する説明変数データx1がc倍(x1→cx1)となった場合は,左辺のdetX⊤Xがc2倍になり,右辺は行列式の性質から
==−(cx1)⊤1x2⊤1(cx1)⊤x2x2⊤x21⊤y+1⊤1x2⊤11⊤x2x2⊤x2(cx1)⊤y−1⊤1(cx1)⊤11⊤x2(cx1)⊤x2x2⊤y−cx1⊤1x2⊤1x1⊤x2x2⊤x21⊤y+c1⊤1x2⊤11⊤x2x2⊤x2x1⊤y−c1⊤1x1⊤11⊤x2x1⊤x2x2⊤yc(−x1⊤1x2⊤1x1⊤x2x2⊤x21⊤y+1⊤1x2⊤11⊤x2x2⊤x2x1⊤y−1⊤1x1⊤11⊤x2x1⊤x2x2⊤y)
とc倍となる.これによってβ^1は1/c倍となる.一方,x2がc倍となった場合は左辺のdetX⊤Xも右辺もc2倍となるため,β^1は不変である.
定数の足し引き
β1について
β1について考察する.
(detX⊤X)β^1=−(i)x1⊤1x2⊤1x1⊤x2x2⊤x21⊤y+(ii)1⊤1x2⊤11⊤x2x2⊤x2x1⊤y−(iii)1⊤1x1⊤11⊤x2x1⊤x2x2⊤y
であり,対応する説明変数データx1が一律bだけ足される(x1→x1+b1)場合,左辺のdetX⊤Xは不変である.右辺の変換の様子を見ていこう.
(i)→(x1+b1)⊤1x2⊤1(x1+b1)⊤x2x2⊤x21⊤y=(i)+b1⊤1x2⊤11⊤x2x2⊤x21⊤y,
(ii)→1⊤1x2⊤11⊤x2x2⊤x2(x1+b1)⊤y=(ii)+b1⊤1x2⊤11⊤x2x2⊤x21⊤y,
(iii)→1⊤1(x1+b1)⊤11⊤x2(x1+b1)⊤x2x2⊤y=(iii)+b=01⊤11⊤11⊤x21⊤x2x2⊤y
より,変換後の右辺は
−{(i)+b1⊤1x2⊤11⊤x2x2⊤x21⊤y}+{(ii)+b1⊤1x2⊤11⊤x2x2⊤x21⊤y}−(iii)
となって結局−(i)+(ii)−(iii)のまま不変である.
対応しない説明変数データx2が一律bだけ足される(x2→x2+b1)場合,
(i)→x1⊤1(x2+b1)⊤1x1⊤(x2+b1)(x2+b1)⊤(x2+b1)1⊤y=x1⊤1x2⊤1x1⊤(x2+b1)x2⊤(x2+b1)1⊤y+bx1⊤11⊤1x1⊤(x2+b1)1⊤(x2+b1)1⊤y=(i)+b=0x1⊤1x2⊤1x1⊤1x2⊤11⊤y+bx1⊤11⊤1x1⊤x21⊤x21⊤y+b2=0x1⊤11⊤1x1⊤11⊤11⊤y,
(ii)→1⊤1(x2+b1)⊤11⊤(x2+b1)(x2+b1)⊤(x2+b1)x1⊤y=1⊤1x2⊤11⊤(x2+b1)x2⊤(x2+b1)x1⊤y+b=01⊤11⊤11⊤(x2+b1)1⊤(x2+b1)x1⊤y=(ii)+b=01⊤1x2⊤11⊤1x2⊤1x1⊤y,
(iii)→1⊤1x1⊤11⊤(x2+b1)x1⊤(x2+b1)(x2+b1)⊤y=1⊤1x1⊤11⊤(x2+b1)x1⊤(x2+b1)x2⊤y+b1⊤1x1⊤11⊤(x2+b1)x1⊤(x2+b1)1⊤y=(iii)+b=01⊤1x1⊤11⊤1x1⊤1x2⊤y+b1⊤1x1⊤11⊤x2x1⊤x21⊤y+b2=01⊤1x1⊤11⊤1x1⊤11⊤y
となる.よって変換後の右辺は
−{(i)+bx1⊤11⊤1x1⊤x21⊤x21⊤y}+(ii)−{(iii)+b1⊤1x1⊤11⊤x2x1⊤x21⊤y}
となりやはり不変である.
切片項
切片β^0について考察する.
(detX⊤X)β^0=(i)′x1⊤x1x2⊤x1x1⊤x2x2⊤x21⊤y−(ii)′1⊤x1x2⊤x11⊤x2x2⊤x2x1⊤y+(iii)′1⊤x1x1⊤x11⊤x2x1⊤x2x2⊤y
であり,説明変数データx1が一律bだけ足される(x1→x1+b1)場合,左辺のdetX⊤Xは不変である.右辺の変換の様子を見ていこう.
(i)′→(x1+b1)⊤(x1+b1)x2⊤(x1+b1)(x1+b1)⊤x2x2⊤x21⊤y=(x1+b1)⊤x1x2⊤x1(x1+b1)⊤x2x2⊤x21⊤y+b(x1+b1)⊤1x2⊤1(x1+b1)⊤x2x2⊤x21⊤y=(i)′+b1⊤x1x2⊤x11⊤x2x2⊤x21⊤y+bx1⊤1x2⊤1x1⊤x2x2⊤x21⊤y+b21⊤1x2⊤11⊤x2x2⊤x21⊤y,
(ii)′→1⊤(x1+b1)x2⊤(x1+b1)1⊤x2x2⊤x2(x1+b1)⊤y=1⊤(x1+b1)x2⊤(x1+b1)1⊤x2x2⊤x2x1⊤y+b1⊤(x1+b1)x2⊤(x1+b1)1⊤x2x2⊤x21⊤y=(ii)′+b1⊤1x2⊤11⊤x2x2⊤x2x1⊤y+b1⊤x1x2⊤x11⊤x2x2⊤x21⊤y+b21⊤1x2⊤11⊤x2x2⊤x21⊤y
(iii)′→1⊤(x1+b1)(x1+b1)⊤(x1+b1)1⊤x2(x1+b1)⊤x2x2⊤y=1⊤(x1+b1)x1⊤(x1+b1)1⊤x2x1⊤x2x2⊤y+b=01⊤(x1+b1)1⊤(x1+b1)1⊤x21⊤x2x2⊤y=(iii)′+b1⊤1x1⊤11⊤x2x1⊤x2x2⊤y
となる.よって右辺は
−+=={(i)′+b1⊤x1x2⊤x11⊤x2x2⊤x21⊤y+bx1⊤1x2⊤1x1⊤x2x2⊤x21⊤y+b21⊤1x2⊤11⊤x2x2⊤x21⊤y}{(ii)′+b1⊤1x2⊤11⊤x2x2⊤x2x1⊤y+b1⊤x1x2⊤x11⊤x2x2⊤x21⊤y+b21⊤1x2⊤11⊤x2x2⊤x21⊤y}{(iii)′+b1⊤1x1⊤11⊤x2x1⊤x2x2⊤y}(i)′−(ii)′+(iii)′+b{x1⊤1x2⊤1x1⊤x2x2⊤x21⊤y−1⊤1x2⊤11⊤x2x2⊤x2x1⊤y+1⊤1x1⊤11⊤x2x1⊤x2x2⊤y}(i)′−(ii)′+(iii)′−b(detXX⊤)β^1
へと変換される.つまり,
β^0→β^0−bβ^1
となる.
行列式の計算による確認まとめ
以上のことを表にまとめると以下の通り.
定数倍
回帰係数 |
1→c1 |
x1→cx1 |
x2→cx2 |
β^0 |
1/c倍 |
不変 |
不変 |
β^1 |
不変 |
1/c倍 |
不変 |
β^2 |
不変 |
不変 |
1/c倍 |
定数の足し引き
回帰係数 |
1→1+b1 |
x1→x1+b1 |
x2→x2+b1 |
β^0 |
1/(1+b)倍 |
β^0−bβ^1 |
β^0−bβ^2 |
β^1 |
不変 |
不変 |
不変 |
β^2 |
不変 |
不変 |
不変 |
複数の説明変数データを同時に足し引きx1→x1+b11,x2→x2+b21する場合でも,説明変数データをひとつずつ動かすことと同じなので,変換は上記の表の通りとなる.
複合的な変換
また,x1→cx1+d1という変換はcx1+d1=c(x1+cd1)とみることで定数の足し引きをした後に定数倍をすることに等しい.よってそれぞれの回帰係数は
β^0→β^0−cdβ^1,β^1→cβ^1,β^2→β^2
と変換される.
つまり,変換前の線形回帰モデルが
fX(t1,t2)=β^0+β^1t1+β^2t2
だとすると,変換後の線形回帰モデルは
fXnew(t1new,t2new)=(β^0−cdβ^1)+cβ^1t1new+β^2t2new=β^0+β^1(ct1new−d)+β^2t2new
となるということである.
Pythonによる確認
以上のことをPythonライブラリstatsmodels
で確かめてみる.
ここでalcohol
について(線形)変換5 * alcohol + 20
を行ってみる.
となる.それぞれの回帰係数は以下のようになる.
となっており,違いがあるのは切片const
と,変換を施した変数alcohol
に対応する回帰係数である.
変換前後の回帰係数の比較
変換前後のalcohol
の回帰係数を比較すると,
となって「変換後のalcohol
の回帰係数 ≒ 変換前のalcohol
の回帰係数の1/c倍」が成り立っている.
変換前後の切片の比較
変換前後の切片const
を比較すると,
となって「変換後の切片const
≒ 変換前の切片const
- 変換前のalcohol
の回帰係数のd/c倍」が成り立っている.
変換前後の予測値と決定係数の比較
変換前後の予測値については
となって等しい.
また,変換前後の決定係数R2についても
と等しいことが確認できる.
まとめ
以上から,記事の最初に触れたように,計画行列X=(1x1x2)から作成した線形回帰モデル
fX(t1,t2)=β^0+β^1t1+β^2t2
があるとしたとき,x1→x1new=c1x1+d11, x2→x2new=c2x2+d21という変換を行ったデータからなる計画行列Xnew=(1x1newx2new)で作られる線形回帰モデルは
fXnew(t1new,t2new)=β^0+β^1(c1t1new−d1)+β^2(c2t1new−d2)
である.
また,入力が学習データである限りにおいては,予測値は学習データの(線形)変換をしても不変である:fXnew(xi1new,xi2new)=fX(xi1,xi2).よって特に決定係数R2も不変である.
Discussion