主成分分析（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
