まずは、相関行列の確認から...。
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%有意な差が認められます。
0 件のコメント:
コメントを投稿