Rで共分散構造分析 {lavaan}パッケージ編 #rstatsj

Rで共分散構造分析を行うにはこれまでは John Fox さんの{sem}パッケージを使うのが普通でした。
私も以前、{sem}パッケージの使い方についてエントリーを書いています。
id:bob3:20091226
id:bob3:20091227
id:bob3:20100530

しかし最近は{OpenMX}や{lavaan}など共分散構造分析のためのパッケージが他にも出てきました。
特に{lavaan}はまだβ版ながら評判が良いようです。
{lavaan}の公式サイトはこちら→ lavaan - latent variable analysis

というわけで、ブルーバックス『原因をさぐる統計学』に載っている事例をいくつが{lavaan}でやってみましょう。

原因をさぐる統計学―共分散構造分析入門 (ブルーバックス)

まずは{lavaan}と関連パッケージのインストール

install.packages(c("lavaan", "psych", "qgraph"))

パッケージの呼び出し。

library(lavaan)
library(psych)
library(qgraph)

まずは『原因を探る統計学』のp129、p237に載っている「ヘッド・スタート計画」の事例から。

まず、分析する相関係数行列を入力します。

# ヘッドスタート計画の相関係数行列
HS <- matrix(c(1.000, 0.484, 0.224, 0.268, 0.230, 0.265,
                0.484, 1.000, 0.342, 0.215, 0.215, 0.295,
                0.224, 0.342, 1.000, 0.387, 0.196, 0.234,
                0.268, 0.215, 0.387, 1.000, 0.115, 0.162,
                0.230, 0.215, 0.196, 0.115, 1.000, 0.635,
                0.265, 0.295, 0.234, 0.162, 0.635, 1.000), 6)
colnames(HS) <- rownames(HS) <- paste("x", 1:6, sep="")

それぞれの変数の意味は、
x1:母親の学歴
x2:父親の学歴
x3:父親の職業
x4:家庭の収入
x5:ＭＲＴ（知能検査）
x6:ＩＴＰＡ（知能検査）

次にモデルの記述です。
「多重指標モデル」と「外生潜在変数が２つのモデル」の２パターンでやってみます。

# 多重指標モデル
HS.mdl.1 <- '
  # 測定方程式の記述
  f1 =~ x1 + x2 + x3 + x4
  f2 =~ x5 + x6
  # 構造方程式の記述
  f2 ~ f1 '

# 外生潜在変数が２つのモデル
# 因子間相関は記述しなくてもOK。
HS.mdl.2 <- '
  # 測定方程式
  f1 =~ x1 + x2
  f2 =~ x3 + x4
  f3 =~ x5 + x6
  # 構造方程式
  f3 ~  f1 + f2 '

{lavaan}パッケージではモデルの記述に以下のような記号の使い方をしています。
=~ 測定方程式
~ 構造方程式（回帰）
~~ 残差の共分散（相関）

それでは分析の実行です。

HS.fit.1 <- sem(HS.mdl.1, sample.cov=HS, sample.nobs=1800)
summary(HS.fit.1, fit.measures=TRUE, standardized=TRUE)

HS.fit.2 <- sem(HS.mdl.2, sample.cov=HS, sample.nobs=1800)
summary(HS.fit.2, fit.measures=TRUE, standardized=TRUE)

出力結果は割愛しますが、CFI、TLI、AIC、BIC、RMSEAといった指標値に続き、推定された係数が出てきます。（なぜかメジャーなGFI、AGFIは出力されない…）
Std.allというのが標準化係数のようです。
また、sem()を実行したときに「相関係数行列からの分析はまだサポートしてないよ！」という警告メッセージが出ます。
このあたりはβ版ゆえ仕方がないでしょう。
それでもテキストと同じ数値が出てきます。

一応、パス図の描画もしてくれます。
ただ、本当に“一応”という感じで、モデルの確認程度にしか使えないと思います。

diagram(HS.fit.1, "多重指標モデル", errors=TRUE, cut=0, lr=FALSE)

qgraph(HS.fit.1, titles=FALSE)

※はてなフォトライフがメンテ中で画像を載せられませんでした…

diagram(HS.fit.2, "潜在変数が２つのモデル", errors=TRUE, cut=0, lr=FALSE)

※はてなフォトライフがメンテ中で画像を載せられませんでした…

qgraph(HS.fit.2, titles=FALSE)

※はてなフォトライフがメンテ中で画像を載せられませんでした…

といった具合で、まだまだ不具合や不満点はあるものの、多母集団同時分析や平均構造モデルにも対応するということで、{lavaan}は今後が楽しみなパッケージです。
時間ができたら他の事例でも試してみます。

あ、そうそう、転職しました。