主成分分析(princomp を援用)     Last modified: Jul 29, 2004

目的

R には,princomp および prcomp という,二種類の関数が用意されている。
しかし,これらが返す「loadings」は固有ベクトルそのものであって,いわゆる負荷量ではない。
ここに示す関数は,princomp を下請け関数として用いて,通常の主成分分析結果として提示するものである。

なお,pca という関数も書いたので,そちらも参照してみるとよい。

使用法

princomp2(dat, pcs=0, cor=TRUE, scores=FALSE, verbose=TRUE)

引数

dat	   データ行列(行がケース,列が変数)
pcs     求める主成分の個数。
        0 を指定する(デフォールト)と固有値が 1 以上の主成分を求める。
cor     TRUE(デフォールト)の場合には相関係数行列を主成分分析する(変数の単位に異なるものがある場合にはこれを選択する以外の手はない)。
        FALSE の場合には,分散・共分散行列を主成分分析する。
scores  TRUEなら主成分得点を計算する。
vervose TRUE(デフォールト) なら,分析結果を print 文によって出力する
        FALSEの場合には,結果はリストによって返される。

戻り値

戻り値のリスト

princomp が返すもの
$sdev       固有値
$loadings   固有ベクトル(いわゆる負荷量ではない)
$center     分析に使用した各変数の平均値
$scale      分析に使用した各変数の標準偏差(変動をデータの個数で割って平方根を取った方)
$n.obs      分析に使用したケース数
$scores     主成分得点(分散が固有値になるように調整されたもの)

付加されたもの
$loadings  整形された主成分分析結果(下記参照)
$x         主成分得点(不偏分散が固有値になるように調整されたもの)

ソース

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

# princomp を援用して,主成分分析の結果をまとめる
princomp2 <- function(       dat,                                                    # データ行列
                        pcs=0,                                                  # 求める主成分の個数
                        cor=TRUE,                                               # TRUE の場合には相関係数行列,FALSE なら分散共分散行列を対象にする
                        scores=FALSE,                                           # 主成分得点を計算するかどうか
                        verbose=TRUE)                                           # 饒舌に結果を表示するかどうか
{
        p <- ncol(dat)                                                               # 変数の個数
        if (is.null(colnames(dat))) colnames(dat) <- paste("Var", 1:p, sep=".")      # 変数名がついていないとき
        n <- nrow(dat)                                                               # サンプルサイズ
        if (is.null(rownames(dat))) rownames(dat) <- paste("Obs", 1:n, sep=".")      # 行名がついていないとき
        dat <- subset(dat, complete.cases(dat))                                      # 欠損値を持つケースを除く
        n <- nrow(dat)                                                               # サンプルサイズは小さくなったかもしれない
        result<-princomp(dat, cor=cor, scores=scores)                                # princomp を呼び出す
        if (pcs == 0) {                                                         # 抽出する主成分数が指定されていないときは,
                pcs <- sum(result$sdev^2 >= 1)                                       # 1 以上の固有値の数
        }
        loadings <- t(result$sde*t(result$loadings))[, 1:pcs, drop=FALSE]    # 主成分負荷量
        Contribution <- rowSums(loadings^2)                                  # 寄与率
        Eigen.values <- c(result$sdev[1:pcs]^2, sum(result$sdev[1:pcs]^2))   # 固有値
        denominator <- if (cor) p else sum(diag(var(dat)*(n-1)/n))           # 累積寄与率を計算する分母(cor=FALSE なら,分散共分散行列の対角成分の和)
        Proportion <- Eigen.values/denominator*100                           # 寄与率
        Cumulative.prop. <- cumsum(Proportion)                                       # 累積寄与率
        Cumulative.prop.[pcs+1] <- NA                                                # 末尾は不要
        result$loadings <- rbind(cbind(loadings, Contribution),                      # 主要分析結果
                                Eigen.values, Proportion, Cumulative.prop.)
        result$x <- result$scores * sqrt((n-1)/n)                            # 主成分得点
        if (verbose) {
                cat(sprintf("%16s", ""), sprintf("%8s", paste("PC", 1:pcs, sep="")), " Contribution\n", sep="")
                for (i in 1:p) {
                        cat(sprintf("%-16s", rownames(result$loadings)[i]), sprintf("%8.3f", loadings[i,]), sprintf("%10.3f\n", Contribution[i]), sep="")
                }
                cat("Eigen.values    ", sprintf("%8.3f", result$loadings[p+1, 1:pcs]), "\n", sep="")
                cat("Proportion      ", sprintf("%8.3f", result$loadings[p+2, 1:pcs]), "\n", sep="")
                cat("Cumulative.prop.", sprintf("%8.3f", result$loadings[p+3, 1:pcs]), "\n", sep="")
                if (scores) {
                        cat(sprintf("\n%16s", ""), sprintf("%8s", paste("PC", 1:pcs, sep="")), "\n", sep="")
                        for (i in 1:n) {
                                cat(sprintf("%-16s", rownames(result$scores)[i]), sprintf("%8.3f", result$scores[i, 1:pcs]), "\n", sep="")
                        }
                }
        }
        invisible(result)
}


