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

「共分散構造分析［Amos編］－構造方程式モデリング－」の63頁図3.18をRで計算します。
セミナー評価という潜在因子は誤差を考慮せずに観測変数に重みをかけて足し上げた合成指標ということです。図3.17のセミナー評価は重回帰分析のように誤差を計算させましたが、図3.18では誤差(因子独自の分散)を0に制約するということになるはずです。
モデルは次のように書き、 semi.3.18.txtとして保存しました。

テキスト -> セミナー評価, p1, NA
プレゼン -> セミナー評価, p2, NA
ペース   -> セミナー評価, p3, NA
講師対処 -> セミナー評価, p4, NA
セミナー評価 -> 満足度,   p5, NA
セミナー評価 -> 理解度,   p6, NA
セミナー評価 -> 目的一致, NA, 1　　　　　 #ここは1に制約
テキスト <-> テキスト, v1, NA
プレゼン <-> プレゼン, v2, NA
ペース <-> ペース,     v3, NA
講師対処 <-> 講師対処, v4, NA
満足度 <-> 満足度,     e1, NA
理解度 <-> 理解度,     e2, NA
目的一致 <-> 目的一致, e3, NA
セミナー評価 <-> セミナー評価, NA, 0 　　    #ここは0に制約
テキスト <-> プレゼン, c1, NA
テキスト <-> ペース,   c2, NA
テキスト <-> 講師対処, c3, NA
プレゼン <-> ペース,   c4, NA
プレゼン <-> 講師対処, c5, NA
ペース <-> 講師対処,   c6, NA

ファイルを読み込んで、実行します。

semi.mdl <- specifyModel("semi.3.18.txt")
semi.sem <- sem(semi.mdl, data=semi.data)
summary(semi.sem)
stdCoef(semi.sem)
round(1-diag(semi.sem$P[-8,-8] / semi.sem$S),2)
round(cov2cor(semi.sem$P),2)

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

3.6 MIMICモデル

「共分散構造分析［Amos編］－構造方程式モデリング－」の62ページです。
ＭＩＭＩＣモデルの記述は長いので、次の別ファイル(semi.3.17.txt)を用意しました。セミナー評価から出ている矢印のうち、1つを制約しないと識別できません。 セミナー評価から目的一致のパスを1に制約しました。

テキスト -> セミナー評価, p1, NA
プレゼン -> セミナー評価, p2, NA
ペース   -> セミナー評価, p3, NA
講師対処 -> セミナー評価, p4, NA
セミナー評価 -> 満足度,   p5, NA
セミナー評価 -> 理解度,   p6, NA
セミナー評価 -> 目的一致, NA, 1　 #このパスを1に制約
テキスト <-> テキスト, v1, NA
プレゼン <-> プレゼン, v2, NA
ペース <-> ペース,     v3, NA
講師対処 <-> 講師対処, v4, NA
満足度 <-> 満足度,     e1, NA
理解度 <-> 理解度,     e2, NA
目的一致 <-> 目的一致, e3, NA
セミナー評価 <-> セミナー評価, e4, NA
テキスト <-> プレゼン, c1, NA
テキスト <-> ペース,   c2, NA
テキスト <-> 講師対処, c3, NA
プレゼン <-> ペース,   c4, NA
プレゼン <-> 講師対処, c5, NA
ペース <-> 講師対処,   c6, NA

これを次のように実行します。

semi.mdl <- specifyModel("semi.3.17.txt")
semi.sem <- sem(semi.mdl, data=semi.data)
summary(semi.sem)
stdCoef(semi.sem)
diag(semi.sem$P[-8,-8] / semi.sem$S)
cov2cor(semi.sem$P)

最後の三行でそれぞれ

• パスの標準化推定値
•  満足度、理解度、目的一致の分散
• テキスト、プレゼン、ペース、講師対処の相関

を確認します。
本に近い値が出ています。

次は、図3.18のモデルに取り組みます。

共分散構造分析（Amos編） 著者：豊田秀樹 価格：3,360円（税込、送料込） 楽天ブックスで詳細を見る

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

3.5 共分散分析

共変量分析というと、以前、日常的に運動をしている人は最大酸素摂取量(VO2MAX)が多くなるか確認できないかと尋ねられて検討したと思います。この時の補助変数(共変量って言ってたような)は年齢でした。群を分けるのも「日常的に運動していますか」という設問への自己申告で心もとない感じだったのですが。
同じようにボディマス指数(BMI)も確認したいということもあったのですが、

