特異値分解・主成分分析・バイプロット

主成分分析と因子分析で説明したことを、より詳しく説明します。

特異値分解

ここでは縦長の行列 $A$ を考えます。$A$ の特異値分解とは、$A$ を次のような形の三つの行列の積に分解することです。


$A$ = $U$ $D$ $V'$

$U$ は直交行列の右側を削って $A$ と同じ形にしたものです。直交行列ですから、各列の2乗和が 1 で、各列は直交しています(内積 0)。$D$ は対角行列で、左上の対角成分ほど大きいように並べられています。$V'$ は直交行列 $V$ の転置です($V'$ も直交行列です)。

主成分分析

主成分分析は、各変数からその平均値を引いた(場合によってはさらに標準偏差で割った)縦長の行列を特異値分解することと同じです。

例えば主成分分析と因子分析での例を考えます:

chu = read.csv("http://okumuralab.org/~okumura/stat/data/atest2014chu.csv",
               fileEncoding="CP932")
row.names(chu) = chu[,1]
chu = chu[,2:5]

これを主成分分析した結果:

> prcomp(chu)
Standard deviations (1, .., p=4):
[1] 4.5647586 1.1615710 0.6880953 0.3647093

Rotation (n x k) = (4 x 4):
            PC1        PC2         PC3        PC4
国語A 0.3236353  0.4347444 -0.50391566 -0.6725522
国語B 0.4296899  0.6692517  0.04669831  0.6043906
数学A 0.5571265 -0.5789391 -0.52250679  0.2853527
数学B 0.6326427 -0.1671196  0.68620328 -0.3177410

> prcomp(chu)$x
                  PC1          PC2          PC3         PC4
北海道    -1.49201813  0.022160600  0.590740297 -0.89633743
青森県     2.58910026 -0.004574465 -0.949666441 -0.17459994
(略)
沖縄県   -15.06049683  1.007170004  0.740618970  0.53323546

同じことを特異値分解でやってみましょう。まずは各科目から平均値を引きます:

> a = t(t(chu) - colMeans(chu))
> a
               国語A       国語B        数学A       数学B
北海道   -0.16808511 -1.14042553 -1.408510638 -0.25744681
青森県    1.43191489  0.95957447  1.891489362  1.04255319
(略)
沖縄県   -5.16808511 -5.44042553 -9.208510638 -9.35744681
> r = svd(a)
> r
$d
[1] 30.959699  7.878158  4.666889  2.473579

$u
               [,1]          [,2]          [,3]         [,4]
 [1,] -0.0481922679  0.0028129164  0.1265811654 -0.362364581
 [2,]  0.0836280814 -0.0005806516 -0.2034902399 -0.070585956
(略)
[47,] -0.4864548788  0.1278433390  0.1586964911  0.215572438

$v
          [,1]       [,2]        [,3]       [,4]
[1,] 0.3236353  0.4347444 -0.50391566 -0.6725522
[2,] 0.4296899  0.6692517  0.04669831  0.6043906
[3,] 0.5571265 -0.5789391 -0.52250679  0.2853527
[4,] 0.6326427 -0.1671196  0.68620328 -0.3177410

この r$dr$ur$v がそれぞれ $D$ の対角要素、$U$、$V'$ です。r$d を $\sqrt{47-1}$ で割ったものは、主成分の Standard deviations として出力された値と一致します(47 は都道府県の数です):

> r$d / sqrt(47 - 1)
[1] 4.5647586 1.1615710 0.6880953 0.3647093

$UD$ が各都道府県の主成分得点(さきほどの prcomp(chu)$x)です:

> r$u %*% diag(r$d)
              [,1]         [,2]         [,3]        [,4]
 [1,]  -1.49201813  0.022160600  0.590740297 -0.89633743
 [2,]   2.58910026 -0.004574465 -0.949666441 -0.17459994
(略)
[47,] -15.06049683  1.007170004  0.740618970  0.53323546

r$v つまり $V'$ は、主成分分析の Rotation として出力されたものとまったく同じです。

バイプロット

バイプロットは特異値分解の第2成分(主成分分析の第2主成分)までを図示したものです。特異値分解の形で書けば、$A = UDV'$ を $UD$ と $V'$ の内積または $U$ と $DV'$ の内積と見て、前半を黒、後半を赤矢印で表したものです。

ヘルプ ?biplot.prcomp を見るとわかるように、デフォルトで biplot(prcomp(chu)) とすれば、$U$ を黒、$DV'$ を赤で表します。また、biplot(prcomp(chu), scale=0) とすれば、$UD$ を黒、$V'$ を赤で表します(こちらのほうが主成分分析と対応づけやすい形です)。