使用例

dat <- matrix(c(	# 10 ケース,5 変数のデータ行列例(ファイルから読んでも良い)
	-1.89, -0.02, 0.42, 1.23, -1.53,
	0.06, 1.81, -0.59, -0.75, -0.12,
	2.58, -0.20, -1.92, -0.49, -0.35,
	0.69, -0.66, -0.77, -1.92, 0.38,
	-1.05, 0.07, -0.37, -0.89, -1.62,
	-2.73, 1.40, 0.18, -0.09, 0.13,
	0.95, 0.84, 1.19, 1.19, 0.10,
	0.93, 1.17, -1.70, 1.46, -0.25,
	-0.04, -0.12, -0.34, -0.24, 1.88,
	0.16, 1.03, -0.05, -0.73, 0.04
), byrow=TRUE, ncol=5)

colnames(dat) <- paste("Var", 1:5, sep=".")     # 変数名
rownames(dat) <- paste("case", 1:10, sep="-")   # ケース名

princomp2(dat)

出力結果例

例 1 固有値が 1 以上の主成分を抽出し,結果を自動出力
> princomp2(dat)
                    PC1    PC2    PC3 Contribution	因子負荷量,寄与率など
Var.1            -0.767  0.448 -0.151     0.812
Var.2             0.530  0.618  0.144     0.684
Var.3             0.750 -0.249  0.387     0.774
Var.4             0.557  0.557 -0.324     0.726
Var.5            -0.357  0.302  0.848     0.937
Eigen.values      1.869  1.047  1.017
Proportion       37.375 20.932 20.342
Cumulative.prop. 37.375 58.306 78.648

例 2 2 つの主成分を抽出し,主成分得点とともに結果を自動出力
> princomp2(dat, pcs=2, scores=TRUE)
                    PC1    PC2 Contribution
Var.1            -0.767  0.448     0.789
Var.2             0.530  0.618     0.663
Var.3             0.750 -0.249     0.624
Var.4             0.557  0.557     0.621
Var.5            -0.357  0.302     0.218
Eigen.values      1.869  1.047
Proportion       37.375 20.932
Cumulative.prop. 37.375 58.306

               PC1     PC2	主成分得点
case-1       1.868  -0.946
case-2       0.231   0.751
case-3      -2.405   0.379
case-4      -1.956  -1.385
case-5       0.293  -1.545
case-6       1.771  -0.196
case-7       1.205   0.863
case-8      -0.209   1.940
case-9      -0.898   0.052
case-10      0.100   0.087

例 3 例 1 と同じ出力をする別法
> ans <- princomp2(dat, pcs=0, scores=TRUE)
> ans$loadings
                        PC1        PC2        PC3 Contribution
Var.1             0.7670549  0.4479912 -0.1511371    0.8119117
Var.2            -0.5301729  0.6180870  0.1441703    0.6838998
Var.3            -0.7496095 -0.2494103  0.3872449    0.7740785
Var.4            -0.5570003  0.5573622 -0.3239124    0.7258212
Var.5             0.3565268  0.3016596  0.8477030    0.9367101
Eigen.values      1.8687315  1.0465842  1.0171057    3.9324214
Proportion       37.3746303 20.9316843 20.3421134   78.6484280
Cumulative.prop. 37.3746303 58.3063146 78.6484280           NA

> ans$scores
               PC1         PC2        PC3         PC4        PC5
case-1   1.8684045 -0.94576668 -1.2332697 -0.93926955 -0.3730075
case-2   0.2309233  0.75142992  0.3473489  1.47623977  0.4528447
case-3  -2.4045883  0.37897746 -1.1474093 -0.15470182  0.1186054
case-4  -1.9560246 -1.38510647  0.5622894 -0.06866507  0.1175167
case-5   0.2926468 -1.54510805 -1.0717479  0.52700054  0.1763542
case-6   1.7712949 -0.19635831  0.9154032  0.86695787 -0.6625894
case-7   1.2045357  0.86272982  0.4474658 -1.37756199  1.0911237
case-8  -0.2093047  1.94011453 -1.1422799  0.07247913 -0.7135016
case-9  -0.8983546  0.05163475  1.7564121 -0.93389874 -0.7875738
case-10  0.1004669  0.08745302  0.5657873  0.53141986  0.5802275

・ 解説ページ


・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI

Made with Macintosh