因子分析     Last modified: Aug 31, 2009

目的

因子分析を行う。
R には factanal という関数が用意されている。最尤法により因子を求めプロマックス回転するならそれを利用すべきである。
factanal を使いやすくする関数(factanal2)を作ったので参照されたい。

使用法

pfa(dat)
pfa(dat, method=c("Varimax", "biQuartimax", "Quartimax", "Equimax", "None"), eps1=1e-5, eps2=1e-5, max1=999, max2=999, factors=0)
print(result, before=FALSE)
plot(result, before=FALSE, fac.no=1:2, scores=FALSE, xlab=NULL, ylab=NULL, axis=TRUE, label.cex=0.7, ...)

引数

dat        データ行列(行がケース,列が変数)。欠損値 NA を含むケース(行)は分析に先だって自動的に削除される。
           または,相関係数行列(正方行列であること)
method     因子軸の回転方法の指定(使用法参照 回転しないときは None)
eps1       共通性の反復推定の終了判定条件
eps2       因子軸の回転の終了判定条件
max1       共通性の反復推定の回数の上限
max2       因子軸の回転の回数の上限
factors    求める因子の数(指定しない場合には,固有値が 1 以上の数になる)

result     pfa が返すオブジェクト
before     因子軸の回転前の結果を使うか回転後の結果を使うかを指定
fac.no     横軸と縦軸にプロットする因子番号
scores     因子得点を描くか因子負荷量を描くかの指定
xlab       横軸の名前
ylab       縦軸の名前
axis       座標軸を描くかどうかの指定
label.cex  描画点につけるラベルの文字の大きさ
...        plot 関数に引き渡す引数

ソース

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

