Rで共分散構造分析をする時の簡単なテンプレート

# 追記

共分散構造分析に関する新しい記事を書きました。この記事にも共分散構造分析をする際のテンプレートになる様なコードを書きました。実際に分析する際に使用していたコードを載っけたので、こちらの方が参考になるかもしれません。
[Rで共分散構造分析をする際に参考になりそうな情報まとめ - Qiita](https://qiita.com/Masahiro_T/items/d03b9a5a86eee75eb22e#%E5%AE%9F%E9%9A%9B%E3%81%AB%E5%88%86%E6%9E%90%E3%81%99%E3%82%8B%E9%9A%9B%E3%81%AE%E3%82%B3%E3%83%BC%E3%83%89%E3%81%AE%E4%BE%8B)
「実際に分析する際のコードの例」という項です。

---

Rで共分散構造分析をする時のマニュアルみたいなものが欲しかったので、簡単なテンプレートを作成してみました。

分析に関する主なパラメータや、実行結果として表示される省略語の意味などもコメントしてあります。

このままでも、スクリプトを順に実行して行けば、``PoliticalDemocracy``データを共分散構造分析する事が出来ます。

```R:R
# ===================================================================
# ライブラリをロード
# ===================================================================

library(psych)
library(semTools)
library(semPlot)
library(gplots)

# ===================================================================
# データをロード
# ===================================================================

X <- PoliticalDemocracy

# ===================================================================
# データの特徴を簡単に見る
# ===================================================================

# 項目ごとのヒストグラムを見る。
# 正規分布から離れていそうな項目を探したり、
# カテゴリ分布の項目を探す。
hist(X) 

# 相関係数を見る。
# 列同士（or行同士）のグラデーションが似ている項目同士は似た様な項目である可能性がある。
# つまり、共通の潜在変数を持っている可能性がある。
# 
# また、多重共線性を探す時にも使用できる。
cor.plot(X)

# ===================================================================
# 尖度・歪度確認
# ===================================================================

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# semの推定法の選択のために正規性を確認しておく。
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# -------------------------------------------------------------------
# 一次元
# -------------------------------------------------------------------

# 歪度
# >= 2 で歪んでいる
lapply(X, skew)

# 尖度
# >= 7 で尖っている
lapply(X, kurtosis)

# -------------------------------------------------------------------
# 多次元
# -------------------------------------------------------------------

# 歪度
# b1d が歪度
# p値で解釈
mardiaSkew(X) 

# 尖度
# b2d が尖度
# p値で解釈
# z > 3 で最尤推定法の結果が不正確になる可能性
mardiaKurtosis(X) 

# ===================================================================
# データ整形
# ===================================================================

# 例えば、これまでのデータチェックから削除した方が良さそうな特徴量を削除する
# X <- X[, !(names(X) %in% c('x1'))]

# ===================================================================
# 因子数の見当をつける
# ===================================================================

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# VSSと平行分析の結果から、大体どのくらいの因子数が良さそうを調べる。
# 評価指標によって提案される因子数は変わるので、大体これくらいの数かな、といった所。
# 参考：https://www.slideshare.net/simizu706/r-42283141
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# VSS（Very Simple Structure）
# 因子数を変えて評価指標を計算して最適と思われる因子数を表示してくれる。
# 評価指標は MAP BIC の提案数を見ると良い。（ RMSEA も良いかもしれない。）
'
[主要なパラメータ]
n
  意味    : 因子の最大候補数
  default : 8
rotate
  意味    : 回転方法
  options : none, varimax, oblimin, promax （斜交回転は promax と oblimin）
  default : varimax
  用法    : 基本的に oblimin で良い
fm
  意味    : 推定法
  options : pa, minres, mle, fa（ mle は ml でも可）
  default : minres
  用法    : 基本的に minres、サンプルサイズが大きければ ml で良い
'
VSS(X, rotate='oblimin')

# 平行分析
# 因子分析での提案数（number of factors）と主成分分析での提案数（number of components）が表示される。
'
[主要なパラメータ]
fm
  意味    : 推定法
  options : minres, ml, uls, wls, gls, pa 
  default : minres
  用法    : 基本的に minres、サンプルサイズが大きければ ml で良い
'
fa.parallel(X)

# ===================================================================
# 因子分析を実行
# ===================================================================

# 因子分析を実行
'
[主要なパラメータ]
nfactors
  意味    : 因子数
  用法    : 見当を付けた因子数を試せば良い
rotate
  意味    : 斜交回転の方法
  options : 色々ある
  default : oblimin
  用法    : 基本的に oblimin で良い
fm 
  意味    : 推定法
  options : 色々ある
  default : minres
  用法    : 基本的に minres、サンプルサイズが大きければ ml で良い
'
Xfa_rotate <- fa(X, nfactors=2)

# 結果を表示
# 表示される省略語の意味（一部）は次の通り。
'
h2  : 共通性。共通因子（潜在変数）による説明の割合。
u2  : 独自性。独自因子（誤差変数）による説明の割合。u2 = 1 - h2
com : 複雑性。単純構造の指標。1つの因子だけから影響を受けていると1、全部の因子からで因子数。

SS loadings           : 寄与。全ての観測変数にどのくらい影響を与えているかを意味する。
proportion var        : 寄与率。寄与率 = 寄与 / 観測変数の総数
cumulative var        : 累積寄与率
proportion explained  : 説明率。説明率 = 寄与 / 寄与の和
cumulative proportoin : 累積説明率
'
Xfa_rotate

# 値が大きな因子負荷だけを表示
# これを参考にしてsemのモデルを作成
Xfa_rotate$loadings

# ===================================================================
# sem実行
# ===================================================================

# モデルを定義
model <- '
  # measurement model
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    dem65 =~ y5 + y6 + y7 + y8
  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60
  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 + y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

# sem実行
'
[主要なパラメータ]
estimator
  意味    : 推定法
  options : 色々
  default : ML
  用法    : 基本的に次で良い : 
    尖度・歪度を調べた結果、
    a) 全部の変数に正規性があるならML
    b) 少し正規性がない変数があるくらいならMLR
    c) それ以外ならWLSMV
ordered
  意味    : カテゴリ変数の列名の配列
  用法    : カテゴリ変数があるときに設定すれば良い
std.lv
  意味    : モデル識別のための母数制約
  options : TRUE （潜在変数の分散・残差分散を1に固定）
            FALSE（潜在変数ごとに、最初のパスの因子負荷量を1に固定）
  default : FALSE
std.ov
  意味    : 観測変数を分析前に標準化する
  options : TRUE, FALSE
  default : FALSE
missing
  意味    : 欠損値の取り扱い
  options : 
    listwise : 欠損を含む行は削除される
    fiml     : 完全情報最尤推定法が適用される（ml, direct でも可。欠損の仕方が MCAR / MAR の場合のみ有効。）
  default : listwise
'
fit <- sem(model, data=X, estimator="WLSMV")

# ===================================================================
# 評価
# ===================================================================

# 結果表示
'
[主要なパラメータ]
standardized
  意味    : 標準化推定値も表示
  options : TRUE, FALSE
  default : FALSE
fit.measures
  意味    : 適合度も表示
  options : TRUE, FALSE
  default : FALSE
'
# また、表示される省略語の意味（一部）は次の通り。
'
<<適合度指標関連>>

[章名]
Model Test User Model     : 仮定したモデルの検定
Model Test Baseline Model : 独立モデル（全ての観測変数間にパスを引かず、分散だけ推定するモデル。最悪のモデル）の検定

[項目名]
Test Statistic       : 検定統計量(χ2値)
Degrees of freedom   : 自由度 = (標本分散の数 + 標本共分散の数) - (母数の数 ( = 係数の数 + 残差分散の数 ))。
                       推定するには自由度が1以上ある必要がある。
P-value (Chi-square) : p値。帰無仮説 = モデルが正しい なので有意水準より大きい方が良い。

[列名（評価指標をロバストに算出する推定法（MLRやWLSMV）を使用した場合）]
Standard : 通常通りに算出した場合の値
Robust   : ロバストに算出した場合の値

<<推定値関連>>

[章名]
Latent Variables : 潜在変数の推定値
Regressions      : 構造方程式の推定値
Covariances      : 共分散の推定値
Intercepts       : 切片の推定値
Variances        : 分散の推定値（外生変数なら分散、内生変数なら残差分散）

[列名]
Estimate : 推定値
Std.Err  : 推定値の標準誤差 
z-value  : = 推定値/標準誤差。帰無仮説 = 推定値が0 の検定に使用する。
P(>|z|)  : 帰無仮説 = 推定値が0 の検定のp値
Std.lv   : 潜在変数だけ標準化（平均0・分散1に変換）してから推定した時の推定値
Std.all  : 潜在変数と観測変数を標準化してから推定した時の推定値（標準化解と呼ばれる）
'
summary(fit, standardized=T, fit.measures=T)

# 適合度を表示
fitMeasures(fit, fit.measures="all")

# 決定係数を表示
# モデル変更の参考に使えるかもしれない。
inspect(fit, what="rsquare")

# ===================================================================
# 相関行列
# ===================================================================

# 推定値を基に相関行列を作成
# 潜在変数ペアのうち相関係数が高いものは
# 似た様な概念を表している
# ので、場合によっては一つの潜在変数にしても良い
cor <- inspect(fit, what="cor.all")
heatmap.2(cor, dendrogram='none', Rowv=F, Colv=F)
cor

# 相関行列の残差行列を計算
# 絶対値が大きい要素がある変数に集中していれば、
# その変数を削除したりその変数へのパスを変更するなどすると良いかもしれない
cor <- resid(fit, type="cor")
cor <- cor$cov
heatmap.2(cor, dendrogram='none', Rowv=F, Colv=F)
cor

# ===================================================================
# 修正指標
# ===================================================================

# 修正指標を計算
MI <- modificationIndices(fit)

# mi昇順で表示
MI[order(MI$mi),]

# ===================================================================
# パス図をプロット
# ===================================================================

'
[主要なパラメータ]
whatLabels
  意味    : 何をパスのラベルにするか。
  options : 色々ある
  default : パラメータ「what」の値に依存
  用法    : 例として、stand で標準化解が表示される。
edge.label.cex
  意味    : 係数のフォントサイズ
  default : 0.6
sizeMan
  意味    : 観測変数のフォントサイズ
  default : 5
sizeLat
  意味    : 潜在変数のフォントサイズ 
  default : 8
rotation
  意味    : 図の向き
  default : 1
nCharNodes
  意味    : 変数名を何文字まで表示するか
  default : 3
optimizeLatRes
  意味    : 潜在変数の残差の角度を最適化して表示する。
  default : False
'
semPaths(object=fit, whatLabels="stand")
```

このまま実行した場合に表示されるパス図：
<img width="545" alt="スクリーンショット 2020-05-11 午前0.12.30.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/178570/8c310f0d-9aa7-cdec-b793-ba15c2c9f007.png">


# 主に参考にしたもの
[The lavaan Project](http://lavaan.ugent.be/tutorial/sem.html)
簡素なlavaanチュートリアルです。
データセットとそのモデルを参考にしました。

[Rで因子分析　商用ソフトで実行できない因子分析のあれこれ（SlideShare）](https://www.slideshare.net/simizu706/r-42283141)
因子分析全般で参考にしました。

[共分散構造分析 R編](https://www.amazon.co.jp/%E5%85%B1%E5%88%86%E6%95%A3%E6%A7%8B%E9%80%A0%E5%88%86%E6%9E%90-R%E7%B7%A8%E2%80%95%E6%A7%8B%E9%80%A0%E6%96%B9%E7%A8%8B%E5%BC%8F%E3%83%A2%E3%83%87%E3%83%AA%E3%83%B3%E3%82%B0-%E8%B1%8A%E7%94%B0-%E7%A7%80%E6%A8%B9/dp/4489021801)
共分散構造分析とは？と、Rで共分散構造分析をする方法が書いてある書籍です。
尖度・歪度全般や、sem全般で参考にしました。