COLOR(red){SIZE(30){Rで共分散構造分析・構造方程式モデル}}
#contents
*パッケージ [#s5504cb5]
-[[sem(構造方程式モデル)パッケージ中のオブジェクト一覧]]
-[[plotSEMM: SEMM の非線形潜在変数相互作用のグラフ化:http://cran.r-project.org/web/packages/plotSEMM/plotSEMM.pdf]]
-[[OpenMx: OpenMx(構造方程式モデル)Mx: 共分散構造分析ソフトのR版:http://openmx.psyc.virginia.edu/]]
--[[OpenMx package for Structural Equation Model in R:http://pairach.com/2012/04/16/openmx-package-for-structural-equation-model-in-r/]]
-[[semtools - Methods for SEMs:https://r-forge.r-project.org/projects/semtools/]] - semtools - Methods for Diagnostic Testing of Structural Equation Models
-[[simsem: SIMulated Structural Equation Modeling:http://cran.md.tsukuba.ac.jp/web/packages/simsem/]]
-[[semPLS: Structural Equation Modeling Using Partial Least Squares:http://cran.md.tsukuba.ac.jp/web/packages/semPLS/]]
*使用方法 [#e6587b3f]
**パッケージのロード [#qd9d2629]
library(sem)
***共分散行列計算 [#ie63acc2]
cov(data)で分散共分散行列を、cor(data)で相関行列を計算できる。
**モデルの指定 [#a62e8f9b]
specifyModel()関数でモデルを指定する方法は次の通り。
1.表現は順に、変数間関係、パラメタの名前、初期値である。
2.矢印による関係の表現:
大なり・小なりの記号とハイフンを使った矢印、すなわち「<-」は二変数間の回帰係数を表現する。A->BあるいはB<-Aという表現は同じ意味である。
また、両方向の矢印、A<->AはAの分散を、A<->BはAとBの共分散(相関係数)を表現する。
ここでのAやBは変数名で、データセットの中に同名の変数があれば顕在変数として、無ければ潜在変数を仮定したものとして扱われる。
ハイフンの数に制限はなく、0本以上、何本でも良い。すなわち、A>B、A->B、A-->Bは全て同じ意味である。
3.パラメタの名前
回帰係数、分散、共分散を特定するために名前をつける。同じ名前が二つ以上の矢印に付されているときは、同値制約をかけたものと見なす。
ここにNAと表記すると、固定パラメタとして扱われる。
4.初期値
自由母数の初期値、あるいは固定母数の値を入れる。ここにNAと表記すると、semは初期値も計算する。
初心者の方へ: モデル式中にはパス係数、分散、共分散の全てを表記しなければならない。パス係数だけ指定してもモデルが推定されないので、要注意である。
**semの実行 [#rf4a0e22]
sem(モデル,データ行列,サンプル数)
で推定する。モデルはspecifyModel()関数で指定したもの、データ行列は分散共分散行列か相関行列を与える。
summary()
で推定値、適合度等の結果を返す。
**潜在変数を仮定しないモデル [#j12acd36]
下のMacintoshでの実行例を参照。
**潜在変数を仮定するモデル [#rb420b01]
[[青木先生の因子分析を行うRコード:http://aoki2.si.gunma-u.ac.jp/R/factanal2.html]]を例に、5変数1因子モデルをかいてみる。
コードは次の通り。
dat <- matrix(c( # 10 ケース，5 変数のデータ行列例（ファイルから読んでも良い）
-1.89, -0.02, 0.42, 1.23, -1.53,
0.06, 1.81, -0.59, -0.75, -0.12,
2.58, -0.20, -1.92, -0.49, -0.35,
0.69, -0.66, -0.77, -1.92, 0.38,
-1.05, 0.07, -0.37, -0.89, -1.62,
-2.73, 1.40, 0.18, -0.09, 0.13,
0.95, 0.84, 1.19, 1.19, 0.10,
0.93, 1.17, -1.70, 1.46, -0.25,
-0.04, -0.12, -0.34, -0.24, 1.88,
0.16, 1.03, -0.05, -0.73, 0.04
), byrow=TRUE, ncol=5)
colnames(dat) <- c("v1","v2","v3","v4","v5")
model <- specifyModel()
v1 < F1,path1,1
v2 < F1,path2,NA
v3 < F1,path3,NA
v4 < F1,path4,NA
v5 < F1,path5,NA
v1 <> v1,s1,NA
v2 <> v2,s2,NA
v3 <> v3,s3,NA
v4 <> v4,s4,NA
v5 <> v5,s5,NA
F1 <> F1,NA,1
ans <- sem(model,cor(dat),10)
summary(ans)
因子F1の分散を1に固定する（因子分析の仮定より）ことを忘れずに。
結果は次の通り。
Model Chisquare = 1.9254 Df = 5 Pr(>Chisq) = 0.85937
Chisquare (null model) = 6.2471 Df = 10
Goodness-of-fit index = 0.93485
Adjusted goodness-of-fit index = 0.80454
RMSEA index = 0 90% CI: (NA, 0.24729)
Bentler-Bonnett NFI = 0.6918
Tucker-Lewis NNFI = -0.63854
Bentler CFI = NaN
SRMR = 0.087271
BIC = -9.5876
Normalized Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.98e-01 -2.98e-06 1.55e-02 1.22e-01 2.43e-01 6.97e-01
Parameter Estimates
Estimate Std Error z value Pr(>|z|)
path1 -0.79356 0.60066 -1.32113 0.186458 v1 <--- F1
path2 0.28438 0.39750 0.71543 0.474344 v2 <--- F1
path3 0.67625 0.51706 1.30789 0.190912 v3 <--- F1
path4 0.24147 0.52052 0.46390 0.642723 v4 <--- F1
path5 -0.22484 0.43513 -0.51671 0.605360 v5 <--- F1
s1 0.37027 0.86496 0.42808 0.668594 v1 <--> v1
s2 0.91913 0.45215 2.03280 0.042073 v2 <--> v2
s3 0.54268 0.63075 0.86037 0.389584 v3 <--> v3
s4 0.94169 0.48496 1.94179 0.052162 v4 <--> v4
s5 0.94945 0.46562 2.03909 0.041441 v5 <--> v5
Iterations = 17
path1〜path5が因子負荷量。s1〜s5が残差の分散であり、1-sxで共通性が得られる。
ついでにパス図も画かせてみた。
pathDiagram(ans, ignore.double=FALSE, edge.labels="values", digits=3, node.font=c("C:\WINDOWS\Fonts\msmincho.ttc", 10))
上のコードでR上に出力されたdotファイルを、graphvizのGVeditで読み込み、描画したのが次の図である。
#ref(FAmodel.jpg)
**パス図 [#rcc6f04d]
pathDiagram
***graphviz による出力 [#g64877f2]
pathDiagramコマンドの出力したGraphviz用シンタックスを読み込んで描画。以下に例あり。
***dynamicGraph [#z62979ff]
Amos のようなことはできるか？
*グラフィカルモデリング [#efb871b8]
GMのソフトMIMは、CRANのパッケージmimRを使って
リモート操作できるようです。
-[[MIMのHP:http://www.hypergraph.dk/]]
-[[mimRのHP:http://gbi.agrsci.dk/~shd/public/mimR/]]
*Macintoshでの実行例 [#q59a14d4]
小島隆矢 「Excel で学ぶ共分散構造分析とグラフィカルモデリング」 Ohmsha の第4，5章で取り上げられているカメラの満足度についてのデータを分析してみた。~
図を日本語で描くことにする。
camera <- structure(list(小型軽量 = c(3, 5, 2, 4, 4, 5, 1,
1, 4, 4, 5, 2, 5, 1, 2, 4, 1, 5, 2, 2, 4, 3, 5, 5, 4, 3, 2, 1,
3, 2, 4, 4, 2, 5, 2, 1, 1, 2, 2, 2, 2, 3, 2, 4, 5, 1, 5, 1, 5,
2, 2, 4, 4, 1, 4, 2, 2, 3, 3, 2, 5, 3, 1, 3, 1, 4, 3, 4, 3, 5,
4, 2, 2, 3, 3, 2, 1, 5, 2, 2, 4, 1, 1, 3, 1, 1, 4, 4, 1, 2, 3,
2, 3, 5, 2, 1, 5, 2, 5, 1), 可搬性 = c(3, 4,
2, 4, 4, 4, 3, 2, 4, 3, 4, 2, 5, 2, 2, 3, 2, 5, 2, 2, 2, 3, 5,
4, 3, 2, 2, 1, 3, 2, 3, 4, 2, 4, 2, 1, 2, 2, 2, 3, 2, 3, 3, 4,
5, 2, 4, 1, 4, 1, 2, 4, 3, 2, 3, 2, 2, 1, 3, 3, 4, 4, 2, 3, 2,
4, 3, 4, 3, 3, 4, 1, 3, 3, 3, 2, 2, 4, 2, 3, 4, 1, 1, 2, 1, 1,
3, 4, 1, 2, 3, 3, 3, 5, 2, 3, 4, 2, 5, 1), 操作性 = c(4,
2, 2, 3, 2, 2, 3, 3, 5, 4, 2, 2, 1, 4, 2, 3, 5, 3, 3, 4, 4, 5,
1, 1, 1, 1, 5, 4, 1, 4, 3, 2, 2, 1, 4, 3, 3, 3, 3, 2, 2, 2, 3,
2, 1, 4, 1, 3, 3, 5, 4, 3, 3, 3, 3, 5, 3, 1, 3, 4, 2, 2, 3, 3,
2, 3, 3, 3, 1, 3, 4, 4, 3, 3, 3, 4, 3, 3, 3, 4, 1, 4, 2, 3, 3,
5, 1, 5, 3, 3, 2, 4, 3, 2, 5, 3, 2, 2, 2, 3), 満足度 = c(3,
2, 1, 2, 3, 1, 3, 1, 5, 3, 3, 1, 4, 1, 1, 2, 4, 5, 2, 2, 2, 4,
3, 1, 2, 1, 2, 2, 2, 3, 3, 3, 2, 2, 5, 1, 3, 1, 2, 4, 1, 2, 2,
2, 2, 5, 2, 1, 3, 3, 1, 4, 1, 2, 1, 3, 3, 1, 4, 4, 1, 3, 2, 2,
1, 4, 1, 2, 2, 3, 2, 3, 4, 3, 1, 2, 4, 4, 2, 4, 3, 1, 1, 3, 1,
3, 3, 4, 1, 1, 2, 3, 1, 4, 2, 2, 3, 1, 4, 2)), .Names = c("小型軽量",
"可搬性", "操作性", "満足度"
), class = "data.frame", row.names = 1:100)
r <- cor(camera)
r
model <- specifyModel()
小型軽量 -> 満足度, gamma, NA
操作性 -> 満足度, delta, NA
可搬性 -> 満足度, kappa, NA
小型軽量 -> 操作性, alpha, NA
小型軽量 -> 可搬性, beta, NA
小型軽量 <-> 小型軽量, theta1, NA
操作性 <-> 操作性, theta2, NA
可搬性 <-> 可搬性, theta3, NA
満足度 <-> 満足度, theta4, NA
ans <- sem(model, r, nrow(camera))
summary(ans)
pathDiagram(ans, "camera", ignore.double=FALSE, edge.labels="values", digits=3, node.font=c("Osaka", 10))
pathDiagram 関数で，ファイル "camera.dot" に書き出され（ファイル名は拡張子抜きで指定する），さらに "camera.pdf" に画像が描き出される。日本語フォントを指定すること。
分析結果は以下のごとし。
> r <- cor(camera)
> r
小型軽量 可搬性 操作性 満足度
小型軽量 1.0000000 0.8421738 -0.4114817 0.2123832
可搬性 0.8421738 1.0000000 -0.3721873 0.3844472
操作性 -0.4114817 -0.3721873 1.0000000 0.2728988
満足度 0.2123832 0.3844472 0.2728988 1.0000000
> model <- specify.model()
1: 小型軽量 -> 満足度, gamma, NA
2: 操作性 -> 満足度, delta, NA
3: 可搬性 -> 満足度, kappa, NA
4: 小型軽量 -> 操作性, alpha, NA
5: 小型軽量 -> 可搬性, beta, NA
6: 小型軽量 <-> 小型軽量, theta1, NA
7: 操作性 <-> 操作性, theta2, NA
8: 可搬性 <-> 可搬性, theta3, NA
9: 満足度 <-> 満足度, theta4, NA
10:
Read 9 records
> ans <- sem(model, r, nrow(camera))
> summary(ans)
Model Chisquare = 0.27002 Df = 1 Pr(>Chisq) = 0.60332
Chisquare (null model) = 185.63 Df = 6
Goodness-of-fit index = 0.99864
Adjusted goodness-of-fit index = 0.9864
RMSEA index = 0 90% CI: (NA, 0.21369)
Bentler-Bonnett NFI = 0.99855
Tucker-Lewis NNFI = 1.0244
Bentler CFI = 1
BIC = -4.3351
Normalized Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.41e-01 -1.36e-01 0.00e+00 -7.39e-02 0.00e+00 1.07e-15
Parameter Estimates
Estimate Std Error z value Pr(>|z|)
gamma -0.22909 0.151662 -1.5105 1.3092e-01 満足度 <--- 小型軽量
delta 0.45680 0.088107 5.1847 2.1641e-07 満足度 <--- 操作性
kappa 0.74739 0.148927 5.0185 5.2069e-07 満足度 <--- 可搬性
alpha -0.41148 0.091601 -4.4921 7.0521e-06 操作性 <--- 小型軽量
beta 0.84217 0.054192 15.5405 0.0000e+00 可搬性 <--- 小型軽量
theta1 1.00000 0.142162 7.0342 2.0037e-12 小型軽量 <--> 小型軽量
theta2 0.83068 0.118097 7.0339 2.0079e-12 操作性 <--> 操作性
theta3 0.29074 0.041353 7.0308 2.0537e-12 可搬性 <--> 可搬性
theta4 0.63666 0.090519 7.0334 2.0155e-12 満足度 <--> 満足度
Iterations = 0
> pathDiagram(ans, "camera", ignore.double=FALSE, edge.labels="values", digits=3, node.font=c("Osaka", 10))
Graphvis による図は以下のごとし（日本語で描けた）。
#ref(camera.png)
*lavaanパッケージの使用例 [#geb21fab]
[[上の例 (潜在変数を仮定するモデル) :http://www.okada.jp.org/RWiki/?R%A4%C7%B6%A6%CA%AC%BB%B6%B9%BD%C2%A4%CA%AC%C0%CF%A1%A6%B9%BD%C2%A4%CA%FD%C4%F8%BC%B0%A5%E2%A5%C7%A5%EB#rb420b01]]を[[lavaanパッケージ:http://lavaan.ugent.be/]]で実行する。
dat <- data.frame(matrix(c( # 上のものと同じデータ。型をdata.frameにする。
-1.89, -0.02, 0.42, 1.23, -1.53,
0.06, 1.81, -0.59, -0.75, -0.12,
2.58, -0.20, -1.92, -0.49, -0.35,
0.69, -0.66, -0.77, -1.92, 0.38,
-1.05, 0.07, -0.37, -0.89, -1.62,
-2.73, 1.40, 0.18, -0.09, 0.13,
0.95, 0.84, 1.19, 1.19, 0.10,
0.93, 1.17, -1.70, 1.46, -0.25,
-0.04, -0.12, -0.34, -0.24, 1.88,
0.16, 1.03, -0.05, -0.73, 0.04
), byrow=TRUE, ncol=5))
colnames(dat) <- c("v1","v2","v3","v4","v5")
library(lavaan)
model.lavaan <- '
F1 =~ v1+v2+v3+v4+v5
'
ans.lavaan <- cfa(model.lavaan, data=dat, std.lv=TRUE, mimic="EQS")
# mimicオプションを"EQS"とするとsemパッケージの結果と一致する。デフォルトでは"Mplus"
summary(ans.lavaan, fit.measures=TRUE, standardize = TRUE)
結果は以下のとおり。semパッケージの結果とほぼ一致
Lavaan (0.4-9) converged normally after 28 iterations
Number of observations 10
Estimator ML
Minimum Function Chi-square 1.925
Degrees of freedom 5
P-value 0.859
Chi-square test baseline model:
Minimum Function Chi-square 6.247
Degrees of freedom 10
P-value 0.794
Full model versus baseline model:
Comparative Fit Index (CFI) NaN
Tucker-Lewis Index (TLI) -0.639
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -70.837
Loglikelihood unrestricted model (H1) -69.768
Number of free parameters 10
Akaike (AIC) 161.675
Bayesian (BIC) 164.701
Root Mean Square Error of Approximation:
RMSEA 0.000
90 Percent Confidence Interval 0.000 0.235
P-value RMSEA <= 0.05 0.866
Standardized Root Mean Square Residual:
SRMR 0.087
Parameter estimates:
Information Expected
Standard Errors Standard
Estimate Std.err Z-value P(>|z|) Std.lv Std.all
Latent variables:
F1 =~
v1 -1.214 0.747 -1.625 0.104 -1.214 -0.794
v2 0.233 0.315 0.739 0.460 0.233 0.284
v3 0.631 0.422 1.494 0.135 0.631 0.676
v4 0.264 0.423 0.625 0.532 0.264 0.241
v5 -0.221 0.380 -0.582 0.561 -0.221 -0.225
Variances:
v1 0.866 1.550 0.559 0.576 0.866 0.370
v2 0.616 0.301 2.049 0.040 0.616 0.919
v3 0.473 0.463 1.020 0.308 0.473 0.543
v4 1.129 0.545 2.071 0.038 1.129 0.942
v5 0.917 0.441 2.079 0.038 0.917 0.949
F1 1.000 1.000 1.000
Machintoshでの実行例 (小島隆矢 「Excel で学ぶ共分散構造分析とグラフィカルモデリング」) も。結果は省略。
camera2 <- data.frame(camera) # データ入力は上のコードで。列名を半角英字にします。
colnames(camera2) <- c("koga", "kaso", "sosa", "mann")
library(lavaan)
model.camera <- '
sosa~koga
kaso~koga
mann~koga+sosa+kaso
'
ans.camera <- cfa(model.camera, data=camera2, std.lv=TRUE, mimic="EQS")
summary(ans.camera, fit.measures=TRUE, standardize = TRUE)
**記号の使い方 [#r4ffb3ee]
|"=~" | 潜在変数の定義 |
|"~" | 回帰 (単方向矢印 <-) |
| "~~" | 共変 (双方向矢印 <-> ) |
|"~ 1" | 切片 |
|"~ 0*" | 直交 (f1 ~ 0*f2 とするとf1とf2が直交すると仮定) |
*Windowsでの日本語フォント使用のTips [#tdaabb78]
Windowsでの日本語の利用方法がうまくいかず、ネットを調べていたら、ずばりこちらで解決されておりましたので、転載させて頂きました。(感謝)
で、Windowsでやる場合、
+フォントファイルを直接指定する(C:\WINDOWS\Fonts\msmincho.ttc)
+out.fileを指定せず、Rに結果を表示させ、エディタにコピペし、UTF-8で保存する。
> pathDiagram(ans, ignore.double=FALSE, edge.labels="values", digits=3, node.font=c("C:\WINDOWS\Fonts\msmincho.ttc", 10))
そしてGraphvizで処理することで、日本語が表示できるようです。(というか、できました。)
-- コメント：もし経常的に日本語表示をするということなら、pathDiagram.sem 関数を開き、その中の最初の方にある
handle <- file(out.file, "w")
を
handle <- file(out.file, "w", encoding="utf-8")
とでもしておけば、出力ファイルは utf-8 になります。と、思いますがいかが？
-[[RでSEM:http://d.hatena.ne.jp/kosugitti/20070528]]
-[[共分散構造分析:http://www.ic.nanzan-u.ac.jp/~kamiya/r/content/sem.html]]
*リンク [#e66c30b0]
-[[R のpackageSEM使用経験:http://www.nuis.ac.jp/~mat/fpr/fpr2003/0103.html]]
-[[Structural Equation Models - Appendix to An R and S-PLUS Companion to Applied Regression:http://socserv.socsci.mcmaster.ca/jfox/Books/Companion/appendix-sems.pdf]] --- sem の作者、John Fox のマニュアル(CRAN のものとは異なる)
-[[Ｒ言語による統計学:http://akimichi.homeunix.net/~emile/aki/medical/Stat/]] --- sem パッケージを使った分析の例があります。
-『Rによるやさしい統計学』 山田 剛史 杉澤 武俊 村井 潤一郎 共著 に解説あり
-[[SEM による確定型因子分析の例 in "Structual Equation Modeling ":http://www.statmethods.net/advstats/factor.html]]
-[[RでSEM:http://sites.google.com/site/officeoga/r/rde-sem]]
-[[RでSEM入門 SEMは（そんなに）怖くない:http://m884.jp/Osaka.R/_media/osaka.r_no.3/r%E3%81%A7sem%E5%85%A5%E9%96%80.pdf]]
-[[lavaan - latent variable analysis:http://lavaan.ugent.be/?q=node/14]]
-[[主成分分析と因子分析:http://eau.uijin.com/advstats/factor.html]]
-[[semdiag - Structural equation modeling diagnostics:http://www.rforge.net/semdiag/]]
----
-お分かりの方、どうぞ、どんどん書いちゃってください -- &new{2004-05-30 (日) 07:18:58};
- dotファイルからjpegを作るには、dot.exeにパスを通してから、system("dot camera.dot -Tjpg -o camera.jpg")とでもすればよいでしょう。 -- &new{2008-07-15 (火) 17:54:05};
- ついでに、system('dot camera.dot -Tjpg -o "camera.jpg" -Glabel="GFI: 0.95"')のようにすると、図の下に適合度とかの統計量を埋め込めます。 -- &new{2008-07-15 (火) 18:05:07};
- BioconductorのパッケージRgraphvizを使えば、Graphvizと協調して、Rの世界だけでパス図が書けますが、できばえはいまひとつでした。 -- &new{2008-08-04 (月) 18:42:56};
- pathDiagramのout.fileの拡張子を.dotとしたら，MacのGraphviz 2.20.2で上手く行きました。.datの拡張子ではGraphvizで参照できません。 -- &new{2008-08-07 (木) 14:57:03};
- ＞ pathDiagramのout.fileの拡張子を.dotとしたら~
ファイルの属性をテキストファイルにするとドラッグ＆ドロップで立ち上がる（クリエーターはなんでも良いが，よく使うエディタなどを指定すると良いだろう）。また，File-->Open で選択することもできる。当然ながらファイルをダブルクリックすると，クリエータに指定したエディタがファイルを開く。~
クリエータに "GRVZ" を指定しても良いが。
ans <- sem(model, r, nrow(camera))
summary(ans)
pathDiagram(ans, out.file="camera.dat", ignore.double=FALSE,
edge.labels="values", digits=3, node.font=c("Osaka", 10))
system("/Developer/Tools/SetFile -t TEXT -c JEDX camera.dat")
こんな風に。 -- &new{2008-08-07 (木) 15:35:08};
- 修正指数の算出の仕方 -- [[覚え書き]] &new{2009-11-16 (月) 23:22:57};
ans <- sem(model, r, nrow(data))
mod.indices(ans)
- OpenMxはどうかな。 -- &new{2009-11-17 (火) 18:40:34};
- [[OpenMx - Advanced Structural Equation Modeling:http://openmx.psyc.virginia.edu/]] -- &new{2009-11-17 (火) 18:41:10};
- pathDiagramで、standardize=TRUE のオプション。標準化係数を出力。 -- &new{2009-12-05 (土) 15:39:47};
- [[「食物とガン」のデータをもとにした３種類の共分散構造分析:http://www.ic.nanzan-u.ac.jp/~kamiya/r/content/sem.html]] -- &new{2009-12-05 (土) 15:44:43};
- lavaanパッケージのlavaan()コマンドで多母集団が分析できます。 -- &new{2011-07-01 (金) 15:10:22};
- 現在のバージョンでは，pathDiagram の第二引数は名前を持たない。したがって，out.file="camera.dot" とすると期待通りの動作（描画）をしない。pathDiagram(ans, "camera", ...) などとすべし（ファイル名は拡張子抜きで指定する。"foo" を指定すると "foo.dot" というファイルができる）。~
Mac版では，"foo.dot" ファイルを出力した後 "foo.pdf" が出力される。 -- [[]] &new{2011-07-05 (火) 16:28:39};
- lavaanパッケージについて書いてみました。ご意見や修正等ありましたらお願いします。 -- [[aor]] &new{2011-07-11 (月) 21:42:32};
// 「可搬性」は，「かはんせい」と読むのですけど？kaso とはどういう省略？
#comment
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
