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

「共分散構造分析［Amos編］－構造方程式モデリング－」の84頁の多母集団因子分析。
semパッケージでのモデルの記述方法は複数あります。今回は因子分析なので確認的因子分析に便利なcfa()関数を使った記述方法を使ってみたいと思います。
まずは、データを読み込んで、"学年"という変数をファクターに変換し、その次にモデルを記述します。

c04fact$学年 <- factor(c04fact$学年, labels = c("小学６年生","中学１年生")) #ファクターに変換
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 &lt- 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.mgmdl1 <- multigroupModel(c04fact.mdlA, c04fact.mdlB,  groups = levels(c04fact$学年))
summary(c04fact.gmsem1 <- sem(c04fact.mgmdl1, data = c04fact, group = "学年",
        formula = ~ 記号操作 + 計算 + 文章問題 + 図形))

測定不変モデルは配置不変モデルと比べると、RMSEAとAICが小さくなっています。
本と比較するとAICは記載どおりの値ですが、RMSEAはやや大きく習慣的基準値を若干上回り、信頼区間はかなり広くなっています。パス係数はほぼ一致しています。本ではパス係数を制約して分散を計算させているのでsem()の出力とそのまま比較できますが、少し大きいようです。

ここまでの多母集団分析でパス係数やGFI、AICはAmosとほぼ一致するのですが、RMSEAは大きくなりました。もしsemパッケージでは必ずそうなるのならば、RMSEAはモデルの採用時にはほとんど常に参照される適合度指標になっているように思うので、少し厄介だなと感じます。

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

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

「共分散構造分析［Amos編］－構造方程式モデリング－」の90頁からの平均共分散構造分析を試してみます。
まず、モデルの記述は次のようにします。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
#

平均構造を計算する場合、sem()関数では必ずrawオプションにTRUEを指定します。デフォルトではFALSEとなっていて分散共分散を計算します。このまま共分散行列を計算してもらうと逆行列が計算できませんということでエラーになります。平均構造を導入する場合は、ここをTRUEに指定して一次のモーメントを計算してもらいます。fixed.xオプションではInterceptを指定します。formulaオプションにはInterceptを指定しません。

summary(semi.sem <- sem(semi.mdl, data = semidata, raw = TRUE, fixed.x = "Intercept",
                        formula = ~ 満足度 + 理解度 + プレゼン + 講師対処 ))

これで、91頁のパス図に書かれている値が推定されていると思います。ちなみに、平均と切片を計算する場合、適合度指標は出力されません。Amosでも。

さて、パス図(左)を作成してみると大変見苦しくなってしまいました。ちょっと手を加えてみたのが右の図です。

共分散構造分析[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
#

この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 = "学校"))

適合度をみると、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 = "学校"))

配置不変モデルと比較すると、適合度指標が悪化しています。これらのグループの集団の異質性を示唆する結果ということでしょうか。

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

 残念ながら、71頁図3.31のモデルで座礁してしまいました。これは宿題として積み残して、「共分散構造分析［Amos編］－構造方程式モデリング－」の74頁から81頁あたりをやってみます。多母集団分析は、昨年(2012)にsemパッケージに追加された機能のようです。
　まずは、相関行列の確認から．．．。

cor(c04school[,-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

次に学校ごと、それぞれ計算します。他母集団分析の時は、集団を区別する変数を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)

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="学校"))

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] #差の検定統計量を計算して、必要な箇所を取り出す

図4.15「パラメータ間の差に対する検定統計量」に該当する表が表示されました。変数の順序は違いますが、小数点3桁表示でほぼ同じ値といっても差し支えないと思います。a1とa3のパスでは、A校とB校の間で5%有意な差が認められます。

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

「共分散構造分析［Amos編］－構造方程式モデリング－」の70頁、１因子モデルによる信頼性、図3.30を書いてみます。

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 )

適合度指標は表示されている少数第三位まで一致しています。パス係数はやや異なっています。
信頼性係数は、上記のように計算して丸めると、表示桁(少数第二位)まで一致します。

さて、 クロンバックのα信頼性係数を出力してくれるRの関数はいくつかあるようですが、psy パッケージの cronbach() という関数を使ってみます。

cronbach( kajyou )

本書に書かれているα係数と値が異なります。