重回帰分析 (Lisp Advent Calendar 2023)

さて、車輪の再発明を続けて、ようやく逆行列の計算までたどりつきました。今回は逆行列を使った応用をしたいと思います。

重回帰分析

逆行列の応用例として、重回帰分析を取り上げます。

重回帰分析は、ひとつの量的な応答変数を、複数の量的な説明変数で予測・説明するための線形モデルです。例えば、アパートの家賃（応答変数、単位：円）を、駅からの徒歩時間（説明変数、単位：分）と延べ床面積（説明変数、単位：m²）から推測するようなときに用いることができます。数式は 𝑦 = 𝛽₀ + 𝛽₁ 𝑥₁ + 𝛽₂ 𝑥₂ + ... + 𝛽ₙ 𝑥ₙ + 𝜖 のように書かれます。家賃の例の場合、家賃 = 係数0 + 係数1×徒歩時間 + 係数2×面積 + 誤差のようになるでしょうか。係数0は賃貸物件としてのベースラインの金額、係数1は徒歩時間が1分増えるごとに増える（あるいは減る）金額、係数2は延べ床面積が1 m²増えるごとに上がる金額を表します。また、誤差には時間と面積だけでは説明しきれない要素（築年数、階数、インターネット回線の品質など）や物件ごとの個体差が含まれます。多くのデータをもとに係数の値を求めるのが重回帰分析です。係数0と係数1だけを求めればいい場合（説明変数がひとつしかないとき）を単回帰分析と区別することもありますが、やっていることに違いはありません。

係数が決まったとして、そこから計算できる予測値は実際の値からは多少ずれます。そのずれの総量を最小化すれば最適な回帰式が得られるはずです。その発想で係数を求めるのが最小二乗法です。イメージとしては、下のグラフのように各データ点から近似直線までの距離（つまり誤差）の二乗和を最小にするように近似直線を定めるかんじです。

誤差の二乗和を最小にするためには微分を使って計算を進めることになります。やたらとめんどくさそうですが、実は簡単な行列計算に帰着します。英語版WikipediaのLinear Regressionの項に詳しく書かれていますが、まず説明変数𝑥₁ ... 𝑥ₙを以下のような行列にまとめます。ここで𝑛は説明変数の数、𝑚はデータの数となっています。行列の1列目に並んでいる1は𝛽₀と掛け算されるために追加されています。このようにすると 𝑦 = 𝛽₀1 + 𝛽₁ 𝑥₁ + 𝛽₂ 𝑥₂ + ... + 𝛽ₙ 𝑥ₙ + 𝜖 を 𝑦 = 𝑋𝛽 + 𝜖 と書け、𝑋𝛽 の部分は行列とベクトルの積で一発計算ができるようになるわけです。

$\displaystyle X = \begin{bmatrix} \mathbf{x}^\mathsf{T}_1 \\ \mathbf{x}^\mathsf{T}_2 \\ \vdots \\ \mathbf{x}^\mathsf{T}_m \end{bmatrix} = \begin{bmatrix} 1 & x_{11} & \cdots & x_{1n} \\ 1 & x_{21} & \cdots & x_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{m1} & \cdots & x_{mn} \end{bmatrix}$

ここで予測誤差である 𝜖 = 𝑦 - 𝑋𝛽 の総和を最小にしたいです。𝜖 の値は正にも負にもなるので、単純に総和を計算するとゼロになってしまいます。そこで、𝜖 の二乗和を使います。つまり、𝑋𝛽 で予測した値と、データで与えられている応答変数 𝑦 との差の二乗和ということです。以下の式を微分して最小になる 𝛽 を探ります。

$\displaystyle \|X\vec{\beta} - Y\|^2$

ここでは導出は省略しますが、実は以下のような行列の転置、逆行列、行列の積という操作だけで 𝛽 の最適値が求まってしまいます。導出については前述のWikipedia記事を参照してください。

$\displaystyle \vec{\hat{\beta}} = \left(X^\textsf{T} X\right)^{-1} X^\textsf{T} Y$

Lisp実装

この式を単純にCommon Lispコードにしたのが以下です。

