相関係数、回帰直線

1 散布図

１個体が複数の変量で表される（データがベクトルで与えられる）場合で、２変量のデータ（平面上のデータ）であれば、

(x_{1}, y_{1}), (x_{2}, y_{2}), ⋯, (x_{n}, y_{n})

平面状にプロットした図を散布図という。また３変量のデータ（空間上のデータ）であれば、つぎのような形がデータである。

(x_{1}, y_{1}, z_{1}), (x_{2}, y_{2}, z_{2}) ⋯, (x_{n}, y_{n}, z_{n})

このようなデータを整理するには、いろいろな方法があるが、 2 変量データについて、つぎのような概念が重要である。

[共分散]

2 変量の間の全体的なばらつき具合を考える。一つの変数の増加について、他方も増加であれば、正の値、逆の場合には、負の値をとる。増減がはっきりしないときには、ゼロに近い値となる。つぎの相関係数と対応する。

s_{XY} = \frac{1}{n} ∑_{n}\{(X_{i} − \overline{X})(Y_{i} − \overline{Y})\}

2 相関係数

相関係数

共分散を各変量の標準偏差で基準化したもの。

ρ = \frac{s_{XY}}{s_{X}s_{Y}} = \frac{s_{XY}}{\sqrt{s_{X}^{2}s_{Y}^{2}}} = \frac{1}{n} ∑_{n} (\frac{X_{i} − \overline{X}}{s_{X}} ) (\frac{Y_{i} − \overline{Y}}{s_{Y}} )

−1 ≤ r ≤ 1 が成り立ち、絶対値が１に近いときには直線傾向が強く、絶対値がゼロに近いときには、直線傾向がみられないことを意味する。

3 曲線の当てはめ

曲線の当てはめ、最小 2 乗法

散布図からある程度の直線的な関係があるとき、これらの 2 変量データに直線をつぎの基準で当てはめる。

∑_{i}(Y_{i} − a − bX_{i})^{2}

を最小にするような a, b を求める。

[回帰直線] この方法を最小 2 乗法といい、求めた直線を Y の X への回帰直線という。最小 2 乗法で定めた直線の係数を回帰係数という。

b = \frac{s_{XY}}{s_{X}^{2}} = r \frac{s_{Y}}{s_{X}} , a = \overline{Y} − r \frac{s_{Y}}{s_{X}} \overline{X}

問 3.1 相関係数は線形変換をしても同じ値になることを示しなさい。つまり (X_{1}, Y_{1}), ⋯, (X_{n}, Y_{n}) と (aX_{1} + b, cY_{1} + d), ⋯, (aX_{n} + b, cY_{n} + d) （ただし ac > 0）とは同じ相関係数をもつ。

問 3.2 与えられた 2 変量データ X, Y が順位データであるとき、すなわちそれぞれが 1, 2, 3, ⋯, n までを順位づけられていて、並び替えたものである。ただし同順位はないとする。このときには相関係数が

r = 1 − \frac{6 ∑_{i} d_{i}^{2}}{n(n^{2} − 1)}, (d_{i} = X_{i} − Y_{i})

で与えられることを示しなさい。これをスピアマンの順位相関係数という。

4 多変量データの解析

4.1 回帰分析

いま２つの変数 (x, y) の間に一次の関数関係 y = α + βx があるものとする。 x の値を与えると、係数 α, β が既知ならば、 y の値を定めることができる。具体的な例として、測定誤差など・・・しかし、いまこの構造式の中に確率変動の要素が加わると考えよう。観測値 y のなかには確率変数が含まれるとする。いま n 個のデータ

(x_{1}, y_{1}), (x_{2}, y_{2}), ⋯, (x_{n}, y_{n})

について

y_{i} = α + βx_{i} + ϵ_{i}, i = 1, 2, ⋯

なる関係があり、この ϵ_{i} は、各 i について独立で同一の分布 N(0, σ^{2}0) にしたがうとする。すなわち各 x_{i} が指定されると、上の式で定まる y_{i} は ϵ_{i} によって確率変動し、平均 μ_{i} は、 α + βx_{i} 、分散 σ_{i}^{2} = σ^{2} (一定) の正規分布にしたがう。

