理学療法士がまとめたノート

統計学備忘録(R言語のメモ)

since2016 ときどきTEXのメモ

共分散構造分析(パス図の描き方)

Rを使ったパス図作成の方法を忘れないうちに簡単に書いておきます
青木先生のデータを借用しまして勉強していきます.
R -- 因子分析(factanal を援用する)

dat <- matrix(c(   
  -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")    #列名を加えます
x <- dat

相関行列、共分散行列

#後でsemパッケージのために必要になります
r <- cor(x)      # 相関行列
cov <- cov(x)    # 共分散行列

ここはモデルをイメージするところで、色々なパターンが考えられます
RjpWikiに習って、下記のようなイメージで作業を進めて行きます.

f:id:yoshida931:20180301180850p:plain:w300

library(sem) 

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

#終了したらCtrl+R

いよいよsemパケージを使用します

相関行列を使用したSEM

library(sem)       # semパッケージ
ans1 <- sem(model, r, N=10)    #相関行列を使用してsemを実行
summary(ans1)

 Model Chisquare =  1.925355   Df =  5 Pr(>Chisq) = 0.8593742
 AIC =  21.92536
 BIC =  -9.58757

 Normalized Residuals
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.3976000 -0.0000002  0.0155300  0.1221000  0.2435000  0.6971000 

 R-square for Endogenous Variables
    v1     v2     v3     v4     v5 
0.6297 0.0809 0.4573 0.0583 0.0506 

 Parameter Estimates
      Estimate   Std Error z value    Pr(>|z|)             
path1  0.7935558 0.4882341  1.6253592 0.10408604 v1 <--- F1
path2 -0.2843858 0.3849962 -0.7386717 0.46010637 v2 <--- F1
path3 -0.6762525 0.4525192 -1.4944171 0.13506663 v3 <--- F1
path4 -0.2414672 0.3860898 -0.6254172 0.53169727 v4 <--- F1
path5  0.2248359 0.3864879  0.5817413 0.56074096 v5 <--- F1
s1     0.3702693 0.6626860  0.5587401 0.57633909 v1 <--> v1
s2     0.9191248 0.4485423  2.0491376 0.04044867 v2 <--> v2
s3     0.5426826 0.5322169  1.0196644 0.30788765 v3 <--> v3
s4     0.9416935 0.4546109  2.0714275 0.03831887 v4 <--> v4
s5     0.9494489 0.4567576  2.0786712 0.03764758 v5 <--> v5

 Iterations =  15 

#Estimate の部分が非標準化係数ですが、相関行列から求めているので不正確

標準化係数

stdCoef(ans1)   # 標準化係数

            Std. Estimate           
path1 path1     0.7935557 v1 <--- F1
path2 path2    -0.2843858 v2 <--- F1
path3 path3    -0.6762525 v3 <--- F1
path4 path4    -0.2414672 v4 <--- F1
path5 path5     0.2248359 v5 <--- F1
s1       s1     0.3702693 v1 <--> v1
s2       s2     0.9191247 v2 <--> v2
s3       s3     0.5426826 v3 <--> v3
s4       s4     0.9416936 v4 <--> v4
s5       s5     0.9494488 v5 <--> v5
                1.0000000 F1 <--> F1

共分散行列を使用したSEM

ans2 <- sem(model, cov, N=284) # 非標準化係数
summary(ans2)

 Model Chisquare =  60.54173   Df =  5 Pr(>Chisq) = 9.392241e-12
 AIC =  80.54173
 BIC =  32.29686

 Normalized Residuals
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-2.230000 -0.000001  0.087080  0.684600  1.365000  3.909000 

 R-square for Endogenous Variables
    v1     v2     v3     v4     v5 
0.6297 0.0809 0.4573 0.0583 0.0506 

 Parameter Estimates
      Estimate   Std Error  z value   Pr(>|z|)               
path1  1.2135217 0.13314543  9.114257 7.921129e-20 v1 <--- F1
path2 -0.2329064 0.05622870 -4.142127 3.440994e-05 v2 <--- F1
path3 -0.6310501 0.07530436 -8.379995 5.293016e-17 v3 <--- F1
path4 -0.2643794 0.07538517 -3.507048 4.531080e-04 v4 <--- F1
path5  0.2209609 0.06773505  3.262135 1.105764e-03 v5 <--- F1
s1     0.8658805 0.27636060  3.133155 1.729383e-03 v1 <--> v1
s2     0.6164834 0.05365106 11.490610 1.470704e-30 v2 <--> v2
s3     0.4725590 0.08264694  5.717805 1.079092e-08 v3 <--> v3
s4     1.1288825 0.09718675 11.615601 3.433641e-31 v4 <--> v4
s5     0.9170030 0.07867069 11.656221 2.133019e-31 v5 <--> v5

#Estimate ここは間違いなく非標準化係数です

パス図の作成

install.packages("DiagrammeR")
install.packages("scales")
pathDiagram(ans1, ignore.double=FALSE, edge.labels="values", digits=3, node.font=c("HGRGE", 10))
pathDiagram(ans2, ignore.double=FALSE, edge.labels="values", digits=3, node.font=c("HGRGE", 10))
f:id:yoshida931:20180301175623p:plain:w250\hspace{2cm} f:id:yoshida931:20180302084342p:plain:w250

lavaanパッケージによる方法

RjpWikiより

installed.packages("lavaan")
library("lavaan")
model2 <- 'f1 =~ v1 + v2 + v3 + v4 + v5
'
ans3 <- cfa(model=model2, data=x, std.lv=TRUE, mimic="EQS")    # mimicオプションを"EQS"とするとsemパッケージの結果と一致する。デフォルトでは"Mplus"
summary(ans3, fit.measures=TRUE, standardize = TRUE)
standardizedSolution(object =ans3 )# 標準化推定値

   lhs op rhs est.std se  z pvalue
1   f1 =~  v1  -0.794 NA NA     NA
2   f1 =~  v2   0.284 NA NA     NA
3   f1 =~  v3   0.676 NA NA     NA
4   f1 =~  v4   0.241 NA NA     NA
5   f1 =~  v5  -0.225 NA NA     NA
6   f1  ~  f1  -3.576 NA NA     NA
7   v1 ~~  v1   0.370 NA NA     NA
8   v2 ~~  v2   0.919 NA NA     NA
9   v3 ~~  v3   0.543 NA NA     NA
10  v4 ~~  v4   0.942 NA NA     NA
11  v5 ~~  v5   0.949 NA NA     NA
12  f1 ~~  f1   1.000 NA NA     NA

参考
Rを使った分析(SEM) | 外国語教育研究ハンドブック