# 因子分析
pfa <- function(dat,                                                 # データ行列または相関係数行列
                method=c("Varimax", "Biquartimax", "Quartimax", "Equimax", "None"),     # 回転法
                eps1=1e-5,                                              # 共通性の収束限界
                eps2=1e-5,                                              # バリマックス回転の収束限界
                max1=999,                                               # 共通性の推定のための収束計算の上限回数
                max2=999,                                               # バリマックス回転を行う上限回数
                factors=0)
{
        method <- match.arg(method)                                  # 回転法を補完
        dat <- subset(dat, complete.cases(dat))                              # 欠損値を持つケースを除く
        nr <- nrow(dat)                                                      # ケース数
        nc <- ncol(dat)                                                      # 変数の個数
        if (is.null(colnames(dat))) {                                   # 変数名が無いときには仮の名前を付ける
                colnames(dat) <- paste("Var", 1:nc, sep=".")
        }
        vnames <- colnames(dat)                                              # 変数名を記録
        if (nr != nc && is.null(rownames(dat))) {                 # ケース名がないときには仮の名前を付ける
                rownames(dat) <- paste("Obs", 1:nr, sep=".")
        }
        cnames <- rownames(dat)                                              # ケース名を記録
        r0 <- r <- if (nr == nc) dat else cor(dat)                        # 与えられたのがデータ行列なら,相関係数行列を計算する
        communality0 <- 1-1/diag(solve(r))                           # 共通性の初期値(SMC)
        diag(r) <- communality0                                              # 相関係数行列の対角成分を共通性で置き換える
        result <- eigen(r)                                           # 固有値・固有ベクトルを求める
        eval <- result$values                                                # 固有値
        evec <- result$vectors                                               # 固有ベクトル
        if (factors == 0) {                                             # 因子数の指定がないときには,
                factors <- sum(eval >= 1)                            # 1 以上の固有値の数とする
        }
        converged <- FALSE                                           # 共通性の収束計算を行う
        for (i in 1:max1) {                                             # 上限回数まで繰り返す
                eval <- eval[1:factors]                                      # 必要な因子数まで,固有値を取る
                evec <- evec[,1:factors]                             # 必要な因子数まで
                evec <- t(sqrt(eval)*t(evec))                                # 因子負荷量を計算
                r <- r0                                                      # 相関係数行列を復帰
                communality <- rowSums(evec^2)                               # 共通性を計算
                if (all(abs(communality-communality0) < eps1)) {     # 共通性の変化が限界以内になれば,
                        converged <- TRUE                            # 収束したとして,
                        break                                           # ループを抜ける
                }
                communality0 <- communality                          # 現在の共通性を保存
                diag(r) <- communality                                       # 相関係数行列の対角成分を置き換える
                result <- eigen(r)                                   # 固有値・固有ベクトルを求める
                eval <- result$values                                        # 固有値
                evec <- result$vectors                                       # 固有ベクトル
        }
        if (converged == FALSE) {                                       # 収束フラッグが FALSE なら,
                warning("Not converged.")                               # 収束しなかったと警告する
        }
        else if (any(communality >= 1)) {                               # 1 を超える共通性があったら,
                warning("Communality >= 1.")                            # 警告する
        }

        if (factors == 1 || method == "None") {                         # 因子数が 1 であるか,因子軸の回転をしないのなら,
                w <- solve(r0)%*%evec                                        # 因子得点係数を計算する
                scores <- (scale(dat)*sqrt(nr/(nr-1)))%*%w           # 因子得点を計算する
                rownames(evec) <- names(communality) <- vnames            # 名前を付けて,結果を返す
                rownames(scores) <- cnames
                colnames(scores) <- colnames(evec) <- names(eval) <- paste("FAC", 1:factors, sep=".")
                return(structure(list(rotation=method, correlation.matrix=r0, communality=communality,
                        before.rotation.fl=evec, before.rotation.eval=eval, scores=scores), class="pfa"))
        }
        else {                                                          # 因子軸を回転するなら,
                fl <- evec/sqrt(communality)                         # 因子負荷量を共通性で基準化
                eig <- numeric(factors)
                ov <- 0
                wt <- switch (method,                                        # 回転方法によって,重みを変える
                        "Varimax" = 1,
                        "biQuartimax" = 0.5,
                        "Quartimax" = 0,
                        "Equimax"       = nc/2)
                fnp <- nc
                for (loop in 1:max2) {                                  # 回転の上限回数まで収束計算
                        for (k in 1:(factors-1)) {
                                for (k1 in (k+1):factors) {
                                        x <- fl[,k]
                                        y <- fl[,k1]
                                        xy <- x^2-y^2
                                        a <- sum(xy)
                                        b <- 2*sum(x*y)
                                        c <- sum(xy^2-4*x^2*y^2)
                                        d <- 4*sum(x*y*xy)
        
                                        dd <- d-2*wt*a*b/fnp
                                        theta <- atan(dd/(c-wt*(a^2-b^2)/fnp))
                                        if(sin(theta)*dd < 0) {
                                                if (theta > 0) {
                                                        theta <- theta-pi
                                                }
                                                else {
                                                        theta <- theta+pi
                                                }
                                        }
                                        theta <- theta/4
                                        cs <- cos(theta)
                                        si <- sin(theta)
                                        fljk <- fl[,k]
                                        fljk1 <- fl[,k1]
                                        fl[,k] <- fljk*cs+fljk1*si
                                        fl[,k1] <- fljk1*cs-fljk*si
                                }
                        }
                        v <- sum((t(fl)^2-colSums(fl^2)*wt/fnp)^2)
        
                        if (abs(v-ov) < eps2) {                              # 収束したら,
                                break                                   # ループから抜ける
                        }
                        ov <- v
                }
                fl <- fl*sqrt(communality)                           # 因子負荷量
                w <- solve(r0)%*%fl                                  # 因子得点係数
                scores <- (scale(dat)*sqrt(nr/(nr-1)))%*%w           # 因子得点
                eval2 <- colSums(fl^2)                                       # 因子負荷量の二乗和
                rownames(evec) <- rownames(fl) <- names(communality) <- vnames
                rownames(scores) <- cnames
                colnames(scores) <- colnames(evec) <- colnames(fl) <- names(eval) <- names(eval2) <- paste("FAC", 1:factors, sep=".")
                return(structure(list(rotation=method, correlation.matrix=r0, communality=communality,
                        before.rotation.fl=evec, before.rotation.eval=eval,
                        after.rotation.fl=fl, after.rotation.eval=eval2, scores=scores), class="pfa"))
        }
}
# print メソッド
print.pfa <- function(       result,                                         # pfa が返すオブジェクト
                        before=FALSE)
{
        communality <- result$communality
        vnames <- sapply(names(communality), function(i) substring(i, 1, min(nchar(i), 7)))
        if (before || is.null(result$after.rotation.fl)) {
                fl <- result$before.rotation.fl
                eval <- result$before.rotation.eval
                label <- "E-value"
                if (result$rotation == "None") {
                        printf("\nResult without rotation\n\n")
                }
                else {
                        printf("\nBefore %s rotation\n\n", result$rotation)
                }
        }
        else {
                fl <- result$after.rotation.fl
                eval <- result$after.rotation.eval
                label <- "Sq.Sum"
                printf("\nAfter %s rotation\n\n", result$rotation)
        }
        nv <- nrow(fl)
        nc <- ncol(fl)
        cat("       ", sprintf(" Fac.%02i", 1:nc), " Communality\n", sep="")
        for (i in 1:nv) {
                cat(sprintf("%-7s", vnames[i]), sprintf("%7.3f", fl[i,]), sprintf("%7.3f\n", communality[i]), sep="")
        }
        cat(sprintf("%-7s", label),   sprintf("%7.3f", eval), "\n", sep="")
        cat(sprintf("%-7s", "Cont."), sprintf("%7.1f", eval/nv*100), "\n", sep="")
        cat(sprintf("%-7s", "Cum."),  sprintf("%7.1f", cumsum(eval/nv*100)), "\n", sep="")
}
# plot メソッド
plot.pfa <- function(        result,                                         # pfa が返すオブジェクト
                        before=FALSE,                                   # 因子軸の回転前の結果を使うか回転後の結果を使うかを指定
                        fac.no=1:2,                                     # 横軸と縦軸にプロットする因子番号
                        scores=FALSE,                                   # 因子得点を描くか因子負荷量を描くかの指定
                        xlab=NULL, ylab=NULL,                           # 横軸,縦軸の名前
                        axis=TRUE,                                      # 座標軸を描くかどうかの指定
                        label.cex=0.7,                                  # 描画点につけるラベルの文字の大きさ
                        ...)                                            # plot 関数に引き渡す引数
{
        fac.name <- names(result$before.rotation.eval)
        if (length(fac.name) > 1) {
                ax1 <- fac.no[1]
                ax2 <- fac.no[2]
                if (is.null(xlab)) {
                        xlab <- fac.name[ax1]
                }
                if (is.null(ylab)) {
                        ylab <- fac.name[ax2]
                }
                if (scores) {
                        x <- result$scores[, ax1]
                        y <- result$scores[, ax2]
                        labels <- 1:length(x)
                }
                else {
                        if (before || is.na(result$after.rotation.fl)) {
                                fl <- result$before.rotation.fl
                        }
                        else {
                                fl <- result$after.rotation.fl
                        }
                        x <- fl[, ax1]
                        y <- fl[, ax2]
                        labels <- names(result$communality)
                }
                plot(x, y, xlab=xlab, ylab=ylab, ...)
                old <- par(xpd=TRUE)
                text(x, y, labels, cex=label.cex, pos=1)
                par(old)
        }
        if (axis) {
                abline(h=0, v=0)
        }
}


