探索的因子分析で完全情報最尤法を実行するRの関数を作りました

 

久々の記事です。※追記あり

タイトルのまんまですが、Rで探索的因子分析(EFA)で完全情報最尤法(FIML)を実行できる方法について質問されたので、そういやそういうパッケージないなーと思ってました。たぶん探せばあるとは思うのですが、探すより作ったほうが早いと思い、作ってみました。

何を思いついたかと言えば、lavaanはCFAでFIMLが実行できます。EFAはCFAのパス図を工夫すれば、推定ができますから、あとはlavaan用のモデルを自動生成する関数を作ればいいだけです。

必要パッケージはlavaanだけですが、推定した因子負荷行列について因子軸の回転がいると思うので、GPArotationパッケージも入れておくといいと思います。あと、psychも入れておくとFIMLじゃない普通の推定結果と比較できますね。

よくよく考えたら、回転したあとの因子得点がないとFIMLの意味がないと思ったので、関数を修正しました。引数に回転方法(faと同じ引数です)もいれるようにしました。

というわけで関数fimlEFAのコードです。そんなに長くないです。

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
  rot.loading <- do.call(rotate,list(raw.loading))
  fscore <- lavPredict(fit, fsm = F)%*%rot.loading$rotmat
  print(rot.loading)
  output <- list(fit=rot.loading,fscore=fscore,model=model)
  return(output)
}

 

引数はデータと因子数だけと回転方法です。

使い方は次のような感じ。

library(lavaan)
library(GPArotation)
library(psych)

dat <- bfi[1:25]
head(dat)

result <- fimlEFA(dat,nfactors=5)
result$fit$loadings #因子負荷行列
result$fscore #因子得点

結果は次のような感じ。

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ができましたね。めでたしめでたし。

 

This entry was posted in 心理統計学. Bookmark the permalink.