このように２つの変数間が一次式で結ばれて変動するような関係を線形回帰モデル (Linear Regression Model) とよぶ。変数 x のことを説明変数（explanatory variable) または独立変数 (independent variable) とよび、これに対する変数 y のことを従属変数（dependent variable) という。一般に、説明変数は非確率的なものとする。係数 α, β を回帰係数（regression coefficient) とよび、これは未知の母数である。 x を与えて観測された y の値をデータとして、この回帰係数を推定することが問題となる。

回帰係数の推定には、ガウスによる最小２乗の原理 (principle of minimum square) に基づく。いま n 個のデータ：

(x_{1}, y_{1}), (x_{2}, y_{2}), ⋯, (x_{n}, y_{n})

が与えられたとする。未知パラメ−タ α, β を

S(α，β) = ∑_{i=1}^{n}(y_{i} − (α + βx_{i}))^{2}

が最小となるように決定しようと考える。このような方法を最小２乗法 (least square method) とよぶ。データは固定されていて、 α と β がどんな値であっても

S(α̂，β̂) ≤ S(α，β)

を満たす最小値 α̂, β̂ を求めようとするわけで、このようにして定まる係数の値を最小２乗推定値とよぶ。 y = α + βx は直線の方程式であるから、係数を変化させることは、いろいろ直線を動かすことで、最小２乗法とは各点での残差 y_{i} − α − βx_{i} の２乗和が小さくなるよう直線をあてはめようとすることである。

実際、 α̂, β̂ を決定するには、微分積分学の極値問題を使って、まず必要条件として S(α, β) に対する２変数の微分係数をゼロとした方程式を解き、最小値の候補を求める。すなわち、

\begin{aligned}\frac{∂S}{∂α} & = & −2 ∑_{i=1}^{n}(y_{i} − α − βx_{i}) = 0 \\ \frac{∂S}{∂β} & = & −2 ∑_{i=1}^{n} x_{i}(y_{i} − α − βx_{i}) = 0 \end{aligned}

これを整理すれば

\{ \begin{matrix}nα + (∑ x_{i})β = ∑ y_{i} & · · · · · · (1) \\ (∑ x_{i})α + (∑ x_{i}^{2})β = ∑ x_{i}y_{i} & · · · · · · (2) \end{matrix}

を得る。これを正規方程式 (normal equation) とよぶ。 α と β に関する２元連立線形方程式である。 (1) より、 α = \overline{y} − β\overline{x} (ただし\overline{x} は, x_{1}, ⋯, x_{n} の平均値、\overline{y} は y_{1}, ⋯, y_{n} の平均値) となるから、 (2) に代入して、 β を得る。さらに α も得られる。こうして得た α, β の値は

\begin{aligned}α̂ & = \overline{y} − β̂\overline{x} = \overline{y} − \frac{σ_{xy}}{σ_{xx}} \overline{x} \\ β̂ & = \frac{_{∑i}(x_{i} − \overline{x})(y_{i} − \overline{y})}{_{∑i}(x_{i} − \overline{x})^{2}} = \frac{σ_{xy}}{σ_{xx}} \end{aligned}

ここで σ_{xx} と σ_{xy} はそれぞれ x, y についての分散、共分散を表すとする。さらに、いま求めた値に対して極値問題の十分条件 S_{αβ}^{2} − S_{αα}S_{ββ} < 0 を調べ、これを満たしているからことがわかる。以上により求める直線の方程式は

y − \overline{y} = \frac{σ_{xy}}{σ_{xx}} (x − \overline{x})

で与えられ、この回帰直線は点 (X, \overline{y}) を通り、傾きが\frac{σ_{xy}}{σ_{xx}} である。つぎに最小２乗法の幾何的意味を考えよう。与えられた構造を表すモデルの式は

y_{i} = α + βx_{i} + e_{i}, i = 1, ⋯, n

と書けていたから、まずこれをベクトルや行列を用いて表す。

[ y_{1} ; y_{2} ; ⋮ ; y_{n} ] = [ \begin{matrix}1 & x_{1} \\ 1 & x_{2} \\ ⋮ & ⋮ \\ 1 & x_{n} \end{matrix} ] [ \begin{matrix}α \\ β \end{matrix} ] + [ e_{1} ; e_{2} ; ⋮ ; e_{n} ]

つまり、行列とベクトルを使って

y = Xβ + e

が線形回帰モデルの表現である。ただし

y = [ y_{1} ; y_{2} ; ⋮ ; y_{n} ] , X = [ \begin{matrix}1 & x_{1} \\ 1 & x_{2} \\ ⋮ & ⋮ \\ 1 & x_{n} \end{matrix} ] , β = [ \begin{matrix}α \\ β \end{matrix} ] , e = [ e_{1} ; e_{2} ; ⋮ ; e_{n} ] .

行列の転置 (transpose) を表すため、左上に t をつけるとしよう。

tX = [ \begin{matrix}1 & 1 & ⋯ & 1 \\ x_{1} & x_{2} & ⋯ & x_{n} \end{matrix} ]

これを用いると、正規方程式は

[^{t}XX]β = ^{t}Xy

と書ける。２次の正方行列

tXX = [ \begin{matrix}n & ∑ x_{i} \\ ∑ x_{i} & ∑ x_{i}^{2} \end{matrix} ]

を考えると、その行列式 (determinant) は n ∑_{i} x_{i}^{2} − (∑ x_{i})^{2} で、これは n^{2} で割ると分散 σ_{xx} であるから、 0 でないとしてよい。したがって行列は正則である。逆行列が存在するから β = [^{t}XX]^{−1}^{t}Xy 回帰直線はしたがって y = X[^{t}XX]^{−1}^{t}Xy + e となる。この式において右辺の y の係数行列を P とおくと、

P = X [^{t}XX]^{−1} ^{t}X

この P はつぎを満たしている。

1) P^{2} = P (ベキ等という)

