R

sem_多重指標モデル

sem,qgraph,Rgraphviz

(参考)
はやしさんのブログ Rev.3
名城大学人間学部神谷研究室

(準備)　graphvizをインストールする

1
2
3
4
5
6
7
8
9
10
library(sem)
#library(semPlot)
library(qgraph)
library(knitr)
#blue backs「原因をさぐる統計学」p.88 データ
data <- read.table("http://wwwhum.meijo-u.ac.jp/labs/hh002/r/rfile/gan.csv", header=TRUE, sep=",", na.strings="NA", dec=".", 
strip.white=TRUE)
names(data)<-c("総熱量","肉類","乳製品","酒類","大腸ガン","直腸ガン")
soukan<-cor(data)
soukan
総熱量 肉類 乳製品 酒類 大腸ガン 直腸ガン
総熱量 1.0000000 0.7731109 0.7155250 0.6341241 0.7394903 0.8100869
肉類 0.7731109 1.0000000 0.7227997 0.4515933 0.8526245 0.6718607
乳製品 0.7155250 0.7227997 1.0000000 0.2321334 0.5787253 0.4785180
酒類 0.6341241 0.4515933 0.2321334 1.0000000 0.5881328 0.6095649
大腸ガン 0.7394903 0.8526245 0.5787253 0.5881328 1.0000000 0.7523879
直腸ガン 0.8100869 0.6718607 0.4785180 0.6095649 0.7523879 1.0000000 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#reference.indicators=T:1つめの測定変数のパラメータを1に固定
model.a1 <- cfa(reference.indicators=T)
洋食傾向: 総熱量, 肉類, 乳製品, 酒類
model.a2 <- cfa(reference.indicators=T)
消化管ガン: 大腸ガン, 直腸ガン
model.a3 <- specifyEquations()
消化管ガン = g21*洋食傾向
#モデルを結合
model.a<- combineModels(model.a1, model.a2, model.a3)
#標準的なsemmodクラスに変換
model <- removeRedundantPaths(model.a)
fit <- sem:::sem(model, soukan, N=nrow(data))
kable(stdCoef(fit))
Std. Estimate
lam[総熱量:洋食傾向] 0.9035308 総熱量 <—- 洋食傾向
lam[肉類:洋食傾向] 0.8886019 肉類 <—- 洋食傾向
lam[乳製品:洋食傾向] 0.7205247 乳製品 <—- 洋食傾向
lam[酒類:洋食傾向] 0.6199373 酒類 <—- 洋食傾向
1.0000000 洋食傾向 <—> 洋食傾向
V[総熱量] 0.1836320 総熱量 <—> 総熱量
V[肉類] 0.2103866 肉類 <—> 肉類
V[乳製品] 0.4808442 乳製品 <—> 乳製品
V[酒類] 0.6156778 酒類 <—> 酒類
0.8959051 大腸ガン <—- 消化管ガン
lam[直腸ガン:消化管ガン] 0.8398075 直腸ガン <—- 消化管ガン
V[消化管ガン] 0.0385689 消化管ガン <—> 消化管ガン
V[大腸ガン] 0.1973540 大腸ガン <—> 大腸ガン
V[直腸ガン] 0.2947234 直腸ガン <—> 直腸ガン
g21 0.9805259 消化管ガン <—- 洋食傾向 
1
summary(fit)

Model Chisquare = 46.91723 Df = 8 Pr(>Chisq) = 1.590119e-07
AIC = 72.91723
BIC = 16.11605

Normalized Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.3290000 -0.3108000 -0.0000003 -0.0440900 0.3667000 0.5985000

R-square for Endogenous Variables
総熱量 肉類 乳製品 酒類 消化管ガン 大腸ガン
0.8164 0.7896 0.5192 0.3843 0.9614 0.8026
直腸ガン
0.7053

Parameter Estimates
Estimate Std Error z value Pr(>|z|)
lam[総熱量:洋食傾向] 0.90353060 0.11570035 7.8092299 5.753844e-15
lam[肉類:洋食傾向] 0.88860195 0.11697802 7.5963158 3.046806e-14
lam[乳製品:洋食傾向] 0.72052466 0.12999089 5.5428856 2.975272e-08
lam[酒類:洋食傾向] 0.61993740 0.13612749 4.5540942 5.261177e-06
V[総熱量] 0.18363196 0.05628021 3.2628159 1.103112e-03
V[肉類] 0.21038662 0.05992060 3.5110903 4.462728e-04
V[乳製品] 0.48084422 0.10824543 4.4421666 8.905755e-06
V[酒類] 0.61567805 0.13405385 4.5927668 4.374078e-06
lam[直腸ガン:消化管ガン] 0.93738433 0.12049701 7.7793162 7.291758e-15
V[消化管ガン] 0.03095721 0.05846916 0.5294622 5.964848e-01
V[大腸ガン] 0.19735405 0.06694604 2.9479572 3.198814e-03
V[直腸ガン] 0.29472338 0.07700310 3.8274224 1.294922e-04
g21 0.87845827 0.11829002 7.4263092 1.116698e-13

lam[総熱量:洋食傾向] 総熱量 <—- 洋食傾向
lam[肉類:洋食傾向] 肉類 <—- 洋食傾向
lam[乳製品:洋食傾向] 乳製品 <—- 洋食傾向
lam[酒類:洋食傾向] 酒類 <—- 洋食傾向
V[総熱量] 総熱量 <—> 総熱量
V[肉類] 肉類 <—> 肉類
V[乳製品] 乳製品 <—> 乳製品
V[酒類] 酒類 <—> 酒類
lam[直腸ガン:消化管ガン] 直腸ガン <—- 消化管ガン
V[消化管ガン] 消化管ガン <—> 消化管ガン
V[大腸ガン] 大腸ガン <—> 大腸ガン
V[直腸ガン] 直腸ガン <—> 直腸ガン
g21 消化管ガン <—- 洋食傾向 

1
qgraph(fit,edge.labels=T,layout = "tree",vsize.man = 6,vsize.lat = 8,fade=F)

1
2
3
library(Rgraphviz)
path.diagram(fit,edge.labels="values",ignore.double=T,standardize=T,digits=2,file="sem1",graphics.fmt ="png",
min.rank="総熱量, 肉類, 乳製品, 酒類",max.rank="大腸ガン, 直腸ガン",same.rank=NULL)

従来の記法

modelA <- specifyModel()

model.c1

洋食傾向 -> 総熱量, b11, NA
洋食傾向 -> 肉類, b21, NA
洋食傾向 -> 乳製品, b31, NA
洋食傾向 -> 酒類, b41, NA
洋食傾向 <-> 洋食傾向, NA,1

model.c2

消化管ガン -> 大腸ガン, NA, 1
消化管ガン -> 直腸ガン, b62, NA
消化管ガン <-> 消化管ガン, delta2, NA

model.c3

洋食傾向 -> 消化管ガン, g21, NA

表記の必要無し

総熱量 <-> 総熱量, e1, NA
肉類 <-> 肉類, e2, NA
乳製品 <-> 乳製品, e3, NA
酒類 <-> 酒類, e4, NA
大腸ガン <-> 大腸ガン, e5, NA
直腸ガン <-> 直腸ガン, e6, NA