使用例 1

> dat <- read.table("pfa.data", header=TRUE)	# ファイルからデータを読んだ
> dat
   var1 var2 var3 var4 var5 var6 var7 var8 var9 var10
1   935  955  926  585 1010  925 1028  807  769   767
2   817  905  901  632 1004  950  957  844  781   738
3   768  825  859  662  893  900  981  759  868   732
4   869  915  856  448  867  874  884  802  804   857
5   787  878  880  592  871  874  884  781  782   807
6   738  848  850  569  814  950  957  700  870   764
7   763  862  839  658  887  900 1005  604  709   753
8   795  890  841  670  853  874  859  701  680   772
9   903  877  919  460  818  849  884  700  718   716
10  761  765  881  485  846  900  981  728  781   714
11  747  792  800  564  796  849  932  682  746   767
12  771  802  840  609  824  874  859  668  704   710
13  785  878  805  527  911  680  884  728  709   747
14  657  773  820  612  810  849  909  698  746   771
15  696  785  791  578  774  725  932  765  706   795
16  724  785  870  509  746  849  807  763  724   760
17  712  829  838  516  875  725  807  754  762   585
18  756  863  815  474  873  725  957  624  655   620
19  622  759  786  619  820  769  807  673  698   695
20  668  753  751  551  834  849  807  601  655   642

