R

共分散行列を利用4

lavaan、semPlot パッケージ　

(参考)blue backs「原因をさぐる統計学」
R-SVG.org
SVG画像

原因を探る統計学

p.210 社会階層モデルの分析 図表 補-9 p.245
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
library(lavaan)
library(semPlot)
mat<-'
1 
0.4703 1
0.4609 0.2655 1
0.2988 0.21 0.4392 1
0.2335 0.2019 0.1917 0.1462 1
0.1172 0.1262 0.0946 0.0564 0.5581 1
0.2782 0.2381 0.2747 0.2359 0.3597 0.325 1
0.1571 0.1455 0.1514 0.0942 0.3352 0.4242 0.4337 1
0.2512 0.2083 0.2517 0.2337 0.2676 0.2723 0.4076 0.3291 1
0.2684 0.2311 0.2657 0.2392 0.4099 0.4406 0.5266 0.4329 0.4761 1
0.2576 0.2123 0.2779 0.2274 0.4294 0.4382 0.5048 0.4089 0.4482 0.7457 1
0.2901 0.2464 0.2842 0.2248 0.4291 0.4812 0.4837 0.3857 0.4342 0.612 0.5634 1
0.2983 0.272 0.2919 0.2205 0.4455 0.5005 0.4592 0.3932 0.4452 0.6455 0.5574 0.845 1
0.2856 0.2136 0.3096 0.2439 0.3851 0.3966 0.4096 0.3426 0.374 0.4822 0.5176 0.6345 0.645 1
0.2521 0.1892 0.2838 0.2275 0.3525 0.343 0.3618 0.297 0.3528 0.4139 0.4674 0.5424 0.5626 0.68 1
'
wheaton.cov <- getCov(mat,names=c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12","x13","x14","x15"))
model<-'
F1=~x1+x2+x3+x4
F2=~x5+x6
F3=~x7+x8+x9
F4=~x10+x11
F5=~x12+x13
F6=~x14+x15
F1~~F2
F3~F1+F2
F4~F2+F3
F5~F4
F6~F5
'
fit <- lavaan:::sem(model, sample.cov=wheaton.cov, sample.nobs=2038)
summary(fit, standardized=TRUE)
manifests<-c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12","x13","x14","x15")
latents<-c("F1","F2","F3","F4","F5","F6")
labels=c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12","x13","x14","x15","親の\n社会経済的\n地位","本人の能力","周囲の励まし","教育的・\n職業的向上心","　学　歴　","職業的地位")
semPaths(fit,"std",fade=F,manifests=manifests,latents=latents,nodeLabels=labels,style="lisrel")

graphvizを使った図

p.186 双方向の因果関係 図表 補-6 p.241
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
library(lavaan)
library(semPlot)
mat<-'
1
0.56 1
0.06 0.05 1
0.16 0.21 0.32 1
0.01 -0.04 0.1 -0.06 1
-0.07 -0.05 0.16 -0.07 0.42 1
-0.02 -0.01 0.14 0.08 0.18 0.31 1
0.17 0.3 0.29 0.4 0.01 0.13 0.21 1
0.16 0.21 0.28 0.19 0.12 0.27 0.27 0.5 1
0.05 0.04 0.08 0.13 0.07 0.15 0.25 0.28 0.24 1
0.1 0.1 0.09 0.17 0.02 0.08 0.08 0.23 0.18 0.59 1
0.1 0.17 0.14 0.17 0.08 0.17 0.33 0.32 0.4 0.55 0.49 1
'
wheaton.cov <- getCov(mat,names=c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12"))
model<-'
F1=~x1+x2
F2=~x3+x4
F3=~x5+x6+x7
F4=~x8+x9
F5=~x10+x11+x12
F1~~F2
F1~~F3
F2~~F3
F4~F1+F2
F5~F3
F5~F4
F4~F5
'
fit <- lavaan:::sem(model, sample.cov=wheaton.cov, sample.nobs=249)
summary(fit, standardized=TRUE)
manifests<-c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12")
latents<-c("F1","F2","F3","F4","F5")
labels=c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10","x11","x12","親の\n社会経済的\n地位","　知　能　","大人からの\n評価","学業成績","友達からの\n評価")
#png("graph.png", family="Takao P明朝", width=1000, height=800)
semPaths(fit,"std",fade=F,manifests=manifests,latents=latents,nodeLabels=labels,style="lisrel",gray=T,edge.label.cex=0.8,esize=TRUE,fixedStyle=c("black",1), freeStyle=c("black",1))
#dev.off()

  • 学業成績と友達からの評価の双方向の因果関係を表すエッジが共分散のような両矢印になる。

グラフをsvg形式で出力し、Inkscapeを使って重なった部分を分離する。

RでグラフをSVG形式で描くには

  • svg()関数（文字も図形になり，ビューアに依存しない）
  • RSVGTipsDevice:::devSVGTips関数（UTF-8の文字列を書き込むので後で変更可）

ここではsvg()関数を使う。

1
2
3
4
#標準（cairo）の svg() を使う(日本語のフォント名を指定する必要があり)
svg("graph.svg", family="Takao P明朝", width=12, height=10)
semPaths(fit,"std",fade=F,manifests=manifests,latents=latents,nodeLabels=labels,style="lisrel",gray=T,edge.label.cex=0.8,esize=TRUE,fixedStyle=c("black",1), freeStyle=c("black",1))
dev.off()

Inkscapeを使って重なった部分を上下に分離する。
-0.06を赤字にする。

graphvizを使った図（手を加えても綺麗な図にならなかった）

##### p.148 因果モデルを作る 図表 補-4 p.238

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
library(lavaan)
library(semPlot)
mat<-'
1
0.36 1 
0.21 0.32 1
0.08 0.15 0.2 1
0.09 0.13 0.18 0.69 1
0.13 0.2 0.23 0.41 0.38 1
0.04 0.06 0.1 0.44 0.39 0.29 1
0.07 0.13 0.15 0.38 0.39 0.3 0.32 1
-0.05 0.02 0.08 0.26 0.22 0.14 0.21 0.19 1
-0.00 0.07 0.13 0.39 0.35 0.24 0.35 0.3 0.53 1
'
wheaton.cov <- getCov(mat,names=c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10"))
model<-'
F1=~x1+x2+x3
F2=~x4+x5+x6+x7+x8
F3=~x9+x10
F2~F1
F3~F2
'
fit <- lavaan:::sem(model, sample.cov=wheaton.cov, sample.nobs=2560)
summary(fit, standardized=TRUE)
manifests<-c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10")
latents<-c("F1","F2","F3")
labels=c("x1","x2","x3","x4","x5","x6","x7","x8","x9","x10"," 学歴志向 ","偏差値に対\nする信頼"," 受験産業 \nへの依存")
#png("graph.png", family="Takao P明朝", width=1000, height=800)
semPaths(fit,"std",fade=F,manifests=manifests,latents=latents,nodeLabels=labels,style="lisrel",gray=T,
#edge.label.cex=0.8,esize=TRUE,fixedStyle=c("black",1), freeStyle=c("black",1))
dev.off()

1
2
3
4
5
6
#latents<-c("F1","F2","F3")
makeawk(latents)
pars <- parameterEstimates(fit, standardized = TRUE) 
g<-data.frame(pars$lhs,pars$op,pars$rhs," [label=",signif(pars$std.all,digits=3),"];")
write.table(g,"out",sep=" ",col.names=FALSE,row.names=FALSE,quote=FALSE, na="")
system( 'gawk -f "lavaan_dot.awk" "out">out.dot')

dotファイル編集

1
system('dot -Tpng out.dot -o out.png')