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

9.4 多変量の場合

(1) 多変量正規分布とマハラノビスの汎距離

次に変数が2つ以上の時の群の判別について考えてみましょう。 変数が1つの時はデータが正規分布すると仮定して尤度を求めました。 それと同様に、変数が2つ以上の時もデータが正規分布すると仮定して尤度を求めることができます。 ただしその場合の正規分布は普通のものではなく、変数が2つ以上あり、しかもその変数間に相関関係があると仮定した拡張正規分布であり、多変量正規分布(multivariate normal distribution)と呼ばれています。 多変量正規分布の式は恐ろしく複雑で、見たとたんに頭が痛くなるので(注1)を見ていただくとして、例えば変数が2つの時の姿は図9.4.1のような感じになります。 (注1)

図9.8 2次元正規分布

変数が1つの時にデータのバラツキ具合を表す指標は偏差(xi - m)でした。 偏差はバラツキの指標であると同時に、見方によっては、あるデータと平均値の距離つまりあるデータが分布の中心からどの程度離れているかを表す指標と解釈することもできます。 しかし偏差が例えば「10」だとしても、データのバラツキが大きくて分布の幅が広い時の「10」と、データのバラツキが小さくて分布の幅が狭い時の「10」では平均値との距離感は違います。

そこで偏差を標準偏差で割って分布の幅に左右されない値にします。 これは元のデータを標準化(規準化、standardization)したことに相当し、平均値が0、標準偏差が1になります。 この標準化した距離は分布の幅に左右されずに距離を表すことができる汎用距離なので、色々な場面で利用されます。 例えば検定の効果量(effect size)はこの標準化距離を利用した指標です。 また距離の正負をなくすために標準化距離を平方した値もよく用いられ、それを平方距離といいます。 (→1.3 データの要約方法 (注1)1.6 統計的仮説検定の考え方 (6)検定の特性式と必要例数の計算式)

偏差:di = xi - m
標準化距離:   平方距離:

この標準化距離を多変量に拡張したものをマハラノビスの汎距離(Mahalanobis's generalized distance)といいDで表します。 そしてそれを平方した値をマハラノビスの平方距離(Mahalanobis's squared distance)といいD2で表します。 図9.4.1から何となくわかるように、一般に多変量正規分布は方向によって分布の幅が違います。 この分布を山に見立てた時、傾斜が急なところを頂上から一定距離だけ山を下った地点の頂上からの水平距離と、傾斜が緩やかなところを頂上から一定距離だけ山を下った地点の頂上からの水平距離を比べると、後者の方が相対的に長くなります。

これを考慮した距離がマハラノビスの汎距離であり、マハラノビスの汎距離が同じ点を結ぶとちょうど地図の等高線のようになります。 この分布の山の表面を空間と考えると、これは曲がった空間になり、マハラノビスの汎距離は曲がった空間における距離になります。 そして曲がった空間は非ユークリッド幾何学におけるリーマン空間(Rieman space)に相当し、マハラノビスの汎距離はリーマン空間における距離つまり線素に相当します。 線素が非ユークリッド幾何学で重要な役目を果たすのと同様に、マハラノビスの汎距離は判別分析で重要な役目を果たします。 (→当館の「ベクトルと行列 第13章 リーマン空間」参照)

図9.4.2 マハラノビスの汎距離(r=0.5)

元のデータがp次元正規分布をする時、D2は自由度pのχ2分布をするという性質があります。 この性質を利用すると、図9.4.2の等高線状の楕円内部に含まれるデータの割合を計算することができます。 例えばD = 1(D2 = 1)の楕円内には約39%、D = 2(D2 = 4)の楕円内には約86%、D = 3(D2 = 9)の楕円内には約99%のデータが含まれます。

逆に95%のデータが含まれる楕円は、D2 = χ2(2,0.05) = 5.991より、分布の中心からの汎距離がD = 2.448になる点を結んだ楕円になります。 この楕円を95%信頼楕円または95%等確率偏差楕円といい、お互いに相関関係がある2つのデータの2次元的な分布状態を表すのに利用されます。 図9.1.1の正常群と動脈硬化症群の楕円は、これを利用して描いたものです。 (注2)