> ( result <- pfa(dat) )	# 因子分析を実行

After Varimax rotation

        Fac.01 Fac.02 Fac.03 Communality
var1    -0.818  0.381 -0.254  0.878
var2    -0.867  0.226 -0.054  0.806
var3    -0.581  0.619 -0.133  0.738
var4     0.041  0.070  0.730  0.539
var5    -0.834  0.026  0.315  0.796
var6    -0.187  0.631  0.371  0.571
var7    -0.412  0.339  0.317  0.385
var8    -0.349  0.560 -0.149  0.457
var9    -0.118  0.734  0.125  0.568
var10   -0.086  0.562  0.055  0.326
Sq.Sum   2.803  2.265  0.997
Cont.     28.0   22.7   10.0
Cum.      28.0   50.7   60.7

その他の結果は,str(result) のようにして要素を確認後,
個別の要素を print 関数によって出力する

> print(result$correlation, digits=3)	# 相関係数行列

        var1    var2   var3    var4    var5  var6  var7    var8   var9   var10
var1   1.000  0.8449  0.775 -0.2257  0.5586 0.375 0.429  0.4447 0.2466  0.3362
var2   0.845  1.0000  0.614 -0.0212  0.7167 0.239 0.390  0.4467 0.2244  0.3136
var3   0.775  0.6141  1.000 -0.1006  0.4470 0.563 0.403  0.5982 0.4815  0.2651
var4  -0.226 -0.0212 -0.101  1.0000  0.2136 0.331 0.189 -0.0544 0.0718  0.1688
var5   0.559  0.7167  0.447  0.2136  1.0000 0.249 0.483  0.3976 0.2189 -0.0407
var6   0.375  0.2391  0.563  0.3311  0.2486 1.000 0.419  0.1989 0.5319  0.3893
var7   0.429  0.3897  0.403  0.1888  0.4827 0.419 1.000  0.1419 0.4076  0.2742
var8   0.445  0.4467  0.598 -0.0544  0.3976 0.199 0.142  1.0000 0.5689  0.4473
var9   0.247  0.2244  0.481  0.0718  0.2189 0.532 0.408  0.5689 1.0000  0.3638
var10  0.336  0.3136  0.265  0.1688 -0.0407 0.389 0.274  0.4473 0.3638  1.0000

あるいは,print.default(result) のようにすればすべての結果を出力することができる。

> print.default(result)

$rotation	# 因子軸の回転法

[1] "Varimax"

$correlation.matrix	# 相関係数行列

            var1        var2       var3        var4       var5      var6      var7        var8
var1   1.0000000  0.84487442  0.7746348 -0.22572126  0.5585596 0.3746964 0.4288587  0.44474253
var2   0.8448744  1.00000000  0.6140948 -0.02116982  0.7166837 0.2391217 0.3897472  0.44666435
var3   0.7746348  0.61409482  1.0000000 -0.10056276  0.4470492 0.5628616 0.4030052  0.59819255
var4  -0.2257213 -0.02116982 -0.1005628  1.00000000  0.2136123 0.3311048 0.1888090 -0.05438703
var5   0.5585596  0.71668368  0.4470492  0.21361227  1.0000000 0.2485824 0.4826623  0.39756234
var6   0.3746964  0.23912166  0.5628616  0.33110479  0.2485824 1.0000000 0.4191012  0.19885885
var7   0.4288587  0.38974717  0.4030052  0.18880904  0.4826623 0.4191012 1.0000000  0.14187410
var8   0.4447425  0.44666435  0.5981926 -0.05438703  0.3975623 0.1988589 0.1418741  1.00000000
var9   0.2466372  0.22442056  0.4814589  0.07176501  0.2189180 0.5319358 0.4075537  0.56891499
var10  0.3362334  0.31363051  0.2651037  0.16878738 -0.0407357 0.3893330 0.2741899  0.44734094
            var9      var10
