図1 第1因子の因子負荷量のヒストグラム
<< 目次へ
最終更新日： 2009.7.30
Bootstarp法による因子負荷量の推定
Davison, A. C. and Hinkley, D. V. (1997). Bootstrap Methods and Their Application. Cambridge University Press.には回帰分析における回帰係数を推定する方法として、観測値(xi, yi)をリサンプリングすることによってブートストラップ標本を作成する観測値ベースの標本抽出(case-based resampling)を紹介している。回帰分析の場合での例示はシミュレーション入門のページで紹介しているので、そちらを参照していただきたい。ここでは因子分析における因子負荷量の推定も、これと同じ方法によって実行できるかどうかを検証してみる。
まずは因子分析モデルの応用のページで用いたデータセットを使い、
という手順でシミュレーションを行うことにする。もし回帰分析と同様に利用できるのであれば、100個の因子負荷量のヒストグラムは正規分布の形状を描くはずであり、その中心に位置する値は、通常の因子分析によって得られる因子負荷量と同等な値になることが期待される。また、これが有効であるならば、推定された因子負荷量の信頼区間を求めることができる。
■ 通常の因子分析の結果
Loadings: # 因子負荷量
Factor1 Factor2
国語 0.470
地理 -0.116 0.614
英語 0.555
数学 0.523 0.104
化学 0.652 -0.209
物理 0.536 0.128
■ ブートストラップ法によるシミュレーションを行う関数を定義する
以下に関数boot.factanal()を示す。引数datには分析するデータセットを、repeat.numはリサンプリングして因子分析を行う回数、s.sizeにはブートストラップ標本の大きさを指定する。なお、この関数は今回のためだけに定義されたもので、因子が2つであると仮定されている場合にのみ適用可能である。例えば、因子数を3つにするならば、関数内にあるfactanal()の引数でfactor=3と指定し、「因子負荷量を取り出す」という部分と「i行目に因子負荷量を代入していく」という部分を書き換える必要がある。
# datはデータ行列
# repeat.numは因子分析を行う数
# s.sizeはリサンプリングするときのサンプルサイズ
boot.factanal <- function(dat, repeat.num, s.size){
# 引数に指定された行列の行数と列数の取得
n <- nrow(dat)
k <- ncol(dat)
# n行k列の空のマトリックスを作成
f1.mat <- matrix(numeric(repeat.num*k), nrow=repeat.num, ncol=k)
f2.mat <- matrix(numeric(repeat.num*k), nrow=repeat.num, ncol=k)
# サンプリングするときに使う
x <- 1:n
for(i in 1:repeat.num){
# datからs.size個のデータをリサンプリング
ana.dat <- dat[sample(x, s.size, replace=TRUE), ]
# ana.datに対して因子分析を行う
res <- factanal(ana.dat, factors=2, rotation="promax", scores="regression")
# 因子負荷量を取り出す
f1 <- res$loadings[ ,1]
f2 <- res$loadings[ ,2]
# i行目に因子負荷量を代入していく
f1.mat[i, ] <- f1
f2.mat[i, ] <- f2
}
list(factor1 = f1.mat, factor2 = f2.mat)
}
■ 実行例
以下に因子分析を行う数を10000回（1万回）、ブートストラップ標本の大きさを100（元の標本と同じ大きさ）として実行した例を示す。図は左上から行方向へ順に国語、地理、英語、数学、化学、物理という並びになっている。また、この順に因子負荷量の95%信頼区間を示す。ちなみに、10000回の繰り返し計算は私のパソコン（環境）で約10分かかった。
図1 第1因子の因子負荷量のヒストグラム
> for(i in 1:6) print(quantile(result$factor1[,i], c(0.025, 0.975)))
2.5% 97.5%
-0.2038923 1.0045372
2.5% 97.5%
-0.2791156 1.0409719
2.5% 97.5%
-0.1535238 1.0023570
2.5% 97.5%
-0.1643438 1.0217891
2.5% 97.5%
-0.3297843 1.0090630
2.5% 97.5%
-0.1362845 1.0291706
図2 第2因子の因子負荷量の因子負荷量
> for(i in 1:6) print(quantile(result$factor2[,i], c(0.025, 0.975)))
2.5% 97.5%
-0.1780796 0.9988578
2.5% 97.5%
-0.2763899 1.0276176
2.5% 97.5%
-0.1420974 0.9969547
2.5% 97.5%
-0.1732718 0.8949310
2.5% 97.5%
-0.3731202 0.9778171
2.5% 97.5%
-0.1502557 0.8606970
■ 結果の解釈
まずfactor1とfactor2のいずれについても、ヒストグラムの右裾を見てみると1.0を超える（不適解が生じている）ケースがある。一般に最尤法は不適解を生じやすいという欠点が指摘されているが、今回の結果もその危険性を表しているのかもしれない。また、信頼区間をみると全ての因子負荷量について、95%の信頼区間に0が含まれていることが分かる。これは5%の有意水準で有意ではないことを意味している。すなわち、推定された因子負荷量は0であるという帰無仮説を採択するということになる。さらに、ヒストグラムの形状が正規分布の形とは若干ずれていることが見てとれる。10000回（1万回）のリサンプリングを行っているが、うまく中心極限定理が生まれていないようである。
以上のことを総合して考えると、どうやら因子負荷量の推定に用いる信頼区間を求める場合、ブートストラップ法を適用するのは適切ではないようである。しかし、もしかしたら、そもそも通常の因子分析によって得られた結果が偶然的に得られた結果なのであって、ブートストラップ法による信頼区間の結果の方が正しいのかもしれない。これについて、ヒストグラムを見ると0.5付近の値が推定される頻度は決して起こりえないほどではないことがわかる。というよりも、0から0.5付近の頻度が高くなっている特長が見られるので、確率的には0から0.5までの値が推定される確率がたかいといえるだろう。この現象は、因子分析は得られる標本の変動によって推定される因子負荷量も影響されるということを示しているのだと考えられる。
不明な点は(1)回帰分析の回帰係数を推定する場合に適用できるこの方法が、因子分析の因子付加量の推定に適用するのが好ましくないという根拠が明らかではない、(2)因子負荷量の推定法として最尤法を用いたが、分布が正規分布でない場合に頑健性の高い重み付け最小自乗法を用いた場合にどのような結果が得られるか、(3)通常の因子分析による結果が正しいのか、それともブートストラップ法による信頼区間から得られる知見が正しいのか判断することができない、ということである。
■ 違うデータセットでの解析
ここで1因子構造を持つようなデータセットを用いて、同じことをやってみる（使用するデータ）。
> result <- factanal(dat, factor=1, rotation="promax", score="regression")
> result
Loadings:
Factor1
[1,] 0.391
[2,] 0.465
[3,] 0.442
[4,] 0.593
[5,] 0.639
boot.factanal()を次のように書き換えて使用する。
boot.factanal <- function(dat, repeat.num, s.size){
# 引数に指定された行列の行数と列数の取得
n <- nrow(dat)
k <- ncol(dat)
# n行k列の空のマトリックスを作成
f1.mat <- matrix(numeric(repeat.num*k), nrow=repeat.num, ncol=k)
# サンプリングするときに使う
x <- 1:n
for(i in 1:repeat.num){
# datからs.size個のデータをリサンプリング
ana.dat <- dat[sample(x, s.size, replace=TRUE), ]
# ana.datに対して因子分析を行う
res <- factanal(ana.dat, factors=1, rotation="promax", scores="regression")
# 因子負荷量を取り出す
f1 <- res$loadings[ ,1]
# i行目に因子負荷量を代入していく
f1.mat[i, ] <- f1
}
list(factor1 = f1.mat)
}
今回は繰り返し数を1000回、サンプルサイズを154（元の標本と同じ大きさ）の場合と、繰り返し数を10000回、サンプルサイズを154とした場合の2種類のやり方で得られる結果を比較してみる。
図3 繰り返しすうが1000回の場合のヒストグラム
> result2 <- boot.factanal(dat, 1000, 154)
> for(i in 1:5) print(quantile(result2$factor1[,i], c(0.025, 0.975)))
2.5% 97.5%
0.1970008 0.5675890
2.5% 97.5%
0.2747490 0.6723675
2.5% 97.5%
0.2270197 0.6515273
2.5% 97.5%
0.3206293 0.7580170
2.5% 97.5%
0.3622818 0.8117094
図4 繰り返し数が10000回の場合のヒストグラム
> result2 <- boot.factanal(dat, 10000, 154)
> for(i in 1:5) print(quantile(result2$factor1[,i], c(0.025, 0.975)))
2.5% 97.5%
0.1830492 0.5701160
2.5% 97.5%
0.2730817 0.6690280
2.5% 97.5%
0.2253723 0.6521873
2.5% 97.5%
0.3499537 0.7687571
2.5% 97.5%
0.3846923 0.8114813
以上の結果より、繰り返し数が10000回のほうが信頼区間は狭くなっていることが分かる。しかし、1000回の繰り返しでも正規分布の形状をしていることから、1000回程度の繰り返しで中心極限定理が生まれるのだと予想できる。少なくとも、この場合において10000回もの繰り返しは不要であろう（私のパソコンで10分ほど要した）。また、最初のデータセット（6教科の成績のデータセット）とは異なり、今度は通常の因子分析の結果と同じ0.5を中心とした正規分布に従っているので、これより正確な95%の信頼区間が推定される。どの因子負荷量についても信頼区間に0が含まれていないことから、5%の有意水準のもとでは有意である（推定された因子負荷量は0ではない）。そして不適解も生じていないようである。
今回の結果を踏まえて、改めて最初のデータセットにおける結果を考察してみよう。まずこのブートストラップ法による信頼区間の推定は因子分析に対しては有効ではない、という主張は正しくはないだろう。なぜなら、今回のように正規分布の形状をしており、正確な信頼区間を推定することもできているからである。だとすると、やはり最初のデータセットにおける通常の因子分析の結果は元々妥当なものではなかったということなのだろうか。
ここで最初のデータセットを2因子構造ではなく、1因子構造と仮定して分析してみる。
> r <- factanal(bsf.data, factors=1, rotation="promax", scores="regression")
> r
Loadings:
Factor1
国語 0.388
地理 0.305
英語 0.375
数学 0.538
化学 0.377
物理 0.602
この結果を、1000回の繰り返し計算結果と照らし合わせてみよう。
> ares <- boot.factanal(bsf.data, 1000, 100) > windows() > par(mfrow=c(2,3)) > for(i in 1:6) hist(ares$factor1[,i])
図5 1因子構造として解析した結果
どうやら1因子構造を仮定した場合は上手くいくようである。なぜ2因子構造、ないしはそれ以上の因子構造を持つ場合には上手くいかないのだろうか。1つ考えられる問題は不適解が生じた場合の因子負荷量が悪い影響を及ぼしているということである。
# 不適解が生じているときの因子構造
Loadings:
Factor1 Factor2
国語 0.348 0.198
地理 1.098 -0.296
英語 0.413 0.169
数学 0.269 0.592
化学 -0.398 0.715
物理 0.506
そこで1000回の繰り返し計算の途中で、不適解が生じている（因子負荷量に1.00を越す値が存在する）結果を削除するようにboot.factanal()を再定義してみる。
boot.factanal <- function(dat, repeat.num, s.size){
# 引数に指定された行列の行数と列数の取得
n <- nrow(dat)
k <- ncol(dat)
# n行k列の空のマトリックスを作成
f1.mat <- matrix(numeric(repeat.num*k), nrow=repeat.num, ncol=k)
f2.mat <- matrix(numeric(repeat.num*k), nrow=repeat.num, ncol=k)
# サンプリングするときに使う
x <- 1:n
for(i in 1:repeat.num){
# datからs.size個のデータをリサンプリング
ana.dat <- dat[sample(x, s.size, replace=TRUE), ]
# ana.datに対して因子分析を行う
res <- factanal(ana.dat, factors=2, rotation="promax", scores="regression")
# 因子負荷量を取り出す
f1 <- res$loadings[ ,1]
f2 <- res$loadings[ ,2]
# 不適解が存在すれば何も代入せずに次のループへ
if(any(f1 > 1 | f1 < -1) | any(f2 > 1 | f2 < -1)) next;
# i行目に因子負荷量を代入していく
f1.mat[i, ] <- f1
f2.mat[i, ] <- f2
}
# 代入されなかった行には0が代入されているので、
# その行を削除したものをf1.matに再代入する
f1.mat <- f1.mat[apply(f1.mat, 1, function(f1.mat) all(f1.mat != 0)), ]
f2.mat <- f2.mat[apply(f2.mat, 1, function(f2.mat) all(f2.mat != 0)), ]
list(factor1 = f1.mat, factor2 = f2.mat)
}
この関数を実行すると、不適解が生じなかった因子分析の結果のみが得られることになる。もし不適解が問題となっているならば、これによりきれいな正規分布の形を描くようになるはずである。
> res <- boot.factanal(bsf.data, 1000, 100) > windows() > par(mfrow=c(2,3)) > for(i in 1:6) hist(res$factor1[,i]) --- > dim(res$factor1) [1] 679 6
図6 不適解が生じた結果を除いた場合のヒストグラム（第1因子負荷量）
図6を見ても、やはり問題は解決されていない。どうやら不適解の問題だけではないようである。このヒストグラムを見ると、0付近と0.6付近とを中心とした2つの正規分布形状があるように見える。これは因子分析を行うたびに、得られる因子構造が変化するからであると考えられる。例えば、1回目の因子分析では第1因子より受ける国語、地理、英語の因子負荷量は0.6くらい、第2因子より受ける因子負荷量が0.0くらいであるような場合がある。この場合は第1因子を文系科目、第2因子を理系科目と名づけられるだろう。対して2回目の分析ではこの結果が全く逆になることもある。つまり、第1因子を理系科目、第2因子を文系科目とラベル付けできるような結果が得られるということである（両者は実質的には同じ結果であるが： 表1と表2）。
おそらく2因子以上の構造を仮定した因子分析で、この観測値ベースの標本抽出によるブートストラップ法が有効ではないのは、こういった理由によるものであろう。すなわち、分析のたびに変わる因子構造の変動が問題といえる。もし同じ因子構造をもつ結果のみを上手く抽出することができれば、この問題は解決するだろう。しかし、そういった作業は分析者にとって都合の良い結果のみを選ぶ出すことになるので、好ましいことではないだろう。極端にいえば、たまたま上手くいった分析結果のみをたくさん集めれば、たとえそれが誤った結果であろうと、もっともらしい結果として主張できることになる。
表1 第1因子：理系科目、第2因子文系科目
|
表2 第1因子文系科目、第2因子理系科目
|
まとめると、1因子構造の場合にはブートストラップ法による因子負荷量の推定値の信頼区間を求める方法は有効である。しかし、2因子以上の構造を持つ場合には不適切である。しかし、探索的因子分析を行うことによって特定の構造が見出せたならば、ある因子から影響を受ける変数のみを1因子構造として解析することによって、因子負荷量の信頼区間を求めることができよう。
# 国語、地理、英語を1因子構造として解析
> res1 <- boot.factanal(bsf.data[,1:3], 1000, 100)
# 95%信頼区間
> for(i in 1:3) print(quantile(res1$factor1[,i], c(0.025, 0.975)))
2.5% 97.5%
0.2509642 0.7379250
2.5% 97.5%
0.2664056 0.7586053
2.5% 97.5%
0.3306858 0.9974971
---
# 数学、化学、物理を1因子構造として解析
> res2 <- boot.factanal(bsf.data[,4:6], 1000, 100)
# 95%信頼区間
> for(i in 1:3) print(quantile(res2$factor1[,i], c(0.025, 0.975)))
2.5% 97.5%
0.3400165 0.9338976
2.5% 97.5%
0.2810358 0.7377767
2.5% 97.5%
0.2994456 0.9469305
図7 それぞれを1因子構造として解析した場合のヒストグラム