(2) マハラノビスの汎距離を利用した群の判別

D2が自由度pのχ2分布をするという性質を利用すると、p次元正規分布をする特定の母集団からp種類の項目値を持つ個体が得られる確率、つまりその母集団の尤度を計算することができます。 そして尤度が計算できれば、前節で説明したように最尤法を利用して群の判別を行うことができます。 しかしマハラノビスの汎距離と尤度の間には反比例的な関係があるので、いちいち尤度を計算しなくても母集団の中心からの汎距離を計算して一番近い群に判別するという方法で良いことになります。

例えば表9.1.1のTCとTGのデータを正常群と動脈硬化症群の母集団からサンプリングした標本集団のデータと考え、それぞれの母集団のTCとTGの母数を次のように推定します。

○正常群:TCの母平均推定値 = 標本平均値 = 207  母標準偏差推定値 = 不偏標準偏差 = 18
TGの母平均推定値 = 標本平均値 = 206  母標準偏差推定値 = 不偏標準偏差 = 59
TCとTGの母相関係数推定値 = 標本相関係数 = 0.79
○動脈硬化症群:TCの母平均推定値 = 標本平均値 = 251  母標準偏差推定値 = 不偏標準偏差 = 19
TGの母平均推定値 = 標本平均値 = 209  母標準偏差推定値 = 不偏標準偏差 = 65
TCとTGの母相関係数推定値 = 標本相関係数 = 0.75

これらの母数推定値とデータが2次元正規分布するという仮定から、あるTCとTGの値を持つ点について、それぞれの母集団の中心からの汎距離を計算することができます。 これにより正常か動脈硬化か不明な被検者についてTCとTGを測定し、その値に対する2つの母集団からの汎距離を比較することによって、どちらの群に属するか判別することが可能になります。 これは2つの母集団からの汎距離の差を計算し、その正負によって群を判別しても同じことであり、その方がより直観的になります。

例えばあるTCとTGの値について正常群の中心からの平方距離がD22になり、動脈硬化症群の中心からの平方距離がD12になったとします。 そしてこれらの差をzとすると、次のようにして群を判別することができます。 汎距離ではなく平方距離を用いるのは、距離の正負をなくして話を単純にするためです。

z = D22 - D12 → z > 0の時、動脈硬化症群と判別  z < 0の時、正常群と判別

さらに話を単純にするために2つの母集団の標準偏差と相関係数が等しいと仮定すると、平方距離の計算が簡単になる上、どちらの母集団から見ても、もう一方の母集団の中心までの平方距離が同じになります。 そしてその母集団間の平方距離の中間点ではzが0になり、2群を判別する境界値になります。 これは2次元正規分布するTCとTGのデータをzという1つのデータに集約し、それを用いて第3節の図9.3.1と同様の判別を行っていることに相当します。

この時、2つの母集団間の平方距離が大きいほど2つの母集団が離れていることになり、判別確率が高くなります。 そのためこの中心間の平方距離のことを判別効率(discrimant efficiency)と呼ぶことがあります。 またこのように平方距離を単純化して計算した時のzが判別スコアになり、zを計算するための式が判別関数になります。

そして動脈硬化症群の判別スコアz1は平均値がD2/2、標準偏差がDの正規分布をし、正常群の判別スコアz2は平均値が-D2/2、標準偏差がDの正規分布をします。 この性質を利用すると、2群の中間点つまり判別スコアが0になる点までの汎距離に基づいて2つの母集団に含まれるデータの割合を計算することができます。 これが理論的な判別確率になります。

図9.3.4 マハラノビスの汎距離を利用した判別

一般に判別関数は次のような式で表されます。 (注3)

z = a0 + a1x1 + … + apxp
z > 0の時、動脈硬化症群と判別  z < 0の時、正常群と判別
判別効率(マハラノビスの平方距離)をD2とすると、D/2をp値変換して誤判別確率 = p (判別確率 = 1 - p)

