![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
Rで共分散構造分析・構造方程式モデル
library(sem)
cov(data)で分散共分散行列を、cor(data)で相関行列を計算できる。
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(モデル,データ行列,サンプル数)
で推定する。モデルはspecifyModel()関数で指定したもの、データ行列は分散共分散行列か相関行列を与える。
summary()
で推定値、適合度等の結果を返す。
下のMacintoshでの実行例を参照。
青木先生の因子分析を行うRコードを例に、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で読み込み、描画したのが次の図である。
pathDiagram
pathDiagramコマンドの出力したGraphviz用シンタックスを読み込んで描画。以下に例あり。
Amos のようなことはできるか？
GMのソフトMIMは、CRANのパッケージmimRを使って リモート操作できるようです。
小島隆矢 「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 による図は以下のごとし（日本語で描けた）。
上の例 (潜在変数を仮定するモデル) をlavaanパッケージで実行する。
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)
| "=~" | 潜在変数の定義 |
| "~" | 回帰 (単方向矢印 <-) |
| "~~" | 共変 (双方向矢印 <-> ) |
| "~ 1" | 切片 |
| "~ 0*" | 直交 (f1 ~ 0*f2 とするとf1とf2が直交すると仮定) |
Windowsでの日本語の利用方法がうまくいかず、ネットを調べていたら、ずばりこちらで解決されておりましたので、転載させて頂きました。(感謝)
で、Windowsでやる場合、
> pathDiagram(ans, ignore.double=FALSE, edge.labels="values", digits=3, node.font=c("C:\WINDOWS\Fonts\msmincho.ttc", 10))
そしてGraphvizで処理することで、日本語が表示できるようです。(というか、できました。)
handle <- file(out.file, "w")を
handle <- file(out.file, "w", encoding="utf-8")とでもしておけば、出力ファイルは utf-8 になります。と、思いますがいかが？
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")
こんな風に。 -- 2008-08-07 (木) 15:35:08ans <- sem(model, r, nrow(data)) mod.indices(ans)
FAmodel.jpg
camera.png