黒も赤も、文字の位置の座標を見てください。赤の矢印は、文字と重ならないように、長さを 0.8 倍して表示してあります。

念のためにバイプロットのソースも引用しておきます(getAnywhere(biplot.prcomp)getAnywhere(biplot.default) で見られます):

biplot.prcomp <- function (x, choices = 1L:2L, scale = 1, pc.biplot = FALSE, ...) 
{
    if (length(choices) != 2L) 
        stop("length of choices must be 2")
    if (!length(scores <- x$x)) 
        stop(gettextf("object '%s' has no scores", deparse1(substitute(x))), domain = NA)
    if (is.complex(scores)) 
        stop("biplots are not defined for complex PCA")
    lam <- x$sdev[choices]
    n <- NROW(scores)
    lam <- lam * sqrt(n)
    if (scale < 0 || scale > 1) 
        warning("'scale' is outside [0, 1]")
    if (scale != 0) 
        lam <- lam^scale
    else lam <- 1
    if (pc.biplot) 
        lam <- lam/sqrt(n)
    biplot.default(t(t(scores[, choices])/lam), t(t(x$rotation[, choices]) * lam), ...)
    invisible()
}

biplot.default <- function (x, y, var.axes = TRUE, col, cex = rep(par("cex"), 2), 
    xlabs = NULL, ylabs = NULL, expand = 1, xlim = NULL, ylim = NULL, 
    arrow.len = 0.1, main = NULL, sub = NULL, xlab = NULL, ylab = NULL, 
    ...) 
{
    n <- nrow(x)
    p <- nrow(y)
    if (missing(xlabs)) {
        xlabs <- dimnames(x)[[1L]] %||% 1L:n
    }
    xlabs <- as.character(xlabs)
    dimnames(x) <- list(xlabs, dimnames(x)[[2L]])
    if (missing(ylabs)) {
        ylabs <- dimnames(y)[[1L]] %||% paste("Var", 1L:p)
    }
    ylabs <- as.character(ylabs)
    dimnames(y) <- list(ylabs, dimnames(y)[[2L]])
    if (length(cex) == 1L) 
        cex <- c(cex, cex)
    if (missing(col)) {
        col <- par("col")
        if (!is.numeric(col)) 
            col <- match(col, palette(), nomatch = 1L)
        col <- c(col, col + 1L)
    }
    else if (length(col) == 1L) 
        col <- c(col, col)
    unsigned.range <- function(x) c(-abs(min(x, na.rm = TRUE)), 
        abs(max(x, na.rm = TRUE)))
    rangx1 <- unsigned.range(x[, 1L])
    rangx2 <- unsigned.range(x[, 2L])
    rangy1 <- unsigned.range(y[, 1L])
    rangy2 <- unsigned.range(y[, 2L])
    if (missing(xlim) && missing(ylim)) 
        xlim <- ylim <- rangx1 <- rangx2 <- range(rangx1, rangx2)
    else if (missing(xlim)) 
        xlim <- rangx1
    else if (missing(ylim)) 
        ylim <- rangx2
    ratio <- max(rangy1/rangx1, rangy2/rangx2)/expand
    on.exit(par(op))
    op <- par(pty = "s")
    if (!is.null(main)) 
        op <- c(op, par(mar = par("mar") + c(0, 0, 1, 0)))
    plot(x, type = "n", xlim = xlim, ylim = ylim, col = col[1L], 
        xlab = xlab, ylab = ylab, sub = sub, main = main, ...)
    text(x, xlabs, cex = cex[1L], col = col[1L], ...)
    par(new = TRUE)
    dev.hold()
    on.exit(dev.flush(), add = TRUE)
    plot(y, axes = FALSE, type = "n", xlim = xlim * ratio, ylim = ylim * 
        ratio, xlab = "", ylab = "", col = col[1L], ...)
    axis(3, col = col[2L], ...)
    axis(4, col = col[2L], ...)
    box(col = col[1L])
    text(y, labels = ylabs, cex = cex[2L], col = col[2L], ...)
    if (var.axes) 
        arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], 
            length = arrow.len)
    invisible()
}

バイプロットを使わない主成分分析の図示

バイプロットの赤の矢印が不要の場合は、次のようにして、主成分分析の結果を直接プロットします:

plot(prcomp(chu)$x[,1:2])

アスペクト比を 1 にする(両軸の目盛を合わせる)なら

plot(prcomp(chu)$x[,1:2], asp=1)

さらに点ではなく県名をプロットするなら

plot(prcomp(chu)$x[,1:2], asp=1, type="n")
text(prcomp(chu)$x[,1:2], row.names(chu))