共分散構造分析[R編]ってないの？[第6章1]

更新の気力を失っていましたが、ついに投稿する気になりました。

「共分散構造分析［Amos編］－構造方程式モデリング－」の102頁の成長曲線モデル(GCM)を計算してみます。
まずは、モデルを記述します。Interceptの制約が厄介です。観測変数の分散は省略してspecifyModel()に補完してもらいます。

mdl1 <- specifyModel()
傾き -> 体重1年, NA, 0
傾き -> 体重2年, NA, 1
傾き -> 体重3年, NA, 2
傾き -> 体重4年, NA, 3
切片 -> 体重1年, NA, 1
切片 -> 体重2年, NA, 1
切片 -> 体重3年, NA, 1
切片 -> 体重4年, NA, 1
傾き <-> 切片, cor1, NA
傾き <-> 傾き, vs, NA
切片 <-> 切片, vi, NA
Intercept -> 傾き, sl, NA
Intercept -> 切片, ic, NA
Intercept -> 体重1年, NA, 0
Intercept -> 体重2年, NA, 0
Intercept -> 体重3年, NA, 0
Intercept -> 体重4年, NA, 0

モデルのパス図をpathDiagrm()でgraphvizに出力したのが右の図です。

このモデルを下のように実行します。オプションで par.size = "ones" を指定しました。というのも、デフォルト(このデータでは多分 par.size = "startvalues")では収束しませんでしたので。収束しないときは、このオプションを変更してみるべしと、説明書きにも書いてあったように思います。



 

weight.sem1 <- sem(mdl1, data = c06Weight, fixed.x = "Intercept", raw = TRUE,
                   formula = ~ 体重1年 + 体重2年 + 体重3年 + 体重4年,
                   par.size = "ones")
summary(weight.sem1)

CFIとRMSEAは表示されません。
計算結果は一致します。
因子の利用方法が書かれているので、因子のデータを計算してもらいましょう。

fs <- fscores(weight.sem1, data=cbind(Intercept=1,c06Weight[,-5]), center=F)

 とすればいきなり因子得点がでてきそうですが、Rの場合は因子得点ウェイトやデータにInterceptが入っているので、求める値が計算されません。面倒ですが、一旦、因子得点を出力してInterceptを除いて因子得点を計算してやりました。

head(fw <- fscores(weight.sem1, data = NULL)) #因子得点ウエイト
ind <- as.matrix(c06Weight)[ , -5] %*% fw[ -1, ] #個人の因子得点

それではggplotなんぞを使って図6.7を描いてみます。

library(ggplot2)
library(reshape2)
ind.path <- ind[c(58,405,496), ] %*% matrix(c(rep( 1, 4), 0:3), 2, , byrow = TRUE) #各年゜の体重計算
dimnames(ind.path) <- list(c(58,405,496), colnames(c06Weight)[-5]) #行名、列名
fig <- ggplot(data = melt(ind.path)) +
  geom_path(aes(x = Var2, y = value, group = Var1 )) +
  geom_point(aes(x = Var2, y= value)) +
  labs(x = "体重", y = "測定時点")
print(fig)

ところで、今年の夏は暑かった。バテました。