2) ^{t}P = P (対称という)

一般にこのような行列を係数として定まる線形写像を射影 (projection) という。 y = [y_{1}y_{2} ⋯ y_{n}] は n 次元ベクトル空間の点で、 P y は２次元ベクトル空間 (平面) の点で、 y を射影したベクトルとなる。

図

y = P y + e より e = (I − P)y I は単位行列。

残差ベクトル e は n 次元ベクトルの点 y より２次元ベクトル空間へおろして垂線を表す。もし

β = [^{t}XX]^{−1} ^{t}Xy

でない一般の β を用いると、 y = Xβ + e であり、このときの残差ベクトルは e = y − Xβ は e = (I − P)y と比べれば、ベクトルの大きさは長くなっている。これが最小２乗法の幾何的意味である。

さて α、 β のそれぞれの推定値 α̂、 β̂ を S(α, β) に代入した最小値は ∑_{i}(y_{i} − ŷ_{i})^{2} となり、これを残差平方和 (sum of square of residue) といい、 n − 2 で割ったものを残差分散 (variance of residue) とよぶ。

S_{e} = \frac{1}{n − 2} ∑_{i=1}^{n}(y_{i} − ŷ_{i})^{2}

この式はつぎのように変形できる。

∑_{i}(y_{i} − ŷ_{i})^{2} = ∑(y_{i} − \overline{y})^{2} − 2β̂ ∑(y_{i} − \overline{y})(x_{i} − \overline{x}) + β̂^{2} ∑(x_{i} − \overline{x})^{2}
= ∑(y_{i} − \overline{y})^{2} − β̂a^{2} ∑(x_{i} − \overline{x})^{2}
= ∑(y_{i} − \overline{y})^{2} − ∑(ŷ_{i} − \overline{y})^{2}

ゆえに

∑(y_{i} − \overline{y})^{2} = ∑(ŷ_{i} − \overline{y})^{2} + ∑(y_{i} − \overline{y})^{2}

が成り立つ。この式は y の全変動 ∑_{i}(y_{i} − \overline{y})^{2} が回帰による変動 ∑_{i}(ŷ_{i} − \overline{y})^{2} と残差変動 ∑_{i}(y_{i} − \overline{y})^{2} に分解されることを示している。変量 y の全変動のうち、回帰による変動の占める割合、

ρ = \frac{∑(ŷ_{i} − \overline{y})^{2}}{∑(y_{i} − \overline{y})^{2}}

を回帰の決定係数 (coefficient of determination) という。この ρ の値が１に近ければ、 y の変動の大部分が x に対する回帰によって説明されることを意味し、もし、 0 に近ければ y の変動の説明に x はほとんど無力であることになる。

また\frac{1}{n} ∑_{i}(ŷ_{i} − \overline{y})^{2} = \frac{2β̂}{n} ∑_{i}(y_{i} − \overline{y})(x_{i} − \overline{x}) − \frac{β̂^{2}}{n} ∑_{i}(x_{i} − \overline{x})^{2} = \frac{σ_{xy}^{2}}{σ_{xx}} より

