第 4 章 多変量解析序説
4.3 重回帰分析
ある変量x1, . . . , xp, yのデータ
(x11, . . . , xp1, y1), . . . ,(x1n, . . . , xpn, yn)
に対し、このデータが表すyのx1, . . . , xpへの従属関係を１次式y=b0+b1x1+
· · ·+bpxpで近似的に表したいとする。このときx1, . . . , xpを説明変量、yを目 的変量という。１次式f(x1, . . . , xp) =b0+b1x1· · ·+bpxpは以下の方法（最小 ２乗法）で決定することができる。
1. 目的変数と説明変数を選択する。
2. 数式のパラメータb0, . . . , bpを決めるために以下の基準を設ける。
S =
∑n i=1
(yi−f(x1i, . . . , xpi))2 を最小にする。Sを偏差２乗和という。
3. 上の基準でパラメータb0, . . . , bpを求めるには、
∂S
∂bi
= 0 (i= 0, . . . , p) の連立方程式を解く。
ここで注意すべきは、一般に最小２乗法において上の連立方程式の解が存在す るとしても、偏微分が0という条件は極値の必要条件に過ぎなく、(b0, . . . , bp)が 上の連立方程式の解であるというだけではこの解におけるSの最小性は保証さ れないことである。しかし以下に示すようにここでのy=b0+b1x1+· · ·+bpxp のあてはめに関しては上の方法からSを最小にする解を構成することが可能で ある。
4.3.1 正規方程式
前節の最小２乗法を実際に遂行してみる。偏差２乗和Sをb0, . . . , bpの関数と 考えれば、求める(b0, . . . , bp)においてはSは極小値をとるから、求めるbkの条 件は
∂S
∂bk =
∑n i=1
2(yi−b0−b1x1i− · · · −bpxpi)(−xki) = 0, k= 0, . . . , p となる。ただし、x0i= 1とする。これを整理すると、
( n
∑
i=1
xki
) b0+
( n
∑
i=1
xkix1i
)
b1+· · ·+ ( n
∑
i=1
xkixpi
) bp=
∑n i=1
xkiyi, k= 0, . . . , p このb0, . . . , bpの連立１次方程式を、正規方程式という。求めるbkはこの方程式 を解くことによって得られる。p= 1, p= 2の場合の正規方程式を書き下ろして みると、
p= 1
nb0 + (∑n
i=1
xi )
b1 =
∑n i=1
yi (∑n
i=1
xi
) b0 +
(∑n
i=1
xi2
) b1 =
∑n i=1
yixi
ns2x=∑
x2i −n¯x2, nsxy=∑
xiyi−n¯x¯y であることに注意すれば、
b0 =
∑x2i∑ yi−∑
xi
∑xiyi
n∑
x2i −(∑ xi)2
= n¯y∑
x2i −n¯x∑
xiyi+n2x¯2y¯−n2x¯2y¯ n2s2x
= rn2ys¯ 2x−n2xs¯ xy
n2s2x
= y¯−x¯sxy s2x
b1 =
∑x2i∑ yi−∑
xi
∑xiyi
n∑
x2i −(∑ xi)2
=
∑(xi−x)(y¯ i−y)¯
∑(xi−x)¯ 2
= rsy sx
p= 2
nb0 +
(∑n
i=1
x1i
) b1 +
(∑n
i=1
x2i
) b2 =
∑n i=1
yi
(∑n
i=1
x1i
) b0 +
(∑n
i=1
x21i )
b1 + (∑n
i=1
x1ix2i
) b2 =
∑n i=1
x1iyi
(∑n
i=1
x2i )
b0 + (∑n
i=1
x2ix1i )
b1 +
(∑n
i=1
x2i2 )
b2 =
∑n i=1
x2iyi こうして求めた一次式y=f(x1, . . . , xp) =b0+b1x1+· · ·+bpxpを重回帰式 という。また、bkを回帰係数という。
4.3.2 幾何学的考察による正規方程式の解の存在
正規方程式の解は一意でないかもしれないが常に存在する。以下にこれを示 す。データ行列Xを
X =
1 x11 · · · xp1 1 x12 · · · xp2
... ... ... 1 x1n · · · xpn
, y=
y1 y2
... yn
, b=
b0 b1
... bp
とおくと、一般次の正規方程式の係数行列は、
X′X =
n ∑
x1i
∑x2i · · · ∑ xpi
∑x1i ∑
x1i2 ∑
x1ix2i · · · ∑ x1ixpi
∑x2i
∑x2ix1i
∑x2i2 ...
... ... . ..
∑xpi ∑
xpix1i · · · ∑ xpi2
で、正規方程式はX′Xb=X′yと書ける。ここに′は行列の転置を表す。
Xの列が張る線形空間をV ={Xb : b ∈Rp+1}とおく。X に列に関する行 列の基本操作を施すことで、V の基底を得ることができ、さらにこれをグラム シュミットの直交化法で正規直交基底にすることができる。よってu0, . . . , urを V の正規直交基底とする(V の正規直交基底は実対称行列X′Xの固有ベクトル からも構成できるが省略する-「多変量解析の徹底研究」現代数学社を参照)。
目的変量ベクトルyのV への射影とは、y−uがすべてのv∈V と直交する ようなu∈V のこととする。このようなuはu=
∑r i=0
(y, ui)uiとして構成でき る。一般に
|v+w|2 = (v+w, v+w)
= |v|2+ 2(v, w) +|w|2 であるから、v∈V に対しv−uとy−uが直交することより、
|y−v|2=|v−u|2+|y−u|2
ここでv=Xb, u=Xˆbと書けば、|y−v|2=S(b), |y−u|2=S(ˆb)である。こ のことより、ˆbがSの最小値を与えることがわかる。
また、xj, j= 0, . . . , pをXの列ベクトルとすると、
(y−Xˆb, v) = 0, ∀v∈V
⇐⇒ (y−Xˆb, xj), ∀j
⇐⇒ X′(y−Xˆb) = 0
であるから、XˆbがyのV への射影であることはˆbが正規方程式を満たすことと 同値である。
よって、正規方程式には解が存在しSの最小値を与えることがわかった。
X′Xの一般化逆行列(X′X)−とはX′X(X′X)−X′X =X′X なるもののこ とである。これは実対称行列X′Xの固有値から構成できる。正規方程式の解は (X′X)−X′yと表すことができる(「多変量解析の徹底研究」現代数学社参照)。
また正規方程式の解が、Sの最小値を与えることは次のように代数的に示すこ とができる。
解をˆb=
(ˆb0,· · ·,ˆbp
)
として、b= (b0,· · · , bp)を任意の組とすると、
S(b) = (y−Xb)′(y−Xb)
= {
y−Xˆb+X(ˆb−b) }′{
y−Xˆb+X(ˆb−b) }
(4.1)
= (y−Xˆb)′(y−Xˆb) + (Xb−Xˆb)′(Xb−Xˆb) (4.2)
≥ S(ˆb)
ここで（１）式から（２）式の導出においては、(Xˆb−Xb)′(y−Xˆb)、(y− Xˆb)′(Xˆb−Xb)の交差積項はともに等しく、X′(y−Xˆb) =0なる関係を使え ば０となることを使う。
4.3.3 重回帰式に関する推測
最小２乗解の推定値としての性質
最小２乗解に関する統計的推測をするために問題を以下のようにとらえる(仮 定を設ける)。
1. 目的変量yと説明変量x1, . . . , xpの間には
y=β0+β1x1+· · ·+βpxp
の関係を想定したいが、実際にはxkたちの値によってyの値は確定しない。
2. そこで目的変量yと説明変量x1, . . . , xpの間には、説明変量の値の組(x1i, . . . , xpi) を与える実体2ごとに
yi=β0+β1x1i+· · ·+βpxpi+εi
の関係があるとする。εiはi番目の実体についての1の関係からの’ずれ’(誤 差)である。
2重回帰分析の実際などあとに出てくる例を見れば、ここでいう実体が何に当たるかは明白であろ う。
3. εiは各々独立に正規分布N(0, σ2)に従うものとする(σは既知とはかぎら ない)。
この仮定のもとでは、各yiはεiを通して確率的に決まるので、それらから計 算される重回帰式
y=b0+b1x1+· · ·+bpxp
における回帰係数bkは確率変数となる。
bkをβkの最小２乗推定値という。
E(bk) =βk
が成立する。すなわち最小２乗推定値bkはβkの不偏推定値である。
同じ一次式でも変量の一次式でなくパラメータb0, . . . , bpの一次式を考える：
l0b0+· · ·+lpbp
ここでl0, . . . , lpは1, x, . . . , xpでも1, x1, . . . , xpでもよい。こういったモデル(現 象を表す数式)を線形モデルという。
正規方程式が解を持ち解がSの最小値を与えることはl0, . . . , lpが1, x, . . . , xp の場合も重回帰とほぼ同様に証明できる。
線形モデルのうち正規方程式の解、最小２乗解には次の統計的性質がある。
定理 4.3.1 (Gauss-Markovの定理の簡単な場合). l1b1+· · ·+lpbpの平均値が
¯
yに等しいもののうちで、分散が最小になるのは最小２乗解l1ˆb1+· · ·+lpˆbpで ある。
重回帰式の有意性の検証(1) - 分散分析による検定
ここで扱うF分布による方法は、回帰係数のうち何個かが0であるという帰 無仮説を検定するものである。したがって次節のt検定による単一の回帰係数が 0であるという帰無仮説を検定する方法を代替できるものである。
ここでは、単純に回帰係数全部が0という帰無仮説：
H0 : β1=· · ·=βp= 0
の検定を考える。もしこの仮定が真なら説明変量はすべて目的変量に影響を与え ないことになる。すなわち、
重回帰式が(全体として)有意であるかどうか を検定することになる。
Se=
∑n i=1
(yi−yˆi)2, SR=
∑n i=1
(ˆyi−y)¯ˆ2 とおく。全変動を
Syy =
∑n i=1
(yi−y)¯ 2 とおくと、
Se=Syy−SR
となっている。
一般に(帰無仮説H0がなくとも)、Se/σ2はN(0,1)に従う互いに独立なνe個 の確率変数の2乗和として表すことができることが知られている(「自然科学の 統計学」p.56)。ここにνe =n−rank(X)(Xはn×(p+ 1)の説明変量のデー タと１からなるいわゆる計画行列)であるが、ここでは計画行列の正則性を仮定 して、
νe=n−p−1 としておく。
よってχ2に関する前述の定理よりSe/σ2は自由度n−p−1のχ2分布に従う ことがわかる。
一方、SR/σ2は帰無仮説H0の下でN(0,1)に従う互いに独立なp個の確率変 数の2乗和として表すことができることが知られている(「自然科学の統計学」
p.56、「多変量解析の徹底研究」など)。
よってχ2分布に関する前述の定理よりSR/σ2は自由度pのχ2分布に従うこ とがわかる。
したがって、F分布に関する前述の定理より、帰無仮説H0の下では、
F =
SR p Se n−p−1
は、自由度(p, n−p−1)のF分布に従う。
これにより、帰無仮説の検定ができる。Fは推定による変動が推定誤差に相対 的にどのぐらい大きいかの指標であるから大きい程推定式は有効であることにな
る。いいかえると、全変動が推定による変動SRで多く説明できる程推定式は有 効であるから、
F = 1/
(Syy/(n−p−1)
SR/p − p
n−p−1 )
が大きい程推定式は有効である。よってF分布のα点をfνA,νe(α)とすると、
F > fνA,νe(α) のときH0を棄却する。
実は、Ve=Se/(n−p−1)は一般に(帰無仮説H0がなくとも)σ2の不偏推定 値であり、VR = SR/pも帰無仮説の下でならσ2の不偏推定値であるので、上 のFは不偏分散比であるとみなすことができる。そこでこの検定を分散分析と いう。
重回帰式の有意性の検証(2) - 重相関係数による方法
次の式は、回帰式による推定誤差（２項目の分子）がデータの全変動（２項目 の分母）に占める割合いの低さを表すものである。
R2= 1−
∑n i=1
(yi−b0−b1x1− · · · −bpxp)2
∑n i=1
(yi−y)¯ 2
すでに述べたように、これはR２乗値（重決定係数）と呼ばれ、その平方根 は重相関係数と呼ばれる。重相関係数は、目的変量y のデータ値yi と推定値 b0+b1x1i+. . .+bpxpiの相関係数に等しい。
さらにつぎがいえる。
定理4.3.2. 目的変量yと最大の相関係数を持つ線形モデルb0+b1l1+· · ·+bplp は重回帰式であり、そのときの最大相関係数は重相関係数である。(「多変量解 析論」北川敏男著、共立出版、p.16参照)
分散分析によるF検定を行って有意であるという結論を得ても、これは「重 回帰式は何らの役にも立たない」という帰無仮説が棄却されたということである から、「何らの役にも立たない」とはいえない、という意味にすぎなく、積極的 に役に立つということではない。普通F値が相当に大きくなければ実際に用い るには有効でないことが多い。
重相関係数は、分散分析による有意の結論を補助するために使うことができ る。すなわち、重相関係数が十分大きければ、回帰式は有効、そうでなければあ まり有効ではないと判断するのである。
回帰係数の区間推定と有意性の検定 定理 4.3.3.
tk= bk−βk
√akkVe, k= 1, . . . , p
t0= vuut b0−β0
(
1 n +
∑p j=1
∑p k=1
¯ xjx¯kajk
) Ve
なる統計量は、自由度n−p−1のt分布に従う。ここに、Veは残差の不偏分散、
すなわち
Ve=
∑n i=1
(yi−f(x1, . . . , xp))2 n−p−1 であり、またajkは、偏差積和行列
A= (ajk) = ( n
∑
i=1
(xji−x¯j)(xki−¯xk) )
の逆行列の(j, k)成分である。
(証明の概略) (「自然科学の統計学」p.60、「多変量解析の徹底研究」など)
これを利用して、回帰係数の区間推定と有意性(βk= 0)の検定ができる。まず、
P(tk> tn−p−1(α)) =P(tk <−tn−p−1(α)) = α 2 となるtn−p−1(α)を求めておく。
すると、信頼係数1−αのβkの信頼区間は [
bk−√
akkVe, bk+√ akkVe
]
, k= 1, . . . , p
b0− vu uu t
1 n+
∑p j=1
∑p k=1
¯ xjx¯kajk
Ve, b0+ vu uu t
1 n+
∑p j=1
∑p k=1
¯ xjx¯kajk
Ve
, k= 0
である。
帰無仮説H0:βk = 0の検定を考える。この仮説は説明変量xkが目的変量の 変化に全く影響しないことを意味するので、これはxkまたはbkの有意性を検定 することになる。
この仮説のものでは、k= 1, . . . , p、k= 0に応じて、
tk= bk
√akkVe
, k= 1, . . . , p
t0= b0
vu ut (
1 n +
∑p j=1
∑p k=1
¯ xjx¯kajk
) Ve
であり、bkが0から遠い程絶対値の大きな値となる。よって、
|tk|> tn−p−1(α) であれば、危険率αでH0を棄却できる。