• 20歳代は運動すると、筋肉がたくさんつくみたいで運動する人の方がBMIが大きい
• 30歳代半ば過ぎると、運動する人のほうが贅肉がつきにくいということなのかBMIが小さい

という関係があるようでした。つまり、(1)の回帰直線の傾きが群間で等しいという仮定が成り立たなかったのです。(2)の補助変数の影響に関しても、年齢と運動強度の間に関係がある？影響がでているのかもしれません。
さらに、年齢とBMIの関係が対象とした年齢(20-60)では、50歳程度でピークとなる非線形の関係があるようでした。人間って50歳過ぎると萎むんですね。私なんか、なんだか縮んじゃって勤め先でも居場所がなくなってきてるんですよ・・・。(; _ ;)

そんなボヤキはさておき、こうしてみました。

english.mdl <- specifyModel()
学力テスト -> 英語の試験, p1, NA
群 -> 英語の試験, p2, NA
英語の試験 <-> 英語の試験, e1, NA
学力テスト <-> 学力テスト, v1, NA
群 <-> 群, v2, NA
1 -> 英語の試験, i0, NA

english.raw <- rawMoments(cbind(english.data,1))
english.sem <- sem(english.mdl, S=english.raw, N=24, raw=TRUE, fixed.x=c("1"))
summary(english.sem)
modIndices(english.sem)

切片を計算するために、

• 元のデータに1ばかりの列を追加して、rawMoments()でsem()用に積和行列を計算
• sem()関数のSに計算した積和行列を指定して、raw=TRUEで積和行列を使用していることを知らせ
• fixed.xオプションに変数"1"を指定

summary()の出力では、p2の標準誤差が本より少し小さく、信頼区間もすこし狭くなります。
修正指数を出力するmodIndices()関数には、しきい値の設定がないようです。出力された指数が無視してよいものなのかどうなのか判断に苦しみます。ただ、本で指定しているしきい値(4)よりはずっと小さいです。

次に、群 -> 英語の試験のパスを0に制約してみます。モデルは次のように書きました。

english.mdl <- specifyModel()
学力テスト -> 英語の試験, p1, NA
群 -> 英語の試験, NA, 0
英語の試験 <-> 英語の試験, e1, NA
学力テスト <-> 学力テスト, v1, NA
群 <-> 群, v2, NA
1 -> 学力テスト, i1, NA

改善指数が一番大きいのは、群 -> 英語の試験にパスを引くことになります。修正指数はAmosのそれに近い値が出ています。

共分散構造分析[R編]ってないの？[第３章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];
 }

完成です。

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

3.3 非逐次モデル

さて、非逐次モデルに進みます。野球のAチームに対する応援態度についての架空データです。

fan.mdl <- specifyModel()
夫の父 -> 夫, p1, NA
妻の父 -> 妻, p2, NA
夫 -> 妻, p3, NA
妻 -> 夫, p4, NA
夫の父 <-> 妻の父, c1, NA
夫の父 <-> 夫の父, v1, NA
妻の父 <-> 妻の父, v2, NA
夫 <-> 夫, e1, NA
妻 <-> 妻, e2, NA
#ここはモデルの最後を示すための空改行
fan.sem <- sem(fan.mdl, data=fan.data)
summary(fan.sem)

• 適合度指標も四捨五入した値でGFI,AGFIが、第三位でAmosの出力と0.01小さくなりました。RMSEAは0.08も大きくなっています。
• 夫の父 -> 父のパス係数は四捨五入した値で、0.01小さくなっています。
• 夫と妻に刺さっている誤差分散の係数がそれぞれ0.03大きくなっています。

標準化推定値を確認してみます。

stdCoef(fan.sem)

• パス係数は一致していますが、
• 夫と妻に刺さっている誤差分散がそれぞれ0.10、0.15大きくなっています。

これはあまり面白くない違いに思われます。モデルを同じように解釈できるでしょうか。
モデルの書き方が正しくないのでしょうか。

これに関連するとおもわれるトピックがRの掲示板にありました。
SEM: question regarding how standard errors are calculated
John Fox先生もびっくりということですが、Amosが何をしているのかよくわからないので、はっきりした理由はわかりません。結局、Dochtermannさんはいろいろやってみたんだろうか？

---

semパッケージのダブルヘッドの出力は分散、共分散ですが、Amosは標準偏差と相関がでているみたいです。本にも、「誤差分散」と繰り返されているので・・・。(^ ^;)
ということで、v_やe_の平方根を計算してみると、Amosの出力と同じになりました。でも、私はパス図では見せるなら標準偏差と相関の出力のほうが解釈しやすいと思います。
しまりのないオチでした。m(_ _)m

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