var1  0.24663719  0.3362334
var2  0.22442056  0.3136305
var3  0.48145888  0.2651037
var4  0.07176501  0.1687874
var5  0.21891798 -0.0407357
var6  0.53193583  0.3893330
var7  0.40755369  0.2741899
var8  0.56891499  0.4473409
var9  1.00000000  0.3637659
var10 0.36376593  1.0000000

$communality	# 共通性

     var1      var2      var3      var4      var5      var6      var7      var8      var9     var10 
0.8782869 0.8058379 0.7376875 0.5392022 0.7962812 0.5710556 0.3852000 0.4574828 0.5682011 0.3258559 

$before.rotation.fl	# 回転前の因子負荷量

            FAC.1       FAC.2       FAC.3
var1  -0.85127134 -0.38093708 -0.09225476
var2  -0.80339537 -0.36668776  0.16104011
var3  -0.83136537 -0.04526261 -0.21088020
var4  -0.06449186  0.52620168  0.50808939
var5  -0.67637810 -0.22526754  0.53670139
var6  -0.57352146  0.49150559  0.02347232
var7  -0.55425633  0.17470446  0.21789501
var8  -0.61545313  0.03446450 -0.27841064
var9  -0.56939070  0.42882823 -0.24515650
var10 -0.42989384  0.30434024 -0.22005500

$before.rotation.eval	# 回転前の固有値

    FAC.1     FAC.2     FAC.3 
4.0468459 1.1590592 0.8591861 

$after.rotation.fl	# 回転後の因子負荷量

            FAC.1      FAC.2       FAC.3
var1  -0.81752240 0.38107025 -0.25441993
var2  -0.86716928 0.22560900 -0.05436886
var3  -0.58096539 0.61850727 -0.13272345
var4   0.04083323 0.06991167  0.72982686
var5  -0.83442587 0.02649545  0.31513908
var6  -0.18699075 0.63098561  0.37141241
var7  -0.41170256 0.33927289  0.31716695
var8  -0.34888652 0.56003956 -0.14871696
var9  -0.11753949 0.73398757  0.12509133
var10 -0.08607321 0.56160143  0.05523714

$after.rotation.eval	# 回転後の因子負荷量の二乗和

    FAC.1     FAC.2     FAC.3 
2.8031906 2.2652820 0.9966186 

$scores	# 因子得点

         FAC.1       FAC.2        FAC.3
1  -2.35863255  0.63067260  0.329220617
2  -1.27981713  0.58243865  1.174903000
3  -0.06492078  1.18107012  0.889716371
4  -0.77918746  0.87233382 -1.166728853
5  -0.34036193  0.73957325  0.183002359
6   0.74062447  1.37316866  0.447702634
7  -0.46817258 -0.47855965  1.607620790
8  -0.30749799 -0.17331089  0.370604867
9  -0.83009214  0.44493091 -1.941757717
10  0.41762185  0.65611503 -0.170693247
11  0.79272411  0.15871956 -0.164168839
12  0.22027197 -0.14173994 -0.006421007
13 -1.08753778 -1.29541749 -0.451412084
14  1.16049740  0.05318377  0.933139700
15  1.31738669 -0.05099253 -0.583232459
16  1.33358487  0.69466121 -1.145936980
17  0.03587380 -0.72694479 -0.881790129
18 -0.43093588 -1.74621258 -0.621098893
19  1.07368915 -1.08529807  0.731121280
20  0.85488189 -1.68839164  0.466208591

attr(,"class")
[1] "pfa"

plot メソッド

> plot(result, axis=FALSE, pch=19, main="因子負荷量")	# 因子負荷量のプロット

graph

> plot(a, pch=19, scores=TRUE)	# 因子得点のプロット

graph

・ 解説ページ


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

Made with Macintosh