因子分析（factanal を援用する） Last modified: Jan 13, 2007
目的
最尤法により因子を求める。因子軸の回転はプロマックス回転，バリマックス回転から選べる。
別に書いた関数pfaは主因子解によるものである。
使用法
factanal2(dat, factors=0, rotation=c("promax", "varimax", "none"),
scores=c("none", "regression", "Bartlett"), verbose=TRUE)
引数
dat データ行列（行がケース，列が変数）。欠損値 NA を含むケース（行）は分析に先だって自動的に削除される。
factors 求める因子の数，省略すると理論的に最も多い因子数が仮定される。
rotation 因子軸の回転方法の指定（使用法参照 回転しないときは none）
scores 因子得点の求め方の指定（求めないときは none）
verbose お手軽モードのときには TRUE（デフォールト），結果を付値して使うときには FALSE
戻り値
戻り値のリスト
factanal が返すもの
$converged 収束したかどうか
$loadings 因子負荷量
$uniquenesses 独自性
$correlation 相関係数
$criteria 最適化の結果
$factors 求めた因子数
$method 因子を求める方法
$scores 因子得点
$STATISTIC 因子数が十分か検定するカイ二乗値
$dof その自由度
$PVAL その P 値
$n.obs ケース数
$call
付加されるもの
$cosmetic 整形した結果（後述）
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/factanal2.R", encoding="euc-jp")
# factanal 関数の結果を整形して表示する
factanal2 <- function( dat, # データ行列
factors=0, # 抽出する因子数
rotation=c("promax", "varimax", "none"), # 因子軸の回転法
scores=c("none", "regression", "Bartlett"), # 因子得点の算出法
verbose=TRUE) # 結果の画面表示をするかどうか
{
p <- ncol(dat) # 変数の個数
n <- nrow(dat) # 行数（欠損値を含むケースも含む）
if (is.null(colnames(dat))) colnames(dat) <- paste("Var", 1:p, sep=".") # 変数名がないときにデフォルト名をつける
if (is.null(rownames(dat))) rownames(dat) <- paste("Case", 1:n, sep="-")# 行名がないときにデフォルト名をつける
dat <- subset(dat, complete.cases(dat)) # 欠損値を持つケースを除く
rotation <- match.arg(rotation) # 引数の補完
scores <- match.arg(scores) # 引数の補完
if (factors == 0) { # 抽出因子数が指定されないときは，
factors <- max(1, floor((2*p+1-sqrt((2*p+1)^2-4*(p^2-p)))/2)) # デフォルトの因子数
}
txt <- sprintf('factanal(dat, factors=factors, rotation="%s", scores=scores)', rotation) # rotation に実際の文字列を渡すためにこのようにする
result <- eval(parse(text=txt)) # 関数呼び出し
Communality <- 1-result$uniquenesses # 共通性は，斜交回転のときには因子負荷量の二乗和ではない
result$cosmetic <- cbind(result$loadings, Communality) # 共通性を付加
if (rotation!="promax") { # 斜交回転でない場合には，
SS.loadings <- colSums(result$loadings^2) # 因子負荷量の二乗和
SS.loadings <- c(SS.loadings, sum(SS.loadings)) # 総和を加える
Proportion <- SS.loadings/p*100 # 寄与率
Cum.Prop. <- cumsum(Proportion) # 累積寄与率
Cum.Prop.[factors+1] <- NA
result$cosmetic <- rbind(result$cosmetic, SS.loadings, Proportion, Cum.Prop.)
}
if (verbose == TRUE) { # 画面表示をするとき
if (result$dof) { # モデル適合度の自由度が 0 でないとき
test <- data.frame(result$STATISTIC, result$dof, result$PVAL)
colnames(test) <- c("Chi sq.", "d.f.", "P value")
rownames(test) <- ""
cat(sprintf("H0: %i factors are sufficient.\n", factors))
print(test)
}
else { # 自由度が 0 になるとき
cat(sprintf("The degrees of freedom for the model is 0 and the fit was %g\n", result$criteria[1]))
}
cat(sprintf("\nFactor loadings(rotation:%s)\n", rotation)) # 因子負荷量
print(result$cosmetic)
if (scores != "none") {
cat(sprintf("\nFactor scores(%s)\n", scores)) # 因子得点
print(result$scores)
}
}
invisible(result) # 明示的に print しないと，何も出力しない
}
使用例
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)
factanal2(dat)
出力結果例
例 1 # 一番お手軽な利用法
> factanal2(dat)
H0: 2 factors are sufficient.
Chi sq. d.f. P value
0.2545753 1 0.6138718
Factor loadings(rotation:promax)
Factor1 Factor2 Communality
Var.1 1.0431755 0.1481079 0.99500000
Var.2 -0.1574652 0.2685831 0.12845163
Var.3 -0.5203595 0.1265378 0.33586077
Var.4 0.1579708 1.0455338 0.99500000
Var.5 0.2153349 -0.1386064 0.08782451
例 2 # 因子数を指定するとき
> factanal2(dat, factors=1)
H0: 1 factors are sufficient.
Chi sq. d.f. P value
1.247915 5 0.9401988
Factor loadings(rotation:promax)
Factor1 Communality
Var.1 -0.7935450 0.62971284
Var.2 0.2843882 0.08087757
Var.3 0.6762596 0.45732658
Var.4 0.2414748 0.05831380
Var.5 -0.2248329 0.05055701
例 3 # バリマックス回転をするとき
> factanal2(dat, rotation="varimax")
H0: 2 factors are sufficient.
Chi sq. d.f. P value
0.2545753 1 0.6138718
Factor loadings(rotation:varimax)
Factor1 Factor2 Communality
Var.1 0.99724059 -0.02264989 0.99500000
Var.2 -0.21428550 0.28728743 0.12845163
Var.3 -0.54136393 0.20684455 0.33586077
Var.4 -0.07317994 0.99480924 0.99500000
Var.5 0.24292539 -0.16974282 0.08782451
SS.loadings 1.39785002 1.14428981 2.54213983
Proportion 27.95700049 22.88579617 50.84279666
Cum.Prop. 27.95700049 50.84279666 NA
例 4 # バートレットの方法による因子得点を求めるとき
> factanal2(dat, scores="Bartlett")
H0: 2 factors are sufficient.
Chi sq. d.f. P value
0.2545753 1 0.6138718
Factor loadings(rotation:promax)
Factor1 Factor2 Communality
Var.1 1.0431755 0.1481079 0.99500000
Var.2 -0.1574652 0.2685831 0.12845163
Var.3 -0.5203595 0.1265378 0.33586077
Var.4 0.1579708 1.0455338 0.99500000
Var.5 0.2153349 -0.1386064 0.08782451
Factor scores(Bartlett)
Factor1 Factor2
Case-1 -1.36002917 1.3862444
Case-2 0.13754415 -0.5658090
Case-3 1.72347381 -0.5820791
Case-4 0.69223394 -1.6756576
Case-5 -0.55423895 -0.5860368
Case-6 -1.72993900 0.2902698
Case-7 0.45674537 1.0805067
Case-8 0.42174371 1.3180138
Case-9 0.01444485 -0.1072308
Case-10 0.19802129 -0.5582213
解説ページ
直前のページへ戻る E-mail to Shigenobu AOKI
