あまり使う場面はないかもしれませんが，回帰分析の目的変数に相当する部分が名義尺度データで，２つ以上のクラスを持つ場合の分析手法に正準判別分析（重判別分析とも呼ばれる）があります。Rで正準判別分析を実行する関数には，MASSパッケージのlda()やcandisc()パッケージのcandiscがあるのですが，lda()は判別係数しか出力してくれず，candisc()は多変量統計量とその検定結果だけしか出力してくれないので，結果オブジェクトの中で必要な値を探し回ったりしなければならなくて面倒でした。
最近また正準判別分析を使う機会があり，やはりそれが面倒だったので，もう少し一度にいろいろ出力できるよう，正準判別分析用の関数を作成してみました。関数の内部にさらにWilksのラムダを算出するための関数と，candisc用のsummary関数を含んでいるので少し長めですが，実際にやっていることは多変量のlmオブジェクトを受け取って，必要な値を算出し，清書しているだけのシンプルなものです。
関数名についてはいい案が浮かばなかったので，関数の名前はcandiscパッケージのcandisc()と被ってしまっていますが，自家製版candisc()を使うときにはcandiscパッケージを使うことはないのでよしとします。
正準判別分析（改変版） candisc()
candisc<-function(object) #object = 多変量 lmオブジェクト lm(cbind(X1,X2) ~ Y, data=data)
{
#### Wilksのラムダの算出関数
Wilks <- function(eig, q, df.res) {
test <- prod(1/(1 + eig))
p <- length(eig)
tmp1 <- df.res - 0.5 * (p - q + 1)
tmp2 <- (p * q - 2)/4
tmp3 <- p^2 + q^2 - 5
tmp3 <- ifelse (tmp3 > 0, sqrt(((p * q)^2 - 4)/tmp3), 1)
c(test, ((test^(-1/tmp3) - 1) * (tmp1 * tmp3 - 2 * tmp2))/p/q, p * q, tmp1 * tmp3 - 2 * tmp2)
}
#### 結果出力用の関数
summary.candisc<-function (object)
{
ans<-list(
"Class Level Information"=object$means,
"Canonical Discriminant Analysis"=object$eigen,
"Statistics"=object$stat,
"Cannonical Variate Coefficients"=object$coeff,
"Class Means on Canonical Variables"=object$means2)
class(ans) <-c("summary.candisc","listof")
ans
}
#### 係数，統計量の算出
m <- min(object$rank-1,ncol(object$coefficients))
w <- resid(object)
W <- t(w)%*%w
t <- resid(update(object,~1))
T <- t(t)%*%t
B <- T-W
We <- eigen(W,symmetric=TRUE)
W1 <- solve(t(We$vectors %*% diag(sqrt(We$values))))
ev <- eigen( t(W1)%*% B %*% W1)
cv <- W1 %*% ev$vectors * sqrt(nrow(object$model)-object$rank)
cv <- cv[,1:m]
e <- ev$values
corsq <- e[1:m]/(1+e[1:m])
prop <- e[1:m]/sum(e[1:m])
cum <- cumsum(prop)
emat <-data.frame(CanCorr = sqrt(corsq), EigenValue=e[1:m], Proportion=prop, Cumulative=cum)
if (m>1){
colnames(cv) <-colnames(cv,do.NULL=FALSE,prefix="CAN")
rownames(cv) <-colnames(object$coefficients,do.NULL=FALSE)
}
STAT<-matrix(c(rep(0,m*5)),ncol=5)
for(i in 1:m){
lambda <- Wilks(e[i:NROW(e)],object$rank-i,object$df.residual)
q <- lambda[2]
df1 <- lambda[3]
df2 <- lambda[4]
pval <- round(pf(q,df1,df2,lower=F),8)
STAT[i,] <- c(lambda[1],q,df1,df2,pval)
}
colnames(STAT) <-c("Lambda","ApproxF","NumDF","DenDF","Prob")
STAT<-as.data.frame(STAT)
#### クラス情報
f_table<-as.vector(table(object$model[,-1]))
df<-data.frame(Class=object$model[,-1],object$fitted.values)
meanTable <- aggregate(.~Class , df,mean)
rownames(meanTable)<-meanTable[,1]
meanTable$Class <- NULL
meanTable$Frequency<-f_table
meanTable$Proportion<-f_table/sum(f_table)
#### 判別関数のクラス別平均
df<-data.frame(Class=object$model[,-1], (object$fitted.values+object$residuals) %*% cv)
canMean <- aggregate(.~ Class, df, mean)
rownames(canMean)<-canMean[,1]
canMean$Class <- NULL
(summary.candisc(structure(list( "eigen"=emat,
"stat"=STAT,
"coeff"=cv,
"means"=meanTable,
"means2"=canMean))))
}
この関数を使って，有名なirisデータの正準判別分析を実行してみます。candisc()関数には，多変量のlmオブジェクトを渡します。多変量lmオブジェクトは，次のような書式で作成します。
model <- lm( cbind(Petal.Length, Sepal.Length, Petal.Width, Sepal.Width) ~ Species, data= iris)
じつを言うと，この書き方では「花弁（Petal）の長さと幅，萼片（Sepal）の長さと幅を品種（Species）で説明する」という意味で，説明変数（独立変数）と目的変数（従属変数）の関係が逆になり，多変量分散分析などに使用されるモデルになってしまうのですが，目的変数の方に数値変数でないものを入れるとうまくいかないのと，多変量分散分析と正準判別分析は面と裏のような関係なので，この形でlmオブジェクトを作成します。なお，candiscパッケージのcandisc()でも，入力にはこれと同じlmオブジェクトをとります。
さて，分析の結果は次のように出力されます。
> candisc(model)
Class Level Information :
Petal.Length Sepal.Length Petal.Width Sepal.Width Frequency Proportion
setosa 1.462 5.006 0.246 3.428 50 0.3333333
versicolor 4.260 5.936 1.326 2.770 50 0.3333333
virginica 5.552 6.588 2.026 2.974 50 0.3333333
Canonical Discriminant Analysis :
CanCorr EigenValue Proportion Cumulative
1 0.9848209 32.191929 0.991212605 0.9912126
2 0.4711970 0.285391 0.008787395 1.0000000
Statistics :
Lambda ApproxF NumDF DenDF Prob
1 0.02343863 199.1453 8 288 0e+00
2 0.77797337 13.7939 3 145 6e-08
Cannonical Variate Coefficients :
CAN1 CAN2
Petal.Length 2.2012117 0.93192121
Sepal.Length -0.8293776 -0.02410215
Petal.Width 2.8104603 -2.83918785
Sepal.Width -1.5344731 -2.16452123
Class Means on Canonical Variables :
CAN1 CAN2
setosa -5.502493 -6.876606
versicolor 3.930156 -5.933573
virginica 7.887657 -7.174239
まず最初に，それぞれのクラス（品種）ごとの測定値の平均値と度数に関する情報が表示されます。
Class Level Information :
Petal.Length Sepal.Length Petal.Width Sepal.Width Frequency Proportion
setosa 1.462 5.006 0.246 3.428 50 0.3333333
versicolor 4.260 5.936 1.326 2.770 50 0.3333333
virginica 5.552 6.588 2.026 2.974 50 0.3333333
次が正準判別分析に関する統計量です。
Canonical Discriminant Analysis :
CanCorr EigenValue Proportion Cumulative
1 0.9848209 32.191929 0.991212605 0.9912126
2 0.4711970 0.285391 0.008787395 1.0000000
Statistics :
Lambda ApproxF NumDF DenDF Prob
1 0.02343863 199.1453 8 288 0e+00
2 0.77797337 13.7939 3 145 6e-08
このブロックの前半には，それぞれの判別関数の正準相関係数と固有値，そしてその固有値が全体に占める比率が，ブロックの後半にはそれぞれの判別関数についてのWilksのラムダ，近似F統計量と有意確率が示されます。
そして最後に，判別関数の係数（非標準化係数）と，それぞれのクラスにおける正準判別関数の平均値です。
Cannonical Variate Coefficients :
CAN1 CAN2
Petal.Length 2.2012117 0.93192121
Sepal.Length -0.8293776 -0.02410215
Petal.Width 2.8104603 -2.83918785
Sepal.Width -1.5344731 -2.16452123
Class Means on Canonical Variables :
CAN1 CAN2
setosa -5.502493 -6.876606
versicolor 3.930156 -5.933573
virginica 7.887657 -7.174239
なお，これらの係数は標準化されていない値ですので，標準化係数が必要な場合は説明変数をあらかじめ標準化したうえで分析します。標準化した場合の結果は次のようになります（必要部分のみ）。
> iris.scaled<-iris
> iris.scaled[,1:4]<-scale(iris[,1:4])
>
> model <- lm( cbind(Petal.Length, Sepal.Length, Petal.Width, Sepal.Width) ~ Species, data= iris.scaled)
>
> candisc(model)
[ 省 略 ]
Cannonical Variate Coefficients :
CAN1 CAN2
Petal.Length 3.8857950 1.64511887
Sepal.Length -0.6867795 -0.01995817
Petal.Width 2.1422387 -2.16413593
Sepal.Width -0.6688251 -0.94344183
Class Means on Canonical Variables :
CAN1 CAN2
setosa -7.607600 -0.2151330
versicolor 1.825049 0.7278996
virginica 5.782550 -0.5127666
どのような目的で判別分析を使用するのかによっても必要な値は異なるのでしょうが，心理データの解析と報告には最低限これぐらいの出力が必要になると思います。なお，判別分析の結果としてどの数値を報告すべきかについては，山際・田中（1997）を参考にさせてもらいました。
参考までに，lda()とcandiscパッケージ版のcandisc()での分析結果を載せておきます。
lda()による判別分析とその結果
> lda(Species~., data=iris)
Call:
lda(Species ~ ., data = iris)
Prior probabilities of groups:
setosa versicolor virginica
0.3333333 0.3333333 0.3333333
Group means:
Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa 5.006 3.428 1.462 0.246
versicolor 5.936 2.770 4.260 1.326
virginica 6.588 2.974 5.552 2.026
Coefficients of linear discriminants:
LD1 LD2
Sepal.Length 0.8293776 0.02410215
Sepal.Width 1.5344731 2.16452123
Petal.Length -2.2012117 -0.93192121
Petal.Width -2.8104603 2.83918785
Proportion of trace:
LD1 LD2
0.9912 0.0088
candisc()による判別分析とその結果
> candisc::candisc(model)
Canonical Discriminant Analysis for Species:
CanRsq Eigenvalue Difference Percent Cumulative
1 0.96987 32.19193 31.907 99.12126 99.121
2 0.22203 0.28539 31.907 0.87874 100.000
Test of H0: The canonical correlations in the
current row and all that follow are zero
LR test stat approx F numDF denDF Pr(> F)
1 0.02344 199.145 8 288 < 2.2e-16 ***
2 0.77797 13.794 3 145 5.794e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
参考図書
- 山際 勇一郎・田中 敏（1997） ユーザーのための心理データの多変量解析法 —— 方法の理解から論文の書き方まで 教育出版