この章では重回帰分析の原理と各種パラメータの意味、そして各種の変数選択手法について解説します。
重回帰分析(multiple regression analysis)は多変量解析の中で最も繁用される手法で、「ミスター多変量解析」とでも呼ぶべきものです。 名前から想像されるように、この手法は回帰分析を多変量に拡張した手法で、あるデータに直線的な影響を与えているデータが多種類ある時に、その影響のしかたを分析するための手法です。
単回帰分析と同様に、この手法も因果関係がはっきりしていて、しかも直線的な影響を与えているデータでなければ適用できません。 その意味から、重回帰分析を用いるのに最適な場合は、例えば各種の項目を総合して概括評価がされており、その評価規準すなわち各項目の重要度を分析したい場合などです。 なお、重回帰分析と区別するために普通の回帰分析を「単回帰分析(simple regression analysis)」と呼ぶことがあります。 (→5.1 相関係数と回帰直線)
では、例として前章の表6.1のデータに重回帰分析を適用してみましょう。
| 患者No. | TC | TG | 重症度 |
|---|---|---|---|
| 1 | 220 | 110 | 0 |
| 2 | 230 | 150 | 1 |
| 3 | 240 | 150 | 2 |
| 4 | 240 | 250 | 1 |
| 5 | 250 | 200 | 3 |
| 6 | 260 | 150 | 3 |
| 7 | 260 | 250 | 2 |
| 8 | 260 | 290 | 1 |
| 9 | 270 | 250 | 4 |
| 10 | 280 | 290 | 4 |
このデータの場合、目的変数が重症度で説明変数がTCとTGになり因果関係は明白です。 しかしこの表を眺めていただけでは、はたしてTCとTGが重症度に直線的な影響を及ぼしているかどうかはっきりしません。
そこで、例によってデータが目に見えるようにグラフ表示することにしましょう。 といっても、残念ながらウェブではホログラムのような3次元表示はまだできないので、表6.1のデータをそのままプロットするわけにはいきません。 そこで仕方がないので、3つの項目から2項目ずつ選んで3枚の散点図を描き、参考までに3次元プロットの見取り図を添えておくことにしましょう。 図7.1の左上が見取り図、右上がTCと重症度の散点図、左下が重症度とTGの散点図、右下がTCとTGの散点図で、TCと重症度の散点図と座標軸を合わせるために他の2枚の散点図は座標軸を反転しているので注意してください。
図7.1には説明変数と目的変数だけでなく説明変数同士の散点図も描いてあるのを御覧になって、不思議に思われる方もいることでしょう。 重回帰分析は説明変数同士の相関関係も計算に入れて分析しますので、その関係が直線的であるかどうかも確認する必要があるのです。 ただしその場合、対目的変数と違って因果関係は明確でなくてもかまいません。 図7.1を見る限り表6.1における3項目間の関係は直線的で、重回帰分析を適用しても問題なさそうです。
ところで賢明な読者は図7.1が3次元プロットの見取り図、正面図、側面図、および立面図になっていることにお気づきでしょう。 そしてもっと賢明な読者は、もし項目が4つ以上あったなら4次元以上のプロットが必要になり、その見取り図を描くことはもちろん、そんなプロットを想像することすら不可能になってしまうことにもお気づきでしょう。 そしてもっともっと賢明な読者は(←しつこい!(^^;))、この方法ではp個の項目についてp(p-1)/2枚の散点図を描く必要があり、もし10項目あったなら、何と45枚もの散点図を描かねばならない事実にもお気づきでしょう。
正直言って、これは恐ろしく面倒な作業です。 項目数が少ないうちはそれでもまだ何とかなるのですが、10項目以上ともなりますとほとんど手に負えません。 しかし正しい解析結果を得るためには、データを目で見て適切かどうかを判断するという、この重要な作業を省くことはできません。
そこで折衷案として、項目数が多い時には目的変数と各説明変数の散点図だけを描き、説明変数同士の散点図は省略することにしても良いでしょう。 これなら説明変数の数だけ散点図を描くだけですし、多変量解析というものは、1変量解析でそれなりの結果が出ている時にそれらを数学的に総合するためのものですから、重回帰分析を適用する前に目的変数と各説明変数の散点図を描き、基本統計量と相関係数ぐらいはちゃんと計算しておくべきなのです。 それすら面倒だから省略したいという御方には、残念ながら多変量解析をお勧めするわけにはいきません。
さて、重回帰分析では目的変数yと説明変数xの間に次のような直線的関係――数学的には「線形関係」といいます――があると仮定します。
偏回帰係数は回帰直線における傾き、すなわち回帰係数に相当するもので、回帰係数と違って「他の説明変数が一定」という条件の下で、ある説明変数だけが「1」変化した時に目的変数がいくつ変化するかを表す値です。 定数は回帰直線と同じく単なる「ゲタ」で、いくら高いゲタをはいても人の身長そのものは変わらないのと同様に、本質的な部分は偏回帰係数にあります。 なお偏回帰係数と区別するために、普通の回帰係数を「単回帰係数」と呼ぶこともあります。
単回帰分析と同様に重回帰分析も最小2乗法の原理を応用しており、誤差εを最小にするパラメーターbを計算し、それを未知パラメーターβの推定値にします。 その計算方法は単回帰分析以上に複雑怪奇ですので、計算法も具体的な解も省略し、その結果が次のような形の式で表されるとだけ述べておきましょう。
というわけで、「ベクトルと行列」を読んでいただいたことにして、重回帰分析をベクトルと行列を用いて記述してみましょう。
| y= | ┌ | | | | | | └ |
y1 : yi : yn |
┐ | | | | | | ┘ |
| X=[1n x1 … xj … xp]= | ┌ | | | | | | | └ |
1 : 1 : 1 |
x11 : xi1 : xn1 |
… … … |
x1j : xij : xnj |
… … … |
x1p : xip : xnp |
┐ | | | | | | | ┘ |
| β= | ┌ | | | | | | | └ |
β0 β1 : βj : βp |
┐ | | | | | | | ┘ |
| ε= | ┌ | | | | | | └ |
ε1 : εi : εn |
┐ | | | | | | ┘ |
この重回帰モデルに最小2乗法を適用してバラメーターβの最良不偏推定値(BLUE解)bを求めると、次のようになります。
この誤差の2乗和Qを最小にするbは、Qをbで偏微分したものを0p(p次元ゼロベクトル)と置いた連立方程式の解ですから、次のようにして解くことができます。
| ∂Q ――― ∂b |
= | ∂(y'y-2y'Xb+b'X'Xb) ―――――――――――――――― ∂b |
これは「正規方程式(normal equations)」と呼ばれる方程式で、行列[X'X]が正則、すなわち行列式|X'X|≠0ならば解を持ちます。 [X'X]は「単純積和行列」と呼ばれることもあり、具体的には次のような内容の(p+1)次対称行列です。
| [X'X]= | ┌ | | | | | | └ |
1n' : x1' : xp' |
┐ | | | | | | ┘ |
[1n x1 … xj … xp] |
| = | ┌ | | | | | └ |
n Σxi1 : Σxip |
Σxi1 Σxi12 : Σxi1xip |
… … … |
Σxip Σxi1xip : Σxip2 |
┐ | | | | | ┘ |
| [X'X]-1[X'X]b=b=[X'X]-1X'y= | ┌ | | | | | | | └ |
b0 b1 : bj : bp |
┐ | | | | | | | ┘ |
ここで、説明変数が1つの場合について具体的に計算してみましょう。
| [X'X]= | ┌ | | | | └ |
n Σxi1 |
Σxi1 Σxi12 |
┐ | | | | ┘ |
| X'y= | ┌ | | | | └ |
Σyi Σxi1yi |
┐ | | | | ┘ |
| ∴[X'X]-1= | ┌ | | | └ |
Σxi12 -Σxi1 |
-Σxi1 n |
┐ | | | ┘ |
1 ―――― n・S11 |
| = | ┌ | | | └ |
Σxi12/n -mx1 |
-mx1 1 |
┐ | | | ┘ |
1 ―― S11 |
| = | ┌ | | | └ |
Σ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 |
┐ | | | ┘ |
偏回帰係数ベクトルbの解を元の重回帰モデルに代入すると次のようになります。
このことから、重回帰モデルをベクトル空間上で幾何学的に解釈すると次のようになります。 n次元ベクトル空間において、1nベクトルとp個の説明変数ベクトルx1〜xpを基底とする(p+1)次元部分空間が考えられ、これはn次元空間中の「(p+1)次元超平面(hyperplan)」を構成します。 この超平面に目的変数ベクトルyを正射影したものが推定値ベクトルy~になり、その超平面の直交補空間すなわち{n-(p+1)}次元部分空間にyを正射影したものが誤差(残差)ベクトルεになります。
また1nベクトルとp個の説明変数ベクトルx1〜xpがn次元ベクトル空間における(p+1)次元部分空間の基底になるためには、nが少なくとも(p+1)以上である必要があります。 このため例数が(項目数+1)よりも少ない時は、重回帰分析の解を求めることができません。
説明変数が1つもなくて、行列Xの成分が1nベクトルだけの場合は図7.3のようになります。 この図7.3は前章第2節(注2)の図6.3と本質的に同じもので、y~がyの平均値を成分とするベクトルmyに、εが偏差ベクトルdyになることがわかります。 このことは、Xの成分が1nベクトルだけの正規方程式を解いてみるとよりはっきりします。
| Z=X[X'X]-1X'= | 1 ―― n |
[1n][1n]' |
| = | 1 ―― n |
┌ | | | └ |
1 : 1 |
… … |
1 : 1 |
┐ | | | ┘ |
また前章第2節(注2)で説明したように、平均値を求めるということはデータベクトルyを平均値ベクトルmyと偏差ベクトルdyとに直交分解していることに相当します。 それと同様に、重回帰分析もデータベクトルyを推定値ベクトルy~と誤差ベクトルεとに直交分解していることに相当します。