2つの母集団間の平方距離が0の時は2つの母集団が重なっていることになり、判別は不可能になります。 そのため平方距離が0かどうか、つまり2つの母集団が重なっているかどうかを検定するホッテリング(Hotelling)のT2検定という手法が考案されています。 これは2標本t検定を多変量に拡張したものに相当し、「T2検定」という名前はそれを表しています。 そしてt検定と同様に、群と汎距離との関連性を表す指標として相関比を計算することができます。 しかし判別分析の場合は群の平均値の比較が目的ではなく個々の個体の判別が目的なので、検定結果や相関比よりも判別確率の方が重要です。 (注4) (→5.3 計数値の相関分析と回帰分析 (3)名義尺度と計量値の場合)

(3) 判別分析結果の解釈

表9.1.1のデータについて実際に計算すると次のようになります。 この時の判別スコアをグラフ表示したものが図9.1.1の一番下にある、2群のプロットが最も分離して見える方向のプロットです。

z = -61.1636 + 0.3352×TC - 0.0749×TG
ホッテリングのT2検定:T2 = 87.813 (p < 0.001)   判別効率(マハラノビスの平方距離):D2 = 14.6355   相関比:η = 0.89019
D/2 = 1.913  誤判別確率:p = 0.0279(判別確率:1 - p = 0.9721)
・この判別関数を用いて個々のデータを判別した場合
No.1:z = -61.1636 + 0.3352×220 - 0.0749×110 = 4.3414 > 0 … 動脈硬化症と判別
:(No.2〜No.24については自分で計算してみてください(^_-))
No.25:z = -61.1636 + 0.3352×240 - 0.0749×320 = -4.6836 < 0 … 正常と判別
∴感度 = 10/10 = 1(100%)  特異度 = 15/15 = 1(100%)  正診率 = 25/25 = 1(100%)

TGの判別係数の符号がTG単独で判別した時と逆転しているのは、重回帰分析の偏回帰係数と同じような理由によるものです。 判別係数の大きさと符号は、一応、群の判別に各変数が寄与している大きさとその方向性を反映しています。 しかし重回帰分析と同様に、各データの単位やバラツキ具合が異なっていると判別係数の値によって寄与の大きさを直接比較するわけにはいきません。 そこで重回帰分析における標準偏回帰係数と同様に、各変数が1標準偏差増加した時に判別スコアがいくつ変化するかを表す標準判別係数を計算し、これを各変数の寄与の大きさの目安にします。 標準判別係数は表9.1.1のデータでは次のようになります。

TCの標準判別係数 = 6.022  TGの標準判別係数 = -4.602

TCとTGの判別係数の絶対値を比べると、TCの方が1桁ほど大きくなっています。 しかし標準判別係数はそれほど大きな違いではなく、TGよりもTCの方が判別にやや大きく寄与している程度であることがわかります。

判別分析は、一応、群が原因で変数が結果という因果関係を想定しているものの、実際には変数としてリスクファクターを含めることもよくあります。 そのため判別分析の結果を解釈する時は重回帰分析以上にきめ細かい注意が必要になります。 判別分析結果を評価する時の注意点は以下のとおりです。 (→7.2 重回帰分析の解釈)

  1. 誤差の少ない信頼のおける多数のデータ——目安として変数の数の10倍以上の例数——に適用したものか?
  2. 判別分析に組み込んだ項目が適当か?
  3. 組み込んだ項目は診断指標なのかリスクファクターなのか?
  4. 高い判別確率が得られているか?
  5. 判別関数が科学的に納得がいくか?

(注1) 多変量正規分布の確率密度関数はベクトルと行列を利用して次のように表されます。 (→付録1 各種の確率分布 (9)多変量正規分布)


変数ベクトル:   母平均ベクトル:   母(分散)共分散行列:
μj:xjの母平均値   σjj = σj2:xjの母分散   σij = σji = ρijσiσj:xiとxjの母共分散
Σ = Σ'(対称行列)   |Σ|:Σの行列式   Σ-1Σの逆行列

p=1の時、当然のことながら、上式は次のように普通の正規分布の確率密度関数になります。

(注2) 普通の正規分布の確率密度関数において、指数項の中の(x - μ)22は平方距離d2に相当します。 したがって確率密度関数を次のように表すことができます。

  

これと同様に多変量正規分布の確率密度関数において、指数項の中の[-μ]'Σ-1[-μ]をD2で表すと確率密度関数を次のように表すことができます。


