sem_OpenMxパッケージ

OpenMxパッケージ

(参考)OpenMx HP

(注)AIC/BIC is NA, but other FIs are computed

1
2
3
4
5
6
library(OpenMx)
library(semPlot)
#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("総熱量","肉類","乳製品","酒類","大腸ガン","直腸ガン")
多重指標モデル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
manifests <- colnames(data)
diet <- colnames(data)[1:4]
cancer <- colnames(data)[5:6]
latents <- c("洋食傾向","消化管ガン")
model <- mxModel("多重指標モデル",
type="RAM",
manifestVars=manifests,
latentVars=latents,
mxPath(from=manifests, arrows=2,free=T,
labels=c("e1","e2","e3","e4","e5","e6"),values=1),
mxPath(from=latents,arrows=2,free=T,
labels=c("e7","e8"),values=1),
mxPath( from= "洋食傾向",to="消化管ガン",arrows=1,
free=T,labels="cor",values=1),
mxPath(from="洋食傾向",to=diet, arrows=1,
free=T, labels=c("l1","l2","l3","l4"),values=1),
mxPath(from="消化管ガン",to=cancer,arrows=1,
free=T,labels=c("l5","l6"),values=1),
mxData(observed=cor(data),type="cor",numObs=nrow(data))
)
MIMICモデル
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
manifests <- colnames(data)
diet <- colnames(data)[1:4]
cancer <- colnames(data)[5:6]
latents <- "消化管ガン"
model <- mxModel("MIMIC",
type="RAM",
manifestVars=manifests,
latentVars=latents,
mxPath(from=cancer, arrows=2,free=T,
labels=c("e1","e2"),values=1),
mxPath( from=c("総熱量","肉類"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF1","cor1","varG1")),
mxPath( from=c("総熱量","乳製品"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF2","cor2","varG2")),
mxPath( from=c("総熱量","酒類"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF3","cor3","varG3")),
mxPath( from=c("肉類","乳製品"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF4","cor4","varG4")),
mxPath( from=c("肉類","酒類"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF5","cor5","varG5")),
mxPath( from=c("乳製品","酒類"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF6","cor6","varG6")),
mxPath(from=latents, arrows=2,free=T,labels="e7",values=1),
mxPath(from=diet, to="消化管ガン", arrows=1,
free=T, labels=c("l1","l2","l3","l4"),values=1),
mxPath(from="消化管ガン",to=cancer,arrows=1,
free=T,labels=c("l5","l6"),values=1),
mxData(observed=cor(data),type="cor",numObs=nrow(data))
)
PLSモデル
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
manifests <- colnames(data)
diet<- colnames(data)[1:4]
cancer <- colnames(data)[5:6]
latents <-c("洋食傾向", "消化管ガン")
model <- mxModel("PLS",
type="RAM",
manifestVars=manifests,
latentVars=latents,
mxPath(from=cancer, arrows=2,free=T,
labels=c("e1","e2"),values=1),
mxPath(from=diet, to="洋食傾向", arrows=1,
free=T, labels=c("l1","l2","l3","l4"),values=1),
mxPath( from=c("総熱量","肉類"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF1","cor1","varG1")),
mxPath( from=c("総熱量","乳製品"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF2","cor2","varG2")),
mxPath( from=c("総熱量","酒類"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF3","cor3","varG3")),
mxPath( from=c("肉類","乳製品"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF4","cor4","varG4")),
mxPath( from=c("肉類","酒類"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF5","cor5","varG5")),
mxPath( from=c("乳製品","酒類"),arrows=2, connect="unique.pairs",
free=c(F,T,F), values=c(1,0,1), labels=c("varF6","cor6","varG6")),
mxPath(from="消化管ガン",to=cancer,arrows=1,
free=T,labels=c("l5","l6"),values=1),
mxPath( from= "洋食傾向",to="消化管ガン",arrows=1,
free=F,labels="cor",values=1),
mxPath(from= latents, arrows=2,free=F,labels=c("e7","e8"),values=c(0,1)),
mxData(observed=cor(data),type="cor",numObs=nrow(data))
)

分析結果とグラフを出力

1
2
3
4
#omxGraphviz(model, "sem.dot")
Fit <- mxRun(model)
summary(Fit)
semPaths(Fit,"std",rotation =2,edge.label.cex =1,exoVar=F,fade=F)