3.1 重回帰モデル

共分散構造分析[R編]の携帯電話の購入意欲データの重回帰分析モデルのスクリプトを作成してみましょう。
式が少ないので別のファイルを作成しないで、specifyModel()の次の行から続けて式を入力します。何も入力しないで改行すると最終行と認識してくれます。
cellphone.corは相関行列なのでいきなり標準化推定値を計算します。
共分散分散行列や相関行列をデータとして使うときは、Nオプションに標本数を設定します。
Amosのパス図の記法とは違って、外生変数の「性能」にも分散を記述しないと、計算できないです。

## 3.1 重回帰分析 #########################
cellphone.mdl &lt;- specifyModel()
性能 -> 値段, a, NA
値段 -> 購買意欲, b, NA
性能 -> 購買意欲, c, NA
性能 <-> 性能, v1, NA
値段 <-> 値段, e1, NA
購買意欲 <-> 購買意欲, e2, NA

cellphone.sem <- sem(cellphone.mdl, S=cellphone.cor, N=250)

計算結果を確認してみましょう。

summary(cellphone.sem)

購買意欲の決定係数はe2を見ればわかると思うんですが、「1－説明できなかった分散／観測された分散」の計算はこれでできるようです。

1 - diag(cellphone.sem$P)/diag(cellphone.sem$S)

3.2 逐次モデル

次に逐次モデルをに進みます。結婚相手評価データは観測データです。分散共分散行列ではなくそのままsem()関数に使用することにして、dataオプションでデータフレーム名を設定してやりました。

marriage.mdl <- specifyModel()
学歴 -> 年収, p1, NA
職業威信 -> 年収, p2, NA
職業威信 -> 評価, p3, NA
年収 -> 評価, p4, NA
学歴 <-> 職業威信, c1, NA
年収 <-> 年収, e1, NA
評価 <-> 評価, e2, NA
学歴 <-> 学歴, v1, NA
職業威信 <-> 職業威信, v2, NA

marriage.sem <- sem(marriage.mdl, data=marriage.data)
summary(marriage.sem)

「かっこ良さ」とか「人柄」なんていう変数は用意されていません。これでデータとの適合度が非常によいなんて救いのない分析ですが、架空のデータということです。

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

パス図は、graphvizでもってできます。次のスクリプトでは標準化係数のパス図を出力します。日本語を使うときは日本語フォントを必ず指定しなくてはなりません。

pathDiagram(semi.sem, 
            edge.font=c("MSGothic",10), node.font=c("MSMincho",10),
            edge.labels="values", standardize=TRUE, digits=3,
            ignore.double=FALSE,
            rank.dir="TB", max.rank = c("講師の質","充実感"),
            min.rank=c("目的一致","理解度","テキスト","ペース","満足度","プレゼン","講師対処"))

出力は次のごとし。

digraph "semi.sem" {
  rankdir=TB;
  size="8,8";
  node [fontname="IPAMincho" fontsize=10 shape=box];
  edge [fontname="IPAGothic" fontsize=10];
  center=1;
  {rank=min "目的一致"}
(中略)
  "理解度" -> "理解度" [label="0.839" dir=both];
  "目的一致" -> "目的一致" [label="0.848" dir=both];
}
 警告メッセージ： 
'std.coef' is deprecated.
Use 'stdCoef' instead.
See help("Deprecated") and help("sem-deprecated"). 

現時点で、標準化係数のパス図を作成するのに古い関数を使っているようで、警告メッセージがでますが当面問題ないでしょう。
なお、graphvizで使うファイルはテキストファイルをutf-8で保存しなくてはならないのです。僕のようにwindowsな人は、 何もしない場合shift-jisで保存されますから、そのままではパス図は大変なことになります。ファイルをutf-8に変換してやりましょう。
dotモードが設定できるEMACSを使うのが便利かもしれませんが、
導入の敷居が低いRStadioの場合でも、こんなに楽ちんでした。
ただし、予めパスを通すって作業は必要。

1. 左上のアイコンから新しいテキスト[New Text]を作って、fileオプションなしのpathDiagram()でコンソールに出力されたコマンド？をコピペ。
2. [file]-[save with encoding...]でutf-8を指定して保存します。(例えば、作業ディレクトリにsemipath.dotで)
3. [tool]-[shell]でDOS窓を開いて、パス図を描いてくださいとおまじないを打ち込むと、音もなくpdfファイルが作成されます。 

※2

→

dot -Tpdf semipath.dot > semipath.pdf

