目的

因子分析の結果（因子数，独自性）に対して適合度検定を行う
（最尤法でなければ適合度検定ができないと思っていたが，
必要なのは抽出した因子数と分析結果のうちの独自性ベクトル
だけであるらしい。
ならば，どんな因子分析の結果についてでも適合度の検定を
行うことができる。）

使用法

fa.fit.test <- function(data, factors, uniquenesses)


引数

data           データ行列（データフレーム）
factors        抽出した因子数
uniquenesses   独自性ベクトル

ソース

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

# 因子分析の適合度検
fa.fit.test <- function(data,                                        # データ行列（データフレーム）
                        factors,                                # 抽出した因子数
                        uniquenesses)                           # 独自性ベクトル
{
        p <- ncol(data)                                              # 変数の個数
        n <- nrow(data)                                              # ケース数
        r <- cor(data)                                               # 相関係数行列
        sc <- diag(1/sqrt(uniquenesses))                     # 
        r <- sc%*%r%*%sc                                     # 
        e <- eigen(r, symmetric=TRUE, only.values=TRUE)              # 固有値だけ求める
        e <- e$values[-(1:factors)]                          # 抽出された因子で説明されない部分
        s <- -sum(log(e)-e)+factors-p                                # 統計量
        chisq <- (n-1-(2*p+5)/6-(2*factors)/3)*s             # カイ二乗統計量に変換
        df <- ((p-factors)^2-p-factors)/2                    # 自由度
        P <- pchisq(chisq, df, lower.tail=FALSE)             # P 値
        return(c(Chi.sq.=chisq, df=df, P=P))                    # 名前付ベクトルで返す
}


使用例

> 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
> ans <- factanal(dat, 4)	# R の関数を使って分析（デフォルトはバリマックス回転）
> ans

Call:
factanal(x = dat, factors = 4)

Uniquenesses:	# 必要な情報の一つ $uniquenesses に入っている
 var1  var2  var3  var4  var5  var6  var7  var8  var9 var10
0.005 0.187 0.236 0.595 0.005 0.005 0.644 0.005 0.467 0.574

Loadings:
      Factor1 Factor2 Factor3 Factor4
var1   0.820   0.200   0.347  -0.402
var2   0.856   0.204          -0.178
var3   0.559   0.454   0.472  -0.151
var4                   0.112   0.627
var5   0.891   0.174  -0.151   0.385
var6   0.231   0.169   0.877   0.379
var7   0.500           0.263   0.190
var8   0.284   0.955
var9   0.105   0.580   0.376   0.211
var10          0.452   0.432  -0.188

               Factor1 Factor2 Factor3 Factor4
SS loadings      2.907   1.799   1.552   1.021
Proportion Var   0.291   0.180   0.155   0.102
Cumulative Var   0.291   0.471   0.626   0.728

Test of the hypothesis that 4 factors are sufficient.
The chi square statistic is 15.16 on 11 degrees of freedom.
The p-value is 0.175
> fa.fit.test(dat, 4, ans$uniquenesses)	# 適合度の検定だけを行う
   Chi.sq.         df          P	# factanal 関数の出力と同じである
15.1635725 11.0000000  0.1751303