Feb 14 2015

sem_MIMICモデル

sem,qgraph,Rgraphviz

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

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)
#covsオプションに"変数1, 変数2"という形で共分散を指定
model.b1 <- specifyEquations(
 covs="総熱量, 肉類,乳製品,酒類")
消化管ガン = gam1*総熱量 + gam2*肉類+ gam3*乳製品 + gam4*酒類
model.b2 <- cfa(reference.indicators=T)
消化管ガン: 大腸ガン, 直腸ガン
#モデルを結合
model.b<- combineModels(model.b1, model.b2)
#標準的なsemmodクラスに変換
model <- removeRedundantPaths(model.b)
fit <- sem:::sem(model, soukan, N=nrow(data))
kable(stdCoef(fit))

		Std. Estimate
gam1	gam1	0.3312691	消化管ガン <—- 総熱量
gam2	gam2	0.6469231	消化管ガン <—- 肉類
gam3	gam3	-0.1336175	消化管ガン <—- 乳製品
gam4	gam4	0.2036032	消化管ガン <—- 酒類
V[総熱量]	V[総熱量]	1.0000000	総熱量 <—> 総熱量
C[総熱量,肉類]	C[総熱量,肉類]	0.7731109	肉類 <—> 総熱量
C[総熱量,乳製品]	C[総熱量,乳製品]	0.7155250	乳製品 <—> 総熱量
C[総熱量,酒類]	C[総熱量,酒類]	0.6341240	酒類 <—> 総熱量
V[肉類]	V[肉類]	1.0000000	肉類 <—> 肉類
C[肉類,乳製品]	C[肉類,乳製品]	0.7227997	乳製品 <—> 肉類
C[肉類,酒類]	C[肉類,酒類]	0.4515934	酒類 <—> 肉類
V[乳製品]	V[乳製品]	1.0000000	乳製品 <—> 乳製品
C[乳製品,酒類]	C[乳製品,酒類]	0.2321334	酒類 <—> 乳製品
V[酒類]	V[酒類]	1.0000000	酒類 <—> 酒類
V[消化管ガン]	V[消化管ガン]	0.0775069	消化管ガン <—> 消化管ガン
		0.9094601	大腸ガン <—- 消化管ガン
lam[直腸ガン:消化管ガン]	lam[直腸ガン:消化管ガン]	0.8272907	直腸ガン <—- 消化管ガン
V[大腸ガン]	V[大腸ガン]	0.1728824	大腸ガン <—> 大腸ガン
V[直腸ガン]	V[直腸ガン]	0.3155901	直腸ガン <—> 直腸ガン

1	summary(fit)

Model Chisquare = 17.34207 Df = 3 Pr(>Chisq) = 0.0006010001
AIC = 53.34207
BIC = 5.791629

Normalized Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.3885000 -0.0000001 0.0000000 0.0068180 0.0000005 0.5215000

R-square for Endogenous Variables
消化管ガン大腸ガン直腸ガン
0.9225 0.8271 0.6844

Parameter Estimates
Estimate Std Error z value
gam1 0.30127602 0.13426489 2.243893
gam2 0.58835070 0.11056024 5.321539
gam3 -0.12151977 0.10741965 -1.131262
gam4 0.18516894 0.09014137 2.054206
V[総熱量] 1.00000001 0.20851442 4.795832
C[総熱量,肉類] 0.77311090 0.18636689 4.148327
C[総熱量,乳製品] 0.71552499 0.18129822 3.946674
C[総熱量,酒類] 0.63412405 0.17458730 3.632132
V[肉類] 1.00000001 0.20851442 4.795832
C[肉類,乳製品] 0.72279968 0.18192446 3.973076
C[肉類,酒類] 0.45159337 0.16177928 2.791417
V[乳製品] 0.99999999 0.20851441 4.795832
C[乳製品,酒類] 0.23213339 0.15136236 1.533627
V[酒類] 1.00000000 0.20851441 4.795832
V[消化管ガン] 0.06410736 0.04949189 1.295310
lam[直腸ガン:消化管ガン] 0.90965028 0.11431250 7.957575
V[大腸ガン] 0.17288239 0.05975161 2.893351
V[直腸ガン] 0.31559004 0.07671434 4.113834
Pr(>|z|)
gam1 2.483929e-02 消化管ガン <—- 総熱量
gam2 1.028928e-07 消化管ガン <—- 肉類
gam3 2.579448e-01 消化管ガン <—- 乳製品
gam4 3.995577e-02 消化管ガン <—- 酒類
V[総熱量] 1.620014e-06 総熱量 <—> 総熱量
C[総熱量,肉類] 3.349131e-05 肉類 <—> 総熱量
C[総熱量,乳製品] 7.924429e-05 乳製品 <—> 総熱量
C[総熱量,酒類] 2.810896e-04 酒類 <—> 総熱量
V[肉類] 1.620014e-06 肉類 <—> 肉類
C[肉類,乳製品] 7.095042e-05 乳製品 <—> 肉類
C[肉類,酒類] 5.247787e-03 酒類 <—> 肉類
V[乳製品] 1.620014e-06 乳製品 <—> 乳製品
C[乳製品,酒類] 1.251215e-01 酒類 <—> 乳製品
V[酒類] 1.620014e-06 酒類 <—> 酒類
V[消化管ガン] 1.952132e-01 消化管ガン <—> 消化管ガン
lam[直腸ガン:消化管ガン] 1.754438e-15 直腸ガン <—- 消化管ガン
V[大腸ガン] 3.811547e-03 大腸ガン <—> 大腸ガン
V[直腸ガン] 3.891416e-05 直腸ガン <—> 直腸ガン

Iterations = 23

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

library(Rgraphviz)
pathDiagram(fit,edge.labels="values",ignore.double=T,standardize=T,digits=2,file="sem2",graphics.fmt ="png",
node.font = c("serif", 14), edge.font = c("sans", 10),
min.rank="総熱量, 肉類, 乳製品, 酒類",max.rank="大腸ガン, 直腸ガン",same.rank=NULL)

従来の記法

modelB <- specifyModel()

model.b1

総熱量 -> 消化管ガン, gam1, NA
肉類 -> 消化管ガン, gam2, NA
乳製品 -> 消化管ガン, gam3, NA
酒類 -> 消化管ガン, gam4, NA
総熱量 <-> 肉類, phi12, NA
総熱量 <-> 乳製品, phi13, NA
総熱量 <-> 酒類, phi14, NA
肉類 <-> 乳製品, phi23, NA
肉類 <-> 酒類, phi24, NA
乳製品 <-> 酒類, phi34, NA

model.b2

消化管ガン -> 大腸ガン, NA, 1
消化管ガン -> 直腸ガン, lam61, NA
消化管ガン <-> 消化管ガン, psi1, NA

表記不要

総熱量 <-> 総熱量, phi11, NA
肉類 <-> 肉類, phi22, NA
乳製品 <-> 乳製品, phi33, NA
酒類 <-> 酒類, phi44, NA
大腸ガン <-> 大腸ガン, theta1, NA
直腸ガン <-> 直腸ガン, theta2, NA