玄関雑学の部屋雑学コーナー統計学入門

第7章 重回帰分析

この章では重回帰分析の原理と各種パラメータの意味、そして各種の変数選択手法について解説します。

7.1 重回帰モデル

(1) 単回帰分析と重回帰分析

重回帰分析(multiple regression analysis)は多変量解析の中で最も繁用される手法で、「ミスター多変量解析」とでも呼ぶべきものです。 名前から想像されるように、この手法は回帰分析を多変量に拡張した手法で、あるデータに直線的な影響を与えているデータが多種類ある時に、その影響のしかたを分析するための手法です。

単回帰分析と同様に、この手法も因果関係がはっきりしていて、しかも直線的な影響を与えているデータでなければ適用できません。 その意味から、重回帰分析を用いるのに最適な場合は、例えば各種の項目を総合して概括評価がされており、その評価規準すなわち各項目の重要度を分析したい場合などです。 なお、重回帰分析と区別するために普通の回帰分析を「単回帰分析(simple regression analysis)」と呼ぶことがあります。 (→5.1 相関係数と回帰直線)

では、例として前章の表6.1のデータに重回帰分析を適用してみましょう。

表6.1 高脂血症患者のTCとTG
患者No.TCTG重症度
12201100
22301501
32401502
42402501
52502003
62601503
72602502
82602901
92702504
102802904

このデータの場合、目的変数が重症度で説明変数がTCとTGになり因果関係は明白です。 しかしこの表を眺めていただけでは、はたしてTCとTGが重症度に直線的な影響を及ぼしているかどうかはっきりしません。

そこで、例によってデータが目に見えるようにグラフ表示することにしましょう。 といっても、残念ながらウェブではホログラムのような3次元表示はまだできないので、表6.1のデータをそのままプロットするわけにはいきません。 そこで仕方がないので、3つの項目から2項目ずつ選んで3枚の散点図を描き、参考までに3次元プロットの見取り図を添えておくことにしましょう。 図7.1の左上が見取り図、右上がTCと重症度の散点図、左下が重症度とTGの散点図、右下がTCとTGの散点図で、TCと重症度の散点図と座標軸を合わせるために他の2枚の散点図は座標軸を反転しているので注意してください。

図7.1 見取り図と3種類の散点図

図7.1には説明変数と目的変数だけでなく説明変数同士の散点図も描いてあるのを御覧になって、不思議に思われる方もいることでしょう。 重回帰分析は説明変数同士の相関関係も計算に入れて分析しますので、その関係が直線的であるかどうかも確認する必要があるのです。 ただしその場合、対目的変数と違って因果関係は明確でなくてもかまいません。 図7.1を見る限り表6.1における3項目間の関係は直線的で、重回帰分析を適用しても問題なさそうです。

ところで賢明な読者は図7.1が3次元プロットの見取り図、正面図、側面図、および立面図になっていることにお気づきでしょう。 そしてもっと賢明な読者は、もし項目が4つ以上あったなら4次元以上のプロットが必要になり、その見取り図を描くことはもちろん、そんなプロットを想像することすら不可能になってしまうことにもお気づきでしょう。 そしてもっともっと賢明な読者は(←しつこい!(^^;))、この方法ではp個の項目についてp(p-1)/2枚の散点図を描く必要があり、もし10項目あったなら、何と45枚もの散点図を描かねばならない事実にもお気づきでしょう。

正直言って、これは恐ろしく面倒な作業です。 項目数が少ないうちはそれでもまだ何とかなるのですが、10項目以上ともなりますとほとんど手に負えません。 しかし正しい解析結果を得るためには、データを目で見て適切かどうかを判断するという、この重要な作業を省くことはできません。

そこで折衷案として、項目数が多い時には目的変数と各説明変数の散点図だけを描き、説明変数同士の散点図は省略することにしても良いでしょう。 これなら説明変数の数だけ散点図を描くだけですし、多変量解析というものは、1変量解析でそれなりの結果が出ている時にそれらを数学的に総合するためのものですから、重回帰分析を適用する前に目的変数と各説明変数の散点図を描き、基本統計量と相関係数ぐらいはちゃんと計算しておくべきなのです。 それすら面倒だから省略したいという御方には、残念ながら多変量解析をお勧めするわけにはいきません。

(2) 重回帰モデル