ρ = σ_{xy}^{2}/\frac{σ_{xx}}{σ_{yy}} = \frac{σ_{xy}^{2}}{σ_{xx}σ_{yy}}

したがって σ_{x}^{2} = σ_{xx}, σ_{y}^{2} = σ_{yy} であるから

\sqrt{ρ} = \frac{σ_{xy}}{σ_{x}σ_{y}}

すなわち回帰の決定係数 ρ の平方根は x と y との相関係数となる．

4.2 主成分分析

主成分分析とは多変量データがもつ情報をなるべく損なうことなく、縮約するための統計的手法である。たとえば、学生の数学試験を前期と後期の２回おこなった成績データが与えられたとする。この相関係数はたかく、一方の成績が良ければ他方もよいという傾向が認められるであろう。単に算術平均をするのではなく、２つの成績の数値を１つの代表値にまとめることを考える。また英語の成績と数学のそれとをまとめる場合も同じ様に考えられる。いま n 個の２次元データ

(x_{1}, y_{1}), (x_{2}, y_{2}), ⋯, (x_{n}, y_{n})

が与えられたとする。これらに

z_{i} = ax_{i} + by_{i}, i = 1, 2, ⋯, n

という線形変換を施して、２次元データを１次元に縮約する。２次元データの情報がなるべく損なわれないよう縮約するには、どのような変換を行えばよいだろうか？まずこの線形変換が幾何的にどういう関係を表しているかを考えよう。原点と点 (a, b) を結んでできる直線を延長して新しい軸を設ける。与えられた２次元の点 (x, y) からこの軸に垂直におろした垂線の足を z とおけば、線形変換は２次元データの点 (x_{i}, y_{i}) を対応する新しい軸上の点 z_{i} に対応、詳しくいうと、平面の点をこの軸のうえに射影させている。

図

この線形変換によって２次元のデータとして散らばり、広がっていたものが、１次元の新しい軸上の点として並び換えられたことになる。したがって問題は、散布していたデータを、なるべくその特徴をよく代表するよう、適切な軸を選ぶことになる。データの特徴を代表するという考えかたを、分散の最大化、すなわち個体の識別効果を大きくすることと意味づけてみる。また係数 a, b の役割は直線の方向を決めている。

点 (a, b) を単位円周上の点、いいかえると、ベクトル [ \begin{matrix}a \\ b \end{matrix} ] の長さが 1 、 a^{2} + b^{2} = 1 を満たすとしてよい。

この立場に立てば、問題は［条件： a^{2} + b^{2} = 1 の下で、 z = ax + by の分散 σ_{z}^{2} = a^{2}σ_{x}^{2} + 2abσ_{xy} + b^{2}σ_{y}^{2} を最大にせよ。］という形に帰着された。さてこの条件つき最大化問題を解かなければならない。線形代数における対称行列の対角化を用いよう。ベクトルと行列をもちいて、目的関数を表現すると、分散と共分

散を要素とする対称行列 (これを分散共分散行列という) σ = [ \begin{matrix}σ_{x}^{2} & σ_{xy} \\ σ_{xy} & σ_{y}^{2} \end{matrix} ] とおいて

σ_{z}^{2} = a^{2}σ_{x}^{2} + 2abσ_{xy} + b^{2}σ_{y}^{2} = (a, b)σ [ \begin{matrix}a \\ b \end{matrix} ]

これはベクトル [ \begin{matrix}a \\ b \end{matrix} ] に関する２次形式と呼ばれる。係数行列は対称であるから、適当な直交行列 P を

用いて、

P^{−1}σP = [ \begin{matrix}λ_{1} & 0 \\ 0 & λ_{2} \end{matrix} ] = Λ

と対角化可能である。 (詳しくは線形代数のテキストを参考にされたい)。ここで λ_{1}, λ_{2} は行列 σ の固有値である。実際、係数行列 σ − λI の行列式が 0 という２次方程式：

λ^{2} − (σ_{x}^{2} + σ_{y}^{2})λ + σ_{x}^{2}σ_{y}^{2} − σ_{x}y^{2} = 0

