RによるSEM
ＳＥＭの共分散構造、つまり観測変数の共分散に関するパラメータ構造 Σ(θ)の定式化には幾つかの方法が知られています。
Bentler and Weeks の ＥＱＳ モデル (EQuationS model) 、McArdle and McDonald の ＲＡＭモデル (Reticular Action Model)、Jöreskog-Keesling-Wiley の LISRELモデル (LInear Structural RELations model)、McDonald の COSANモデル (COvariance Structure ANalysis model)、 一般化COSANモデル (generalized COSAN model)など。ＲのパッケージではこのうちＲＡＭが利用出来ます。
さて、それではＳＥＭパッケージをダウンロードしましょう。但し、ダウンロードした zip ファイルを解凍ツールで解凍しないで下さい。どうやらコードが文字化けするようでＳＥＭを走らせることが出来ません。
メニューバーの Packages の Install package(s) from local zip files からインストールすると上手くいきます。
まず、観測変数の共分散行列、例でいうとR.DHPですが、これを入力します。
>R.DHP <- matrix(c(
1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
.6247, 1, 0, 0, 0, 0, 0, 0, 0, 0,
.3269, .3669, 1, 0, 0, 0, 0, 0, 0, 0,
.4216, .3275, .6404, 1, 0, 0, 0, 0, 0, 0,
.2137, .2742, .1124, .0839, 1, 0, 0, 0, 0, 0,
.4105, .4043, .2903, .2598, .1839, 1, 0, 0, 0, 0,
.3240, .4047, .3054, .2786, .0489, .2220, 1, 0, 0, 0,
.2930, .2407, .4105, .3607, .0186, .1861, .2707, 1, 0, 0,
.2995, .2863, .5191, .5007, .0782, .3355, .2302, .2950, 1, 0,
.0760, .0702, .2784, .1988, .1147, .1021, .0931, -.0438, .2087, 1
), ncol=10, byrow=TRUE)
次に、モデルをＲＡＭの特定化に沿って入力します。
model.dhp <- matrix(c(
'RParAsp -> RGenAsp', 'gam11', NA,
'RIQ -> RGenAsp', 'gam12', NA,
'RSES -> RGenAsp', 'gam13', NA,
'FSES -> RGenAsp', 'gam14', NA,
'RSES -> FGenAsp', 'gam23', NA,
'FSES -> FGenAsp', 'gam24', NA,
'FIQ -> FGenAsp', 'gam25', NA,
'FParAsp -> FGenAsp', 'gam26', NA,
'FGenAsp -> RGenAsp', 'beta12', NA,
'RGenAsp -> FGenAsp', 'beta21', NA,
'RGenAsp -> ROccAsp', NA, 1,
'RGenAsp -> REdAsp', 'lam21', NA,
'FGenAsp -> FOccAsp', NA, 1,
'FGenAsp -> FEdAsp', 'lam42', NA,
'RGenAsp <-> RGenAsp', 'ps11', NA,
'FGenAsp <-> FGenAsp', 'ps22', NA,
'RGenAsp <-> FGenAsp', 'ps12', NA,
'ROccAsp <-> ROccAsp', 'theta1', NA,
'REdAsp <-> REdAsp', 'theta2', NA,
'FOccAsp <-> FOccAsp', 'theta3', NA,
'FEdAsp <-> FEdAsp', 'theta4', NA),
ncol=3, byrow=TRUE)
'RParAsp -> RGenAsp', 'gam11', NA,を例にとると、RParAsp から RGenAspへの因果関係を仮定し、そのパラメータを gam11と名付け、NA を初期値として推計する。
ということになります。
このパッケージでは最尤法が用いられています。
次に、観測変数がどれなのかということを指定します。何も指定しない場合には共分散行列の行の名称が設定されます。
obs.vars.dhp <- c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp', 'RParAsp', 'RIQ','RSES', 'FSES', 'FIQ', 'FParAsp')
ここまでで準備はお終い。
いよいよ、モデルを推定します。
329 というのは共分散行列を推計するもとになった実際の観測数です。
John Fox氏の"Structural Equation Models / Appendix to An R and S-PLUS Companion to Applied Regression"によると、サンプルで用いられているDuncan, Haller, and Portesのデータセットに基づいています。
[参考]Duncan, O. D., A. O. Haller & A. Portes. 1968. "Peer Influences on Aspirations: A Reinterpretation." American Journal of Sociology 74:119 - 37.
> sem.dhp.1 <- sem(model.dhp, R.DHP, 329, obs.vars.dhp,fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp'))
> summary(sem.dhp.1)
とすると、推計結果が表示されます。
データに対するモデル Σ(θ) の適合度関数 f に対して、次の量
χ2 = ( N - 1 ) f
は、N が大きい時近似的に、m を自由母数の数 + 制約母数の数 - 制約数として
自由度 = (1/2)p(p+1) - m
という χ2 分布に従うことが知られています。
Ｒのアウトプットにある
Model Chisquare = 26.697 Df = 15 Pr(>Chisq) = 0.031302
というのは、
モデルが正しいとする帰無仮説をモデルは正しくないとする対立仮説に対してχ2 分布を用いて検定した結果を示しています。
この場合、
Pr(>Chisq) = 0.031302
ですから、Chisquare = 26.697 よりも大きな値となる確率は約３％ということになり帰無仮説は棄却される。モデルは正しくないということになる。
しかし、ここで注意しなければならないのはカイ自乗検定は観測変数Nに影響受けやすいという点。標本数が多くなるとモデルが棄却されやすくなることが知られています。
次に、Goodness-of-fit index(GFI)は、カイ自乗検定とは異なり観測変数Nに影響を受けにくいモデル評価の指標とされる。
通常は０から１までの値をとり、0.9以上あればかなり良いモデルということになります。
この例でいうと、
Goodness-of-fit index = 0.98439
なので非常に当て嵌まりが良いと言っても良いでしょう。
ここまで来たらもう少し。
> path.diagram(sem.dhp.1, min.rank='RIQ, RSES, RParAsp, FParAsp, FSES, FIQ',
max.rank='ROccAsp, REdAsp, FEdAsp, FOccAsp')
と入力するとパス図の「素」が出力されます。
digraph "sem.dhp.1" {
rankdir=LR;
size="8,8";
node [fontname="Helvetica" fontsize=14 shape=box];
edge [fontname="Helvetica" fontsize=10];
center=1;
{rank=min "RIQ" "RSES" "RParAsp" "FParAsp" "FSES" "FIQ"}
{rank=max "ROccAsp" "REdAsp" "FEdAsp" "FOccAsp"}
"RGenAsp" [shape=ellipse]
"FGenAsp" [shape=ellipse]
"RParAsp" -> "RGenAsp" [label="gam11"];
"RIQ" -> "RGenAsp" [label="gam12"];
"RSES" -> "RGenAsp" [label="gam13"];
"FSES" -> "RGenAsp" [label="gam14"];
"RSES" -> "FGenAsp" [label="gam23"];
"FSES" -> "FGenAsp" [label="gam24"];
"FIQ" -> "FGenAsp" [label="gam25"];
"FParAsp" -> "FGenAsp" [label="gam26"];
"FGenAsp" -> "RGenAsp" [label="beta12"];
"RGenAsp" -> "FGenAsp" [label="beta21"];
"RGenAsp" -> "ROccAsp" [label=""];
"RGenAsp" -> "REdAsp" [label="lam21"];
"FGenAsp" -> "FOccAsp" [label=""];
"FGenAsp" -> "FEdAsp" [label="lam42"];
}
これを、sem.dhp.1.dot として EUC 形式でファイルに保存します。
そして、graphvizのWindows 版をダウンロード。
dot ファイルを gif 形式の画像ファイルに変換します。
Bentler and Weeks の ＥＱＳ モデル (EQuationS model) 、McArdle and McDonald の ＲＡＭモデル (Reticular Action Model)、Jöreskog-Keesling-Wiley の LISRELモデル (LInear Structural RELations model)、McDonald の COSANモデル (COvariance Structure ANalysis model)、 一般化COSANモデル (generalized COSAN model)など。ＲのパッケージではこのうちＲＡＭが利用出来ます。
さて、それではＳＥＭパッケージをダウンロードしましょう。但し、ダウンロードした zip ファイルを解凍ツールで解凍しないで下さい。どうやらコードが文字化けするようでＳＥＭを走らせることが出来ません。
メニューバーの Packages の Install package(s) from local zip files からインストールすると上手くいきます。
まず、観測変数の共分散行列、例でいうとR.DHPですが、これを入力します。
>R.DHP <- matrix(c(
1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
.6247, 1, 0, 0, 0, 0, 0, 0, 0, 0,
.3269, .3669, 1, 0, 0, 0, 0, 0, 0, 0,
.4216, .3275, .6404, 1, 0, 0, 0, 0, 0, 0,
.2137, .2742, .1124, .0839, 1, 0, 0, 0, 0, 0,
.4105, .4043, .2903, .2598, .1839, 1, 0, 0, 0, 0,
.3240, .4047, .3054, .2786, .0489, .2220, 1, 0, 0, 0,
.2930, .2407, .4105, .3607, .0186, .1861, .2707, 1, 0, 0,
.2995, .2863, .5191, .5007, .0782, .3355, .2302, .2950, 1, 0,
.0760, .0702, .2784, .1988, .1147, .1021, .0931, -.0438, .2087, 1
), ncol=10, byrow=TRUE)
次に、モデルをＲＡＭの特定化に沿って入力します。
model.dhp <- matrix(c(
'RParAsp -> RGenAsp', 'gam11', NA,
'RIQ -> RGenAsp', 'gam12', NA,
'RSES -> RGenAsp', 'gam13', NA,
'FSES -> RGenAsp', 'gam14', NA,
'RSES -> FGenAsp', 'gam23', NA,
'FSES -> FGenAsp', 'gam24', NA,
'FIQ -> FGenAsp', 'gam25', NA,
'FParAsp -> FGenAsp', 'gam26', NA,
'FGenAsp -> RGenAsp', 'beta12', NA,
'RGenAsp -> FGenAsp', 'beta21', NA,
'RGenAsp -> ROccAsp', NA, 1,
'RGenAsp -> REdAsp', 'lam21', NA,
'FGenAsp -> FOccAsp', NA, 1,
'FGenAsp -> FEdAsp', 'lam42', NA,
'RGenAsp <-> RGenAsp', 'ps11', NA,
'FGenAsp <-> FGenAsp', 'ps22', NA,
'RGenAsp <-> FGenAsp', 'ps12', NA,
'ROccAsp <-> ROccAsp', 'theta1', NA,
'REdAsp <-> REdAsp', 'theta2', NA,
'FOccAsp <-> FOccAsp', 'theta3', NA,
'FEdAsp <-> FEdAsp', 'theta4', NA),
ncol=3, byrow=TRUE)
'RParAsp -> RGenAsp', 'gam11', NA,を例にとると、RParAsp から RGenAspへの因果関係を仮定し、そのパラメータを gam11と名付け、NA を初期値として推計する。
ということになります。
このパッケージでは最尤法が用いられています。
次に、観測変数がどれなのかということを指定します。何も指定しない場合には共分散行列の行の名称が設定されます。
obs.vars.dhp <- c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp', 'RParAsp', 'RIQ','RSES', 'FSES', 'FIQ', 'FParAsp')
ここまでで準備はお終い。
いよいよ、モデルを推定します。
329 というのは共分散行列を推計するもとになった実際の観測数です。
John Fox氏の"Structural Equation Models / Appendix to An R and S-PLUS Companion to Applied Regression"によると、サンプルで用いられているDuncan, Haller, and Portesのデータセットに基づいています。
[参考]Duncan, O. D., A. O. Haller & A. Portes. 1968. "Peer Influences on Aspirations: A Reinterpretation." American Journal of Sociology 74:119 - 37.
> sem.dhp.1 <- sem(model.dhp, R.DHP, 329, obs.vars.dhp,fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp'))
> summary(sem.dhp.1)
とすると、推計結果が表示されます。
データに対するモデル Σ(θ) の適合度関数 f に対して、次の量
χ2 = ( N - 1 ) f
は、N が大きい時近似的に、m を自由母数の数 + 制約母数の数 - 制約数として
自由度 = (1/2)p(p+1) - m
という χ2 分布に従うことが知られています。
Ｒのアウトプットにある
Model Chisquare = 26.697 Df = 15 Pr(>Chisq) = 0.031302
というのは、
モデルが正しいとする帰無仮説をモデルは正しくないとする対立仮説に対してχ2 分布を用いて検定した結果を示しています。
この場合、
Pr(>Chisq) = 0.031302
ですから、Chisquare = 26.697 よりも大きな値となる確率は約３％ということになり帰無仮説は棄却される。モデルは正しくないということになる。
しかし、ここで注意しなければならないのはカイ自乗検定は観測変数Nに影響受けやすいという点。標本数が多くなるとモデルが棄却されやすくなることが知られています。
次に、Goodness-of-fit index(GFI)は、カイ自乗検定とは異なり観測変数Nに影響を受けにくいモデル評価の指標とされる。
通常は０から１までの値をとり、0.9以上あればかなり良いモデルということになります。
この例でいうと、
Goodness-of-fit index = 0.98439
なので非常に当て嵌まりが良いと言っても良いでしょう。
ここまで来たらもう少し。
> path.diagram(sem.dhp.1, min.rank='RIQ, RSES, RParAsp, FParAsp, FSES, FIQ',
max.rank='ROccAsp, REdAsp, FEdAsp, FOccAsp')
と入力するとパス図の「素」が出力されます。
digraph "sem.dhp.1" {
rankdir=LR;
size="8,8";
node [fontname="Helvetica" fontsize=14 shape=box];
edge [fontname="Helvetica" fontsize=10];
center=1;
{rank=min "RIQ" "RSES" "RParAsp" "FParAsp" "FSES" "FIQ"}
{rank=max "ROccAsp" "REdAsp" "FEdAsp" "FOccAsp"}
"RGenAsp" [shape=ellipse]
"FGenAsp" [shape=ellipse]
"RParAsp" -> "RGenAsp" [label="gam11"];
"RIQ" -> "RGenAsp" [label="gam12"];
"RSES" -> "RGenAsp" [label="gam13"];
"FSES" -> "RGenAsp" [label="gam14"];
"RSES" -> "FGenAsp" [label="gam23"];
"FSES" -> "FGenAsp" [label="gam24"];
"FIQ" -> "FGenAsp" [label="gam25"];
"FParAsp" -> "FGenAsp" [label="gam26"];
"FGenAsp" -> "RGenAsp" [label="beta12"];
"RGenAsp" -> "FGenAsp" [label="beta21"];
"RGenAsp" -> "ROccAsp" [label=""];
"RGenAsp" -> "REdAsp" [label="lam21"];
"FGenAsp" -> "FOccAsp" [label=""];
"FGenAsp" -> "FEdAsp" [label="lam42"];
}
これを、sem.dhp.1.dot として EUC 形式でファイルに保存します。
そして、graphvizのWindows 版をダウンロード。
dot ファイルを gif 形式の画像ファイルに変換します。
<< Home