さて、重回帰分析では目的変数yと説明変数xの間に次のような直線的関係――数学的には「線形関係」といいます――があると仮定します。

y=β01x1+…+βjxj+…+βpxp+ε  (j=1,…,p)
β0:定数(y切片)  βj:偏回帰係数  ε:誤差

これを「(線形)重回帰モデル」といい、β0を「定数(constant、y切片)」、βjを「偏回帰係数(partial regression coefficient)」と呼びます。 例えば、表6.1のデータに当てはめる重回帰モデルは次のようになります。

y(重症度)=β01x1(TC)+β2x2(TG)+ε  (p=2)

偏回帰係数は回帰直線における傾き、すなわち回帰係数に相当するもので、回帰係数と違って「他の説明変数が一定」という条件の下で、ある説明変数だけが「1」変化した時に目的変数がいくつ変化するかを表す値です。 定数は回帰直線と同じく単なる「ゲタ」で、いくら高いゲタをはいても人の身長そのものは変わらないのと同様に、本質的な部分は偏回帰係数にあります。 なお偏回帰係数と区別するために、普通の回帰係数を「単回帰係数」と呼ぶこともあります。

単回帰分析と同様に重回帰分析も最小2乗法の原理を応用しており、誤差εを最小にするパラメーターbを計算し、それを未知パラメーターβの推定値にします。 その計算方法は単回帰分析以上に複雑怪奇ですので、計算法も具体的な解も省略し、その結果が次のような形の式で表されるとだけ述べておきましょう。

y~=b0+b1x1+…+bjxj+…+bpxp  (j=1,…,p)
y~:yの推定値  推定誤差:ε=y−y~  bj:偏回帰係数(推定値)

この式を「重回帰式」といい、直線回帰式を多変量に拡張したものに相当します。 (注1)


(注1) またまた頭の痛い話で恐縮ですが、重回帰分析の計算方法を理解するためにはベクトルと行列の知識の中で、特に逆行列と微分の知識が必要になります。 それらについての解説は「ベクトルと行列」の第3章、第7章、第8章あたりをご覧ください。

というわけで、「ベクトルと行列」を読んでいただいたことにして、重回帰分析をベクトルと行列を用いて記述してみましょう。

・重回帰モデル
yi=β01xi1+…+βjxij+…+βpxipi
=β0n11+…+βjj+…+βpp+εXβ+ε~+ε
~=Xβ







y1

yi

yn







=[1n 1jp]=







1


1


1
x11

xi1

xn1






x1j

xij

xnj






x1p

xip

xnp








β







β0
β1

βj

βp








ε






ε1

εi

εn







この重回帰モデルに最小2乗法を適用してバラメーターβの最良不偏推定値(BLUE解)を求めると、次のようになります。

~=Xb
ε~=Xb
Q=Σεi2ε'ε=[Xb]'[Xb]
 =''[Xb]−[Xb]'+''Xb
 ='−2'Xb+''Xb → 最小

この誤差の2乗和Qを最小にするは、Qをで偏微分したものをp(p次元ゼロベクトル)と置いた連立方程式の解ですから、次のようにして解くことができます。

∂Q
―――
∂('-2'Xb+''Xb)
――――――――――――――――
  =-2'+2[']'p
