【R・Lavaan ( sem )・共分散構造分析】色んなパターンでパスを一つだけ追加した場合のsemの実行結果の適合度を比べるコード

表題の通りです。

用途としては、「パスを追加した方が良いモデルが出来そうだけど、どのパスを追加すれば良いのかが分からない時」などを想定しています。

R

'
<<以下のコードで行う事の概要>>

色んなパターンでパスを一つだけ追加した場合のsem実行結果の適合度を比べる。

比較するパターンは次の通り。

・ 　　　　    潜在変数 =~ 観測変数
・ 　　　　    潜在変数 ~  観測変数 or 潜在変数
・ 観測変数    　　　　 ~  観測変数
・ 観測変数 or 潜在変数 ~~ 観測変数 or 潜在変数

それぞれのパターンでsemを実行した結果を記録し、最後に cfi.scaled で降順ソートしたものを表示する。

（sem実行時に warning や error が出た場合はメッセージだけ表示して次のパターンに進む。）



<<下のコードを実行する前に既に実行済みと想定している事>>

1) 使用するデータセットの定義
# 例
X <- read.csv("data.csv")

2) ベースとなるモデルの定義
# 例
model <- "
  f1 =~ x1 + x2 + x3
  f2 =~ x4 + x5 + x6 + x7
  f2 ~ f1
" 

3) これらを用いた sem の実行結果
# 例
fit <- sem(model, data=X, estimator="WLSMV", std.lv=T)
'


# withTimeout用
library('R.utils')

# モデル内で定義されている変数名の一覧を取得
ov_ary  <- colnames(resid(fit, type="raw")$cov)   # 観測変数
lv_ary  <- colnames(predict(object=fit))          # 潜在変数
all_ary <- colnames(inspect(fit, what="cor.all")) # 両方

# 追加する式のペアを用意
pair_ary <- list(
  
  # c(左辺の値, 右辺の値, 関係性)
  list(lv_ary,  ov_ary,  ' =~ '), # f      =~ x
  list(lv_ary,  all_ary, ' ~ ' ), # f      ~  x or f
  list(ov_ary,  ov_ary,  ' ~ ' ), # x      ~  x
  list(all_ary, all_ary, ' ~~ ')  # f or x ~~ f or x
  
)

# ここに記録していく
result_df_all <- data.frame()

for (pair in pair_ary) {
  
  left_v_ary  <- pair[[1]]
  right_v_ary <- pair[[2]]
  sep         <- pair[[3]]
  
  for (left_v in left_v_ary) {
    for (right_v in right_v_ary) {
      
      # 追加する式を作成
      s <- paste(left_v, right_v, sep=sep)
      print(s)
      
      # 式を追加したモデル定義を作成
      model_modified <- paste(model, s, sep='\n')
      
      # tryCatch する事で、
      # 1) for 文の中でも warning が表示される様にする
      # 2) error handling する事で error で for 文が止まらない様にする
      tryCatch({
        
          # タイムアウト付きでsem実行
          # タイムアウトを付けないと、収束がとても遅い様なパターンで時間がかかってしまう
          # そしてその様なパターンでは（経験的に）良い解が得られない
          # そのため、収束が遅い場合はタイムアウトする様にする
          fit <- withTimeout({
            sem(model_modified, data=X, estimator="WLSMV", std.lv=T)
          }, timeout = 3.0, onTimeout = "warning")
          
          # 結果取得
          result <- fitMeasures(fit, fit.measures=c('rmsea.scaled', 'cfi.scaled', 'tli.scaled', 'agfi', 'srmr'))
          
          # hoge$foo ができる様にlist型にする
          result <- as.list(result)
          
          # 追加した変数も記録する様にする
          result$left_v  <- left_v
          result$sep     <- sep
          result$right_v <- right_v

          # 推定値が > 1、 < -1 の数も記録
          # warning や error が出なくても、1、-1と比べて大きすぎる、小さすぎる推定値があったりする
          # そのため、この様な推定値の有無が分かる様にしておくと良い
          p <- parameterEstimates(fit, standardized=T)
          result$n_std.all_too_small          <- nrow(p[p$std.all < -1,])
          result$n_std.all_too_large          <- nrow(p[p$std.all >  1,])
          result$n_std.all_too_small_or_large <- nrow(p[p$std.all < -1,]) + nrow(p[p$std.all >  1,])
          
          # 追記
          result_df <- data.frame(result)
          result_df_all <- rbind(result_df_all, result_df)
        }, 
        
        # warningが起きたらメッセージを表示する
        warning = function(e) {
          message(e)
          cat('\n') # 改行
        },
        
        # errorが起きたらメッセージを表示する
        error=function(e) {
          message(e)
          cat('\n') # 改行
        }
        
      )
    }
  }
}

# cfi.scaled が降順になる様にする
result_df_all <- result_df_all[order(result_df_all$cfi.scaled, decreasing=T),]
result_df_all

# 保存
readr::write_excel_csv(result_df_all, 'result_df_all.csv')

結果の例（最初の数行）：