• 本と同じ向きの出力にするには、2行目のrankdir=TBをrankdir=BTに書き換えます。現時点では、rank.dirオプションを"BT"とするとエラーになりました。

適合度指標などは、summary()で出力されてました。
次に、直接効果、間接効果、総合効果も計算してもらいました。

effects(semi.sem)

開発環境もRStadioのような立派なものができて、僕のようないんちきユーザーでもかなり便利になりました。
ここまでは、本に書かれているものとほとんど同じに見えます。

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

「共分散構造分析[R編]--構造方程式モデリング」(豊田秀樹編著)が、2014年4月に出版されました。RでSEMを実行したい方は、そちらを購読されることをおすすめします。

共分散構造分析［Amos編］を購入しましたが、［Ｒ編］というのは出ないのでしょうね。個人でAmosなんてなかなか買えないです。
さて、Rのsemパッケージが3.0になり、多母集団分析 もできるようになったようです。これでもって、Amos編に掲載されている事例をそのまま(^^;)なぞってみようではないかと思ったものであります。備忘録も兼ねて・・・。

でも、公開されている本書のデータ[02008DL.exe]を使って分析した成果をwebなんぞに公開するには、著作権者(豊田先生？)に許可を貰わなくてはなりませんと明記されています。こんなつまんないことで私なんぞが豊田先生に許可を求めるなんてあまりに畏れ多いので、結果は表示しません。

まず、第一章の8ページのパス図。
データはsemi.tabにタブ区切りで保存してあります。
このモデルを次のように記述して保存しました。
semi.txt

講師の質 <-> 充実感, sf, NA
講師の質 -> テキスト, p1, NA
講師の質 -> プレゼン, p2, NA
講師の質 -> ペース, p3, NA
講師の質 -> 講師対処, p0, NA
充実感  -> 満足度, NA, 1
充実感  -> 理解度, p5, NA
充実感  -> 目的一致, p6, NA
目的一致 <-> 目的一致, e1, NA
理解度 <-> 理解度  , e2, NA
テキスト <-> テキスト, e3, NA
ペース  <-> ペース  , e4, NA
満足度  <-> 満足度  , e5, NA
プレゼン <-> プレゼン, e6, NA
講師対処 <-> 講師対処, e7, NA
講師の質 <-> 講師の質, NA, 1
充実感   <-> 充実感, sat, NA 

• なんと、例えばeをつけた誤差分散は記述しなくてもspecifyModel()関数が補ってくれます。その時は、こんな風に補完されています。

> semi.mdl <- specifyModel("model/semi.txt")

Read 10 records
NOTE: adding 7 variances to the model
> semi.mdl
   Path                  Parameter   StartValue
1  講師の質 -> 充実感    sf
　　：
11 テキスト <-> テキスト V[テキスト]
12 プレゼン <-> プレゼン V[プレゼン]
13 ペース <-> ペース     V[ペース]
14 講師対処 <-> 講師対処 V[講師対処]
15 満足度 <-> 満足度     V[満足度]
16 理解度 <-> 理解度     V[理解度]
17 目的一致 <-> 目的一致 V[目的一致]

• Amosだと、講師の質を測定するパスのうち一つには1に制約するけれど、semではその制約を入れると期待した結果は求められませんでした。

Rで次のように実行します。

semi <- read.delim("data/semi.tab")
semimdl <- specifyModel("model/semi.txt") #モデルのオブジェクトを作成
semi.sem <- sem(semi.mdl, data = semi[,-1], N = nrow(semi)) #分析を実行
summary(semi.sem)
stdCoef(semi.sem) #標準化係数

標準化パス係数は小数点2桁まで同じ

• GFI, AGFI, CFI, RMSEAは小数点3桁まで同じ

 つぎに、構成概念スコアを出力してみます。fscores()で直接求めることができるようです。 dataオプションにNULLを指定すると因子得点が出力されるのですね。

fscores(semi.sem, data=NULL)                    　#因子得点ウェイト
semi.fs <- fscores(semi.sem, scale=TRUE)   #構成概念スコア
fs.plot <- ggplot(data.frame(semi.fs)) +
  geom_point(aes(x=講師の質, y=充実感))       #ggplot2で散布図を描いた
fs.plot                              

• 本の出力とはかなり違います。ただし、それぞれの因子の中では本の出力と比例しているようです。

次は、第２章の事例をやってみよう。

共分散構造分析（Amos編） 著者：豊田秀樹 価格：3,360円（税込、送料込） 楽天ブックスで詳細を見る