正準判別分析 Last modified: Aug 25, 2009
目的
正準判別分析を行う。
MASS パッケージにある lda は，3群以上の判別の場合には，正準判別分析を行っている。
そして，このプログラムも lda も，全く同じ答えを出す。
使用法
candis(data, group)
print.candis(obj)
summary.candis(obj)
plot.candis(obj, pch=as.integer(obj$group), col=as.integer(obj$group))
引数
data データ行列（行がケース，列が変数）
group 各ケースがどの群であるかを表す変数
obj candis が返すオブジェクト
3 群以上の判別図
pch プロット記号（判別する群の数と同じ長さの整数ベクトル）
col プロット色（判別する群の数と同じ長さの整数ベクトル）
xpos, ypos legend 関数の x, y 引数（座標値の他，xpos="topright" など）
... plot への引数
2 群の判別図
which 図の種類 boxplot か barplot
nclass barplot のときのおよその階級数
... boxplot, barplot への引数
ソース
以下のほかに，一般化固有値問題を解く関数が必要。
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/candis.R", encoding="euc-jp")
# 正準判別分析
candis <- function( data, # 説明変数データ行列
group) # 群変数
{
vnames <- colnames(data)
group <- as.factor(as.matrix(group)) # 群を factor に変換
OK <- complete.cases(data, group) # 欠損値を持つケースを除く
data <- as.matrix(data[OK,])
colnames(data) <- vnames
group <- group[OK]
p <- ncol(data) # 変数の個数
n <- nrow(data) # ケース数
n.i <- table(group) # 各群のケース数
k <- length(n.i) # 群の個数
group.means <- matrix(unlist(by(data, group, mean)), p) # 各変数の群ごとの平均値
grand.means <- colMeans(data) # 各変数の全体の平均値
means <- cbind(grand.means, group.means) # 結果表示のためにまとめる
vars <- by(data, group, function(x) var(x)*(nrow(x)-1)) # 各群の変動・共変動行列
w <- matrix(rowSums(matrix(mapply("+", vars), ncol=k)), p) # 群内変動・共変動行列
b <- (var(data)*(n-1))-w # 群間変動・共変動行列
pooled.cov <- w/(n-k) # プールされた分散・共分散行列
result <- geneig(b, w) # 一般化固有値問題を解く
eigen.values <- result$values # 固有値
nax <- sum(eigen.values > 1e-10) # 解の個数
eigen.values <- eigen.values[1:nax] # 必要な個数の固有値を取り出す
eigen.vectors <- result$vectors[,1:nax,drop=FALSE] # 固有ベクトルも取り出す
canonical.corr.coef <- sqrt(eigen.values/(1+eigen.values)) # 正準相関係数
temp <- diag(t(eigen.vectors) %*% w %*% eigen.vectors)
temp <- 1/sqrt(temp/(n-k))
coeff <- eigen.vectors*temp # 正準判別係数
const <- as.vector(-grand.means %*% coeff) # 定数項
coefficient <- rbind(coeff,const) # 両方併せる
std.coefficient <- sqrt((n-1)/n)*apply(data, 2, sd) * coeff # 標準化正準判別係数
centroids <- t(t(coeff) %*% group.means+const) # 重心
can.score <- t(t(data %*% coeff)+const) # 正準得点
d <- sapply(1:k, function(i) colSums((t(can.score)-centroids[i,])^2)) # 二乗距離
p <- pchisq(d, p, lower.tail=FALSE) # ケースがそれぞれの群に属するとしたとき，その判別値を取る確率
p.Bayes <- t(t(exp(-d/2))*as.numeric(n.i)) # 各ケースが各群に所属するベイズ確率
p.Bayes <- p.Bayes/rowSums(p.Bayes)
gname <- levels(group)
classification <- factor(gname[max.col(p.Bayes)], levels=gname)
colnames(means) <- c("grand mean", gname)
colnames(w) <- rownames(w) <- colnames(pooled.cov) <- rownames(pooled.cov) <- rownames(std.coefficient) <- rownames(means)
colnames(coefficient) <- colnames(std.coefficient) <- colnames(centroids) <- colnames(can.score) <- paste("axis", 1:nax)
rownames(coefficient) <- c(rownames(means), "constant")
rownames(centroids) <- colnames(p.Bayes) <- colnames(p) <- gname
rownames(can.score) <- rownames(p.Bayes) <- rownames(p) <- paste("case", 1:n)
return(structure(list(means=means, between.cov=b, within.cov=w, pooled.cov=pooled.cov,
eigen.values=eigen.values, canonical.corr.coef=canonical.corr.coef,
std.coefficient=std.coefficient, coeff=coefficient, centroids=centroids,
can.score=can.score, p.Bayes=p.Bayes, p=p, classification=classification, group=group, ngroup=k), class="candis"))
}
# print メソッド
print.candis <- function( obj, # candis が返すオブジェクト
digits=5) # 結果の表示桁数
{
cat("\n判別係数\n\n")
print(round(obj$coeff, digits=digits))
cat("\n判別結果\n\n")
print(xtabs(~obj$group+obj$classification))
}
# summary メソッド
summary.candis <- function( obj, # candis が返すオブジェクト
digits=5) # 結果の表示桁数
{
print.default(obj, digits=digits)
}
# plot メソッド
plot.candis <- function(obj, # candis が返すオブジェクト
pch=1:obj$ngroup, # 3 群以上の判別時の記号
col=1:obj$ngroup, # 3 群以上の判別時の色
xpos="topright", ypos=NULL, # 3 群以上の判別時の legend 関数の x, y 引数
which=c("boxplot", "barplot"), # 2 群判別のときのグラフの種類
nclass=20, # barplot のおよその階級数
...) # plot, boxplot, barplot への引数
{
score <- obj$can.score
group <- obj$group
if (ncol(score) >= 2) {
int.group <- as.integer(group)
plot(score, pch=pch[int.group], col=col[int.group], ...)
legend(x=xpos, y=ypos, legend=levels(group), pch=pch, col=col)
}
else {
which <- match.arg(which)
if (which == "boxplot") { # boxplot
plot(score ~ group, xlab="群", ylab="判別値", ...)
}
else { # barplot
tbl <- table(group, cut(score, breaks=pretty(score, n=nclass)))
barplot(tbl, beside=TRUE, legend=TRUE, ...)
}
}
}
使用例
> data(iris)
> (result <- candis(iris[1:4], iris[5]))
print メソッドでは以下の結果だけを表示する
判別係数
axis 1 axis 2
Sepal.Length 0.82938 0.02410
Sepal.Width 1.53447 2.16452
Petal.Length -2.20121 -0.93192
Petal.Width -2.81046 2.83919
constant 2.10511 -6.66147
判別結果
classification
group setosa versicolor virginica
setosa 50 0 0
versicolor 0 48 2
virginica 0 1 49
> summary(result)
sumamry メソッドでは，すべての要素を表示することができる。
$means # 全体および各群の平均値
grand mean setosa versicolor virginica
Sepal.Length 5.8433 5.006 5.936 6.588
Sepal.Width 3.0573 3.428 2.770 2.974
Petal.Length 3.7580 1.462 4.260 5.552
Petal.Width 1.1993 0.246 1.326 2.026
$between.cov # 群間平方和・積和行列
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 63.212 -19.953 165.25 71.279
Sepal.Width -19.953 11.345 -57.24 -22.933
Petal.Length 165.248 -57.240 437.10 186.774
Petal.Width 71.279 -22.933 186.77 80.413
$within.cov # 群内平方和・積和行列
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 38.956 13.6300 24.6246 5.6450
Sepal.Width 13.630 16.9620 8.1208 4.8084
Petal.Length 24.625 8.1208 27.2226 6.2718
Petal.Width 5.645 4.8084 6.2718 6.1566
$pooled.cov # プールした分散・共分散行列
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 0.265008 0.092721 0.167514 0.038401
Sepal.Width 0.092721 0.115388 0.055244 0.032710
Petal.Length 0.167514 0.055244 0.185188 0.042665
Petal.Width 0.038401 0.032710 0.042665 0.041882
$eigen.values # 固有値
[1] 32.19193 0.28539
$canonical.corr.coef # 正準相関係数
[1] 0.98482 0.47120
$std.coefficient # 標準化判別係数
axis 1 axis 2
Sepal.Length 0.68449 0.019892
Sepal.Width 0.66659 0.940292
Petal.Length -3.87282 -1.639626
Petal.Width -2.13509 2.156910
$coeff # 判別係数（生データに対応するもの）
axis 1 axis 2
Sepal.Length 0.82938 0.024102
Sepal.Width 1.53447 2.164521
Petal.Length -2.20121 -0.931921
Petal.Width -2.81046 2.839188
constant 2.10511 -6.661473
$centroids # 各群の重心
axis 1 axis 2
setosa 7.6076 0.21513
versicolor -1.8250 -0.72790
virginica -5.7826 0.51277
$can.score # 各ケースの判別スコア
axis 1 axis 2
case 1 8.06180 0.3004206
case 2 7.12869 -0.7866604
case 3 7.48983 -0.2653845
case 4 6.81320 -0.6706311
case 5 8.13231 0.5144625
【中略】
case 147 -5.17956 -0.3634750
case 148 -4.96774 0.8211405
case 149 -5.88615 2.3450905
case 150 -4.68315 0.3320338
$p.Bayes # 各ケースが各群に所属するベイズ確率
setosa versicolor virginica
case 1 1.0000e+00 3.8964e-22 2.6112e-42
case 2 1.0000e+00 7.2180e-18 5.0421e-37
case 3 1.0000e+00 1.4638e-19 4.6759e-39
case 4 1.0000e+00 1.2685e-16 3.5666e-35
case 5 1.0000e+00 1.6374e-22 1.0826e-42
【中略】
case 147 4.6166e-36 5.8988e-03 9.9410e-01
case 148 5.5490e-35 3.1459e-03 9.9685e-01
case 149 1.6137e-40 1.2575e-05 9.9999e-01
case 150 2.8580e-33 1.7542e-02 9.8246e-01
$p # ケースがそれぞれの群に属するとしたとき，その判別値をとる確率
setosa versicolor virginica
case 1 9.9469e-01 1.7650e-20 2.2729e-40
case 2 8.7264e-01 1.6010e-16 2.3190e-35
case 3 9.9309e-01 5.7625e-18 3.6980e-37
case 4 8.4147e-01 2.3946e-15 1.4239e-33
case 5 9.8525e-01 7.0052e-21 8.8231e-41
【中略】
case 147 2.1872e-34 2.2556e-02 8.8926e-01
case 148 3.0565e-33 1.5412e-02 9.4386e-01
case 149 2.8248e-39 3.2602e-05 4.9821e-01
case 150 1.1969e-31 5.4196e-02 8.7125e-01
> plot(result, xpos="top")
plot メソッドにより，判別図を描くことができる
（主成分分析の結果と比較してみよ）
二群判別の場合
> ( result2 <- candis(iris[51:150, 1:4], iris[51:150, 5]) )
判別係数
axis 1
Sepal.Length 0.94312
Sepal.Width 1.47943
Petal.Length -1.84845
Petal.Width -3.28473
constant 4.41899
判別結果
classification
group versicolor virginica
versicolor 48 2
virginica 1 49
> plot(result2)
> plot(result2, which="barplot", nclass=40, args.legend=list(x="top"))
解説ページ
直前のページへ戻る E-mail to Shigenobu AOKI
