myscratchpad

# 心理統計学の基礎　共分散構造分析　p351-363 の例を実行
# ...してみようと思ったら以下のURLにそのものズバリが載っていた。感謝。自分用にコピペしておく
 # http://home.hiroshima-u.ac.jp/ksatoh/documents/Rsem.txt
 # http://home.hiroshima-u.ac.jp/ksatoh/download.htm

library(sem)

# 表10-8　観測変数間の相関係数 p354
mat <- matrix(c(
1.000,0.160,0.302,0.461,0.299,0.152,0.134,0.182,0.251,0.372,0.157,0.203,
0.160,1.000,0.341,0.400,0.404,0.320,0.403,0.374,0.285,0.100,0.291,-0.014,
0.302,0.341,1.000,0.372,0.552,0.476,0.467,0.572,0.316,0.408,0.393,0.369,
0.461,0.400,0.372,1.000,0.302,0.225,0.256,0.255,0.164,0.236,0.229,0.224,
0.299,0.404,0.552,0.302,1.000,0.708,0.623,0.776,0.361,0.294,0.472,0.342,
0.152,0.320,0.476,0.225,0.708,1.000,0.324,0.769,0.295,0.206,0.351,0.202,
0.134,0.403,0.467,0.256,0.623,0.324,1.000,0.724,0.260,0.071,0.204,0.152,
0.182,0.374,0.572,0.255,0.776,0.769,0.724,1.000,0.284,0.142,0.320,0.189,
0.251,0.285,0.316,0.164,0.361,0.295,0.260,0.284,1.000,0.295,0.290,0.418,
0.372,0.100,0.408,0.236,0.294,0.206,0.071,0.142,0.295,1.000,0.468,0.351,
0.157,0.291,0.393,0.229,0.472,0.351,0.204,0.320,0.290,0.468,1.000,0.385,
0.203,-0.014,0.369,0.224,0.342,0.202,0.152,0.189,0.418,0.351,0.385,1.000),
nr=12,
dimnames = list(paste("y", 1:12, sep=""), paste("y", 1:12, sep=""))
)
mat

# 日本語じゃなくて英語にしておく
model.agg <- specify.model()
mot -> y1, b11, NA
mot -> y2, b12, NA
mot -> y3, b13, NA
mot -> y4, b14, NA
int -> y5,  NA, 1
int -> y6, b22, NA
int -> y7, b23, NA
int -> y8, b24, NA
agg ->  y9,  NA, 1
agg -> y10, b32, NA
agg -> y11, b33, NA
agg -> y12, b34, NA
mot -> int, g11, NA
mot ->       agg, g12, NA
 y1 <->  y1,  e1, NA
 y2 <->  y2,  e2, NA
 y3 <->  y3,  e3, NA
 y4 <->  y4,  e4, NA
 y5 <->  y5,  e5, NA
 y6 <->  y6,  e6, NA
 y7 <->  y7,  e7, NA
 y8 <->  y8,  e8, NA
 y9 <->  y9,  e9, NA
y10 <-> y10, e10, NA
y11 <-> y11, e11, NA
y12 <-> y12, e12, NA
mot <-> mot, NA, 1
int <-> int, delta2, NA
agg <-> agg, delta3, NA

sem.agg <- sem(model.agg, mat, N=50)
summary(sem.agg)
std.coef(sem.agg,digit=4) # 標準解 図10-6と一致するのはこっち

# 図10-6　くらいは自分で書いてみよう
library(Rgraphviz)
path.diagram(sem.agg, ignore.double=FALSE, edge.labels="values", digits=3, file="out2")
# やはり、日本語フォントでは表示できない。Graphvizのバージョンが古いのかもしれない
# しかたないので協調性はagg, 母親価値はmot, 相互作用経験はintとした。でも文字の一部が消えたり、どうもうまくいかない。
# そのうえ、sem.aggの係数で描画するため、標準解を出力できない