久々の記事です。※追記あり
タイトルのまんまですが、Rで探索的因子分析（EFA)で完全情報最尤法（FIML)を実行できる方法について質問されたので、そういやそういうパッケージないなーと思ってました。たぶん探せばあるとは思うのですが、探すより作ったほうが早いと思い、作ってみました。
何を思いついたかと言えば、lavaanはCFAでFIMLが実行できます。EFAはCFAのパス図を工夫すれば、推定ができますから、あとはlavaan用のモデルを自動生成する関数を作ればいいだけです。
必要パッケージはlavaanだけですが、推定した因子負荷行列について因子軸の回転がいると思うので、GPArotationパッケージも入れておくといいと思います。あと、psychも入れておくとFIMLじゃない普通の推定結果と比較できますね。
よくよく考えたら、回転したあとの因子得点がないとFIMLの意味がないと思ったので、関数を修正しました。引数に回転方法（faと同じ引数です）もいれるようにしました。
というわけで関数fimlEFAのコードです。そんなに長くないです。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 |
fimlEFA <- function(data,nfactors=1,rotate="Promax"){ A <- nfactors P <- ncol(data) model <- '' for(a in 1:A){ if(a==1){ model <- paste0(model,"F",a,"=~NA*",colnames(data)[1]) }else{ model <- paste0(model,"F",a,"=~") for(a2 in 1:(a-1)){ if(a2==1){ model <- paste0(model,"0*",colnames(data)[a2]) }else{ model <- paste0(model,"+0*",colnames(data)[a2]) } } model <- paste0(model,"+",colnames(data)[a]) } for(p in (a+1):P){ model <- paste0(model,"+",colnames(data)[p]) } model <- paste0(model,"\n") if(a < A){ for(a2 in (a+1):A){ if(a2==(a+1)){ model <- paste0(model,"F",a,"~~0*F",a2) }else{ model <- paste0(model,"+0*F",a2) } } model <- paste0(model,"\n") } model <- paste0(model,"F",a,"~~1*F",a,"\n") } require(lavaan) fit <- cfa(model,data=data,missing="fiml") raw.loading <- inspect(fit,what="std")$lambda if(nfactors>1){ rot.loading <- do.call(rotate,list(raw.loading)) if(is.null(rot.loading$rotmat)==FALSE){ fscore <- lavPredict(fit, fsm = F)%*%rot.loading$rotmat }else{ fscore <- lavPredict(fit, fsm = F)%*%rot.loading$Th } }else{ rot.loading <- raw.loading fscore <- lavPredict(fit, fsm = F) } print(rot.loading) output <- list(loading=rot.loading,fscore=fscore,model=model) return(output) } |
引数はデータと因子数だけと回転方法です。
使い方は次のような感じ。
|
1 2 3 4 5 6 7 8 9 10 |
library(lavaan) library(GPArotation) library(psych) dat <- bfi[1:25] head(dat) result <- fimlEFA(dat,nfactors=5) result$loading #因子負荷行列 result$fscore #因子得点 |
結果は次のような感じ。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
Call: NULL Standardized loadings (pattern matrix) based upon correlation matrix F1 F2 F3 F4 F5 h2 u2 A1 0.37 -0.12 0.05 0.23 0.04 0.15 0.85 A2 -0.59 -0.07 0.07 -0.03 -0.02 0.40 0.60 A3 -0.64 -0.17 0.00 -0.04 -0.02 0.51 0.49 A4 -0.44 -0.07 0.19 -0.06 0.16 0.28 0.72 A5 -0.56 -0.26 -0.04 -0.14 -0.04 0.48 0.52 C1 0.01 0.06 0.54 0.07 -0.15 0.32 0.68 C2 -0.09 0.13 0.66 0.13 -0.05 0.43 0.57 C3 -0.09 0.10 0.59 0.04 0.07 0.32 0.68 C4 -0.07 -0.01 -0.66 0.11 0.02 0.46 0.54 C5 -0.02 0.12 -0.57 0.13 -0.11 0.43 0.57 E1 0.07 0.64 0.15 -0.13 0.09 0.37 0.63 E2 0.06 0.72 0.03 0.02 0.04 0.55 0.45 E3 -0.27 -0.44 -0.07 0.08 -0.29 0.44 0.56 E4 -0.32 -0.60 -0.04 0.03 0.07 0.52 0.48 E5 -0.04 -0.46 0.25 0.22 -0.19 0.41 0.59 N1 0.16 -0.18 0.03 0.90 0.06 0.71 0.29 N2 0.15 -0.11 0.04 0.86 0.00 0.66 0.34 N3 -0.06 0.06 -0.03 0.68 -0.02 0.52 0.48 N4 -0.08 0.37 -0.12 0.40 -0.09 0.48 0.52 N5 -0.19 0.20 0.00 0.44 0.14 0.34 0.66 O1 -0.02 -0.11 0.03 0.01 -0.52 0.32 0.68 O2 -0.19 -0.03 -0.09 0.16 0.44 0.24 0.76 O3 -0.08 -0.19 -0.03 0.02 -0.62 0.47 0.53 O4 -0.16 0.31 -0.03 0.07 -0.39 0.26 0.74 O5 -0.08 -0.05 -0.03 0.11 0.52 0.27 0.73 F1 F2 F3 F4 F5 SS loadings 1.92 2.32 1.99 2.55 1.56 Proportion Var 0.08 0.09 0.08 0.10 0.06 Cumulative Var 0.08 0.17 0.25 0.35 0.41 Proportion Explained 0.19 0.22 0.19 0.25 0.15 Cumulative Proportion 0.19 0.41 0.60 0.85 1.00 F1 F2 F3 F4 F5 F1 1.00 0.24 -0.23 -0.07 0.21 F2 0.24 1.00 -0.36 0.37 0.15 F3 -0.23 -0.36 1.00 -0.26 -0.24 F4 -0.07 0.37 -0.26 1.00 -0.01 F5 0.21 0.15 -0.24 -0.01 1.00 |
bfiはデータサイズが大きいですが、それでも数十秒で分析が終わりました。
HADにもFIMLでEFAができますが、断然こっちのほうが早いですね。もちろん結果は一致してました。
というわけで、Rでとてもお手軽にFIMLでEFAができましたね。めでたしめでたし。