∴[']'

これは「正規方程式(normal equations)」と呼ばれる方程式で、行列[']が正則、すなわち行列式|'|≠0ならば解を持ちます。 [']は「単純積和行列」と呼ばれることもあり、具体的には次のような内容の(p+1)次対称行列です。

[']=






1n'

1'

p'







[1n 1jp]
  =





n

Σxi1

Σxip
Σxi1
Σxi12

Σxi1xip




Σxip
Σxi1xip

Σxip2






この行列が正則ならばその逆行列を正規方程式の両辺に左から掛けて、解は次のようになります。

[']-1[']=[']-1'







b0
b1

bj

bp








b0:定数項  bj:偏回帰係数(推定値、j=1,…,p)

ここで、説明変数が1つの場合について具体的に計算してみましょう。

[']=




n


Σxi1
Σxi1

Σxi12





'




Σyi

Σxi1yi





余因子行列を用いて単純積和行列の逆行列を求めますと、次のようになります。

|'|=nΣxi12-(Σxi1)2=n・S11
X11=Σxi12  X12=-Σxi1
X21=-Σxi1  X22=n
∴[']-1



Σxi12

-Σxi1
-Σxi1

n




1
――――
n・S11
  =



Σxi12/n

-mx1
 -mx1

1




1
――
S11
=[']-1'
 =



Σxi12/n

-mx1
 -mx1

1








Σyi

Σxi1yi




1
――
S11
 =



myΣxi12-mx1Σxi1yi

-mx1Σyi+Σxi1yi




1
――
S11
 =



my(Σxi12-n・mx12)-mx1(Σxi1yi-n・mx1my)

Σxi1yi-n・mx1my




1
――
S11
 =



myS11-mx1S1y

S1y




1
――
S11
 =



my-mx1(S1y/S11)

S1y/S11








b0

b1




以上のように、説明変数が1つの場合は単回帰分析の結果と一致します。 (→5.1 相関係数と回帰直線)

偏回帰係数ベクトルの解を元の重回帰モデルに代入すると次のようになります。

~=Xb[']-1'
[']-1'

は次のような性質を持ち、n次元ベクトル空間上のを(p+1)次元部分空間に正射影して、~に変換する射影子(projection)になります。 (→ベクトルと行列 6.ベクトルの直交分解と直交変換 (注1))

'=([']-1')'=([']-1')''=[']-1' (n次元対称行列)
2=([']-1')([']-1')
  =[']-1['][']-1'=[']-1' (ベキ等行列、行列のベキが元の行列に等しい行列)

このことから、重回帰モデルをベクトル空間上で幾何学的に解釈すると次のようになります。 n次元ベクトル空間において、nベクトルとp個の説明変数ベクトル1pを基底とする(p+1)次元部分空間が考えられ、これはn次元空間中の「(p+1)次元超平面(hyperplan)」を構成します。 この超平面に目的変数ベクトルを正射影したものが推定値ベクトル~になり、その超平面の直交補空間すなわち{n-(p+1)}次元部分空間にを正射影したものが誤差(残差)ベクトルεになります。

図7.2 重回帰モデルの幾何学的解釈 図7.3 説明変数がない時の重回帰モデル

なぜならば~はn1pの1次結合として、

~=β0n11+…+βpp

と表されるため、これらを基底とする超平面に含まれることになります。 そして~は、を射影子によってその超平面に正射影したものになります。 ~が正射影であることは、εの大きさを最小にする射影、すなわち、

ε2=Σεi2 → 最小

という最小2乗条件を満足する射影は必然的に正射影になり、はその条件を満足する解から作った射影子であることからわかります。 そしてこのことは、図7.2から直観的に理解できると思います。

またnベクトルとp個の説明変数ベクトル1pがn次元ベクトル空間における(p+1)次元部分空間の基底になるためには、nが少なくとも(p+1)以上である必要があります。 このため例数が(項目数+1)よりも少ない時は、重回帰分析の解を求めることができません。

説明変数が1つもなくて、行列の成分がnベクトルだけの場合は図7.3のようになります。 この図7.3は前章第2節(注2)の図6.3と本質的に同じもので、~がの平均値を成分とするベクトルyに、εが偏差ベクトルyになることがわかります。 このことは、の成分がnベクトルだけの正規方程式を解いてみるとよりはっきりします。

[']=[n]  [']-1=[1/n]  '=[Σyi]
=[']-1'=[Σyi/n]=[my]=[b0]
[']-1'= 1
――
n
[1n][1n]'
  = 1
――
n




1

1


1

1




また前章第2節(注2)で説明したように、平均値を求めるということはデータベクトルを平均値ベクトルyと偏差ベクトルyとに直交分解していることに相当します。 それと同様に、重回帰分析もデータベクトルを推定値ベクトル~と誤差ベクトルεとに直交分解していることに相当します。

~+ε
~'ε~'(-~)=~'-~'~
  =()'-()'()=''-''
  ='-'=0
~⊥ε (~とεは直交する)

したがってデータを最小2乗法を用いて単純に要約したものが平均値mであり、1つの説明変数との関係を考慮して要約したものが回帰直線、p個の説明変数との関係を考慮して要約したものが重回帰式である、ということになります。 このことから統計学とはデータを要約することであり、多変量解析がその特別 な場合として1変量解析を包含しているということが理解できると思います。