(defun mra (y X)
  (let ((XX (make-array `(,(array-dimension X 0) ,(1+ (array-dimension X 1))) :initial-element 1)))
    (loop for i from 0 below (array-dimension X 0)
          do (loop for j from 0 below (array-dimension X 1)
                   do (setf (aref XX i (1+ j)) (aref X i j))))
    (matrix-multiply (matrix-multiply (inv (matrix-multiply (matrix-transpose XX) XX)) (matrix-transpose XX)) y)))

なんの工夫もない、単純なコードですね。まず、第1列に1が並び、それ以降にXの内容が入っているXX行列を作り、それに上記式を適用した結果を返すだけです。matrix-multiplyとmatrix-transposeはこのブログでは紹介したことはありませんが、二次元配列で与えられた行列についての積と転置の単純な実装になっています。

実行例

南風原朝和『心理統計学の基礎 https://amzn.to/48K5fKY』（有斐閣アルマ，2002年）の第8章に掲載されている例を使って実行してみましょう。この本はt検定、分散分析、回帰分析、因子分析などの背景について、数式を使いつつも直感的に理解できるように書かれている良書です。

p.226に掲載されている「母親の協調性価値得点 (𝑥₁)」「通年年数 (𝑥₂)」「子どもの協調性得点 (𝑦)」を抜き出しました。まずは、Rで実行してみます。重回帰分析はモデル式をlmで評価させるだけです。

kachi <- c(12, 12, 7, 17, 14, 9, 10, 13, 15, 12, 12, 15, 11, 14, 17, 17, 16, 15, 15, 10, 12, 9, 12, 12, 19, 11, 14, 15, 15, 15, 16, 15, 12, 10, 11, 12, 15, 13, 15, 12, 12, 12, 13, 17, 13, 11, 14, 16, 12, 12)
nensu <- c(2, 2, 2, 3, 2, 2, 3, 3, 3, 1, 3, 3, 2, 2, 4, 2, 4, 3, 4, 2, 2, 1, 2, 2, 4, 2, 3, 2, 3, 3, 2, 3, 2, 2, 3, 1, 2, 3, 2, 2, 2, 3, 3, 3, 2, 3, 2, 4, 2, 2)
kyouchou <- c(6, 11, 11, 13, 13, 10, 10, 15, 11, 11, 16, 14, 10, 13, 12, 15, 16, 14, 14, 8, 13, 12, 12, 11, 16, 9, 12, 13, 13, 14, 12, 15, 8, 12, 11, 6, 12, 15, 9, 13, 9, 11, 14, 12, 13, 9, 11, 14, 16, 8)
res <- lm(kyouchou ~ kachi + nensu)
summary(res)

以下のように、子どもの協調性得点 = 5.1716 + 0.3027×母親の協調性価値得点 + 1.1259×通園年数という結果になりました。

Call:
lm(formula = kyouchou ~ kachi + nensu)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0563 -1.4780 -0.0563  1.4751  4.9437 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   5.1716     1.6647   3.107   0.0032 **
kachi         0.3027     0.1467   2.064   0.0446 * 
nensu         1.1259     0.4713   2.389   0.0210 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.12 on 47 degrees of freedom
Multiple R-squared:  0.3137,    Adjusted R-squared:  0.2845 
F-statistic: 10.74 on 2 and 47 DF,  p-value: 0.000144

今回のLisp実装でこの3つの係数が得られれば成功ということにします。

(defparameter *kachi* #(12 12 7 17 14 9 10 13 15 12 12 15 11 14 17 17 16 15 15 10 12 9 12 12 19 11 14 15 15 15 16 15 12 10 11 12 15 13 15 12 12 12 13 17 13 11 14 16 12 12))
(defparameter *nensu* #(2 2 2 3 2 2 3 3 3 1 3 3 2 2 4 2 4 3 4 2 2 1 2 2 4 2 3 2 3 3 2 3 2 2 3 1 2 3 2 2 2 3 3 3 2 3 2 4 2 2))
(defparameter *kyouchou* #(6 11 11 13 13 10 10 15 11 11 16 14 10 13 12 15 16 14 14 8 13 12 12 11 16 9 12 13 13 14 12 15 8 12 11 6 12 15 9 13 9 11 14 12 13 9 11 14 16 8))

(mra (matrix-transpose *kyouchou*) (cbind (matrix-transpose *kachi*) (matrix-transpose *nensu*)))

matrix-transposeは行列の転置ですが、入力がsimple vectorの場合はn×1のarrayにします。cbindは二つの2次元arrayを横方向につなげる関数です。この結果、#2A((64133/12401) (11263/37203) (13962/12401))が得られました。64133÷12401 = 5.171599、11263÷37203 = 0.3027444、13962÷12401 = 1.1258769なので、Rでの計算結果と同じになりました。

ここでは決定係数（R²）の計算はしていませんが、𝑋𝛽の分散を𝑦の分散で割ると得られます。今回のデータでは0.31367722です。