2012年9月15日土曜日

共分散構造分析[R編]ってないの?[第3章4]

3.4 偏相関と多変量回帰

生活環境と体力の架空データから、偏相関係数を求めるための共分散構造分析です。
strength.mdl <- specifyModel()
年齢 -> 結婚年数, p1, NA
年齢 -> 50m走, p2, NA
年齢 <-> 年齢, v1, NA
結婚年数 <-> 50m走, pc, NA
結婚年数 <-> 結婚年数, v2, NA
50m走 <-> 50m走, v3, NA

strength.sem <- sem(strength.mdl, S=strength.cor, N=50)
summary(strength.sem)
stdCoef(strength.sem)

これでパス係数は正しく計算されました。しかし、本のAmosが誤差の相関を出力している一方で、semパッケージは誤差共分散を出力しています。結婚年数50m走の共分散を結婚年数の分散と50m走の分散の積で除しての平方根を計算すれば良いはずです。
でも、面倒だし、たくさん誤差相関を計算したいときもあるので、
cov2cor(strength.sem$P)

で計算して、目的のセルを参照したら良いのではないでしょうか。
『手計算では面倒だった計算もAmosでは簡単』ということですが、semパッケージ×覚束ないタッチタイピングではそれなりに面倒。

-----
次に、多変量回帰分析を実行します。

study.mdl <- specifyModel()
知能 -> 試験, p1, NA
知能 -> 睡眠時間, p2, NA
勉強時間 -> 試験, p3, NA
勉強時間 -> 睡眠時間, p4, NA
知能 <-> 勉強時間, c1, NA
試験 <-> 睡眠時間, c2, NA
知能 <-> 知能, v1, NA
勉強時間 <-> 勉強時間, v2, NA
試験 <-> 試験, v3, NA
睡眠時間 <-> 睡眠時間, v4, NA

study.sem <- sem(study.mdl, S=study.cor, N=study.n)
summary(study.sem)

これもパス図の値は四捨五入すれば同じです。Amosと同じように相関係数を計算してみます。
cov2cor(study.sem$P)

これもAmos同様の出力が得られました。
パス図も描いてみましょう。edge.labelsオプションを"value"にすると推計値が出力されます。
pathDiagram(study.sem,
            min.rank=c("試験","睡眠時間"), max.rank=c("知能","勉強時間"),
            ignore.double=FALSE, edge.labels="names",
            edge.font=c("IPAMincho",8), node.font=c("IPAGothic",10), rank.dir="LR")

共分散のパスはheadportとtailportのオプションを書き加えました。向きもLRに書き換えて、本のパス図の向きと同じにします。
digraph "study.sem" {
  rankdir=RL;
  size="8,8";
  node [fontname="IPAGothic" fontsize=10 shape=box];
  edge [fontname="IPAMincho" fontsize=8];
  center=1;
  {rank=min "試験"}
   {rank=min "睡眠時間"}
  {rank=max "知能"}
   {rank=max "勉強時間"}
  "知能" -> "試験" [label="p1"];
  "知能" -> "睡眠時間" [label="p2"];
  "勉強時間" -> "試験" [label="p3"];
  "勉強時間" -> "睡眠時間" [label="p4"];
  "知能" -> "勉強時間" [label="c1" dir=both headport="e" tailport="e"];
  "試験" -> "睡眠時間" [label="c2" dir=both headport="w" tailport="w"];
  "知能" -> "知能" [label="v1" dir=both];
  "勉強時間" -> "勉強時間" [label="v2" dir=both];
  "試験" -> "試験" [label="v3" dir=both];
  "睡眠時間" -> "睡眠時間" [label="v4" dir=both];
 }

完成です。

0 件のコメント:

コメントを投稿