Rで共分散構造分析・構造方程式モデル

パッケージ

  • sem(構造方程式モデル)パッケージ中のオブジェクト一覧?

使用方法

パッケージのロード

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の実行

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で読み込み、描画したのが次の図である。

FAmodel.jpg

パス図

pathDiagram

graphviz による出力

pathDiagramコマンドの出力したGraphviz用シンタックスを読み込んで描画。以下に例あり。

dynamicGraph

Amos のようなことはできるか?

グラフィカルモデリング

GMのソフトMIMは、CRANのパッケージmimRを使って リモート操作できるようです。

Macintoshでの実行例

小島隆矢 「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 による図は以下のごとし(日本語で描けた)。

camera.png

lavaanパッケージの使用例

上の例 (潜在変数を仮定するモデル) 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での日本語フォント使用のTips

Windowsでの日本語の利用方法がうまくいかず、ネットを調べていたら、ずばりこちらで解決されておりましたので、転載させて頂きました。(感謝)

で、Windowsでやる場合、

  1. フォントファイルを直接指定する(C:\WINDOWS\Fonts\msmincho.ttc)
  2. 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によるやさしい統計学』 山田 剛史 杉澤 武俊 村井 潤一郎 共著 に解説あり

  • お分かりの方、どうぞ、どんどん書いちゃってください -- 2004-05-30 (日) 07:18:58
  • dotファイルからjpegを作るには、dot.exeにパスを通してから、system("dot camera.dot -Tjpg -o camera.jpg")とでもすればよいでしょう。 -- 2008-07-15 (火) 17:54:05
  • ついでに、system('dot camera.dot -Tjpg -o "camera.jpg" -Glabel="GFI: 0.95"')のようにすると、図の下に適合度とかの統計量を埋め込めます。 -- 2008-07-15 (火) 18:05:07
  • BioconductorのパッケージRgraphvizを使えば、Graphvizと協調して、Rの世界だけでパス図が書けますが、できばえはいまひとつでした。 -- 2008-08-04 (月) 18:42:56
  • pathDiagramのout.fileの拡張子を.dotとしたら,MacのGraphviz 2.20.2で上手く行きました。.datの拡張子ではGraphvizで参照できません。 -- 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")
    こんな風に。 -- 2008-08-07 (木) 15:35:08
  • 修正指数の算出の仕方 -- 覚え書き? 2009-11-16 (月) 23:22:57
    ans <- sem(model, r, nrow(data))
    mod.indices(ans)
  • OpenMx?はどうかな。 -- 2009-11-17 (火) 18:40:34
  • OpenMx - Advanced Structural Equation Modeling -- 2009-11-17 (火) 18:41:10
  • pathDiagramで、standardize=TRUE のオプション。標準化係数を出力。 -- 2009-12-05 (土) 15:39:47
  • 「食物とガン」のデータをもとにした3種類の共分散構造分析 -- 2009-12-05 (土) 15:44:43
  • lavaanパッケージのlavaan()コマンドで多母集団が分析できます。 -- 2011-07-01 (金) 15:10:22
  • 現在のバージョンでは,pathDiagram の第二引数は名前を持たない。したがって,out.file="camera.dot" とすると期待通りの動作(描画)をしない。pathDiagram(ans, "camera", ...) などとすべし(ファイル名は拡張子抜きで指定する。"foo" を指定すると "foo.dot" というファイルができる)。
    Mac版では,"foo.dot" ファイルを出力した後 "foo.pdf" が出力される。 -- [[]] 2011-07-05 (火) 16:28:39
  • lavaanパッケージについて書いてみました。ご意見や修正等ありましたらお願いします。 -- aor? 2011-07-11 (月) 21:42:32


添付ファイル: fileFAmodel.jpg 1884件 [詳細] filecamera.png 2197件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (904d)