因子分析(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

Made with Macintosh