あまり使う場面はないかもしれませんが，回帰分析の目的変数に相当する部分が名義尺度データで，２つ以上のクラスを持つ場合の分析手法に正準判別分析（重判別分析とも呼ばれる）があります。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）ユーザーのための心理データの多変量解析法 —— 方法の理解から論文の書き方まで　教育出版