図1 線形判別関数による判別
<< 目次へ
最終更新日： 2008.11.27
判別分析モデルの応用
■ 2群の判別分析
判別分析とは既存のデータを基に未知のサンプル（被験体）がどの群に属するかを予測するための分析手法である。未知のサンプルがどの群に属するかを予測するための判別ルールを作るわけであるが、ここでは線形判別分析といって直線によって判別する方法を主に紹介する。
まずはもっとも単純な2群の判別分析の例を示す。表1のデータは男女の身長と体重のデータであり、仮に未知のサンプル（被験体）の身長と体重のデータだけが得られていたとする。このとき、No. 11のサンプル（被験体）は性別が不明であるので、判別分析によってNo. 11のサンプルの性別を予測することを考える。
| No. | 身長（Height） | 体重（Weight） | 性別（Sex） |
| 1 | 177 | 75 | male |
| 2 | 180 | 73 | male |
| 3 | 175 | 70 | male |
| 4 | 182 | 85 | male |
| 5 | 170 | 69 | male |
| 6 | 166 | 65 | female |
| 7 | 159 | 60 | female |
| 8 | 160 | 61 | female |
| 9 | 155 | 58 | female |
| 10 | 163 | 68 | female |
| 11 | 170 | 67 | ? |
判別分析のモデル式は次のように表すことができる。
Sex = α + β1Height + β2Weight
ここで従属変数Sexは2つのカテゴリをもつカテゴリカル型の変数であり、独立変数HeightとWeightは連続型の変数である。勘のいい人は気づいたかもしれないが、これは2値ロジットモデルとして解析されるべきモデルである。実際、表1のデータセットに対してロジスティック回帰分析を行っても問題なく、得られる結論も判別分析と同等であるといえる。これについては後述することにして、とりあえず判別分析モデルとして解析してみよう。
Rで判別分析を行うにはMASSパッケージに含まれているlda()という関数を用いる。
# MASSパッケージの呼び出し
> library(MASS)
# 身長・体重・性別のデータを用意する
> Height <- c(177, 180, 175, 182, 170, 166, 159, 160, 155, 163)
> Weight <- c(75, 73, 70, 85, 69, 65, 60, 61, 58, 68)
> Sex <- rep(c("male", "female"), c(5, 5))
> Sex <- factor(Sex, levels=c("male", "female"))
# lda()を用いて判別分析を行う
> disc.model <- lda(Sex ~ Height + Weight)
# 結果の出力
> disc.model
Call:
lda(Sex ~ Height + Weight)
Prior probabilities of groups: # 事前確率
male female
0.5 0.5
Group means: # 群平均
Height Weight
male 176.8 74.4
female 160.6 62.4
Coefficients of linear discriminants: # 線形係数
LD1
Height -0.3114300
Weight 0.0986707
Prior probabilities of groups: という部分は事前確率、つまりmaleとfemaleに属する確率は等しく0.5であると仮定した、ということである。これは引数priorで指定できるが、省略するとR側が勝手に計算してくれる。今回の場合だと、事前に分かっているmaleの割合は5 / 10 = 0.5であり、femaleの割合は5 / 10 = 0.5であるから、その値がPrior probabilities of group: として出力されているわけである。
Group means: は群平均であり、maleについて身長の平均値は176.8、体重は74.4である。femaleについて身長の平均値は160.6、体重は62.4であることを示している。
Coefficients of liner discriminants: は線形係数、要するに回帰分析でいうところの偏回帰係数のことである。判別分析のモデル式
Sex = α + β1Height + β2Weight
のβ1とβ2の値のことである。Heightの係数はβ1 = -0.311であり、Weightの係数はβ2 = 0.099であることが分かる。
定数項αは次のようにして計算することができる（少々面倒だが、自分で計算しなければならない）。
# 群平均のデータ行列をgroup.meanに代入
> (group.mean <- disc.model$mean)
Height Weight
male 176.8 74.4
female 160.6 62.4
# 線形係数のデータ行列をcoefに代入
> (coef <- disc.model$scaling)
LD1
Height -0.3114300
Weight 0.0986707
# 行列のかけ算は%*%で行える
# その結果をaに代入
> (a <- group.mean %*% coef)
LD1
male -47.71973
female -43.85861
# 行列aの列平均を計算する
# 引数に2ではなく1を指定すれば行に対して計算される
# 仮にmeanをsumにすると、平均ではなく合計が計算される
> apply(a, 2, mean)
LD1
-45.78917
# 以上の手順は1回で済ますこともできる
> apply(disc.model$means %*% disc.model$scaling, 2, mean)
LD1
-45.78917
以上から、次式のような線形判別式が得られる。
Sex = -(- 45.789) - 0.311*Height + 0.099*Weight
この判別式のHeightとWeightそれぞれにあるサンプル（被験体）のデータを代入することによって判別得点が得られる。そして、その判別得点がプラスの値かマイナスの値かによってそのサンプル（被験体）がどちらのグループに属するかを判別するというわけである。
# 線形判別関数を定義する > f <- function(H, W) -(- 45.789) - 0.311*H + 0.099*W # 引数にHeightとWeightを指定すれば、各サンプルの判別得点が計算される > f(Height, Weight) [1] -1.833 -2.964 -1.706 -2.398 -0.250 0.598 2.280 2.068 3.326 1.828
これはpredict()を利用して簡単に求めることができる（丸め誤差の関係で上記の値とは若干異なる）。
# 判別得点の出力
> predict(disc.model)$x
LD1
1 -1.9336427
2 -3.0652743
3 -1.8041362
4 -2.5040859
5 -0.3456566
6 0.5053807
7 2.1920375
8 1.9792782
9 3.2404163
10 1.7356830
表1のデータセットをプロットしたものに先ほど求めた判別式が描く直線を加えると図1のようになる。これを見ればよく分かるが、あるサンプル（被験体）のデータポイントが直線よりも上側にあれば（あるサンプル（被験体）の判別得点がプラスであれば）maleであると判別される。逆に下側（判別得点がマイナス）であればfemaleであると判別される。
# 図1を描くためのコード
plot(Weight, Height, col=c("blue", "red")[Sex], cex=1.5)
f2 <- function(x) (45.789/0.311) + (0.099/0.311) * x
x <- seq(55, 85, 1)
points(x, f2(x), type="l")
図1 線形判別関数による判別
ここでNo. 11のサンプルのデータ（Height = 170, Weight = 67）を判別式に代入して、このサンプルがどちらのグループに属するかを予測してみる。また、No. 11のサンプルを図示するにはpoints()を使えばよい。
# 先ほど定義した判別関数の中身を確認 > f function(H, W) -(- 45.789) - 0.311*H + 0.099*W # 身長170、体重67として判別得点を求める > f(170, 67) [1] -0.448 # No.11のサンプルを新たに書き加える > points(67, 170, pch=17, col="gray", cex=1.5)
図2 No.11のサンプルを書き加えたグラフ
結局、未知のサンプルはmaleであることが予測された（これは図2を見るとよく分かる）。しかし、このように図示できるのは独立変数の数が2つの場合のみである。独立変数が3つ以上あるときは頭の中で想像するしかないのである。もっとも3つならば3次元プロットができるが、見やすさの観点からいってもあまりおススメはできるものではない。
ところで、図1で描いた直線の方程式はどのようにして導かれたものであるかを説明しておこう。そもそもパラメータβ1、β2は次の連立方程式を解くことによって得られるものである。
α + β1Height + β2Weight = 0
Sex = α + β1Height + β2Weight
前者は2つのグループを分割する直線の方程式で、後者は線形判別関数である。そして定数項であるαは、群平均の中点を通るように求められることになる。
前述した線形判別関数
Sex = -(- 45.789) - 0.311*Height + 0.099*Weight
から、2つのグループを分割する方程式が
-(- 45.789) - 0.311*Height + 0.099*Weight = 0
であることが分かるだろう。この式において、定数項とWeightを左辺に移動させると、
- 0.311*Height = - 45.789 - 0.099*Weight
Height = (-45.789 / -0.311) + (-0.099 / -0.311)*Weight
となる。これが図1に描かれている直線の方程式である。
【補足1】 重回帰分析と判別分析
2群の判別分析と重回帰分析には関連がある。試しにlm()を使って解析して、そこから得られた直線の方程式がどのような直線であるかを図示してみよう。
# 図1のグラフを新しく描き直す
> plot(Weight, Height, col=c("blue", "red")[Sex], cex=1.5)
> points(x, f2(x), type="l")
# 性別を表す変数を連続型として用意する
> Sex2 <- rep(c(1, 0), c(5, 5))
> Sex2
[1] 1 1 1 1 1 0 0 0 0 0
# lm()を使って分析する
> lm.model <- lm(Sex2 ~ Height + Weight)
> lm.model
Call:
lm(formula = Sex2 ~ Height + Weight)
Coefficients:
(Intercept) Height Weight
-9.26337 0.06640 -0.02104
# 式の変形は上でやったとおり
> lm.f <- function(x) (-9.263 / -0.066) + (-0.021 / -0.066)*x
> points(x, lm.f(x), type="l", col="pink", lwd=2)
図3 重回帰によって得られた直線（ピンク色の直線）
図3を見れば明らかだが、2つの直線は切片（つまり定数項）が異なるだけで、傾きは同じである。つまり、2つの平行な直線が描かれるかたちになっている。これはどういうことかというと、重回帰によって得られた直線を、maleとfemaleの群平均の中点を通るように切片を調整しているということである。
具体的に例を示そう。
> disc.model$mean
Height Weight
male 176.8 74.4
female 160.6 62.4
これはmaleとfemaleの群平均である。maleについては(Weight, Height) = (74.4, 176.8)という座標に示される。femaleについては(Weight, Height) = (62.4, 160.6)という座標に示される。この2つの座標を図示し、さらにその2点を直線でつないでみると図4のようになる。
# 図4を描くためのコード： # 男性の座標 青いダイヤで示す points(74.4, 176.8, pch=18, col="blue", cex=1.5) # 女性の座標 赤いダイヤで示す points(62.4, 160.6, pch=18, col="red", cex=1.5) # 座標(74.4, 176.8)と(62.4, 160.6)を結ぶ直線を描く segments(74.4, 176.8, 62.4, 160.6, col="green")
図4 男性と女性の群平均の座標を直線で結ぶ
このように、2群の判別分析とは重回帰分析によって得られた直線の方程式を各群の群平均の中点を通るように定数項を調節したものである。つまり、重回帰分析のプログラムがあれば2群の判別分析が可能であるということである。
【補足2】 判別分析とロジスティック回帰分析
判別分析といっても分割直線で判別する方法とマハラノビスの距離で判別する方法の2つがある。今までにやってきたのは線形判別関数による判別分析で、これを拡張した非線形判別関数による判別分析（これを特に2次の判別分析という）もある。いずれにしても、ある情報を基にあるサンプルがどのグループに属するかを予測することが目的であることに変わりはない。
ここで改めて表1を見てみよう。このページの冒頭でも述べたように、これは2値ロジットモデルとして扱えるデータセットである。では実際に2値ロジットモデルとして解析した場合はどのような結果が得られるのかを試してみよう。
表1 男女別身長・体重データ
| No. | 身長（Height） | 体重（Weight） | 性別（Sex） |
| 1 | 177 | 75 | male |
| 2 | 180 | 73 | male |
| 3 | 175 | 70 | male |
| 4 | 182 | 85 | male |
| 5 | 170 | 69 | male |
| 6 | 166 | 65 | female |
| 7 | 159 | 60 | female |
| 8 | 160 | 61 | female |
| 9 | 155 | 58 | female |
| 10 | 163 | 68 | female |
| 11 | 170 | 67 | ? |
2値ロジットモデルを扱うにはRのglm()を用いるが、表1のデータセットは使うことができない。なぜかというと、glm()ではある独立変数が完全に従属変数を説明してしまうようなデータセットは受けつけないからである。モノは試しでやってみると、
> logit.model <- glm(Sex ~ Height + Weight, binomial) Warning messages: 1: 'pdf" ファイル引数 'C:\Documents and Settings\Owner\My Documents\R\workdir\a.pdf' を開くことができません 2: In glm.fit(x = X, y = Y, weights = weights, start = start, etastart = etastart, : 数値的に 0 か 1 である確率が生じました
といったエラーメッセージが出てしまうのが確認できるだろう。
それなので、少しデータセットを改変して性別を表す新たな変数new.Sexを用意することにしよう。
> new.Sex <- c(1, 1, 0, 1, 1, 1, 0, 0, 0, 0)
> new.Sex <- factor(new.Sex, levels=c(1, 0), labels=c("male", "female"))
> new.Sex
[1] male male female male male male female female female female
Levels: male female
これについて、glm()を用いたロジットモデルとlda()を用いた線形判別分析モデルとの結果を比べてみよう。
# ロジットモデルとして解析 logit.model <- glm(new.Sex ~ Height + Weight, binomial) # 線形判別分析モデルとして解析 d.model <- lda(new.Sex ~ Height + Weight)
ロジットモデルについて、関数fitted()を用いることによって各サンプルが「男性であるか否かという確率」が得られる。判別分析モデルについても、predict()をによって各サンプルが「男性のグループに属する確率と女性のグループに属する確率」が得られる。
> fitted(logit.model)
1 2 3 4 5 6 7 8
0.10345534 0.09279978 0.23665446 0.01376522 0.41361913 0.68295867 0.91941496 0.89616468
9 10
0.96317745 0.67799031
> predict(d.model)$posterior
male female
1 0.88262112 0.11737888
2 0.93518131 0.06481869
3 0.81321691 0.18678309
4 0.96575521 0.03424479
5 0.57647561 0.42352439
6 0.33660726 0.66339274
7 0.08549437 0.91450563
8 0.10685845 0.89314155
9 0.03484643 0.96515357
10 0.21208160 0.78791840
ロジットモデルの解析の結果より、例えば、サンプル1が男性でない確率は0.10（10%）であることがわかる。それに対して判別分析モデルの解析結果からは、サンプル1が男性でない確率（女性のグループに属する確率）は0.11（11%）と両者の結果はそれなりに一致している（少なくとも今回のケースにおいては）。
要するに（大雑把にいえば）
> lg.prob <- fitted(logit.model)
> lg.prob
1 2 3 4 5 6 7 8
0.10345534 0.09279978 0.23665446 0.01376522 0.41361913 0.68295867 0.91941496 0.89616468
9 10
0.96317745 0.67799031
と、
> d.prob <- predict(d.model)$posterior[,2]
> d.prob
1 2 3 4 5 6 7 8
0.11737888 0.06481869 0.18678309 0.03424479 0.42352439 0.66339274 0.91450563 0.89314155
9 10
0.96515357 0.78791840
は同じ事実を示しているのである。
それではロジスティック回帰分析と判別分析は同じ手法として、同様のデータセットに適用しても問題ないのかというと、そういうものでもない。この2つの手法による違いについて、どのような場面でどのような問題が生じうるのか、それは私にも分からない。もし興味があるならば、そういった論文がないかを探してみるとよいかもしれない。
■ 3群以上の判別分析
2群の判別分析の節においては、直線によって2群を分割して判別するという、線形判別関数による判別の方法を示した。これを拡張して、3群以上の場合も複数の直線によって判別を行うことができる。この方法を特に正準判別分析というが、一般に3群以上の判別を行う場合はマハラノビスの距離による判別法を用いた方がよい。この方法は線形モデルの範疇ではないので、ここでの解説は省略することにする。