D2 = [ - μ]'Σ-1[ - μ]

このD2がマハラノビスの平方距離です。 この式からマハラノビスの汎距離Dがリーマン空間の線素に相当し、母共分散行列の逆行列Σ-1が計量行列に相当することがわかると思います。 そして変数xjがお互いに独立で分散が1の時、多変量正規分布は一様分布つまりユークリッド空間に相当し、マハラノビスの汎距離Dは通常のユークリッド距離に相当します。 (→当館の「ベクトルと行列 第13章 リーマン空間」参照)

また上式から変数がp次元正規分布をする時、平方距離が同じ値なら確率密度は等しくなることもわかると思います。 この性質を2次元でグラフ表示したものが図9.4.2の信頼楕円または等確率偏差楕円です。 そして通常の平方距離が自由度1のχ2分布をするのと同様に、マハラノビスの平方距離は自由度pのχ2分布をします。

そこでこれらの性質を利用して、品質管理における2σ法を多変量に拡張した多変量管理を行なうことができます。 2σ法は各データの偏差を調べ、それが標準偏差の2倍以上ある時は「何らかの異常によって得られた可能性の高いデータ」として検出する方法です。 これは標準化距離が標準正規分布をし、それが2(正確には1.96)以上になる確率は5%しかないことを利用しています。

|x - μ| > 2σの時警告 つまり の時警告

マハラノビスの平方距離を利用すれば、1変量の場合も多変量の場合も同じように「何らかの異常によって得られた可能性の高いデータ」を検出することができます。

D2 > χ2(p,0.05)の時警告
p = 1の時: の時警告
p = 2の時:D2 > χ2(2,0.05) = 5.991の時警告

(注3) p次元正規分布に従う2つの母集団g1とg2があり、それぞれの母平均値ベクトルと母共分散行列をμ1Σ1μ2Σ2とします。 今、どちらの群に属すのかわからない個体のデータベクトルをoとすると、最尤法の原理に従ってこの個体がどちらの群に属すのか判別することができます。

           
D12 = [o - μ1]'Σ1-1[o - μ1]   D22 = [o - μ2]'Σ2-1[o - μ2]
  
の時、g1に属すと判別
の時、g1に属すと判別

判別をもう少し簡単にするために2つの尤度を対数にし、その差をzとすると次のようになります。

さらに2群の母共分散行列が等しく、Σ1 = Σ2 = Σとすると次のようになります。

ここでΣ-1[μ1 - μ2]をと置くと次のようになります。

  
z = a0 + o' = a0 + a1xo1 + … + ajxoj + … + apxop

