目的

正準相関分析を行う。
R にも cancor 関数があるが，出力されるのは標準化されていない係数だけで，
しかもその係数は「データの個数から1引いたものの平方根」で割り算されているものである。

使用法

my.cancor(x, gr1, gr2)

引数

x     データ行列
gr1   第一群の変数がある列位置のベクトル
gr2   第二群の変数がある列位置のベクトル

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/my_cancor.R", encoding="euc-jp")

# 正準相関分析
my.cancor <- function(       x,                                              # データ行列
                        gr1,                                            # 第一群の変数がある列位置のベクトル
                        gr2)                                            # 第二群の変数がある列位置のベクトル
{
        geneig2 <- function(a, b, k, sd)                             # 一般化固有値問題を解く関数
        {
            a <- as.matrix(a)
            b <- as.matrix(b)
            if (nrow(a) == 1) {
                res <- list(values=a/b, vectors=as.matrix(1))
            }
            else {
                res <- eigen(b)
                g <- diag(1/sqrt(res$values))
                v <- res$vectors
                res <- eigen(g %*% t(v) %*% a %*% v %*% g)
                res$vectors <-v %*% g %*% res$vectors
            }
            std.vectors <- res$vectors[,1:k]
            unstd.vectors <- std.vectors/sd
            list(values=res$values[1:k], std.vectors=std.vectors, unstd.vectors=unstd.vectors)
        }

        k <- min(length(gr1), length(gr2))                           # 第一変数群と第二変数群の個数の少ない方の個数
        x <- subset(x, complete.cases(x))                            # 欠損値を持つケースを除く
        r <- cor(x)                                                  # 相関係数行列
        S11 <- r[gr1, gr1, drop=FALSE]                                       # 第一変数群の相関係数
        S22 <- r[gr2, gr2, drop=FALSE]                                       # 第二変数群の相関係数
        S12 <- r[gr1, gr2, drop=FALSE]                                       # 第一変数群と第二変数群の相関係数
        x1 <- as.matrix(x[, gr1, drop=FALSE])                                # 第一変数群のデータ行列
        x2 <- as.matrix(x[, gr2, drop=FALSE])                                # 第二変数群のデータ行列
        sd1 <- apply(x1, 2, sd)                                              # 第一変数群の標準偏差
        sd2 <- apply(x2, 2, sd)                                              # 第二変数群の標準偏差
        res1 <- geneig2(S12 %*% solve(S22) %*% t(S12), S11, k, sd1)  # 第一変数群に対する解
        res2 <- geneig2(t(S12) %*% solve(S11) %*% S12, S22, k, sd2)  # 第二変数群に対する解
        score1 <- scale(x1 %*% res1[[3]])                            # 第一変数群に対する正準得点
        score2 <- scale(x2 %*% res2[[3]])                            # 第二変数群に対する正準得点
        list(canonical.correlation.coefficients=sqrt(res1[[1]]),
                standardized.coefficients=list(group1=res1[[2]], group2=res2[[2]]),
                coefficients=list(group1=res1[[3]], group2=res2[[3]]),
                canonical.scores=list(group1=score1, group2=score2))
}


使用例

奥野忠一ほか「多変量解析法」日科技連，379 ページからの数値例
x <- matrix(c(
2, 1, 2, 2,
1, 2, -1, -1,
0, 0, 0, 0,
-1, -2, -2, 1,
-2, -1, 1, -2), byrow=TRUE, ncol=4)
my.cancor(x, 1:2, 3:4)

出力結果例

> my.cancor(x, 1:2, 3:4)
$canonical.correlation.coefficients	# 正準相関係数
[1] 1.0000000 0.3015113

$standardized.coefficients	# 標準化された係数
$standardized.coefficients$group1	# 第一群の変数に対して
          [,1]          [,2]
[1,]  1.666667 -1.685062e-16
[2,] -1.333333  1.000000e+00

$standardized.coefficients$group2	# 第二群の変数に対して
             [,1]       [,2]
[1,] 1.064368e-18  1.0050378
[2,] 1.000000e+00 -0.1005038


$coefficients	# 係数
$coefficients$group1	# 第一群の変数に対して
          [,1]          [,2]
[1,]  1.054093 -1.065727e-16
[2,] -0.843274  6.324555e-01

$coefficients$group2	# 第二群の変数に対して
             [,1]        [,2]
[1,] 6.731652e-19  0.63564173
[2,] 6.324555e-01 -0.06356417


$canonical.scores	# 正準得点
$canonical.scores$group1	# 第一群の変数に対して
           [,1]          [,2]
[1,]  1.2649111  6.324555e-01
[2,] -0.6324555  1.264911e+00
[3,]  0.0000000 -2.220446e-17
[4,]  0.6324555 -1.264911e+00
[5,] -1.2649111 -6.324555e-01
attr(,"scaled:center")
[1] 0.000000e+00 2.220446e-17
attr(,"scaled:scale")
[1] 1 1

$canonical.scores$group2	# 第二群の変数に対して
           [,1]          [,2]
[1,]  1.2649111  1.144155e+00
[2,] -0.6324555 -5.720776e-01
[3,]  0.0000000 -2.220446e-17
[4,]  0.6324555 -1.334848e+00
[5,] -1.2649111  7.627701e-01
attr(,"scaled:center")
[1] 0.000000e+00 2.220446e-17
attr(,"scaled:scale")
[1] 1 1

 解説ページ