ついに、待望の「共分散構造分析[R編]―構造方程式モデリング」が出版されました。AmosやOpenMxといった商用ソフトが容易に使える環境にない私には良い報せです。幾分天邪鬼な私は、`楽天ブックス`に注文していましたが、この金曜日に配送されてきました。手に取ると[Amos編]より幾分重い感じがしました。ページ数は1割増で、内容もぎっしりです。やはり'lavaan'を使用した内容になっていました。当分の間、本物の「共分散構造分析[R編]」を読みつつ勉強したいと思います。
[Amos編]の再現も、7章まではやってみようと思っていましたが、ここでこのプロジェクトも一区切りという感じです。答えがわかっていてやっているのだから気楽なものですが、それでも収束させることができなかったモデルも残ってしまいました。残念。気が向いたらまた試してみることにしよう。
2014年4月13日日曜日
2013年9月7日土曜日
共分散構造分析[R編]ってないの？[第6章4]
「共分散構造分析［Amos編］－構造方程式モデリング－」の106頁の成長曲線モデルと他の変数の組み合わせを計算させてみましょう。
運動志向から傾きや切片に矢印を入れてやります。運動志向という変数は外生変数でありかつ観測変数なので、fixed.xオプションに追加します。
パス係数は出力は一致します。分散はおかしな値になっています。平均構造を導入したモデルでは、運動志向の双方向矢印には何が出力されているのか、よくわかりません。(変数[運動志向]の二乗和の平均sum(c06Weight$運動志向^2)/500) = 100.8883が出力されているようです。)
ところで、semパッケージで平均構造のあるモデルを計算するとAICやBIC以外の適合度指数は出力されません。本では、AMOSだとCFIやRMSEAが出力されるように読めます。
いずれにしても、適合度を見ずに済ますと言うわけにもいかないです。semパッケージで出力されないなら、自分で何とかしなくてはね、そのうち．．．。
運動志向から傾きや切片に矢印を入れてやります。運動志向という変数は外生変数でありかつ観測変数なので、fixed.xオプションに追加します。
mdl3 <- 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
傾き <-> 切片, cov1, 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
運動志向 -> 傾き, exsl, NA
運動志向 -> 切片, exic, NA
運動志向 <-> 運動志向, vex, NA
#
weight.sem3 <- sem(mdl3, data = c06Weight, fixed.x = "Intercept", raw = TRUE,
formula = ~ 体重1年 + 体重2年 + 体重3年 + 体重4年 + 運動志向,
par.size = "ones")
summary(weight.sem3)
傾き -> 体重1年, NA, 0
傾き -> 体重2年, NA, 1
傾き -> 体重3年, NA, 2
傾き -> 体重4年, NA, 3
切片 -> 体重1年, NA, 1
切片 -> 体重2年, NA, 1
切片 -> 体重3年, NA, 1
切片 -> 体重4年, NA, 1
傾き <-> 切片, cov1, 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
運動志向 -> 傾き, exsl, NA
運動志向 -> 切片, exic, NA
運動志向 <-> 運動志向, vex, NA
#
weight.sem3 <- sem(mdl3, data = c06Weight, fixed.x = "Intercept", raw = TRUE,
formula = ~ 体重1年 + 体重2年 + 体重3年 + 体重4年 + 運動志向,
par.size = "ones")
summary(weight.sem3)
パス係数は出力は一致します。分散はおかしな値になっています。平均構造を導入したモデルでは、運動志向の双方向矢印には何が出力されているのか、よくわかりません。(変数[運動志向]の二乗和の平均sum(c06Weight$運動志向^2)/500) = 100.8883が出力されているようです。)
ところで、semパッケージで平均構造のあるモデルを計算するとAICやBIC以外の適合度指数は出力されません。本では、AMOSだとCFIやRMSEAが出力されるように読めます。
いずれにしても、適合度を見ずに済ますと言うわけにもいかないです。semパッケージで出力されないなら、自分で何とかしなくてはね、そのうち．．．。
2013年9月1日日曜日
共分散構造分析[R編]ってないの？[第6章3.2]
引き続いて、「共分散構造分析［Amos編］－構造方程式モデリング－」の105頁の曲線的な成長の表現。
ここは、 6.1からモデルの記述の変更はわずかです。なるほど、こうすれば測定時点ごとに傾きが自由になるんだなぁ。
par.size = "ones"を設定して実行しました。
計算結果は、一致しているようです。一応、図6.8や個人の図も描いてみます。
あまりはっきりしませんが、線形ではないことは確認できます。
元のデータは架空データなので現実の大学生のデータとは全然違うのでしょうが、私個人は、4年間で7kgくらい痩せました。当時、ちゃんとご飯を食べてなかったのでしょう。元々痩せっぽちだったので、就活の時は「病気しただろ」と疑われました。最近の大学生の中には親のすねが細って、経済的に苦しい人も少なくないとか聞きます。大学生の方が見ていたら老婆(爺?)心ながら申しますと、食費は削りたくなるものですが、ご飯はしっかり食べましょうね。
ここは、 6.1からモデルの記述の変更はわずかです。なるほど、こうすれば測定時点ごとに傾きが自由になるんだなぁ。
par.size = "ones"を設定して実行しました。
mdl2 <- specifymodel()
傾き -> 体重1年, NA, 0
傾き -> 体重2年, NA, 1
傾き -> 体重3年, w3, NA #
傾き -> 体重4年, w4, NA #
切片 -> 体重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
#
weight.sem2 <- sem(mdl2, data = c06Weight, fixed.x = "Intercept", raw = TRUE,
formula = ~ 体重1年 + 体重2年 + 体重3年 + 体重4年,
par.size = "ones")
summary(weight.sem2)
傾き -> 体重1年, NA, 0
傾き -> 体重2年, NA, 1
傾き -> 体重3年, w3, NA #
傾き -> 体重4年, w4, NA #
切片 -> 体重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
#
weight.sem2 <- sem(mdl2, data = c06Weight, fixed.x = "Intercept", raw = TRUE,
formula = ~ 体重1年 + 体重2年 + 体重3年 + 体重4年,
par.size = "ones")
summary(weight.sem2)
計算結果は、一致しているようです。一応、図6.8や個人の図も描いてみます。
 |
 |
