Help us understand the problem. What is going on with this article?

重回帰分析における回帰係数の計算方法

はじめに

 今回は, 重回帰分析における偏回帰係数(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$の計算方法を確認する.

参考文献

neuman71
Data Analyst: Qiitaでの投稿内容は個人の見解の表明であり, 所属組織の見解を代表しない.
lifull
日本最大級の不動産・住宅情報サイト「LIFULL HOME'S」を始め、人々の生活に寄り添う様々な情報サービス事業を展開しています。
https://lifull.com/
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした