はじめに
今回は, 重回帰分析における偏回帰係数(partial regression coefficient)の最小二乗法(ordinary least squares)による計算過程を見ていく. 振り返るまでもないが, 重回帰分析とは, 1つの目的変数$Y$を複数の説明変数$X_i$を用いて予測することである.
Y=\alpha+\beta_1X_1+\beta_2X_2 ... +\beta_iX_i\\
また, 最小二乗法による各偏回帰係数$\beta_i$の計算を理解する上で必要となる線形代数の前提知識は以下となる.
- 行列の写像・変換
- 余因子行列
- 逆行列
それでは早速, 例題とともに今回の主題について確認していく.
最小二乗法で切片αと回帰係数βを求める式
まず, 複数の説明変数を持つ回帰式における切片$\alpha_i$と偏回帰係数$\beta_i$を求める際, 最小二乗法に従って計算すると以下2つの式が導出できる(※途中式を省略).
\mu_Y=\alpha+\beta_1\mu_{X_1}+\beta_2\mu_{X_2} ... +\beta_i\mu_{X_i}\\
A\vec{W}=\vec{C}\\
このとき, $A$は, 各説明変数$X_i$の共分散$Cov(X_i, X_i)$と分散$Var(X_i)$が並ぶ「分散共分散行列」を表している.
A = \left(
\begin{array}{ccc}
Var(X_1) & Cov(X_1, X_2) & Cov(X_1, X_3) & ... & Cov(X_1, X_i) \\
Cov(X_2, X_1) & Var(X_2) & Cov(X_2, X_3) & ... & Cov(X_2, X_i) \\
Cov(X_3, X_1) & Cov(X_3, X_2) & Var(X_3) & ... & Cov(X_3, X_i) \\
... & ... & ... & ... & ... \\
Cov(X_i, X_1) & Cov(X_i, X_2) & Cov(X_i, X_3) & ... & Var(X_i)
\end{array}
\right)
また, $\vec{W}$は偏回帰係数が並んだ縦ベクトルを表している.
\vec{W} = \left(
\begin{array}{ccc}
\beta_1 \\
\beta_2 \\
\beta_3 \\
... \\
\beta_i
\end{array}
\right)
そして, $\vec{C}$は目的変数$Y_i$と説明変数$X_i$の共分散$Cov(Y, X_i)$が並んだ縦ベクトルを表している.
\vec{C} = \left(
\begin{array}{ccc}
Cov(Y, X_1) \\
Cov(Y, X_2) \\
Cov(Y, X_3) \\
... \\
Cov(Y, X_i)
\end{array}
\right)
いま求めたいのは, 偏回帰係数の縦ベクトル$\vec{W}$と切片のスカラー$\alpha$の2つなので, それらを左辺に持つ様に式変形をする.
\vec{W}=\vec{C} A^{-1}\\
\alpha=\mu_Y-(\beta_1\mu_{X_1}+\beta_2\mu_{X_2} ... +\beta_i\mu_{X_i})\\
以上で, 偏回帰係数$\beta_i$と切片$\alpha$の計算の基となる式が明らかになったので, 以降は実際に行列を使った計算に移る.
A. Rの重回帰モデルのメソッドで計算する
今回は例題として, 説明変数が2つの重回帰モデルを想定する. 一応ではあるが, 説明変数の$x_1$は身長, $x_2$は体重, 目的変数の$y$は年収を意識した数値となっている. なお, 各変数の値は正規分布に従う乱数を適当に生成した.
# サンプルデータを作成する
x1<-rnorm(10,160,5)
x2<-rnorm(10,60,3)
y<-rnorm(10,40000000,500000)
df<-data.frame("x1"=x1, "x2"=x2, "y"=y)
df
>>> output:
| $x_1$ | $x_2$ | $y$ |
|---|---|---|
| 154.9222 | 61.04467 | 39060658 |
| 153.4911 | 62.06968 | 40616943 |
| 165.6492 | 60.43379 | 39436220 |
| 164.7504 | 62.62436 | 39665936 |
| 168.0310 | 62.25878 | 40059877 |
| 169.3117 | 57.76524 | 40382747 |
| 163.9869 | 59.86037 | 40399481 |
| 163.6964 | 61.54883 | 39450105 |
| 167.7211 | 61.21164 | 40374236 |
| 159.6782 | 59.34171 | 40028748 |
以降で公式を用いて手計算する前に, Rの線形回帰メソッドlmで, 上記のデータに当てはまる切片(Intercept)と説明変数の偏回帰係数xiを最小二乗法によって計算する.
# 線形回帰モデル: Linear Model
model<-lm(formula=y~x1+x2, data=df)
model
>>> output:
//モデル式と切片, 回帰係数
>>> Call:
>>> lm(formula = y ~ x1 + x2, data = df)
>>>
>>> Coefficients:
>>> (Intercept) x1 x2
>>> 25389813 68090 65026
B. Rで公式を基に行列を使って計算する
改めて, 定式化すると以下によって切片$\alpha_i$と回帰係数$\beta_i$が求まる.
\vec{W}=\vec{C} A^{-1}\\
\alpha=\mu_Y-(\beta_1\mu_{X_1}+\beta_2\mu_{X_2} ... +\beta_i\mu_{X_i})\\
まず, 説明変数$X_i$同士の分散共分散行列(variance-covariance matrix)となる$A$を求める.
# 説明変数の分散共分散行列(variance-covariance matrix)
vcm=var(df)[1:2,]
a<-vcm[1:2,1:2]
a
>>> output:
| $A$ | $x_1$ | $x_2$ |
|---|---|---|
| $x_1$ | 30.925886 | -7.657957 |
| $x_2$ | -7.657957 | 4.320784 |
次に, 説明変数$X_i$と目的変数$Y$の共分散行列の縦ベクトルとなる$\vec{C}$を求める.
# 説明変数と目的変数の共分散行列
vcm=var(df)[1:2,]
c<-vcm[1:2,3]
c
>>> output:
| $\vec{C}$ | $y$ |
|---|---|
| $x_1$ | 1607784.88015929 |
| $x_2$ | -240467.744191122 |
さらに, 切片を除く回帰係数を並べた縦ベクトル$\vec{W}$を求めたい. 今回はたまたま説明変数が2つ(2×2行列)だったため, 手計算してみる. 逆行列については, ヨビノリのたくみ先生の説明がわかりやすい.
ある正方行列$A$($|A|\neq0$)において, 余因子行列を用いて逆行列を求める公式は以下となる. なお, $|A|$は行列式を, $\tilde{A}$は余因子行列を指す.
A^{-1}=\frac{1}{|A|}\tilde{A}
今回は説明変数$X_i$が2つなので, $X_i$による分散共分散行列$A$は2×2の正方行列となる. ちなみに2×2行列$A$を以下で表した場合, 行列式$|A|$と余因子行列$\tilde{A}$は同様に以下となる.
A=
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\\
|A|=
\begin{vmatrix}
a & b \\
c & d
\end{vmatrix}
=ad-bc\\
\tilde{A}=
\begin{pmatrix}
d & -b \\
-c & a
\end{pmatrix}
それでは早速, 行列式$|A|$と余因子行列$\tilde{A}$から逆行列$A^{-1}$を求める.
# 行列式を求める
det_a<-a[1,1]*a[2,2]-a[1,2]*a[2,1]
# 余因子行列を求める
til_a<-cbind(c(a[2,2],-a[2,1]), c(-a[1,2],a[1,1]))
#余因子行列を1/行列式の解でスカラー倍することで分散共分散行列aの逆行列を求める
inv_a<-(1/det_a)*til_a
inv_a
>>> output:
| $A^{-1}$ | $x_1$ | $x_2$ |
|---|---|---|
| $x_1$ | 0.0576260 | 0.1021337 |
| $x_2$ | 0.1021337 | 0.4124564 |
なお, 2×2行列以上の行列における逆行列の計算など少し面倒な場合には, Rの標準関数を使うことをお勧めする. Rでは逆行列を求めるための関数solve()が存在する.
#分散共分散行列sigmaの逆行列を求める
c%*%solve(a)
>>> output:
| $A^{-1}$ | $x_1$ | $x_2$ |
|---|---|---|
| $x_1$ | 0.0576260 | 0.1021337 |
| $x_2$ | 0.1021337 | 0.4124564 |
ようやく, 回帰係数ベクトル$\vec{W}$が求められる. 上述の通り$\vec{W}$は, 説明変数$X$と目的変数$Y$の共分散行列$\vec{C}$と分散共分散行列の逆行列$A^{-1}$の内積となる. なお, Rで内積を求める関数は%*%である.
w<-c%*%inv_a
w
>>> output:
| $\beta_i$ | $coef$ |
|---|---|
| $w_1$ | 68090.37 |
| $w_2$ | 65026.48 |
上記で, $\vec{W}$の各要素, つまり各説明変数$x_1,x_2$に対する回帰係数が求められた. 最後に切片$\alpha$を求める.
alpha<-mean(df$y)-(w[1]*mean(df$x1)+w[2]*mean(df$x2))
alpha
>>> output:
>>> 25389812.79
よって, 改めて公式から求められた切片$\alpha$と回帰係数$\beta_i$を表にまとめると以下となる.
| $coef$ | |
|---|---|
| $\alpha$ | 25389812.79 |
| $\beta_1$ | 68090.37 |
| $\beta_2$ | 65026.48 |
確認のために, 最初にlmメソッドで求めた回帰係数の値と比較してみる.
# モデルから回帰係数のみを抜き出す.
model$coefficients
>>> output:
//Coefficients:
>>> (Intercept) x1 x2
>>> 25389813 68090 65026
公式から求めた値の四捨五入と, lmメソッドで求めた値が一致していることがわかる.
さいごに
今回は, 最小二乗法を使った偏回帰係数の計算過程を見てきた. これらを理解するには線形代数の基礎知識となる, 余因子展開や余因子行列や, 逆行列や行列の内積などの理解が必要となった. また, 共分散行列(および共分散分散行列)など統計学の概念も必要となる. 次回は最尤推定法による偏回帰係数$\beta$の計算方法を確認する.