RでSEM（共分散構造分析／構造方程式モデリング）3

半年ぐらい前にR言語のsemパッケージで共分散構造分析の練習をしました。
id:bob3:20091226#p1
id:bob3:20091227#p1
このたび、次々回のTokyo.R#7で『Ｒによるやさしい統計学』の第17章「共分散構造分析」パートを担当することになったので復習中。

Rによるやさしい統計学

例題として『Excelで学ぶ共分散構造分析とグラフィカルモデリング』p22の「住意識」データを使ってみます。

Excelで学ぶ共分散構造分析とグラフィカルモデリング

以下の相関行列をCSVファイルとして保存。ここでは「cor.csv」としておきます。

x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13
1,0.423491329,0.403553642,0.127291949,0.179424357,0.132039721,0.148471449,0.213510548,0.224842949,0.25177753,0.125013841,0.122467773,0.191899015
0.423491329,1,0.75117574,0.187228911,0.267487322,0.115166239,0.188329495,0.22264277,0.264941282,0.300143211,0.143769574,0.096972479,0.192815288
0.403553642,0.75117574,1,0.238411139,0.275971308,0.178031892,0.239260433,0.24933465,0.280945687,0.325355384,0.174399645,0.146946528,0.220178984
0.127291949,0.187228911,0.238411139,1,0.411340708,0.399252508,0.349942359,0.286442891,0.273901189,0.272518259,0.149680419,0.19198312,0.176169959
0.179424357,0.267487322,0.275971308,0.411340708,1,0.54530188,0.329966496,0.302289993,0.296037879,0.332500945,0.234793591,0.248269672,0.224471563
0.132039721,0.115166239,0.178031892,0.399252508,0.54530188,1,0.298241169,0.327162049,0.228691718,0.362303562,0.173482684,0.242372426,0.143496196
0.148471449,0.188329495,0.239260433,0.349942359,0.329966496,0.298241169,1,0.265943207,0.205411512,0.239727945,0.140376895,0.152658635,0.16090992
0.213510548,0.22264277,0.24933465,0.286442891,0.302289993,0.327162049,0.265943207,1,0.416647368,0.553811105,0.338026113,0.292691735,0.324305158
0.224842949,0.264941282,0.280945687,0.273901189,0.296037879,0.228691718,0.205411512,0.416647368,1,0.371119334,0.222578881,0.170531959,0.221756995
0.25177753,0.300143211,0.325355384,0.272518259,0.332500945,0.362303562,0.239727945,0.553811105,0.371119334,1,0.286346812,0.291387323,0.283200002
0.125013841,0.143769574,0.174399645,0.149680419,0.234793591,0.173482684,0.140376895,0.338026113,0.222578881,0.286346812,1,0.334934818,0.452021148
0.122467773,0.096972479,0.146946528,0.19198312,0.248269672,0.242372426,0.152658635,0.292691735,0.170531959,0.291387323,0.334934818,1,0.332915322
0.191899015,0.192815288,0.220178984,0.176169959,0.224471563,0.143496196,0.16090992,0.324305158,0.221756995,0.283200002,0.452021148,0.332915322,1

それぞれの変数の意味は以下の通り。
x1 「給食サービス」
x2 「デイサービス」
x3 「老人ホーム」
x4 「自然の緑」
x5 「ウォーキング、ジョギングしたい道」
x6 「大きな公園」
x7 「河川・湖」
x8 「コミュニティ活発イベント」
x9 「多世代が住む」
x10「活発な地域コミュニティ」
x11「趣味の合う交流」
x12「ホームパーティー」
x13「自分の経験、能力を還元」

モデルもテキストファイルで用意。

モデルの表現は、左から変数間の関係、係数の名前、初期値。
ここでF1〜F4は潜在変数（因子）で、F1「福祉の充実」、F2「自然環境」、F3「地域コミュニティ」、F4「手作りの交流」。

まずは検証的因子分析モデル(p125)。「model1.txt」として保存。