あまりはっきりしませんが、線形ではないことは確認できます。
元のデータは架空データなので現実の大学生のデータとは全然違うのでしょうが、私個人は、4年間で7kgくらい痩せました。当時、ちゃんとご飯を食べてなかったのでしょう。元々痩せっぽちだったので、就活の時は「病気しただろ」と疑われました。最近の大学生の中には親のすねが細って、経済的に苦しい人も少なくないとか聞きます。大学生の方が見ていたら老婆(爺?)心ながら申しますと、食費は削りたくなるものですが、ご飯はしっかり食べましょうね。
2013年8月31日土曜日
共分散構造分析[R編]ってないの？[第6章1]
更新の気力を失っていましたが、ついに投稿する気になりました。
「共分散構造分析［Amos編］－構造方程式モデリング－」の102頁の成長曲線モデル(GCM)を計算してみます。
まずは、モデルを記述します。Interceptの制約が厄介です。観測変数の分散は省略してspecifyModel()に補完してもらいます。
モデルのパス図をpathDiagrm()でgraphvizに出力したのが右の図です。
このモデルを下のように実行します。オプションで par.size = "ones" を指定しました。というのも、デフォルト(このデータでは多分 par.size = "startvalues")では収束しませんでしたので。収束しないときは、このオプションを変更してみるべしと、説明書きにも書いてあったように思います。
CFIとRMSEAは表示されません。
計算結果は一致します。
因子の利用方法が書かれているので、因子のデータを計算してもらいましょう。
とすればいきなり因子得点がでてきそうですが、Rの場合は因子得点ウェイトやデータにInterceptが入っているので、求める値が計算されません。面倒ですが、一旦、因子得点を出力してInterceptを除いて因子得点を計算してやりました。
それではggplotなんぞを使って図6.7を描いてみます。
ところで、今年の夏は暑かった。バテました。
「共分散構造分析［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
傾き -> 体重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)
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)
head(fw <- fscores(weight.sem1, data = NULL)) #因子得点ウエイト
ind <- as.matrix(c06Weight)[ , -5] %*% fw[ -1, ] #個人の因子得点
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)
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)
ところで、今年の夏は暑かった。バテました。
2013年4月17日水曜日
共分散構造分析[R編]ってないの？[第5章4]
しばらく、更新を中断していました。第5章の2はうまく収束せず、またまた宿題とします。ここでは、「共分散構造分析［Amos編］－構造方程式モデリング－」の96頁の平均構造を含めた構成概念間のパス解析をやってみます。
このモデルでは、構成概念間に平均構造を含めるために、観測変数の平均は0に制約するということのようです。構成概念が観測変数の平均を決めるので、観測変数独自の平均は0ということです。そこで、以下のようなファイルを用意しました。
講師の質 -> 講師対処, NA, 1
講師の質 -> ペース , path1, NA
講師の質 -> プレゼン, path2, NA
講師の質 -> テキスト, path3, NA
充実感 -> 目的一致, NA, 1
充実感 -> 理解度, path5, NA
充実感 -> 満足度, path6, NA
講師の質 -> 充実感, path7, NA
講師の質 <-> 講師の質, vq, NA
充実感 <-> 充実感, vf, NA
Intercept -> 講師の質, iq, NA
Intercept -> 充実感, if, NA
Intercept -> 講師対処 , NA, 0
Intercept -> ペース,NA, 0
Intercept -> プレゼン, NA, 0
Intercept -> テキスト, NA, 0
Intercept -> 目的一致, NA, 0
Intercept -> 理解度, NA, 0
Intercept -> 満足度, NA, 0
講師の質 -> ペース , path1, NA
講師の質 -> プレゼン, path2, NA
講師の質 -> テキスト, path3, NA
充実感 -> 目的一致, NA, 1
充実感 -> 理解度, path5, NA
充実感 -> 満足度, path6, NA
講師の質 -> 充実感, path7, NA
講師の質 <-> 講師の質, vq, NA
充実感 <-> 充実感, vf, NA
Intercept -> 講師の質, iq, NA
Intercept -> 充実感, if, NA
Intercept -> 講師対処 , NA, 0
Intercept -> ペース,NA, 0
Intercept -> プレゼン, NA, 0
Intercept -> テキスト, NA, 0
Intercept -> 目的一致, NA, 0
Intercept -> 理解度, NA, 0
Intercept -> 満足度, NA, 0
次のスクリプトを実行すれば、収束すると思ったのですが。
summary(c5semi.sem <- sem(c5semi.mdl, data=semi, fixed.x="Intercept", raw=TRUE,
debug=TRUE))
debug=TRUE))
ダメですぅ。
「収束で困ったらpar.sizeオプションの値を変えてみれば」というふうにも読めなくもなかったので、 par.size = "startvalues"を追加して実行すると収束しました。
デフォルトのまま、モデル記述の初期値を変えても収束しました。(講師の質 <-> 講師の質, vq, 1)
CFI,RSMEAの出力はありません。
係数、平均値、切片、分散については、小数3桁まで一致します。
標準化推定値については、stdCoef(c5semi.sem)で出力した値は、全く違います。明らかに不自然な値てす。semパッケージでは、平均構造を導入したsemオブジェクトにstdCoef()を使って標準化係数を求めることができないのではないでしょうか。
＜ 下線部は、2013/8/31修正 ＞
2013年2月25日月曜日
共分散構造分析[R編]ってないの？[第４章6]
「共分散構造分析［Amos編］－構造方程式モデリング－」の84頁の多母集団因子分析。
semパッケージでのモデルの記述方法は複数あります。今回は因子分析なので確認的因子分析に便利なcfa()関数を使った記述方法を使ってみたいと思います。
まずは、データを読み込んで、"学年"という変数をファクターに変換し、その次にモデルを記述します。
これは楽ちんですね。"記号操作 -> 数的処理"のパスにはちゃんと制約がついていて、誤差分散もきっちり入っています。このc04fact.mdlAをつかって、集団ごとの分析と配置不変の検討に進みたいと思います。
次の計算では、sem()の中で共分散を計算する変数を指定するformulaオプションを書いています。デフォルトでは formula = ~. ということです。これまでの例では指定する必要がなかったのですが、このモデルで指定しないと「res_u %*% w_adf : 適切な引数ではありません」と叱られてしまいました。
適合度について、集団ごとの評価は記載のとおりとなりました。配置不変性を検討すると、GFIとCFIは記載どおり記載どおりの値です。RMSEAは記載よりも大きく、習慣的基準よりもすこし大きくなりました(習慣的基準とは0.05以下と理解していますが)。パス解析でもそうでしたが、semパッケージで多母集団分析でモデルから推定されるRMSEAはAmosのそれより大きくなる傾向があるのかもしれません。インチキな私は推計方法を理解していないので、詳しくはわかりません。
4.7の測定不変モデルを推計してみます。
ここで、cfa()でモデルを記述したことの失敗が判明しました。どうやら、semパッケージが自動で付与してくれたパラメータの名称は変更できないようなのです。削除して追加しなおすのもわかりにくいので、やむなく矢印式で書きなおし、update()を使って修正したモデルをもうひとつ記述することにします。
測定不変モデルは配置不変モデルと比べると、RMSEAとAICが小さくなっています。
本と比較するとAICは記載どおりの値ですが、RMSEAはやや大きく習慣的基準値を若干上回り、信頼区間はかなり広くなっています。パス係数はほぼ一致しています。本ではパス係数を制約して分散を計算させているのでsem()の出力とそのまま比較できますが、少し大きいようです。
ここまでの多母集団分析でパス係数やGFI、AICはAmosとほぼ一致するのですが、RMSEAは大きくなりました。もしsemパッケージでは必ずそうなるのならば、RMSEAはモデルの採用時にはほとんど常に参照される適合度指標になっているように思うので、少し厄介だなと感じます。
semパッケージでのモデルの記述方法は複数あります。今回は因子分析なので確認的因子分析に便利なcfa()関数を使った記述方法を使ってみたいと思います。
まずは、データを読み込んで、"学年"という変数をファクターに変換し、その次にモデルを記述します。
c04fact$学年 <- factor(c04fact$学年, labels = c("小学６年生","中学１年生")) #ファクターに変換
c04fact.mdlA <- cfa()
数的処理 : 記号操作, 図形, 計算, 文章問題
#
(c04fact.mdlA) #モデル記述を表示
c04fact.mdlA <- cfa()
数的処理 : 記号操作, 図形, 計算, 文章問題
#
(c04fact.mdlA) #モデル記述を表示
これは楽ちんですね。"記号操作 -> 数的処理"のパスにはちゃんと制約がついていて、誤差分散もきっちり入っています。このc04fact.mdlAをつかって、集団ごとの分析と配置不変の検討に進みたいと思います。
次の計算では、sem()の中で共分散を計算する変数を指定するformulaオプションを書いています。デフォルトでは formula = ~. ということです。これまでの例では指定する必要がなかったのですが、このモデルで指定しないと「res_u %*% w_adf : 適切な引数ではありません」と叱られてしまいました。
opt <- options(fit.indices = c("GFI","AGFI","RMSEA","CFI","AIC"))
c04fact.mgmdl0 <- multigroupModel(c04fact.mdlA, groups=levels(c04fact$学年))
summary(c04fact.gmsem0 <-
sem(c04fact.mgmdl0, data = c04fact, group = "学年",
formula = ~ 記号操作 + 計算 + 文章問題 + 図形))
c04fact.mgmdl0 <- multigroupModel(c04fact.mdlA, groups=levels(c04fact$学年))
summary(c04fact.gmsem0 <-
sem(c04fact.mgmdl0, data = c04fact, group = "学年",
formula = ~ 記号操作 + 計算 + 文章問題 + 図形))
適合度について、集団ごとの評価は記載のとおりとなりました。配置不変性を検討すると、GFIとCFIは記載どおり記載どおりの値です。RMSEAは記載よりも大きく、習慣的基準よりもすこし大きくなりました(習慣的基準とは0.05以下と理解していますが)。パス解析でもそうでしたが、semパッケージで多母集団分析でモデルから推定されるRMSEAはAmosのそれより大きくなる傾向があるのかもしれません。インチキな私は推計方法を理解していないので、詳しくはわかりません。
4.7の測定不変モデルを推計してみます。
ここで、cfa()でモデルを記述したことの失敗が判明しました。どうやら、semパッケージが自動で付与してくれたパラメータの名称は変更できないようなのです。削除して追加しなおすのもわかりにくいので、やむなく矢印式で書きなおし、update()を使って修正したモデルをもうひとつ記述することにします。
##測定不変モデル
c04fact.mdlA <- specifyModel()
数的処理 -> 記号操作, NA,1
数的処理 -> 図形, path1, NA
数的処理 -> 計算, path2, NA
数的処理 -> 文章問題, Path3, NA
数的処理 <-> 数的処理, V1, NA
記号操作 <-> 記号操作, V2, NA
図形 <-> 図形, V3, NA
計算 <-> 計算, V4, NA
文章問題 <-> 文章問題 ,V5, NA
#
c04fact.mdlB <- update(c04fact.mdlA)
replace, V1, V1B
replace, V2, V2B
replace, V3, V3B
replace, V4, V4B
replace, V5, V5B
#
c04fact.mdlA <- specifyModel()
数的処理 -> 記号操作, NA,1
数的処理 -> 図形, path1, NA
数的処理 -> 計算, path2, NA
数的処理 -> 文章問題, Path3, NA
数的処理 <-> 数的処理, V1, NA
記号操作 <-> 記号操作, V2, NA
図形 <-> 図形, V3, NA
計算 <-> 計算, V4, NA
文章問題 <-> 文章問題 ,V5, NA
#
c04fact.mdlB <- update(c04fact.mdlA)
replace, V1, V1B
replace, V2, V2B
replace, V3, V3B
replace, V4, V4B
replace, V5, V5B
#
c04fact.mgmdl1 <- multigroupModel(c04fact.mdlA, c04fact.mdlB, groups = levels(c04fact$学年))
summary(c04fact.gmsem1 <- sem(c04fact.mgmdl1, data = c04fact, group = "学年",
formula = ~ 記号操作 + 計算 + 文章問題 + 図形))
summary(c04fact.gmsem1 <- sem(c04fact.mgmdl1, data = c04fact, group = "学年",
formula = ~ 記号操作 + 計算 + 文章問題 + 図形))
測定不変モデルは配置不変モデルと比べると、RMSEAとAICが小さくなっています。
本と比較するとAICは記載どおりの値ですが、RMSEAはやや大きく習慣的基準値を若干上回り、信頼区間はかなり広くなっています。パス係数はほぼ一致しています。本ではパス係数を制約して分散を計算させているのでsem()の出力とそのまま比較できますが、少し大きいようです。
ここまでの多母集団分析でパス係数やGFI、AICはAmosとほぼ一致するのですが、RMSEAは大きくなりました。もしsemパッケージでは必ずそうなるのならば、RMSEAはモデルの採用時にはほとんど常に参照される適合度指標になっているように思うので、少し厄介だなと感じます。
2013年2月23日土曜日
共分散構造分析[R編]ってないの？[第5章1]
「共分散構造分析［Amos編］－構造方程式モデリング－」の90頁からの平均共分散構造分析を試してみます。
まず、モデルの記述は次のようにします。Interceptから平均や切片を計算する観測変数に向けて矢印を書きます。dataを入力する場合は、Interceptという変数名にします。sem()関数が自動でIntercept込のraw momentを計算してくれるようです。raw momentを別に計算してSに代入する場合は例えば rawMoment( cbind( semidata, 切片=1 ))などと計算すれば良いのではないかと思います。
平均構造を計算する場合、sem()関数では必ずrawオプションにTRUEを指定します。デフォルトではFALSEとなっていて分散共分散を計算します。このまま共分散行列を計算してもらうと逆行列が計算できませんということでエラーになります。平均構造を導入する場合は、ここをTRUEに指定して一次のモーメントを計算してもらいます。fixed.xオプションではInterceptを指定します。formulaオプションにはInterceptを指定しません。
これで、91頁のパス図に書かれている値が推定されていると思います。ちなみに、平均と切片を計算する場合、適合度指標は出力されません。Amosでも。
さて、パス図(左)を作成してみると大変見苦しくなってしまいました。ちょっと手を加えてみたのが右の図です。
まず、モデルの記述は次のようにします。Interceptから平均や切片を計算する観測変数に向けて矢印を書きます。dataを入力する場合は、Interceptという変数名にします。sem()関数が自動でIntercept込のraw momentを計算してくれるようです。raw momentを別に計算してSに代入する場合は例えば rawMoment( cbind( semidata, 切片=1 ))などと計算すれば良いのではないかと思います。
semi.mdl <- specifyModel()
理解度 -> 満足度, path1, NA
プレゼン -> 満足度, path2, NA
講師対処 -> 満足度, path3, NA
理解度 <-> プレゼン, c1, NA
理解度 <-> 講師対処, c2, NA
プレゼン <-> 講師対処, c3, NA
Intercept -> 満足度, i0, NA
Intercept -> 理解度, i1, NA
Intercept -> プレゼン, i2, NA
Intercept -> 講師対処, i3, NA
#
理解度 -> 満足度, path1, NA
プレゼン -> 満足度, path2, NA
講師対処 -> 満足度, path3, NA
理解度 <-> プレゼン, c1, NA
理解度 <-> 講師対処, c2, NA
プレゼン <-> 講師対処, c3, NA
Intercept -> 満足度, i0, NA
Intercept -> 理解度, i1, NA
Intercept -> プレゼン, i2, NA
Intercept -> 講師対処, i3, NA
#
平均構造を計算する場合、sem()関数では必ずrawオプションにTRUEを指定します。デフォルトではFALSEとなっていて分散共分散を計算します。このまま共分散行列を計算してもらうと逆行列が計算できませんということでエラーになります。平均構造を導入する場合は、ここをTRUEに指定して一次のモーメントを計算してもらいます。fixed.xオプションではInterceptを指定します。formulaオプションにはInterceptを指定しません。
summary(semi.sem <- sem(semi.mdl, data = semidata, raw = TRUE, fixed.x = "Intercept",
formula = ~ 満足度 + 理解度 + プレゼン + 講師対処 ))
formula = ~ 満足度 + 理解度 + プレゼン + 講師対処 ))
これで、91頁のパス図に書かれている値が推定されていると思います。ちなみに、平均と切片を計算する場合、適合度指標は出力されません。Amosでも。
さて、パス図(左)を作成してみると大変見苦しくなってしまいました。ちょっと手を加えてみたのが右の図です。
2013年2月16日土曜日
共分散構造分析[R編]ってないの？[第４章5]
「共分散構造分析［Amos編］－構造方程式モデリング－」の82頁の等値制約をおいた多母集団パス解析を行います。
モデルの一部に等値制約を付す場合は、モデルを複数記述して等値制約をパラメータの名前を同じにすると良いらしいので、ここでは2つのモデルを記述してみます。等値制約を置かないパラメータの名称には学校の名前(A,B)をつけて、名前を区別しました。
c04school.mdlA <- specifyModel() #A校
入学時勉強時間 -> 入学時学力テスト, a1, NA #等値制約
入学時勉強時間 -> 卒業時学力テスト, a2, NA #等値制約
入学後趣味活動 -> 卒業時学力テスト, a3, NA #等値制約
入学時学力テスト -> 卒業時学力テスト, Aa4, NA
入学時学力テスト <->入学時学力テスト, e1, NA #等値制約
卒業時学力テスト <-> 卒業時学力テスト, Ae2, NA
入学時勉強時間 <-> 入学時勉強時間, Avar1, NA
入学後趣味活動 <-> 入学後趣味活動, Avar2, NA
#
c04school.mdlB <- specifyModel() #B校
入学時勉強時間 -> 入学時学力テスト, a1, NA #等値制約
入学時勉強時間 -> 卒業時学力テスト, a2, NA #等値制約
入学後趣味活動 -> 卒業時学力テスト, a3, NA #等値制約
入学時学力テスト -> 卒業時学力テスト, Ba4, NA
入学時学力テスト <->入学時学力テスト, e1, NA #等値制約
卒業時学力テスト <-> 卒業時学力テスト, Be2, NA
入学時勉強時間 <-> 入学時勉強時間, Bvar1, NA
入学後趣味活動 <-> 入学後趣味活動, Bvar2, NA
#
入学時勉強時間 -> 入学時学力テスト, a1, NA #等値制約
入学時勉強時間 -> 卒業時学力テスト, a2, NA #等値制約
入学後趣味活動 -> 卒業時学力テスト, a3, NA #等値制約
入学時学力テスト -> 卒業時学力テスト, Aa4, NA
入学時学力テスト <->入学時学力テスト, e1, NA #等値制約
卒業時学力テスト <-> 卒業時学力テスト, Ae2, NA
入学時勉強時間 <-> 入学時勉強時間, Avar1, NA
入学後趣味活動 <-> 入学後趣味活動, Avar2, NA
#
c04school.mdlB <- specifyModel() #B校
入学時勉強時間 -> 入学時学力テスト, a1, NA #等値制約
入学時勉強時間 -> 卒業時学力テスト, a2, NA #等値制約
入学後趣味活動 -> 卒業時学力テスト, a3, NA #等値制約
入学時学力テスト -> 卒業時学力テスト, Ba4, NA
入学時学力テスト <->入学時学力テスト, e1, NA #等値制約
卒業時学力テスト <-> 卒業時学力テスト, Be2, NA
入学時勉強時間 <-> 入学時勉強時間, Bvar1, NA
入学後趣味活動 <-> 入学後趣味活動, Bvar2, NA
#
この2つのモデルの記述をmultigroupModel() で分析用のsemmodlistオブジェクトにし、sem()で計算させます。なお、今回は警告が出ないように、学校変数はA校、B校のラベルをつけてファクターに変換しておきます。
c04school$学校 <- factor(c04school$学校, labels = c("A校","B校")) #ファクター化
opt <- options(fit.indices = c("GFI","AGFI","RMSEA","CFI","AIC"))
#2つのモデルを読み込む
c04school.mgmdl2 <- multigroupModel( c04school.mdlA, c04school.mdlB, groups = c("A校","B校"))
summary(c04school.mgsem2 <- sem( c04school.mgmdl2, data = c04school, group = "学校"))
opt <- options(fit.indices = c("GFI","AGFI","RMSEA","CFI","AIC"))
#2つのモデルを読み込む
c04school.mgmdl2 <- multigroupModel( c04school.mdlA, c04school.mdlB, groups = c("A校","B校"))
summary(c04school.mgsem2 <- sem( c04school.mgmdl2, data = c04school, group = "学校"))
適合度をみると、RMSEAはAmosよりすこし大きく、もちろん「習慣的基準」を超えています。その他の適合度指標も悪化しました。
graphvizによるパス係数などの出力は下図のとおりです。両矢印は誤差の"分散"です。Amosでは、誤差分散を1に制約してパス係数を出力するのが普通のようです。Amos式だと、標準偏差相当の値が出力されます。
4.5.2その他の等値制約の置き方の節に、 「構造モデルのウェイト」ですべてのパス係数に等値制約を置くというのがあります。multigroupModel()のallEqualオプションをFALSEに設定することで可能です。
c04school.mgmdl3 <- multigroupModel(c04school.mdlA, groups = c("A校","B校"), allEqual=TRUE)
summary(c04school.mgsem2 <- sem(c04school.mgmdl3, datac= c04school, group = "学校"))
summary(c04school.mgsem2 <- sem(c04school.mgmdl3, datac= c04school, group = "学校"))
配置不変モデルと比較すると、適合度指標が悪化しています。これらのグループの集団の異質性を示唆する結果ということでしょうか。
2013年2月15日金曜日
共分散構造分析[R編]ってないの？[第４章1-4]
残念ながら、71頁図3.31のモデルで座礁してしまいました。これは宿題として積み残して、「共分散構造分析［Amos編］－構造方程式モデリング－」の74頁から81頁あたりをやってみます。多母集団分析は、昨年(2012)にsemパッケージに追加された機能のようです。
まずは、相関行列の確認から．．．。
検討を終えて、早速、モデルを記述します。
次に学校ごと、それぞれ計算します。他母集団分析の時は、集団を区別する変数をfactorにしておくことになっているようなので、予め"学校"という変数をfactorに変換します。学校ごとの分析では、"学校"変数はデータフレームから除いています。
AMOSと比べるとRMSEAが幾分異なります。
次に、配置不変のモデルを計算します。多母集団分析のモデルは multigroupModel()を使って指定します。配置不変の場合は、パス図が共通なのでパス図の記述は一つで、allEqual オプションをFALSEに指定しておけばOKのようです。allEqual=FALSEはデフォルト設定です。モデル記述は、上記のc04school.mdlaをそのまま用います。groupsオプションには、A校、B校と指定しましたが、 上記で因子に変換しておいた"学校"という変数の水準の順に出力が割り当てられるらしいです。
sem()の中では、dataとして因子化した"学校"変数を含むデータフレームを指定します。そして、groupオプションで集団を区別する変数名"学校"を指定します。
groupsの内容と学校のラベルが異なっているので警告が出ます。問題ないと思いますので無視します。最初に学校変数を因子化するときにラベルを指定するか、groups=c("1","2")としておけば警告は出ないでしょう。
適合度を確認すると、やはりRMSEAの違いが目立ちます。GFIとCFIは表示桁までは同じ値となっています。
さて、パラメータの差の検定をやってみます。たいした計算をしているわけでもないのに、ちょっと強引で汚い感じがします。エレガントなスクリプトかけませんねえ。
図4.15「パラメータ間の差に対する検定統計量」に該当する表が表示されました。変数の順序は違いますが、小数点3桁表示でほぼ同じ値といっても差し支えないと思います。a1とa3のパスでは、A校とB校の間で5%有意な差が認められます。
まずは、相関行列の確認から．．．。
cor(c04school[,-5])#併合データ
lapply(1:2, function(x) cor(c04school[c04school$学校==x,-5]))#学校別
lapply(1:2, function(x) cor(c04school[c04school$学校==x,-5]))#学校別
検討を終えて、早速、モデルを記述します。
c04school.mdla <- specifyModel()
入学時勉強時間 -> 入学時学力テスト, Aa1, NA
入学時勉強時間 -> 卒業時学力テスト, Aa2, NA
入学時学力テスト -> 卒業時学力テスト, Aa4, NA
入学後趣味活動 -> 卒業時学力テスト, Aa3, NA
入学時学力テスト <->入学時学力テスト, Ae1, NA
卒業時学力テスト <-> 卒業時学力テスト, Ae2, NA
入学時勉強時間 <-> 入学時勉強時間, Avar1, NA
入学後趣味活動 <-> 入学後趣味活動, Avar2, NA
入学時勉強時間 -> 入学時学力テスト, Aa1, NA
入学時勉強時間 -> 卒業時学力テスト, Aa2, NA
入学時学力テスト -> 卒業時学力テスト, Aa4, NA
入学後趣味活動 -> 卒業時学力テスト, Aa3, NA
入学時学力テスト <->入学時学力テスト, Ae1, NA
卒業時学力テスト <-> 卒業時学力テスト, Ae2, NA
入学時勉強時間 <-> 入学時勉強時間, Avar1, NA
入学後趣味活動 <-> 入学後趣味活動, Avar2, NA
次に学校ごと、それぞれ計算します。他母集団分析の時は、集団を区別する変数をfactorにしておくことになっているようなので、予め"学校"という変数をfactorに変換します。学校ごとの分析では、"学校"変数はデータフレームから除いています。
opt <- options(fit.indices = c("GFI","AGFI","RMSEA","CFI"))
c04school$学校 <- factor(c04school$学校) #学校区分を因子に変換
#A校
summary(c04schoolA.sem <-
sem(c04school.mdla, data = c04school[c04school$学校 == 1,-5]))
stdCoef(c04schoolA.sem)
#B校
summary(c04schoolB.sem <-
sem(c04school.mdla, data = c04school[c04school$学校 == 2,-5]))
stdCoef(c04schoolB.sem)
c04school$学校 <- factor(c04school$学校) #学校区分を因子に変換
#A校
summary(c04schoolA.sem <-
sem(c04school.mdla, data = c04school[c04school$学校 == 1,-5]))
stdCoef(c04schoolA.sem)
#B校
summary(c04schoolB.sem <-
sem(c04school.mdla, data = c04school[c04school$学校 == 2,-5]))
stdCoef(c04schoolB.sem)
AMOSと比べるとRMSEAが幾分異なります。
次に、配置不変のモデルを計算します。多母集団分析のモデルは multigroupModel()を使って指定します。配置不変の場合は、パス図が共通なのでパス図の記述は一つで、allEqual オプションをFALSEに指定しておけばOKのようです。allEqual=FALSEはデフォルト設定です。モデル記述は、上記のc04school.mdlaをそのまま用います。groupsオプションには、A校、B校と指定しましたが、 上記で因子に変換しておいた"学校"という変数の水準の順に出力が割り当てられるらしいです。
sem()の中では、dataとして因子化した"学校"変数を含むデータフレームを指定します。そして、groupオプションで集団を区別する変数名"学校"を指定します。
c04school.mgmdl <-
multigroupModel( c04school.mdla, groups = c("A校","B校"))
summary(c04school.mgsem <-
sem(c04school.mgmdl, data = c04school, group="学校"))
multigroupModel( c04school.mdla, groups = c("A校","B校"))
summary(c04school.mgsem <-
sem(c04school.mgmdl, data = c04school, group="学校"))
groupsの内容と学校のラベルが異なっているので警告が出ます。問題ないと思いますので無視します。最初に学校変数を因子化するときにラベルを指定するか、groups=c("1","2")としておけば警告は出ないでしょう。
適合度を確認すると、やはりRMSEAの違いが目立ちます。GFIとCFIは表示桁までは同じ値となっています。
さて、パラメータの差の検定をやってみます。たいした計算をしているわけでもないのに、ちょっと強引で汚い感じがします。エレガントなスクリプトかけませんねえ。
ecov <- vcov(c04school.mgsem)#推定された分散共分散行列
sde <- sqrt(sweep(sweep(-2 * ecov, 2, diag(ecov), "+"), 1, diag(ecov), "+")) #差の標準誤差
pdif <- matrix(coef(c04school.mgsem), 16, 16) -
matrix(coef(c04school.mgsem), 16, 16, byrow=TRUE) #パラメータの差
( pdif / sde )[-1:-8, 1:8] #差の検定統計量を計算して、必要な箇所を取り出す
sde <- sqrt(sweep(sweep(-2 * ecov, 2, diag(ecov), "+"), 1, diag(ecov), "+")) #差の標準誤差
pdif <- matrix(coef(c04school.mgsem), 16, 16) -
matrix(coef(c04school.mgsem), 16, 16, byrow=TRUE) #パラメータの差
( pdif / sde )[-1:-8, 1:8] #差の検定統計量を計算して、必要な箇所を取り出す
図4.15「パラメータ間の差に対する検定統計量」に該当する表が表示されました。変数の順序は違いますが、小数点3桁表示でほぼ同じ値といっても差し支えないと思います。a1とa3のパスでは、A校とB校の間で5%有意な差が認められます。
2013年2月4日月曜日
共分散構造分析[R編]ってないの？[第３章10]
「共分散構造分析［Amos編］－構造方程式モデリング－」の70頁、１因子モデルによる信頼性、図3.30を書いてみます。
適合度指標は表示されている少数第三位まで一致しています。パス係数はやや異なっています。
信頼性係数は、上記のように計算して丸めると、表示桁(少数第二位)まで一致します。
さて、 クロンバックのα信頼性係数を出力してくれるRの関数はいくつかあるようですが、psy パッケージの cronbach() という関数を使ってみます。
本書に書かれているα係数と値が異なります。
opt <- br="" fit.indices="c(" options="">
kajyou330.mdl <- specifyModel()
v1 <- 過剰適応傾向, a1, NA
v2 <- 過剰適応傾向, a2, NA
v3 <- 過剰適応傾向, a3, NA
v4 <- 過剰適応傾向, a4, NA
v5 <- 過剰適応傾向, a5, NA
v6 <- 過剰適応傾向, a6, NA
過剰適応傾向 <-> 過剰適応傾向, NA, 1
#
summary( kajyou330.sem <- sem( kajyou330.mdl, data = kajyou))
a <- sum( coef( kajyou330.sem )[1:6] )^2
b <- sum( coef( kajyou330.sem )[7:12] )
a / ( a + b )
kajyou330.mdl <- specifyModel()
v1 <- 過剰適応傾向, a1, NA
v2 <- 過剰適応傾向, a2, NA
v3 <- 過剰適応傾向, a3, NA
v4 <- 過剰適応傾向, a4, NA
v5 <- 過剰適応傾向, a5, NA
v6 <- 過剰適応傾向, a6, NA
過剰適応傾向 <-> 過剰適応傾向, NA, 1
#
summary( kajyou330.sem <- sem( kajyou330.mdl, data = kajyou))
a <- sum( coef( kajyou330.sem )[1:6] )^2
b <- sum( coef( kajyou330.sem )[7:12] )
a / ( a + b )
適合度指標は表示されている少数第三位まで一致しています。パス係数はやや異なっています。
信頼性係数は、上記のように計算して丸めると、表示桁(少数第二位)まで一致します。
さて、 クロンバックのα信頼性係数を出力してくれるRの関数はいくつかあるようですが、psy パッケージの cronbach() という関数を使ってみます。
cronbach( kajyou )
本書に書かれているα係数と値が異なります。
2013年1月31日木曜日
共分散構造分析[R編]ってないの？[第３章9]その2
「共分散構造分析［Amos編］－構造方程式モデリング－」の69頁、図3.27のファントム変数を含む因果モデルを計算するために、次のテキストファイル company327.txt を用意しました。
ファントム変数は母親の統制的態度因子の符号が反転しただけのものなので、ファントム変数へのパスを-1に固定するだけでなく、ファントム変数そのものの誤差分散を0に固定しています。
そして、次のスクリプトを実行しました。
と、最近のsemパッケージのアップデートで、適合度指標の出力にはおまじないが一行、またはオプションの設定が必要になったようです。一行目がそれです。sem のヘルプのサンプルコードやsummry.semのヘルプを参照してください。
今回はshift-jisが使われているWindowsではなくて、utf-8が標準エンコードのubuntuで作業するので、pathDiagram()もいきなり.dotファイル -> .pdfファイルに出力です。
さて、path1とpath2の平均の差の統計量は、推計されたパラメータの分散共分散行列を利用して、
z = (path1 - path2 ) / sqrt ( path1の誤差分散 - 2 * path1とpath2の誤差共分散 + path2の誤差分散)
で計算するらしい。
semパッケージでは推計されたパラメータの分散共分散行列は、
vcov(c3mother.sem) で得られる行列だと思うので、この行列を使って統計量を計算してみましょう。
本に書かれていた出力と随分違う。んんん。
もちろん、sqrt(diag(vc)[8:9] はsemパッケージのpath1とpath2の標準誤差の出力と一致している。何か大きな勘違いをしているのでしょうか。AMOSの非標準化解はどうなっているんでしょう。いかん！おいらがインチキだってことがバレバレです。semパッケージの標準誤差の計算方法は今後の宿題にしておこう。
図3.28。等値制約を置いたモデルを試してみます。パスに同じ名前をつければ等値制約がかかるはずなので、c3mother327のモデルのpath2をpath1 に書き換えただけのファイルをsem()に計算してもらいました。
非標準化解では、「母親の絶対視 <--- F1」の値だけが表示されています。
標準化解を出力すると、「母親の絶対視 <--- F1」と「母親の絶対視 <--- 母親の同一化傾向」は一致していません。他の推計値も「Amos編」とは微妙に違う。
ファントム変数は母親の統制的態度因子の符号が反転しただけのものなので、ファントム変数へのパスを-1に固定するだけでなく、ファントム変数そのものの誤差分散を0に固定しています。
#観測方程式
母親の統制的態度 -> v1, NA, 1
母親の統制的態度 -> v2, p2, NA
母親の統制的態度 -> v3, p3, NA
母親の同一化傾向 -> v4, NA, 1
母親の同一化傾向 -> v5, p5, NA
母親の同一化傾向 -> v6, p6, NA
母親の絶対視 -> v7, NA, 1
母親の絶対視 -> v8, p8, NA
母親の絶対視 -> v9, p9, NA
#構成概念間の関係
母親の統制的態度 <-> 母親の同一化傾向, cor1, NA
母親の統制的態度 -> F1, NA, -1 #ファントム変数へのパスを-1に制約
F1 -> 母親の絶対視, path1, NA #ファントム変数からのパスは自由に推計
母親の同一化傾向 -> 母親の絶対視, path2, NA
#因子の分散
母親の統制的態度 <-> 母親の統制的態度, var1, NA
母親の同一化傾向 <-> 母親の同一化傾向, var2, NA
母親の絶対視 <-> 母親の絶対視, var3, NA
F1 <-> F1, NA, 0 #ファントム変数の誤差分散は0
母親の統制的態度 -> v2, p2, NA
母親の統制的態度 -> v3, p3, NA
母親の同一化傾向 -> v4, NA, 1
母親の同一化傾向 -> v5, p5, NA
母親の同一化傾向 -> v6, p6, NA
母親の絶対視 -> v7, NA, 1
母親の絶対視 -> v8, p8, NA
母親の絶対視 -> v9, p9, NA
#構成概念間の関係
母親の統制的態度 <-> 母親の同一化傾向, cor1, NA
母親の統制的態度 -> F1, NA, -1 #ファントム変数へのパスを-1に制約
F1 -> 母親の絶対視, path1, NA #ファントム変数からのパスは自由に推計
母親の同一化傾向 -> 母親の絶対視, path2, NA
#因子の分散
母親の統制的態度 <-> 母親の統制的態度, var1, NA
母親の同一化傾向 <-> 母親の同一化傾向, var2, NA
母親の絶対視 <-> 母親の絶対視, var3, NA
F1 <-> F1, NA, 0 #ファントム変数の誤差分散は0
そして、次のスクリプトを実行しました。
と、最近のsemパッケージのアップデートで、適合度指標の出力にはおまじないが一行、またはオプションの設定が必要になったようです。一行目がそれです。sem のヘルプのサンプルコードやsummry.semのヘルプを参照してください。
今回はshift-jisが使われているWindowsではなくて、utf-8が標準エンコードのubuntuで作業するので、pathDiagram()もいきなり.dotファイル -> .pdfファイルに出力です。
opt <- options( fit.indices = c("GFI","AGFI","CFI","RMSEA"))
c3mother327.mdl <- specifyModel("model/c3mother327.txt")
summary(c3mother327.sem <- sem(c3mother327.mdl, data=c3mother))
stdCoef(c3mother327.sem)
pathDiagram(c3mother327.sem, edge.labels="values", standardize=TRUE,
node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
rank.direction="LR", ignore.double=FALSE,
same.rank=c("母親の統制的態度,母親の同一化傾向"), #並べるノード
min.rank=c("v1,v2,v3,v4,v5,v6"), #最も左に配置するノード
file="graphviz/c3mother327") #Windowsではshift-jifで処理される
c3mother327.mdl <- specifyModel("model/c3mother327.txt")
summary(c3mother327.sem <- sem(c3mother327.mdl, data=c3mother))
stdCoef(c3mother327.sem)
pathDiagram(c3mother327.sem, edge.labels="values", standardize=TRUE,
node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
rank.direction="LR", ignore.double=FALSE,
same.rank=c("母親の統制的態度,母親の同一化傾向"), #並べるノード
min.rank=c("v1,v2,v3,v4,v5,v6"), #最も左に配置するノード
file="graphviz/c3mother327") #Windowsではshift-jifで処理される
さて、path1とpath2の平均の差の統計量は、推計されたパラメータの分散共分散行列を利用して、
z = (path1 - path2 ) / sqrt ( path1の誤差分散 - 2 * path1とpath2の誤差共分散 + path2の誤差分散)
で計算するらしい。
semパッケージでは推計されたパラメータの分散共分散行列は、
vcov(c3mother.sem) で得られる行列だと思うので、この行列を使って統計量を計算してみましょう。
vc <- vcov(c3mother327.sem)
coeffs <- c3mother327.sem$coef
(coeffs[8] - coeffs[9]) / sqrt(vc[8,8] + vc[9,9] - 2 * vc[8,9])
path1
-1.276916
coeffs <- c3mother327.sem$coef
(coeffs[8] - coeffs[9]) / sqrt(vc[8,8] + vc[9,9] - 2 * vc[8,9])
path1
-1.276916
本に書かれていた出力と随分違う。んんん。
もちろん、sqrt(diag(vc)[8:9] はsemパッケージのpath1とpath2の標準誤差の出力と一致している。何か大きな勘違いをしているのでしょうか。AMOSの非標準化解はどうなっているんでしょう。いかん！おいらがインチキだってことがバレバレです。semパッケージの標準誤差の計算方法は今後の宿題にしておこう。
図3.28。等値制約を置いたモデルを試してみます。パスに同じ名前をつければ等値制約がかかるはずなので、c3mother327のモデルのpath2をpath1 に書き換えただけのファイルをsem()に計算してもらいました。
非標準化解では、「母親の絶対視 <--- F1」の値だけが表示されています。
標準化解を出力すると、「母親の絶対視 <--- F1」と「母親の絶対視 <--- 母親の同一化傾向」は一致していません。他の推計値も「Amos編」とは微妙に違う。
2013年1月27日日曜日
共分散構造分析[R編]ってないの？[第３章9]その1
「共分散構造分析［Amos編］－構造方程式モデリング－」の68頁、図3.26の因果モデルを計算してみます。
次のテキストファイル company326.txt を用意しました。
そして、次のスクリプトを実行しました。
company326.dotファイルは、見栄えのために赤い箇所を書き換えました。
digraph "c3mother1.sem" {
rankdir=LR;
size="8,8";
node [fontname="IPAGothic" fontsize=12 shape=box];
edge [fontname="IPAMincho" fontsize=10];
center=1;
{rank=same "母親の統制的態度" "母親の同一化傾向"}
"母親の同一化傾向" [shape=ellipse]
"母親の絶対視" [shape=ellipse]
"母親の統制的態度" [shape=ellipse]
"v1" -> "母親の統制的態度" [label="0.8" dir=back];
"v2" -> "母親の統制的態度" [label="0.68" dir=back];
"v3" -> "母親の統制的態度" [label="0.6" dir=back];
"v4" -> "母親の同一化傾向" [label="0.73" dir=back];
"v5" -> "母親の同一化傾向" [label="0.87" dir=back];
"v6" -> "母親の同一化傾向" [label="0.7" dir=back];
"母親の絶対視" -> "v7" [label="0.88"];
"母親の絶対視" -> "v8" [label="0.7"];
"母親の絶対視" -> "v9" [label="0.92"];
"母親の統制的態度" -> "母親の同一化傾向" [label="0.31" dir=both headport=nw tailport=sw];
"母親の統制的態度" -> "母親の絶対視" [label="-0.29"];
"母親の同一化傾向" -> "母親の絶対視" [label="0.37"];
"母親の統制的態度" -> "母親の統制的態度" [label="1" dir=both];
"母親の同一化傾向" -> "母親の同一化傾向" [label="1" dir=both];
"母親の絶対視" -> "母親の絶対視" [label="0.85" dir=both];
"v1" -> "v1" [label="0.36" dir=both];
"v2" -> "v2" [label="0.54" dir=both];
"v3" -> "v3" [label="0.64" dir=both];
"v4" -> "v4" [label="0.47" dir=both];
"v5" -> "v5" [label="0.24" dir=both];
"v6" -> "v6" [label="0.51" dir=both];
"v7" -> "v7" [label="0.23" dir=both];
"v8" -> "v8" [label="0.51" dir=both];
"v9" -> "v9" [label="0.15" dir=both];
fontname="IPAGothic" ;
fontsize=15;
label="構成概念間の因果モデル";
labelloc="t";
}
以下を実行して、パス図を作成します。
次のテキストファイル company326.txt を用意しました。
#観測方程式
母親の統制的態度 -> v1, NA, 1
母親の統制的態度 -> v2, p2, NA
母親の統制的態度 -> v3, p3, NA
母親の同一化傾向 -> v4, NA, 1
母親の同一化傾向 -> v5, p5, NA
母親の同一化傾向 -> v6, p6, NA
母親の絶対視 -> v7, NA, 1
母親の絶対視 -> v8, p8, NA
母親の絶対視 -> v9, p9, NA
#構成概念間の関係
母親の統制的態度 <-> 母親の同一化傾向, cor1, NA
母親の統制的態度 -> 母親の絶対視, path1, NA
母親の同一化傾向 -> 母親の絶対視, path2, NA
#因子の分散
母親の統制的態度 <-> 母親の統制的態度, v1, NA
母親の同一化傾向 <-> 母親の同一化傾向, v2, NA
母親の絶対視 <-> 母親の絶対視, v3, NA
母親の統制的態度 -> v2, p2, NA
母親の統制的態度 -> v3, p3, NA
母親の同一化傾向 -> v4, NA, 1
母親の同一化傾向 -> v5, p5, NA
母親の同一化傾向 -> v6, p6, NA
母親の絶対視 -> v7, NA, 1
母親の絶対視 -> v8, p8, NA
母親の絶対視 -> v9, p9, NA
#構成概念間の関係
母親の統制的態度 <-> 母親の同一化傾向, cor1, NA
母親の統制的態度 -> 母親の絶対視, path1, NA
母親の同一化傾向 -> 母親の絶対視, path2, NA
#因子の分散
母親の統制的態度 <-> 母親の統制的態度, v1, NA
母親の同一化傾向 <-> 母親の同一化傾向, v2, NA
母親の絶対視 <-> 母親の絶対視, v3, NA
そして、次のスクリプトを実行しました。
c3mother1.mdl <- specifyModel("model/c3mother1.txt")
summary(c3mother1.sem <- sem(c3mother1.mdl, data=c3mother))
stdCoef(c3mother1.sem)
pathDiagram(c3mother1.sem, edge.labels="values", standardize=TRUE,
node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
rank.direction="LR", ignore.double=FALSE,
same.rank=c("母親の統制的態度,母親の同一化傾向"),
min.rank=c("v1,v2,v3,v4,v5,v6") )
summary(c3mother1.sem <- sem(c3mother1.mdl, data=c3mother))
stdCoef(c3mother1.sem)
pathDiagram(c3mother1.sem, edge.labels="values", standardize=TRUE,
node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
rank.direction="LR", ignore.double=FALSE,
same.rank=c("母親の統制的態度,母親の同一化傾向"),
min.rank=c("v1,v2,v3,v4,v5,v6") )
company326.dotファイルは、見栄えのために赤い箇所を書き換えました。
- headport=nw tailport=sw は矢印が因子に刺さる位置を指定して矢印を曲線にする
- 終わりの部分はタイトルを表示する
digraph "c3mother1.sem" {
rankdir=LR;
size="8,8";
node [fontname="IPAGothic" fontsize=12 shape=box];
edge [fontname="IPAMincho" fontsize=10];
center=1;
{rank=same "母親の統制的態度" "母親の同一化傾向"}
"母親の同一化傾向" [shape=ellipse]
"母親の絶対視" [shape=ellipse]
"母親の統制的態度" [shape=ellipse]
"v1" -> "母親の統制的態度" [label="0.8" dir=back];
"v2" -> "母親の統制的態度" [label="0.68" dir=back];
"v3" -> "母親の統制的態度" [label="0.6" dir=back];
"v4" -> "母親の同一化傾向" [label="0.73" dir=back];
"v5" -> "母親の同一化傾向" [label="0.87" dir=back];
"v6" -> "母親の同一化傾向" [label="0.7" dir=back];
"母親の絶対視" -> "v7" [label="0.88"];
"母親の絶対視" -> "v8" [label="0.7"];
"母親の絶対視" -> "v9" [label="0.92"];
"母親の統制的態度" -> "母親の同一化傾向" [label="0.31" dir=both headport=nw tailport=sw];
"母親の統制的態度" -> "母親の絶対視" [label="-0.29"];
"母親の同一化傾向" -> "母親の絶対視" [label="0.37"];
"母親の統制的態度" -> "母親の統制的態度" [label="1" dir=both];
"母親の同一化傾向" -> "母親の同一化傾向" [label="1" dir=both];
"母親の絶対視" -> "母親の絶対視" [label="0.85" dir=both];
"v1" -> "v1" [label="0.36" dir=both];
"v2" -> "v2" [label="0.54" dir=both];
"v3" -> "v3" [label="0.64" dir=both];
"v4" -> "v4" [label="0.47" dir=both];
"v5" -> "v5" [label="0.24" dir=both];
"v6" -> "v6" [label="0.51" dir=both];
"v7" -> "v7" [label="0.23" dir=both];
"v8" -> "v8" [label="0.51" dir=both];
"v9" -> "v9" [label="0.15" dir=both];
fontname="IPAGothic" ;
fontsize=15;
label="構成概念間の因果モデル";
labelloc="t";
}
以下を実行して、パス図を作成します。
dot -Tpdf < company325.dot > company326.pdf
2013年1月19日土曜日
共分散構造分析[R編]ってないの？[第３章8]その2
「共分散構造分析［Amos編］－構造方程式モデリング－」の67頁、図3.25階層因子分析です。
次のテキストファイル company325.txt を用意しました。
そして、次のスクリプトを実行しました。
company325.dotファイルは出力のママだと、下位因子が一般因子が同じ並びになるので、赤い箇所の矢印の向きを反転し、dir=back を書き加えてみました。ついでにラベルを書き加えてタイトルまでおごってしまいます。
で、こんなふうな絵が描けました。
次のテキストファイル company325.txt を用意しました。
#下位因子
V11 <- 自信, NA ,1
V10 <- 自信, s10, NA
V9 <- 自信, s9, NA
V8 <- 自信, s8, NA
V7 <- 自信, s7, NA
V6 <- 自信, s6, NA
V5 <- 自信, s5, NA
V4 <- 自信, s4, NA
V3 <- 自信, s3, NA
V2 <- 自信, s2, NA
V1 <- 自信, s1, NA
#一般因子
V11 <- 仲間関係, NA, 1
V10 <- 仲間関係, r1, NA
V9 <- 仲間関係, r2, NA
V8 <- 身体外見, NA, 1
V7 <- 身体外見, l1, NA
V6 <- 身体外見, l2, NA
V5 <- 身体外見, l3, NA
V4 <- 身体能力, NA, 1
V3 <- 身体能力, a1, NA
V2 <- 身体能力, a2, NA
V1 <- 身体能力, a3, NA
#因子の分散
仲間関係 <-> 仲間関係, d1, NA
身体外見 <-> 身体外見, d2, NA
身体能力 <-> 身体能力, d3, NA
自信 <-> 自信, d4, NA
#因子の間の相関
仲間関係 <-> 身体外見, cor1
身体能力 <-> 仲間関係, cor2
身体外見 <-> 身体能力, cor3
V11 <- 自信, NA ,1
V10 <- 自信, s10, NA
V9 <- 自信, s9, NA
V8 <- 自信, s8, NA
V7 <- 自信, s7, NA
V6 <- 自信, s6, NA
V5 <- 自信, s5, NA
V4 <- 自信, s4, NA
V3 <- 自信, s3, NA
V2 <- 自信, s2, NA
V1 <- 自信, s1, NA
#一般因子
V11 <- 仲間関係, NA, 1
V10 <- 仲間関係, r1, NA
V9 <- 仲間関係, r2, NA
V8 <- 身体外見, NA, 1
V7 <- 身体外見, l1, NA
V6 <- 身体外見, l2, NA
V5 <- 身体外見, l3, NA
V4 <- 身体能力, NA, 1
V3 <- 身体能力, a1, NA
V2 <- 身体能力, a2, NA
V1 <- 身体能力, a3, NA
#因子の分散
仲間関係 <-> 仲間関係, d1, NA
身体外見 <-> 身体外見, d2, NA
身体能力 <-> 身体能力, d3, NA
自信 <-> 自信, d4, NA
#因子の間の相関
仲間関係 <-> 身体外見, cor1
身体能力 <-> 仲間関係, cor2
身体外見 <-> 身体能力, cor3
そして、次のスクリプトを実行しました。
ccompany325.mdl <- specifyModel("model/company325.txt")
summary(company325.sem <- sem(company325.mdl, S = cor.company, N = 100, standardized = TRUE))
stdCoef(company325.sem)
pathDiagram(company325.sem, edge.labels = "values",standardize = TRUE,
node.font = c("IPAGothic",12), edge.font = c("IPAMincho",10),
rank.direction = "TB", ignore.double = FALSE,
same.rank = c("仲間関係,身体外見,身体能力", "自信"))
summary(company325.sem <- sem(company325.mdl, S = cor.company, N = 100, standardized = TRUE))
stdCoef(company325.sem)
pathDiagram(company325.sem, edge.labels = "values",standardize = TRUE,
node.font = c("IPAGothic",12), edge.font = c("IPAMincho",10),
rank.direction = "TB", ignore.double = FALSE,
same.rank = c("仲間関係,身体外見,身体能力", "自信"))
company325.dotファイルは出力のママだと、下位因子が一般因子が同じ並びになるので、赤い箇所の矢印の向きを反転し、dir=back を書き加えてみました。ついでにラベルを書き加えてタイトルまでおごってしまいます。
digraph "company325.sem" {
rankdir=TB;
size="8,8";
node [fontname="IPAGothic" fontsize=12 shape=box];
edge [fontname="IPAMincho" fontsize=10];
center=1;
{rank=same "仲間関係" "身体外見" "身体能力"}
"仲間関係" [shape=ellipse]
"身体外見" [shape=ellipse]
"身体能力" [shape=ellipse]
"自信" [shape=ellipse]
"V11" -> "自信" [label="0.63" dir=back];
"V10" -> "自信" [label="0.29" dir=back];
"V9" -> "自信" [label="0.21" dir=back];
"V8" -> "自信" [label="0.51" dir=back];
"V7" -> "自信" [label="0.6" dir=back];
"V6" -> "自信" [label="-0.03" dir=back];
"V5" -> "自信" [label="0.02" dir=back];
"V4" -> "自信" [label="0.38" dir=back];
"V3" -> "自信" [label="0.64" dir=back];
"V2" -> "自信" [label="0.5" dir=back];
"V1" -> "自信" [label="0.2" dir=back];
"仲間関係" -> "V11" [label="0.55"];
"仲間関係" -> "V10" [label="0.82"];
"仲間関係" -> "V9" [label="0.61"];
"身体外見" -> "V8" [label="0.72"];
"身体外見" -> "V7" [label="0.54"];
"身体外見" -> "V6" [label="0.72"];
"身体外見" -> "V5" [label="0.89"];
"身体能力" -> "V4" [label="0.78"];
"身体能力" -> "V3" [label="0.45"];
"身体能力" -> "V2" [label="0.61"];
"身体能力" -> "V1" [label="0.66"];
"仲間関係" -> "仲間関係" [label="1" dir=both];
"身体外見" -> "身体外見" [label="1" dir=both];
"身体能力" -> "身体能力" [label="1" dir=both];
"自信" -> "自信" [label="1" dir=both];
"仲間関係" -> "身体外見" [label="0.69" dir=both];
"身体能力" -> "仲間関係" [label="0.4" dir=both];
"身体外見" -> "身体能力" [label="0.32" dir=both];
"V11" -> "V11" [label="0.3" dir=both];
"V10" -> "V10" [label="0.24" dir=both];
"V9" -> "V9" [label="0.58" dir=both];
"V8" -> "V8" [label="0.22" dir=both];
"V7" -> "V7" [label="0.35" dir=both];
"V6" -> "V6" [label="0.48" dir=both];
"V5" -> "V5" [label="0.21" dir=both];
"V4" -> "V4" [label="0.24" dir=both];
"V3" -> "V3" [label="0.38" dir=both];
"V2" -> "V2" [label="0.38" dir=both];
"V1" -> "V1" [label="0.52" dir=both];
fontname = "IPAGothic" ;
fontsize = 15;
label = "階層因子分析による結果";
labelloc = "t"
}
rankdir=TB;
size="8,8";
node [fontname="IPAGothic" fontsize=12 shape=box];
edge [fontname="IPAMincho" fontsize=10];
center=1;
{rank=same "仲間関係" "身体外見" "身体能力"}
"仲間関係" [shape=ellipse]
"身体外見" [shape=ellipse]
"身体能力" [shape=ellipse]
"自信" [shape=ellipse]
"V11" -> "自信" [label="0.63" dir=back];
"V10" -> "自信" [label="0.29" dir=back];
"V9" -> "自信" [label="0.21" dir=back];
"V8" -> "自信" [label="0.51" dir=back];
"V7" -> "自信" [label="0.6" dir=back];
"V6" -> "自信" [label="-0.03" dir=back];
"V5" -> "自信" [label="0.02" dir=back];
"V4" -> "自信" [label="0.38" dir=back];
"V3" -> "自信" [label="0.64" dir=back];
"V2" -> "自信" [label="0.5" dir=back];
"V1" -> "自信" [label="0.2" dir=back];
"仲間関係" -> "V11" [label="0.55"];
"仲間関係" -> "V10" [label="0.82"];
"仲間関係" -> "V9" [label="0.61"];
"身体外見" -> "V8" [label="0.72"];
"身体外見" -> "V7" [label="0.54"];
"身体外見" -> "V6" [label="0.72"];
"身体外見" -> "V5" [label="0.89"];
"身体能力" -> "V4" [label="0.78"];
"身体能力" -> "V3" [label="0.45"];
"身体能力" -> "V2" [label="0.61"];
"身体能力" -> "V1" [label="0.66"];
"仲間関係" -> "仲間関係" [label="1" dir=both];
"身体外見" -> "身体外見" [label="1" dir=both];
"身体能力" -> "身体能力" [label="1" dir=both];
"自信" -> "自信" [label="1" dir=both];
"仲間関係" -> "身体外見" [label="0.69" dir=both];
"身体能力" -> "仲間関係" [label="0.4" dir=both];
"身体外見" -> "身体能力" [label="0.32" dir=both];
"V11" -> "V11" [label="0.3" dir=both];
"V10" -> "V10" [label="0.24" dir=both];
"V9" -> "V9" [label="0.58" dir=both];
"V8" -> "V8" [label="0.22" dir=both];
"V7" -> "V7" [label="0.35" dir=both];
"V6" -> "V6" [label="0.48" dir=both];
"V5" -> "V5" [label="0.21" dir=both];
"V4" -> "V4" [label="0.24" dir=both];
"V3" -> "V3" [label="0.38" dir=both];
"V2" -> "V2" [label="0.38" dir=both];
"V1" -> "V1" [label="0.52" dir=both];
fontname = "IPAGothic" ;
fontsize = 15;
label = "階層因子分析による結果";
labelloc = "t"
}
dot -Tpdf < company325.dot > company325.pdf
で、こんなふうな絵が描けました。
2013年1月16日水曜日
共分散構造分析[R編]ってないの？[第３章8]その1
「共分散構造分析［Amos編］－構造方程式モデリング－」の66頁、２次因子分析です。
モデルがやや複雑になってきたので、モデルファイルを別に作成します。
観測変数の誤差分散は specifyModel() が自動的に補完してくれるので、省略します。
用意したファイル company324.txt の中身は次のとおりです。
specifyModel() でファイルを読み込むときには、ファイルに空白行があっても無視されるし、
#でコメントアウトすることもできるみたい。
次のスクリプトで計算してもらいます。
graphvizで絵を描いてみましょう。誤差分散も出力します。英語以外のフォントを使うときには、フォントの指定をします。
こんなの出たら、company324.dot に utf-8 のエンコードで保存します。
そして、コマンドプロンプトから以下を実行。
サイズの小さいラスター画像は貧相ですが、pdfはうつくしーですね。
モデルがやや複雑になってきたので、モデルファイルを別に作成します。
観測変数の誤差分散は specifyModel() が自動的に補完してくれるので、省略します。
用意したファイル company324.txt の中身は次のとおりです。
specifyModel() でファイルを読み込むときには、ファイルに空白行があっても無視されるし、
#でコメントアウトすることもできるみたい。
#２次因子
仲間関係 <- 仲間評価, NA, 1
身体外見 <- 仲間評価, p1, NA
身体能力 <- 仲間評価, p2, NA
#一次因子
V11 <- 仲間関係, NA, 1
V10 <- 仲間関係, r1, NA
V9<- 仲間関係, r2, NA
V8 <- 身体外見, NA, 1
V7 <- 身体外見, l1, NA
V6 <- 身体外見, l2, NA
V5 <- 身体外見, l3, NA
V4 <- 身体能力, NA, 1
V3 <- 身体能力, a1, NA
V2 <- 身体能力, a2, NA
V1 <- 身体能力, a3, NA
仲間関係 <-> 仲間関係, d3, NA
身体外見 <-> 身体外見, d2, NA
身体能力 <-> 身体能力, d1, NA
仲間評価 <-> 仲間評価, d0, 1
仲間関係 <- 仲間評価, NA, 1
身体外見 <- 仲間評価, p1, NA
身体能力 <- 仲間評価, p2, NA
#一次因子
V11 <- 仲間関係, NA, 1
V10 <- 仲間関係, r1, NA
V9<- 仲間関係, r2, NA
V8 <- 身体外見, NA, 1
V7 <- 身体外見, l1, NA
V6 <- 身体外見, l2, NA
V5 <- 身体外見, l3, NA
V4 <- 身体能力, NA, 1
V3 <- 身体能力, a1, NA
V2 <- 身体能力, a2, NA
V1 <- 身体能力, a3, NA
仲間関係 <-> 仲間関係, d3, NA
身体外見 <-> 身体外見, d2, NA
身体能力 <-> 身体能力, d1, NA
仲間評価 <-> 仲間評価, d0, 1
次のスクリプトで計算してもらいます。
company324.mdl <- specifyModel("model/company324.txt")
summary(company324.sem <- sem(company324.mdl, S=cor.company, N=100))
stdCoef(company324.sem)
summary(company324.sem <- sem(company324.mdl, S=cor.company, N=100))
stdCoef(company324.sem)
graphvizで絵を描いてみましょう。誤差分散も出力します。英語以外のフォントを使うときには、フォントの指定をします。
pathDiagram(company324.sem, edge.labels="values",standardize=TRUE,
node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
rank.direction="TB", ignore.double=FALSE,
same.rank=c("仲間関係,身体外見,身体能力"))
node.font=c("IPAGothic",12), edge.font=c("IPAMincho",10),
rank.direction="TB", ignore.double=FALSE,
same.rank=c("仲間関係,身体外見,身体能力"))
こんなの出たら、company324.dot に utf-8 のエンコードで保存します。
digraph "company324.sem" {
rankdir=TB;
size="8,8";
node [fontname="IPAGothic" fontsize=12 shape=box];
edge [fontname="IPAMincho" fontsize=10];
center=1;
{rank=same "仲間関係" "身体外見" "身体能力"}
"仲間関係" [shape=ellipse]
"身体外見" [shape=ellipse]
"身体能力" [shape=ellipse]
"仲間評価" [shape=ellipse]
"仲間評価" -> "仲間関係" [label="0.97"];
"仲間評価" -> "身体外見" [label="0.83"];
"仲間評価" -> "身体能力" [label="0.6"];
"仲間関係" -> "V11" [label="0.8"];
"仲間関係" -> "V10" [label="0.81"];
"仲間関係" -> "V9" [label="0.63"];
"身体外見" -> "V8" [label="0.91"];
"身体外見" -> "V7" [label="0.76"];
"身体外見" -> "V6" [label="0.58"];
"身体外見" -> "V5" [label="0.73"];
"身体能力" -> "V4" [label="0.84"];
"身体能力" -> "V3" [label="0.74"];
"身体能力" -> "V2" [label="0.8"];
"身体能力" -> "V1" [label="0.64"];
"仲間関係" -> "仲間関係" [label="0.07" dir=both];
"身体外見" -> "身体外見" [label="0.31" dir=both];
"身体能力" -> "身体能力" [label="0.64" dir=both];
"仲間評価" -> "仲間評価" [label="1" dir=both];
"V11" -> "V11" [label="0.37" dir=both];
"V10" -> "V10" [label="0.34" dir=both];
"V9" -> "V9" [label="0.6" dir=both];
"V8" -> "V8" [label="0.18" dir=both];
"V7" -> "V7" [label="0.43" dir=both];
"V6" -> "V6" [label="0.67" dir=both];
"V5" -> "V5" [label="0.47" dir=both];
"V4" -> "V4" [label="0.3" dir=both];
"V3" -> "V3" [label="0.45" dir=both];
"V2" -> "V2" [label="0.36" dir=both];
"V1" -> "V1" [label="0.59" dir=both];
}
rankdir=TB;
size="8,8";
node [fontname="IPAGothic" fontsize=12 shape=box];
edge [fontname="IPAMincho" fontsize=10];
center=1;
{rank=same "仲間関係" "身体外見" "身体能力"}
"仲間関係" [shape=ellipse]
"身体外見" [shape=ellipse]
"身体能力" [shape=ellipse]
"仲間評価" [shape=ellipse]
"仲間評価" -> "仲間関係" [label="0.97"];
"仲間評価" -> "身体外見" [label="0.83"];
"仲間評価" -> "身体能力" [label="0.6"];
"仲間関係" -> "V11" [label="0.8"];
"仲間関係" -> "V10" [label="0.81"];
"仲間関係" -> "V9" [label="0.63"];
"身体外見" -> "V8" [label="0.91"];
"身体外見" -> "V7" [label="0.76"];
"身体外見" -> "V6" [label="0.58"];
"身体外見" -> "V5" [label="0.73"];
"身体能力" -> "V4" [label="0.84"];
"身体能力" -> "V3" [label="0.74"];
"身体能力" -> "V2" [label="0.8"];
"身体能力" -> "V1" [label="0.64"];
"仲間関係" -> "仲間関係" [label="0.07" dir=both];
"身体外見" -> "身体外見" [label="0.31" dir=both];
"身体能力" -> "身体能力" [label="0.64" dir=both];
"仲間評価" -> "仲間評価" [label="1" dir=both];
"V11" -> "V11" [label="0.37" dir=both];
"V10" -> "V10" [label="0.34" dir=both];
"V9" -> "V9" [label="0.6" dir=both];
"V8" -> "V8" [label="0.18" dir=both];
"V7" -> "V7" [label="0.43" dir=both];
"V6" -> "V6" [label="0.67" dir=both];
"V5" -> "V5" [label="0.47" dir=both];
"V4" -> "V4" [label="0.3" dir=both];
"V3" -> "V3" [label="0.45" dir=both];
"V2" -> "V2" [label="0.36" dir=both];
"V1" -> "V1" [label="0.59" dir=both];
}
そして、コマンドプロンプトから以下を実行。
dot -Tpdf company324.pdf
2013年1月14日月曜日
共分散構造分析[R編]ってないの？[第３章7]その2
しばらく放置していましたが、約３ヶ月ぶりに再開です。
「共分散構造分析［Amos編］－構造方程式モデリング－」の65頁、図3.22の確認的因子分析です。
因子からのパスを1、誤差分散を0に制約する方法という前者の方法です。
チラと適合度指標を表示すると
以前あるセミナーで、因子分析も適合度を確認するので、学生さんが論文を書けなくて大変なんだということを聞きました。実際、アンケートで集めたデータで因子分析してみると、まずいなって思うことがあります。ていうか、ほとんどいつも悲惨。コストかけて、質の悪いデータを集めたと告白すべきでしょうか。
さて、信頼度係数ρ=0.81ということにして、sqrt(ρ)*Sx=9.981をパス、(1-ρ)*Sx^2=23.37を誤差分散として用いて計算してみます。図3.23のパス係数を出すためには、sem()のSオプションには共分散を用いる必要があるようです。(または、data=testとする。)
反対に図3.22の場合は、Sオプションに共分散行列(cov(test))や生データを与えると、 標準化解は図3.22と一致しません。なぜだかはわかりません。教えてください。
「共分散構造分析［Amos編］－構造方程式モデリング－」の65頁、図3.22の確認的因子分析です。
因子からのパスを1、誤差分散を0に制約する方法という前者の方法です。
test2.mdl <- specifyModel()
英語学力 -> 英語, NA, 1
学力 -> 国語, p5, NA
学力 -> 地理, p4, NA
学力 -> 数学, p3, NA
学力 -> 化学, p2, NA
学力 -> 物理, p1, NA
学力 <-> 学力, NA, 1
英語学力 <-> 英語学力, NA, 1
英語学力 <-> 学力, c1, NA
国語 <-> 国語, e1, NA
地理 <-> 地理, e2, NA
英語 <-> 英語, NA, 0
数学 <-> 数学, e4, NA
化学 <-> 化学, e5, NA
物理 <-> 物理, e6, NA
summary( test2.sem <- sem( test2.mdl, S = cor( test ), N = nrow( test )))
英語学力 -> 英語, NA, 1
学力 -> 国語, p5, NA
学力 -> 地理, p4, NA
学力 -> 数学, p3, NA
学力 -> 化学, p2, NA
学力 -> 物理, p1, NA
学力 <-> 学力, NA, 1
英語学力 <-> 英語学力, NA, 1
英語学力 <-> 学力, c1, NA
国語 <-> 国語, e1, NA
地理 <-> 地理, e2, NA
英語 <-> 英語, NA, 0
数学 <-> 数学, e4, NA
化学 <-> 化学, e5, NA
物理 <-> 物理, e6, NA
summary( test2.sem <- sem( test2.mdl, S = cor( test ), N = nrow( test )))
チラと適合度指標を表示すると
Model Chisquare = 26.604 Df = 10 Pr(>Chisq) = 0.003
GFI = 0.912
AGFI = 0.816
RMSEA = 0.130 90% CI: (0.071, 0.190)
CFI = 0.681
SEMではAGFIを0.9以上とか、RMSEAを0.05以下という基準が受け入れられているようですが、この因子分析は、少し当てはまりが良くないようにみえます。GFI = 0.912
AGFI = 0.816
RMSEA = 0.130 90% CI: (0.071, 0.190)
CFI = 0.681
以前あるセミナーで、因子分析も適合度を確認するので、学生さんが論文を書けなくて大変なんだということを聞きました。実際、アンケートで集めたデータで因子分析してみると、まずいなって思うことがあります。ていうか、ほとんどいつも悲惨。コストかけて、質の悪いデータを集めたと告白すべきでしょうか。
さて、信頼度係数ρ=0.81ということにして、sqrt(ρ)*Sx=9.981をパス、(1-ρ)*Sx^2=23.37を誤差分散として用いて計算してみます。図3.23のパス係数を出すためには、sem()のSオプションには共分散を用いる必要があるようです。(または、data=testとする。)
反対に図3.22の場合は、Sオプションに共分散行列(cov(test))や生データを与えると、 標準化解は図3.22と一致しません。なぜだかはわかりません。教えてください。
test3.mdl <- specifyModel()
英語学力 -> 英語, NA, 9.981
学力 -> 国語, p5, NA
学力 -> 地理, p4, NA
学力 -> 数学, p3, NA
学力 -> 化学, p2, NA
学力 -> 物理, p1, NA
学力 <-> 学力, NA, 1
英語学力 <-> 英語学力, NA, 1
英語学力 <-> 学力, c1, NA
国語 <-> 国語, e1, NA
地理 <-> 地理, e2, NA
英語 <-> 英語, NA, 23.37
数学 <-> 数学, e4, NA
化学 <-> 化学, e5, NA
物理 <-> 物理, e6, NA
summary( test2.sem <- sem( test2.mdl, S = cov( test ), N = nrow( test )))
stdCoef( test3.sem )
英語学力 -> 英語, NA, 9.981
学力 -> 国語, p5, NA
学力 -> 地理, p4, NA
学力 -> 数学, p3, NA
学力 -> 化学, p2, NA
学力 -> 物理, p1, NA
学力 <-> 学力, NA, 1
英語学力 <-> 英語学力, NA, 1
英語学力 <-> 学力, c1, NA
国語 <-> 国語, e1, NA
地理 <-> 地理, e2, NA
英語 <-> 英語, NA, 23.37
数学 <-> 数学, e4, NA
化学 <-> 化学, e5, NA
物理 <-> 物理, e6, NA
summary( test2.sem <- sem( test2.mdl, S = cov( test ), N = nrow( test )))
stdCoef( test3.sem )
2012年10月4日木曜日
共分散構造分析[R編]ってないの？[第３章7]その1
「共分散構造分析［Amos編］－構造方程式モデリング－」の64頁、図3.21の確認的因子分析です。
今回はsem()関数には相関行列を使いました。
パス図も描いてみました。edge.labels="names"のdotファイルを示します。
test1.mdl <- specifyModel()
文系学力 -> 国語, p1, NA
文系学力 -> 地理, p2, NA
文系学力 -> 英語, p3, NA
理系学力 -> 数学, p4, NA
理系学力 -> 化学, p5, NA
理系学力 -> 物理, p6, NA
文系学力 <-> 文系学力, NA, 1
理系学力 <-> 理系学力, NA, 1
文系学力 <-> 理系学力, c1, NA
国語 <-> 国語, e1, NA
地理 <-> 地理, e2, NA
英語 <-> 英語, e3, NA
数学 <-> 数学, e4, NA
化学 <-> 化学, e5, NA
物理 <-> 物理, e6, NA
summary(test1.sem <- sem(test1.mdl, S=cor(test), N=nrow(test)))
pathDiagram(test1.sem,ignore.double=FALSE, rank.dir="TB",
same.rank="文系学力,理系学力",
standardized = TRUE,
node.font=c("IPAGothic",12), edge.labels="values")
文系学力 -> 国語, p1, NA
文系学力 -> 地理, p2, NA
文系学力 -> 英語, p3, NA
理系学力 -> 数学, p4, NA
理系学力 -> 化学, p5, NA
理系学力 -> 物理, p6, NA
文系学力 <-> 文系学力, NA, 1
理系学力 <-> 理系学力, NA, 1
文系学力 <-> 理系学力, c1, NA
国語 <-> 国語, e1, NA
地理 <-> 地理, e2, NA
英語 <-> 英語, e3, NA
数学 <-> 数学, e4, NA
化学 <-> 化学, e5, NA
物理 <-> 物理, e6, NA
summary(test1.sem <- sem(test1.mdl, S=cor(test), N=nrow(test)))
pathDiagram(test1.sem,ignore.double=FALSE, rank.dir="TB",
same.rank="文系学力,理系学力",
standardized = TRUE,
node.font=c("IPAGothic",12), edge.labels="values")
今回はsem()関数には相関行列を使いました。
パス図も描いてみました。edge.labels="names"のdotファイルを示します。
digraph "test1.sem" {
rankdir=TB;
size="8,8";
node [fontname="IPAGothic" fontsize=12 shape=box];
edge [fontname="Helvetica" fontsize=10];
center=1;
{rank=same "文系学力" "理系学力"}
"文系学力" [shape=ellipse]
"理系学力" [shape=ellipse]
"文系学力" -> "国語" [label="p1"];
"文系学力" -> "地理" [label="p2"];
"文系学力" -> "英語" [label="p3"];
"理系学力" -> "数学" [label="p4"];
"理系学力" -> "化学" [label="p5"];
"理系学力" -> "物理" [label="p6"];
"文系学力" -> "文系学力" [label="" dir=both];
"理系学力" -> "理系学力" [label="" dir=both];
"文系学力" -> "理系学力" [label="c1" dir=both];
"国語" -> "国語" [label="e1" dir=both];
"地理" -> "地理" [label="e2" dir=both];
"英語" -> "英語" [label="e3" dir=both];
"数学" -> "数学" [label="e4" dir=both];
"化学" -> "化学" [label="e5" dir=both];
"物理" -> "物理" [label="e6" dir=both];
}
rankdir=TB;
size="8,8";
node [fontname="IPAGothic" fontsize=12 shape=box];
edge [fontname="Helvetica" fontsize=10];
center=1;
{rank=same "文系学力" "理系学力"}
"文系学力" [shape=ellipse]
"理系学力" [shape=ellipse]
"文系学力" -> "国語" [label="p1"];
"文系学力" -> "地理" [label="p2"];
"文系学力" -> "英語" [label="p3"];
"理系学力" -> "数学" [label="p4"];
"理系学力" -> "化学" [label="p5"];
"理系学力" -> "物理" [label="p6"];
"文系学力" -> "文系学力" [label="" dir=both];
"理系学力" -> "理系学力" [label="" dir=both];
"文系学力" -> "理系学力" [label="c1" dir=both];
"国語" -> "国語" [label="e1" dir=both];
"地理" -> "地理" [label="e2" dir=both];
"英語" -> "英語" [label="e3" dir=both];
"数学" -> "数学" [label="e4" dir=both];
"化学" -> "化学" [label="e5" dir=both];
"物理" -> "物理" [label="e6" dir=both];
}
共分散構造分析[R編]ってないの？[第３章6]その3 PLSモデル
図3.19は端折って、「共分散構造分析［Amos編］－構造方程式モデリング－」の63頁にあるPLSモデル図3.20を計算してみます。
以下のようにモデルを書いたファイルをsemi.3.20.txtとして保存しました。
実行します。
講師の質の分散を1に制約すると、充実感の分散が負になってしまいます。
淡々とやりましたが、少し気になることがあります。名城大学人間学部人間学科神谷研究室ホームページの共分散構造分析でも、「原因をさぐる統計学」(講談社ブルーバックス)に記載されているデータを基にRでPLSモデルを計算しておられます。ここでは、「洋食傾向」という指標化された変数があって、その分散の制約は1になっています。上で書いたこととは違うじゃん。実際、制約を1にすると本に記載されている係数が求められます。「洋食傾向」の制約を0にすると識別できません。
ついでに試してみると、制約を1にした場合でも、相関係数ではなく元データをそのものをデータフレームでsem()関数に与えると識別できません。 どういうことなんだ！わけわかんない。(@o@)そう、僕は嘘つきというより、ばかなんですよね。
以下のようにモデルを書いたファイルをsemi.3.20.txtとして保存しました。
テキスト -> 講師の質, p1, NA
プレゼン -> 講師の質, p2, NA
ペース -> 講師の質, p3, NA
講師対処 -> 講師の質, p4, NA
充実感 -> 満足度, p5, NA
充実感 -> 理解度, p6, NA
充実感 -> 目的一致, NA, 1 #1に制約
講師の質 -> 充実感, NA, 1 #1に制約
テキスト <-> テキスト, v1, NA
プレゼン <-> プレゼン, v2, NA
ペース <-> ペース, v3, NA
講師対処 <-> 講師対処, v4, NA
満足度 <-> 満足度, e1, NA
理解度 <-> 理解度, e2, NA
目的一致 <-> 目的一致, e3, NA
充実感 <-> 充実感, v5, NA
講師の質 <-> 講師の質, NA, 0 #0に制約
テキスト <-> プレゼン, c1, NA
テキスト <-> ペース, c2, NA
テキスト <-> 講師対処, c3, NA
プレゼン <-> ペース, c4, NA
プレゼン <-> 講師対処, c5, NA
ペース <-> 講師対処, c6, NA
</-></-></-></-></-></-></-></-></-></-></-></-></-></-></->
プレゼン -> 講師の質, p2, NA
ペース -> 講師の質, p3, NA
講師対処 -> 講師の質, p4, NA
充実感 -> 満足度, p5, NA
充実感 -> 理解度, p6, NA
充実感 -> 目的一致, NA, 1 #1に制約
講師の質 -> 充実感, NA, 1 #1に制約
テキスト <-> テキスト, v1, NA
プレゼン <-> プレゼン, v2, NA
ペース <-> ペース, v3, NA
講師対処 <-> 講師対処, v4, NA
満足度 <-> 満足度, e1, NA
理解度 <-> 理解度, e2, NA
目的一致 <-> 目的一致, e3, NA
充実感 <-> 充実感, v5, NA
講師の質 <-> 講師の質, NA, 0 #0に制約
テキスト <-> プレゼン, c1, NA
テキスト <-> ペース, c2, NA
テキスト <-> 講師対処, c3, NA
プレゼン <-> ペース, c4, NA
プレゼン <-> 講師対処, c5, NA
ペース <-> 講師対処, c6, NA
</-></-></-></-></-></-></-></-></-></-></-></-></-></-></->
実行します。
semi.mdl <- specifyModel("semi.3.20.txt")
(semi.sem <- sem(semi.mdl, data=semi.data, debug=TRUE))
summary(semi.sem)
stdCoef(semi.sem)
round(1-diag(semi.sem$P[-8,-8] / semi.sem$S),2)
round(cov2cor(semi.sem$P),2)
(semi.sem <- sem(semi.mdl, data=semi.data, debug=TRUE))
summary(semi.sem)
stdCoef(semi.sem)
round(1-diag(semi.sem$P[-8,-8] / semi.sem$S),2)
round(cov2cor(semi.sem$P),2)
講師の質の分散を1に制約すると、充実感の分散が負になってしまいます。
淡々とやりましたが、少し気になることがあります。名城大学人間学部人間学科神谷研究室ホームページの共分散構造分析でも、「原因をさぐる統計学」(講談社ブルーバックス)に記載されているデータを基にRでPLSモデルを計算しておられます。ここでは、「洋食傾向」という指標化された変数があって、その分散の制約は1になっています。上で書いたこととは違うじゃん。実際、制約を1にすると本に記載されている係数が求められます。「洋食傾向」の制約を0にすると識別できません。
ついでに試してみると、制約を1にした場合でも、相関係数ではなく元データをそのものをデータフレームでsem()関数に与えると識別できません。 どういうことなんだ！わけわかんない。(@o@)そう、僕は嘘つきというより、ばかなんですよね。
2012年9月26日水曜日
共分散構造分析[R編]ってないの？[第３章6]その2
「共分散構造分析［Amos編］－構造方程式モデリング－」の63頁図3.18をRで計算します。
セミナー評価という潜在因子は誤差を考慮せずに観測変数に重みをかけて足し上げた合成指標ということです。図3.17のセミナー評価は重回帰分析のように誤差を計算させましたが、図3.18では誤差(因子独自の分散)を0に制約するということになるはずです。
モデルは次のように書き、 semi.3.18.txtとして保存しました。
ファイルを読み込んで、実行します。
セミナー評価という潜在因子は誤差を考慮せずに観測変数に重みをかけて足し上げた合成指標ということです。図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
プレゼン -> セミナー評価, 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)
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)
2012年9月25日火曜日
共分散構造分析[R編]ってないの？[第３章6]その1
3.6 MIMICモデル
「共分散構造分析［Amos編］－構造方程式モデリング－」の62ページです。
ＭＩＭＩＣモデルの記述は長いので、次の別ファイル(semi.3.17.txt)を用意しました。セミナー評価から出ている矢印のうち、1つを制約しないと識別できません。 セミナー評価から目的一致のパスを1に制約しました。
これを次のように実行します。
最後の三行でそれぞれ
本に近い値が出ています。
次は、図3.18のモデルに取り組みます。
「共分散構造分析［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
プレゼン -> セミナー評価, 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)
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編） |
2012年9月17日月曜日
共分散構造分析[R編]ってないの？[第３章5]
3.5 共分散分析
共変量分析というと、以前、日常的に運動をしている人は最大酸素摂取量(VO2MAX)が多くなるか確認できないかと尋ねられて検討したと思います。この時の補助変数(共変量って言ってたような)は年齢でした。群を分けるのも「日常的に運動していますか」という設問への自己申告で心もとない感じだったのですが。
同じようにボディマス指数(BMI)も確認したいということもあったのですが、
さらに、年齢とBMIの関係が対象とした年齢(20-60)では、50歳程度でピークとなる非線形の関係があるようでした。人間って50歳過ぎると萎むんですね。私なんか、なんだか縮んじゃって勤め先でも居場所がなくなってきてるんですよ・・・。(; _ ;)
そんなボヤキはさておき、こうしてみました。
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)
切片を計算するために、
summary()の出力では、p2の標準誤差が本より少し小さく、信頼区間もすこし狭くなります。
修正指数を出力するmodIndices()関数には、しきい値の設定がないようです。出力された指数が無視してよいものなのかどうなのか判断に苦しみます。ただ、本で指定しているしきい値(4)よりはずっと小さいです。
次に、群 -> 英語の試験のパスを0に制約してみます。モデルは次のように書きました。
english.mdl <- specifyModel()
学力テスト -> 英語の試験, p1, NA
群 -> 英語の試験, NA, 0
英語の試験 <-> 英語の試験, e1, NA
学力テスト <-> 学力テスト, v1, NA
群 <-> 群, v2, NA
1 -> 学力テスト, i1, NA
改善指数が一番大きいのは、群 -> 英語の試験にパスを引くことになります。修正指数はAmosのそれに近い値が出ています。
共変量分析というと、以前、日常的に運動をしている人は最大酸素摂取量(VO2MAX)が多くなるか確認できないかと尋ねられて検討したと思います。この時の補助変数(共変量って言ってたような)は年齢でした。群を分けるのも「日常的に運動していますか」という設問への自己申告で心もとない感じだったのですが。
同じようにボディマス指数(BMI)も確認したいということもあったのですが、
- 20歳代は運動すると、筋肉がたくさんつくみたいで運動する人の方がBMIが大きい
- 30歳代半ば過ぎると、運動する人のほうが贅肉がつきにくいということなのかBMIが小さい
さらに、年齢とBMIの関係が対象とした年齢(20-60)では、50歳程度でピークとなる非線形の関係があるようでした。人間って50歳過ぎると萎むんですね。私なんか、なんだか縮んじゃって勤め先でも居場所がなくなってきてるんですよ・・・。(; _ ;)
そんなボヤキはさておき、こうしてみました。
english.mdl <- specifyModel()
学力テスト -> 英語の試験, p1, NA
群 -> 英語の試験, p2, NA
英語の試験 <-> 英語の試験, e1, NA
学力テスト <-> 学力テスト, v1, NA
群 <-> 群, v2, NA
1 -> 英語の試験, i0, NA
学力テスト -> 英語の試験, 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のそれに近い値が出ています。
登録:
投稿 (Atom)