この時のzを判別スコア判別係数、式全体を(線形)判別関数と呼びます。 そして切片a0を除いたものを、特にフィッシャーの線形判別関数(Fisher's linear discriminant function)と呼ぶことがあります。 この判別関数を利用して個体を次のように判別することができます。

z > 0なら となり、g1に属すと判別
z < 0なら となり、g2に属すと判別

各変数が1標準偏差だけ変化した時の判別係数は、次のように各判別係数に変数の標準偏差をかけた値になります。

標準判別係数 (j = 1,…,p)
σj2 = σjjΣのj番目の対角要素

判別スコアが0になる点つまり境界値は、次のように2群の中心からの平方距離が等しくなる点になり、平方距離スケールにおける2群の中間点になります。 変数が1つの時の境界値は2群の母平均値の中点になり「大山鳴動して鼠一匹」という感じでしたが、変数がp個の時はさすがに多少は物々しく、猫10匹ぐらいの感じになります。

  ∴D12 = D22 = D2

zの期待値と分散を求めると次のようになります。


群1の場合:
群2の場合:

したがって2群のzの期待値の差δzは、次のように2群の中心間の平方距離になります。

群内分散に対するδzの平方の比Δ2判別効率と呼び、この値が大きいほど2群の判別が効率良く行なえることを表します。 この値を求めると、次のようにやはり2群の中心間の平方距離になります。

がp次元正規分布に従う時、D2は自由度pのχ2分布に従い、両群の判別スコアz1、z2は次のような正規分布に従います。

z1 〜 N(D2/2,D2)   z2 〜 N(-D2/2,D2)

このことと2群の境界値は判別スコアが0になる点であるという性質を利用して、誤判別確率と判別確率を求めることができます。 群g1に属す個体を群g2と誤って判別する確率をp12とし、反対に群g2に属す個体を群g1と誤って判別する確率をp21とすると次のようになります。



と標準化すると、z=0の時
∴誤判別確率:  判別確率=1 - p

このように誤判別確率は標準正規分布におけるD/2という値の片側確率になり、判別確率は1から誤判別確率を引いた値になります。

図9.4.4 誤判別確率と判別確率

実際のデータでは母平均値ベクトルと母共分散行列はほとんど不明ですから、それぞれの群に属すことがはっきりしている個体のデータを測定し、それに基いて母平均値ベクトルと母共分散行列を推定します。

群1のデータ行列:   平均値ベクトル:
積和行列:
S1jj = S1xjxj:x1jの平方和   S1ij = S1xixj:x1iとx1jの積和
群2のデータ行列:   平均値ベクトル:
積和行列:
S2jj = S2xjxj:x2jの平方和   S2ij = S2xixj:x2iとx2jの積和
=1 + 2   
μ11   μ22   Σ
-1[1 - 2]   D2 ≒ [1 - 2]'-1[1 - 2]

このような推定値を用いると、D2は自由度pのχ2分布ではなく第1自由度p、第2自由度(n+h-2)のF分布をします。 ただしχ2(p) = p・F(p,∞)という関係から、実際にはD2/pが第1自由度p、第2自由度(n+h-2)のF分布をします。

(注4) 1変量の場合、母平均値がある基準値と異なっているかどうかを検定する「1標本t検定」と、2群の母平均値が異なっているかどうかを検定する「2標本t検定」がありました。 これらを多変量に拡張した手法がホッテリングのT2検定です。 まず1群の時の1標本t検定は次のようになります。 (→3.1 1標本の計量値)

データベクトル:   平均値ベクトル:
  
H0:μ = μ0 に関して
> t(n-1,α)の時、有意水準100α%で有意

後の説明のためにt値の式を次のようにも書いておきましょう。


Fo の時、有意水準100α%で有意

母平均μの100(1 - α)%信頼区間は次のようになります。

  

これらを多変量に拡張したホッテリングのT2検定は次のようになります。 そしてこの検定における汎距離は効果量を多変量に拡張した指標になります。

データ行列:   平均値ベクトル:
  
H0μ = μ0 つまり D2 = 0 に関して

To2の時、有意水準100α%で有意
 または
> F(p,n-p,α)の時、有意水準100α%で有意

母平均値ベクトルμと、任意の変換η = 'μの100(1 - α)%信頼区間は次のようになります。

nD2 = n[-μ]'-1[-μ] ≦ Tα2
  

次に2群の場合の2標本t検定は次のようになります。 (→3.3 2標本の計量値5.1 相関係数と回帰直線 (注4))


H0:μ1 = μ2に関して
> t(n+h-2,α)の時、有意水準100α%で有意

Fo の時、有意水準100α%で有意
相関比:

母平均値の差δ = μ12の100(1-α)%信頼区間は次のようになります。

  

これらを多変量に拡張したホッテリングのT2検定は次のようになります。


H0μ1 = μ2 つまり D2 = 0 に関して

To2の時、有意水準100α%で有意
 または
> F(p,n+h-1-p,α)の時、有意水準100α%で有意
相関比:

母平均の差ベクトルδ = μ1 - μ2と、任意の変換δ = 'δの100(1-α)%信頼区間は次のようになります。


  

多変量データは例数が多いので、ホッテリングのT2検定はたいてい有意になります。 これは群の中心が基準位置とズレている、あるいは2群の中心位置がズレているということを意味しているだけで、群の判別について具体的な情報を提供してくれるわけではありません。 そしてそのズレつまり汎距離の科学的な許容範囲を決めるのは難しいので、ホッテリングのT2検定はたいてい有意性検定になります。 したがってこの検定は実質的にはほとんど無意味であり、判別確率の方が重要な意味を持っています。

今までしつこく述べてきたように、このことはどの検定にもあてはまることであり、実際の要約値とそれに対する科学的考察が最も大切です。