# 潜在変数から観測変数へのパス
F1 -> x1,  p01, NA
F1 -> x2,  p02, NA
F1 -> x3,  p03, NA
F2 -> x4,  p04, NA
F2 -> x5,  p05, NA
F2 -> x6,  p06, NA
F2 -> x7,  p07, NA
F3 -> x8,  p08, NA
F3 -> x9,  p09, NA
F3 -> x10, p10, NA
F4 -> x11, p11, NA
F4 -> x12, p12, NA
F4 -> x13, p13, NA
# 因子間相関係数（共分散）
F1  <-> F2, F1F2, NA
F1  <-> F3, F1F3, NA
F1  <-> F4, F1F4, NA
F2  <-> F3, F2F3, NA
F2  <-> F4, F2F4, NA
F3  <-> F4, F3F4, NA
# 観測変数の残差（分散）
x1  <-> x1,  e01, NA
x2  <-> x2,  e02, NA
x3  <-> x3,  e03, NA
x4  <-> x4,  e04, NA
x5  <-> x5,  e05, NA
x6  <-> x6,  e06, NA
x7  <-> x7,  e07, NA
x8  <-> x8,  e08, NA
x9  <-> x9,  e09, NA
x10 <-> x10, e10, NA
x11 <-> x11, e11, NA
x12 <-> x12, e12, NA
x13 <-> x13, e13, NA
# 潜在変数の分散
F1  <-> F1,  NA, 1
F2  <-> F2,  NA, 1
F3  <-> F3,  NA, 1
F4  <-> F4,  NA, 1

次に多重指標モデル(p127)。「model2.txt」として保存。

# 潜在変数から観測変数へのパス
F1 -> x1,  p01, NA
F1 -> x2,  p02, NA
F1 -> x3,  p03, NA
F2 -> x4,  p04, NA
F2 -> x5,  p05, NA
F2 -> x6,  p06, NA
F2 -> x7,  p07, NA
F3 -> x8,  p08, NA
F3 -> x9,  p09, NA
F3 -> x10, p10, NA
F4 -> x11, p11, NA
F4 -> x12, p12, NA
F4 -> x13, p13, NA
# 因子間のパス
F4 -> F3,  p14, NA
F3 -> F2,  p15, NA
F3 -> F1,  p16, NA
# 観測変数の残差（分散）
x1  <-> x1,  e01, NA
x2  <-> x2,  e02, NA
x3  <-> x3,  e03, NA
x4  <-> x4,  e04, NA
x5  <-> x5,  e05, NA
x6  <-> x6,  e06, NA
x7  <-> x7,  e07, NA
x8  <-> x8,  e08, NA
x9  <-> x9,  e09, NA
x10 <-> x10, e10, NA
x11 <-> x11, e11, NA
x12 <-> x12, e12, NA
x13 <-> x13, e13, NA
# 潜在変数の分散
F1  <-> F1,  NA, 1
F2  <-> F2,  NA, 1
F3  <-> F3,  NA, 1
F4  <-> F4,  NA, 1

さて、それではRを起動しましょう。

install.packages("sem") # semパッケージのインストール
library(sem)            # semパッケージの呼び出し

# データの読み込み
DATA <- as.matrix(read.csv("cor.csv"))
rownames(DATA) <- colnames(DATA)

# モデルの読み込み
MODEL1 <- specify.model("model1.txt")
MODEL2 <- specify.model("model2.txt")

# 検証的因子分析モデル
MODEL1.sem <- sem(MODEL1, DATA, N=1120)

# 結果の表示
# 適合度としてGFI、AGFI、RMSEA、NFI、NNFI、CFI、SRMR、BIC
# 推定されたパラメータなどが出力されます
summary(MODEL1.sem)

# 標準化解の表示
std.coef(MODEL1.sem)

# 修正指数の算出
mod.indices(MODEL1.sem)

# dot言語によるパス図の出力。
path.diagram(MODEL1.sem, ignore.double=FALSE, edge.labels="values", digits=2,
             standardize=TRUE)

# 多重指標モデル
MODEL2.sem <- sem(MODEL2, DATA, N=1120)
summary(MODEL2.sem)
path.diagram(MODEL2.sem, ignore.double=FALSE, edge.labels="values", digits=2,
             standardize=TRUE)

path.diagram()の出力を EasyGraphViz に食わせるとこんな感じになります。

……これ、説明すんの難しいなぁ。
特にモデル識別のための分散の固定。