から、２解 λ_{1}, λ_{2} が固有値であった。方程式はとくに固有方程式とよばれている。上の式における σ を σ = P ΛP^{−1} として目的関数の式に代入すると、 σ_{z}^{2} = λ_{1}ã^{2} + λ_{2}b̃^{2} に等しい。ただしベクトル

[ ã ; b̃ ] = P^{−1} [ \begin{matrix}a \\ b \end{matrix} ]. この固有値の性質をみると、固有方程式の判別式より２つの解は実数解で、 σ_{x}^{2} +σ_{y}^{2} ≥

0, σ_{x}^{2}σ_{y}^{2}−σ_{xy}^{2} ≥ 0 より２解とも非負でとなっている。したがって、条件を用いて λ_{M} = max(λ_{1}, λ_{2}), λ_{m} = min(λ_{1}, λ_{2}) とおくと, λ_{M} = λ_{M}ã^{2} + λ_{M}b̃^{2} ≥ λ_{1}ã^{2} + λ_{2}b̃^{2} ≥ λ_{m}ã^{2} + λ_{m}b̃^{2} = λ_{m} から、結局、 z の分散 σ_{z}^{2} は最大値が行列 σ の最大固有値となり、最小値は最小固有値となる。このようにして問題の最大値は分散共分散行列の最大固有値であることが分かる。また σ = P ΛP^{−1} よりトレースの値は変らない、 tr(Λ) = tr(P ΛP)^{−1} = tr(Λ) 。すなわち σ_{x}^{2} + σ_{y}^{2} = λ_{1} + λ_{2} という関係が得られる。この式から、もとのデータと線形変換されたデータの２つには、分散の合計が同じという意味で等しい情報をもつととらえることができる。

特に大切なポイントとしては、得られた固有値の大きさが問題である。いま対角化するための直交行列をつくるベクトルの順序を換えて、固有値が大きいものから小さいものへとならぶようにする。 λ_{1} ≥ λ_{2} このようにすると最大固有値に対する第一成分が分散の和 λ_{1} + λ_{2} のうち、最も大きな情報を含む。\frac{λ_{1}}{λ_{1} + λ_{2}} この値が大きければ大きいほど、第１成分だけで多次元のデータを縮約している、第１成分で２次元データを説明できていることになる。すべての固有値の総和の中で１つの固有値の占める割合を寄与率とよび、この値の最大値に対する成分が第１主成分、つぎのものが第２主成分とよばれる。最大固有値にたいする

固有ベクトルを [ \begin{matrix}a^{∗} \\ b^{∗} \end{matrix} ] とすると、線形変換 z = a^{∗}x + b^{∗}y が、変数 z の分散を最もおおきくなるような

変換であり、２次元データ (x, y) を１次元データ z に寄与率\frac{λ_{1}}{λ_{1} + λ_{2}} で縮小した第１主成分である。

変量が (x, y, z) の３個であっても、このデータを w = ax + by + cz と１次元に縮約することが考えられる。３次元空間に原点から延びる直線を引き、それを新しい軸として、空間の点をこの軸の上に射影して並べる。軸上に並んだ点の分散がなるべく大きくなるよう方向を決める。このときも、３次の対称行列は直交行列で対角化され、分散共分散行列に対する最大固有値がその最大分散で、対応する固有ベクトルがその方向を決定するのである。さらに大きさの順に固有値をならべ、第１主成分、第２主成分としていく。このように、与えられた高次元のデータをなるべく少ない個数の主成分で説明しようということが主成分分析である。

例題．

変量 x, y の平均、分散、共分散が\overline{x} = k, \overline{y} = l, σ_{x}^{2} = A, σ_{y}^{2} = B, σ_{xy} = C と与えられたとき、 z = ax + by の分散が大きくなるよう a, b を定めよ。解．

分散共分散行列は [ \begin{matrix}A & C \\ C & B \end{matrix} ] 。固有方程式は λ^{2} − (A + B)λ + AB − C^{2} = 0 であるから、解を λ_{1} ≥ λ_{2}

として a と b にかんする連立方程式

\{ \begin{matrix}(A − λ_{1})a + Cb & = 0 \\ Ca + (B − λ_{1})b & = 0 \end{matrix}

を求め、 a^{∗} = \frac{a}{\sqrt{a^{2} + b^{2}}} , b^{∗} = \frac{b}{\sqrt{a^{2} + b^{2}}} とおき、 z = a^{∗}x + b^{∗}y が分散最大となる代表値である。