統計数理（2010）
第 58 巻 第 2 号 207–221
c 2010 統計数理研究所
［研究ノート］

2 次と 3 次の積率を用いたマルチレベルモデルの提案
1

2

尾崎 幸謙 ・中村 健太郎 ・室橋 弘人

3（受付 2010 年 5 月 17 日；改訂 8 月 25 日；採択 8 月 27 日）

要旨マルチレベルモデルは多段抽出データに対する適切な分析手法であり，社会調査や教育調
査を含む社会科学研究においてしばしば使用されている．一方，構造方程式モデリング（SEM;
も社会科学の分野で広く使用されている分析手法であるが，SEM
Structural Equation Modeling）
は多母集団解析を応用することでマルチレベルモデルを表現することが可能である．SEM では
これまで（平均と）分散・共分散の情報を利用したモデリングがなされていたが，近年になって歪度や尖度を利用した nnSEM（non-normal Structural Equation Modeling）と呼ばれる手法が開発された．この手法は，非正規分布に従う変数を分析可能なだけでなく，変数間の因果の方向性を 2 変数の横断データから統計的に判断可能であるという利点がある．本論文では，nnSEMの枠組みで 2 段抽出の場合のマルチレベルモデルを表現し，各階層（抽出単位レベル）において因果の方向性を判断可能な統計モデルの開発を行い，シミュレーション研究でモデルの性質を調べた．
キーワード： マルチレベルモデル，2 段抽出，SEM，nnSEM.
1.

はじめに

1.1 多段抽出データとマルチレベルモデル
社会調査や教育調査等の調査を行う目的は，母集団についての知見を得ることにあり
（土屋，
2009）
，そのためのデータ収集法として単純無作為抽出法がある．しかし，単純無作為抽出法は費用や手間の面から鑑みて現実的に困難であり，しばしば多段抽出法（豊田, 1998; 土屋, 2009）
によってデータが収集される．多段抽出法とは，標本を抽出する前に，まずその上位の 1 次抽
出単位を抽出し，1 次抽出単位から 2 次抽出単位を抽出し，さらに 2 次抽出単位から 3 次抽出
単位を抽出し…と抽出を繰り返し，最終的な抽出単位を標本としてデータを収集する方法である．例えば，日本人の国民性調査や，社会階層と社会移動（SSM）調査などの社会調査では市区町村を 1 次抽出単位とした 2 段抽出法によってデータが収集されている（実際には層化抽出法
を組み合わせた層化 2 段抽出法が用いられることも多い）．他にも，教育調査において，生徒を標本としてデータを抽出する場合，そのデータは，その生徒の所属する学校を 1 次抽出単位とみなした 2 段抽出データとして扱われる．
多段抽出によって収集されたデータは，標本よりも上位の抽出単位を無視して単純無作為
1
2
3

統計数理研究所：〒190–8562 東京都立川市緑町 10–3
埼玉学園大学 経営学部：〒333–0831 埼玉県川口市木曽呂 1510 番地
お茶の水女子大学 人間発達教育研究センター：〒112–8610 東京都文京区大塚 2–1–1

208

統計数理

第 58 巻

第 2 号 2010

抽出によって収集されたデータと同じように分析すると，母数の推定量にバイアスが生じる（Raudenbush and Bryk, 2002）．これは，例えば学校から生徒を抽出した場合，異なる学校に所属する生徒間の違いには学校間の違いが混在しているにも関わらず，後者の違いを無視していることから生じる．逆に，同じ学校に所属するという意味で同一の学校内の生徒は独立な標本ではないが，分析においてそのことが無視されているという問題もある．そのため，多段抽出データに対する適切な分析手法として，マルチレベルモデル（Goldstein, 2002）あるいは階層線形モデル（Raudenbush and Bryk, 2002）と呼ばれる手法が開発されている．
マルチレベルモデルには，MLwiN version 2.1（Rasbash et al., 2009）
や HLM version 6（Rauなどの専門のソフトウェアがある．あるいは SPSS やフリーの統計ソフト
denbush et al., 2004）
ウェア R でも分析を行うことが可能である．また，社会科学研究で欠かすことのできない統計手法となりつつある構造方程式モデリング（SEM; Bollen, 1989）
も多母集団解析を応用することでマルチレベルモデルが表現可能である（Muthén, 1994）．SEM によるアプローチは，データ全体の共分散行列を，各抽出単位ごとの共分散行列に分解し，分解された共分散行列に対して，各抽出単位ごとにモデルをあてはめるものである．このとき，各抽出単位を多母集団と捉えて分析を行う．現在の SEM のソフトウェアでは， Mplus version 6（Muthén and Muthén, 2010）
や EQS 6.1（Bentler, 2006）に実装されている．
1.2 高次積率を使った SEM
SEM では母数推定に最尤法が用いられることが多く，その場合，観測変数に対して多変量正
規分布の仮定が置かれる．しかし，観測変数は必ずしも正規分布に従っているとは限らないため，3 次以上の高次積率を用いることで非正規分布の変数に対応した分析手法が Bentler（1983）
によって提案された．その後，Mooijaart（1985）
では，3 次の積率までを使用して因子分析における回転の不定性問題を解決する方法が提案された． Mooijaart（1985）は，高次積率を用いれば，単に非正規分布の変数が分析可能になるだけでなく，非正規分布であることを利用して 2次の積率までを利用した場合には克服が困難であった問題に対する解決策を提案した研究と位置付けることができる．さらに，Shimizu and Kano（2008）
では，高次積率を利用することで変数間の因果の方向性を 2 変数の横断データから統計的に判断可能であることが示された．これまでは，道具的変数モデルやラグつき変数モデルなど，2 変数以外の変数を収集しなければ適合度の観点から因果の方向が推測できなかったことを考えると，Shimizu and Kano（2008）の方法は非常に大きな進歩であるといえる．しかし，多くの社会科学研究のように質問紙調査がデータ収集法となっている場合には，未観測の交絡変数の存在は因果の方向性の検討に悪影響がある．そこで，Kano and Shimizu（2003）
では，2 次・3 次・4 次の積率を使用することで，未観測の交絡変数の影響を考慮した上で 2 変数間の因果の方向性を検証する方法が提案された．
また，豊田（2007）では，2 変数 x, y 間に相関がある場合，1）x → y ，2）y → x のいずれの適合度が良いのかを判断する方法に加えて，3）擬似相関のモデルが分析可能であることも実例を挙げて示されている．本論文では Shimizu and Kano（2008）にならい，高次積率を用いた構造方程式モデリングを nnSEM（non-normal Structural Equation Modeling）と呼ぶ．
変数間の因果の方向性の検討は， ICA（Independent Component Analysis，Comon, 1994）の枠組みでも行われている． ICA は非正規性をモデルの識別に利用することを特徴としており，
線形逐次モデルにおける変数間の順序性同定の方法開発が近年行われている． Shimizu et al.
（2006）は，データ生成過程が線形であり，交絡変数が存在せず，残差変数に非正規を仮定する
ことで，3 変数以上の場合において因果の方向性に関して最適な逐次モデルが識別可能である
ことを示した．彼らのモデルは LiNGAM（Linear Non-Gaussian Acyclic Models）と呼ばれ，推定のためのアルゴリズムが開発されている．そのアルゴリズムでは，復元行列の行の順序の不2 次と 3 次の積率を用いたマルチレベルモデルの提案

209

定性を，対角行列に 0 が含まれないような復元行列に対する置換行列を定めることが重要な鍵となっている．また，Hoyer et al.（2008）
は，未観測の交絡変数が存在していたとしても，デー
タに対して同等に当てはまる他の逐次モデルと区別はできないものの，逐次モデルのパラメタ推定が可能であることを示した．さらに，Lacerda et al.（2008）
は，複数の同分布に従うモデルを区別できないものの，LiNGAM を非逐次の場合に対するモデルへ拡張した．また，Shimizu
は平均，回帰係数，因果の方向性に潜在クラスを仮定したモデルを提案
and Hyvärinen（2008）
は潜在変数間のモデルに拡張するなど，ICA の文脈で変数間の因果のし，Shimizu et al.（2009）
方向性に関する研究が盛んである．
nnSEM を用いた実データ分析は Shimizu and Kano（2008）や豊田（2007）
で行われている．それらでは，縦断的に収集された 2 変数間の因果の方向性を上記の 1）x → y モデルと 2）y → x モ
デルの適合度比較によって判断し，時系列的に正しい方向性が示されるか否かを調べ，因果方向探索の方法としての nnSEM の妥当性が経験的に検討されている．例えば，Shimizu andKano（2008）
では大学生を対象とした調査における「高校生時の犯罪行動の頻度を尋ねる項目の合計得点」と「昨年の犯罪行動の頻度を尋ねる項目の合計得点」を 2 変数とし，高校生時を独立変数とし，昨年を従属変数としたモデルの適合度が，因果の方向を逆転させたモデルよりも高いことが示されている．あるいは，豊田（2007）では Golton による親子の身長データに対して nnSEM を適用し，親の身長が子の身長を説明するモデルの適合度が，因果の方向を逆転させたモデルよりも高いことを示した．また，豊田（2007）では自動車のエンジンの「排気量」と「最大出力」に対する適用例も示されており，「排気量」の大きな車ほど「最大出力」も大きくなるという物理的に妥当な結果が導かれている．あるいは ICA の枠組みにおける Shimizu
et al.（2006）
は，定常性条件が満たされるような自己回帰過程が成立し，ノイズが非正規の成り立つような時系列データに対しては，因果の方向性が正しく検出され，そうでないデータに対しては誤って検出されるという予想のもと 22 の時系列データセットに対して LiNGAM を適用し，予想通りの結果を得た．
しかし，Shimizu and Kano（2008）
による nnSEM は，データが多段抽出によって収集されていることは想定しておらず，したがって，多段抽出用の方法ではない．そこで，本論文では，
2 次と 3 次の積率を使った nnSEM の枠組みで 2 段抽出データを分析する方法を，2 次抽出単位の変数が 2 つで，1 次抽出単位の変数がその 2 変数の平均となっている場合において提案する．この方法は，上に挙げた 1）x → y ，2）y → x，3）
擬似相関の 3 つのモデルを各抽出単位レベ
ルで分析可能とするものである．例えば，勉強時間と成績の 2 変数を，学校を 1 次抽出単位，
生徒を 2 次抽出単位として収集した場合を考えてみる．勉強時間と成績の関係は，勉強時間が長いほど成績が良いと考えるのが妥当であるが，これは生徒レベルの現象であり，学校レベル
では，成績平均が高い学校は勉強時間の学校平均が長い，となっているかもしれない．このような分析を可能とする方法を提案し，nnSEM の利用範囲を広げることが本論文の目的である．
本論文ではこれ以降，2 次の積率を用いた場合の SEM の 2 段抽出モデルを解説し，3 次の積率を用いた場合の積率構造の期待値を各抽出単位で求め，その後シミュレーション研究を行い，モデルの性質を検討し，最後に今後の研究に向けての展望を述べる．
2.

SEM の 2 段抽出モデル

を参考に，2 次の積率を用いた場合の SEM による 2 段抽出モデルを記述
ここでは，豊田（2000）
する．なお，SEM の 2 段抽出モデルは，任意の数の観測変数 p で記述する．xcj を c(= 1, . . . , C)
番目の 1 次抽出単位の j(= 1, . . . , Nc ) 番目の標本の観測変数ベクトル（サイズは p × 1）とし，以

統計数理

210

第 58 巻下で表わされるとする．
(2.1)

xcj = x∗c + v ∗cj =第 2 号 2010

xc
yc





0
+
v cj




全標本数は N (= C（2.1）式の x∗c は 1 次抽出単位の特徴を表わす変数であり，c=1 Nc ) である．
一方 v ∗cj は 1 次抽出単位 c 内の標本
（2 次抽出単位）の特徴を表わす変数である． 1 次抽出単位 c
内の期待値は E[v ∗cj ] = 0 とする．したがって，µ = E[xcj ] = E[x∗c ] + 0(2.2)となる．µ は一般平均である．
（2.1）式最右辺第 1 項は，x∗c が xc（サイズ pb × 1）と y c（サイズ pw × 1）から構成されていることを表わしている．ここで，xc は 1 次抽出単位に関わる変数であり，例えば学校が 1 次抽出
単位の場合には，「生徒と教員の比率」や「授業時間数」など学校を直接記述する変数が含まれる．一方，yc は（2.2）式の仮定から，1 次抽出単位 c の母平均ベクトルであり，「家庭での勉
強時間の学校 c 内平均」など，学校 c に所属する生徒に関わる変数の平均が含まれる．なお，
pb + pw = p である．
v cj は標本（2 次抽出単位）に関わる変数であり，学校 c に所属する生徒 j の「勉強時間」や「成績」などの変数の c 内偏差が含まれる．（2.1）式最右辺第 2 項は，v ∗cj がサイズ pb × 1 のゼロベクトルと，v cj（サイズ pw × 1）から構成されていることを表わしている．サイズ pb × 1 のゼロベクトルは， 1 次抽出単位 c 内における母平均からの偏差を表わす変数には， 1 次抽出単位の特徴を表わす情報が（当然）含まれないことを表わしている．
標本共分散行列は，1 次抽出単位間と 1 次抽出単位内（標本間）でそれぞれ求められる．学校から生徒を抽出する例でいえば，前者は学校間，後者は学校内（生徒間）の標本共分散行列を指す．1 次抽出単位間の標本共分散行列は，サイズ p × p の
(2.3)

Sb =

C
1 
Nc (x̄c − x̄)(x̄c − x̄)C − 1 c=1と記述される．ここで，
x̄ =

1 
Nc x̄c
N c=1

(2.4)

x̄c =

1 
xcj = x∗c + v̄ ∗c
Nc j=1(2.5)

v̄ ∗c =1  ∗v
Nc j=1 cj

C

Nc

Nc

である．
一方，1 次抽出単位内の標本共分散行列は，サイズ pw × pw の(2.6)

と記述される．

Sw =

C Nc
1 
(xcj − x̄c )(xcj − x̄c )
N − C c=1 j=1

2 次と 3 次の積率を用いたマルチレベルモデルの提案

211

母数で構造化された共分散行列は，1 次抽出単位間については Σb (θ)（サイズ p × p）とする．
1 次抽出単位内については，サイズを p × p として，00Σw (θ) =0 Σv (θ)とする．Σv (θ) はサイズ pw × pw である．Σw (θ) は 1 次抽出単位内の共分散行列であるから，
2 × 2 要素以外はゼロ行列である．θ は 2 次の積率のみを考えた場合の母数ベクトルである．
ここで，（2.3）式で表わされる 1 次抽出単位間の標本共分散行列 Sb の期待値は，母数で構造化された 1 次抽出単位間の共分散行列 Σb (θ) とはならず，
E[Sb ] = Σw (θ) + ωΣb (θ)
N 2 −ΣCN2

c=1 cとなることが知られている．ここで，ω = N(C−1)となる．一方，（2.6）
式で表わされる 1 次
抽出単位内の標本共分散行列 Sw の期待値は，

E[Sw ] = Σv (θ)となることが知られている．
以上から，データが 2 段抽出によって収集されており，2 次の積率を用いた場合の適合度関
数 f は以下となる（Muthén, 1994）．
f =(N − C)(sw − σ(θw )) Ww−1 (sw − σ(θw ))+ C(sb − (σ(θw ) + ωσ(θb ))) Wb−1 (sb − (σ(θ w ) + ωσ(θb )))b は 1 次抽出単位，w は 2 次抽出単位を表わす添え字である．σ(θw ) = vec [Σw (θ)] と σ(θb ) =vec [Σb (θ)] はそれぞれ，母数で表現された 1 次抽出単位内と 1 次抽出単位間における 2 次の積率ベクトルである．また，sw = vec [Sw ]，sb = vec [Sb ] である．vec はベック操作を表わす．
W −1 は重みであり，W の違いによって最小 2 乗法，一般化最小 2 乗法，最尤推定法，ADF 法（Asymptotically Distribution Free method; Browne, 1982, 1984）などが表現される．
2 次に加え，3 次の積率を使用する場合には，1 次抽出単位内と 1 次抽出単位間の 3 次の標
本積率ベクトル sw3 ，sb3 の期待値 E[sw3 ]，E[sb3 ] が必要となる．次章では，データが 2 段抽出で収集されている場合の，3 次の積率構造の期待値が導出される．
3.

nnSEM の 2 段抽出モデル

1 次抽出単位内の 3 次の標本積率ベクトル sw3 と 1 次抽出単位間の 3 次の標本積率ベクトル
sb3 の期待値は以下となる．
C
1  (Nc − 1)(Nc − 2)σ 3 (θw3 )N − C c=1
Nc

(3.1)

E[sw3 ] =(3.2)

E[sb3 ] =ωw σ 3 (θ w3 ) + ωb σ 3 (θb3 )ここで，σ 3 (θ w3 ) と σ 3 (θb3 ) は，それぞれ，母数で表現された 1 次抽出単位内と 1 次抽出単位間における 3 次の積率ベクトルである．（3.2）式の ωw と ωb はそれぞれ，
 C

 1
1
3C
2
ωw =
−
(3.3)
+
C − 1 c=1 Nc
N
N


C
C
1
3  2
2  3
N −
N
(3.4)
ωb =
N+ 2
C−1
N c=1 c
N c=1 c

統計数理

212

第 58 巻

第 2 号 2010

である．以下において，（3.1）式から（3.4）式の具体的な導出を行う．
3.1 1 次抽出単位内の 3 次の標本積率の期待値
2 変数 x と y に関して，1 次抽出単位内の 3 次の標本積率を，sw3 = (sxw3 , sxw2 yw , sxwyw2 , syw3 )とする．まず，sxw3 の期待値を求める．
（2.4）式と（2.5）式を利用すると，

 C N
c

1
3
E[sxw3 ] =(xcj − x̄c )E
N −C
c=1 j=1

 C N
c

1
∗
(x∗c + vcj− x∗c − v̄c∗ )3=E
N −C
c=1 j=1
N

C
c

1 
∗
∗ 3
E
(vcj − v̄c )=N − C c=1
j=1となる．ここで，3 次のキュムラントの不偏推定量 κ3 が k 統計量 k3 の期待値として与えられ，nk3 = (n−1)(n−2)
m3 となる（Stuart et al., 1994, p. 422）こと（ただし， m3 は 3 次の標本積率であり，n は標本数である），および，κ3 は平均値周りの 3 次の母集団積率 µ3 と等しい
（Stuart etことを利用すると，上式は，al., 1994, p. 90）
E[sxw3 ] =C
1  (Nc − 1)(Nc − 2)
σxw3
N − C c=1
Ncとなり，（3.1）式に一致する．ここで，σxw3 は変数 x の 1 次抽出単位内の 3 次の母集団積率である．
次に，sxw2 yw の期待値を求める．まず，k21 = (n−1)(n−2)
m21 を示しておく．ここで，k21 はそのn2期待値をとると 2 次 ×1 次の 2 変量キュムラントの不偏推定量 κ21 を与える k 統計量であり，m21 =
n
n
2
2
2
2
1
1
j=1 (xj − x̄) (yj − ȳ) である．m21 を展開すると nj=1 (xj yj − xj ȳ − 2xj yj x̄ + 2xj x̄ȳ + yj x̄ −nnx̄2 ȳ) となる．ここで，sts = j=1 xtj yjs とおくと，m21 は n1 s21 − n22 s10 s11 n12 s20 s01 + n23 s210 s01
1と表わされる．一方，Stuart et al.（1994, p. 455）より，k21 = n(n−1)(n−2)(n2 s21 − 2ns10 s11 −2
n20 s01 + 2s10 s01 ) であり，以上から，
(3.5)

k21 =

(n − 1)(n − 2)
m21
n2

となる．
sxw2 yw の期待値は，

 C N

c

1
2
E
E[sxw2 yw ] =(xcj − x̄c ) (ycj − ȳc )N −C
c=1 j=1
N

C
c

1 
∗
∗ 2
∗
∗
=
E
(vcj − v̄c ) (wcj − w̄c )N − C c=1
j=1∗となる．wcj，w̄c∗ はそれぞれ，変数 y における 1 次抽出単位 c 内の標本の特徴を表わす変数とその c 内平均である．ここで，（3.5）
式および，κ21 は平均値周りの 2 次 ×1 次の母集団積率 µ21と等しい（Stuart et al., 1994, p. 107）ことを利用すると，上式は，

E[sxw2 yw ] =C
1  (Nc − 1)(Nc − 2)
σxw2 yw
N − C c=1Nc

2 次と 3 次の積率を用いたマルチレベルモデルの提案

213となり，やはり（3.1）
式に一致する．ここで，σxw2 yw は変数 x と y の 1 次抽出単位内の 2 次 ×1
次の母集団積率である．
同様にして，
E[syw3 ] =E[sxwyw2 ] =C
1  (Nc − 1)(Nc − 2) 3
σyw
N − C c=1NcC
1  (Nc − 1)(Nc − 2)
σxwyw2
N − C c=1
Ncとなる．したがって，（3.1）式が成立する．
3.2 1 次抽出単位間の 3 次の標本積率の期待値
2 変数 x̄ と ȳ に関して，1 次抽出単位間の 3 次の標本積率を，sb3 = (sxb3 , sxb2 yb , sxbyb2 , syb3 )とする．E[sb3 ] の計算過程で必要となる確率変数 x に関する期待値 E[(x̄ − µx )3 ] と確率変数 xと y に関する期待値 E[(x̄ − µx )2 (ȳ − µy )] は付録 A で求められている．ここで，µx と µy はそれぞれ x と y の母平均である．
sxb3 の期待値は，

 C

1
3
E[sxb3 ] =Nc (x̄c − x̄)
E
C −1
c=1
⎤
⎡
3
C
C
C



1
1
1
Nc (x∗c − µ) + v̄c∗ −Nc (x∗c − µ) −Nc v̄c∗ ⎦
=
E⎣
C −1
N 
N 
c=1
c =1

c =1

C
1

となる．µ は x∗c の母平均である．ここで，(x∗c − µ) = a1 ，v̄c∗ = a2 ， NC∗1c =1 Nc v̄c = a4 とおくと，上式は
N(3.6)

E[sxb3 ] =∗c =1 Nc (xc − µ) = a3 ，

C


1 
Nc E a31 + a32 − a33 − a34 − 3a21 a3 + 3a1 a23 − 3a22 a4 + 3a2 a24
C − 1 c=133となる．これらの期待値は，E[a31 ] = σxb，E[a32 ] = N12 σxw，E[a33 ] = N13
Nc2
N2

c

C

3 3
3
c=1 Nc σxb ，E[a4 ] =

3
3
3
1
1
σ 3 ，E[a21 a3 ] = NNc σxb，E[a1 a23 ] =σxb，E[a22 a4 ] = NN
σ 3 ，E[a2 a24 ] = N12 σxwとなる．こN 2 xwc xw
3
3
2
2
こで，E[a2 ]，E[a4 ]，E[a2 a4 ]，E[a2 a4 ] の計算には付録（A.2）
式が利用されている．これらを（3.6）

式に代入して整理すると，（3.3）
式で示される ωw と（3.4）
式で示される ωb が求まり，（3.2）
式が成立する．これは E[syb3 ] についても同様である．
付録（A.3）
式を利用することで，同様にして，E[sxb2 yb ] および E[sxbyb2 ] を計算すると，（3.3）
式で示される ωw と（3.4）式で示される ωb が求まり，（3.2）式がやはり成立する．（3.3）式から，1 次抽出単位間の 3 次の標本積率の期待値は，母数で構造化された 1 次抽出単位内の 3 次の積率の影響をほぼ受けないことが分かる．
3.3 適合度関数
sxw3 の計算に含まれる (xcj − x̄c )3 は 1 次抽出単位 c における c 内平均からの個人得点の偏差の 3 乗であり，sxb3 の計算に含まれる (x̄c − x̄)3 は全平均からの c 内平均の偏差の 3 乗である．個人得点の偏差と c 内平均の偏差は独立であるから（学校内での相対的な学力の高低と，その学校の学力レベルの高低には関係がない），1 次抽出単位内のモデルと単位間のモデルを多母集団として捉えることが可能である．
以上から，データが 2 段抽出によって収集されており，2 次と 3 次の積率を用いた場合の適

統計数理

214

第 58 巻

第 2 号 2010

合度関数 f は以下となる．
−1
f =(N − C)(sw23 − σ(θw23 )) Ww23
(sw23 − σ(θw23 ))−1(sb23 − (σ(θw23 ) + ωσ(θb23 )))+ C(sb23 − (σ(θw23 ) + ωσ(θ b23 ))) Wb23

下付き添え字 23 は標本積率ベクトルや母集団積率ベクトルに 2 次と 3 次の積率が含まれていることを表す．
4.

シミュレーション研究

1 次抽出単位内における 2 変数 x，y の因果の方向と，1 次抽出単位間における 2 変数 x̄，ȳの因果の方向のこのモデルによる正判別確率を調べるためにシミュレーション研究を行った．
分析に使用したデータは 2 変数 x と y であり，1 次抽出単位内で x → y ，y → x という 2 通りの単回帰モデルおよび，1 次抽出単位間で x̄ → ȳ ，ȳ → x̄ という 2 通りの単回帰モデルを考えた．x̄ および ȳ は x と y の 1 次抽出単位内平均である．したがって，合計 4（= 2 × 2）通りのモデルを真のモデルとした．ただし，x → y かつ x̄ → ȳ と y → x かつ ȳ → x̄ は 1 次抽出単位内と単位間で方向が同じという意味で， y → x かつ x̄ → ȳ と x → y かつ ȳ → x̄ は 1 次抽出単位内と単位間で方向が逆という意味で同じであるから，4 通りのうち，x → y かつ x̄ → ȳ と，y → x
かつ x̄ → ȳ の 2 通りのみを使用した．2 通りのうちのいずれかを真のモデルとして発生させたシミュレーションデータに対して 4 通りのモデルで分析を行い，適合度関数 f の値の比較によってそれぞれのモデルが採択される回数を計算した．使用する変数は x，y ，x̄，ȳ であるから，pw = 2，pb = 0 である．
平均構造を含めないモデルとしたため，推定すべき母数は，1 次抽出単位内と単位間におけ232る回帰係数 αw ，αb ，独立変数の分散 σw，σb2 ，独立変数の 3 次の積率 σw，σb3 ，誤差の分散 σew，
233σeb ，誤差の歪度 σew ，σeb の合計 10 個である．したがって，1 次抽出単位内における 2 次と 3次の積率構造は以下のようになる．1 次抽出単位間については，添え字の w を b に変更すればよい．
2
E[Sxw2 ] = σw
2
E[Sxwyw ] = αw σw
2
2
E[Syw2 ] = α2w σw+ σew
3E[Sxw3 ] = σw
3
E[Sxw2 yw ] = αw σw
3
E[Sxwyw2 ] = α2w σw
3
3
E[Syw3 ] = α3w σw+ σew

回帰係数は 1 次抽出単位内と単位間においてともに真値を 1 とした．独立変数と誤差は 1 次
抽出単位内・単位間ともに χ2 (df ) から発生させた． df = 8, 32 の 2 通りを真の構造としたため，独立変数の分散と誤差の分散については 16 と 64，独立変数の 3 次の積率と誤差の 3 次の積率については 64 と 256 が真値である．なお，自由度が 8 と 32 いずれの場合も，2 変数間の相関は 0.707 となるため，変数間の関係の強さは同じである． 1 次抽出単位内と単位間の独立変数と誤差の 4 変数をそれぞれ異なる分布から発生させることも可能であるが，ここでは同じ χ2分布から発生させた．従属変数は，発生させた独立変数と誤差を使って 1× 独立変数 + 誤差とした．また，1 次抽出単位の数 C を，25，50，100，200，400 の 5 通り考え，1 次抽出単位 c 内のオブザベーション数 Nc は 10–20 の一様分布および， 20–40 の一様分布の 2 通りから発生さ2 次と 3 次の積率を用いたマルチレベルモデルの提案

215

せた．
以上から，シミュレーション研究では，χ2 分布の自由度（2 通り），1 次抽出単位の数 C（5 通り），一様分布の範囲（2 通り）の 3 つをパラメタとして，2 通りの真のモデルからデータを発生させ，各 4 つのモデルが採択される回数をそれぞれ求めた．シミュレーションの組み合わせは40(= 2 × 5 × 2 × 2) 通りであり，それぞれにおいて 100 のデータセットを発生させた．なお，推定には最小 2 乗法を用い，プログラムは R（version 2.10.0）で記述した．最適化計算には R の
optim 関数に含まれる BFGS 法（準ニュートン法）を使用した．
4.1 結果
結果は表 1 から表 8 の通りである．1 次抽出単位内（x，y ）と 1 次抽出単位間（x̄ ，ȳ）で区別して結果の解釈を行う．1 次抽出単位内については，C = 25，Nc = 10 − 20，df = 32 という最も正判
別確率が低いケースであっても，その確率は 1 次抽出単位内・間で方向が同じ場合で 91％，逆
表 1. x → y ，x̄ → ȳ ，10–20， df = 8 の場合のシミュレーション結果．

表 2. x → y ，x̄ → ȳ ，20–40， df = 8 の場合のシミュレーション結果．

表 3.

x → y ，x̄ → ȳ ，10–20，df = 32 の場合のシミュレーション結果．

表 4.

x → y ，x̄ → ȳ ，20–40，df = 32 の場合のシミュレーション結果．

統計数理

216

第 58 巻

第 2 号 2010

表 5. y → x，x̄ → ȳ ，10–20， df = 8 の場合のシミュレーション結果．

表 6. y → x，x̄ → ȳ ，20–40， df = 8 の場合のシミュレーション結果．

表 7.

y → x，x̄ → ȳ ，10–20，df = 32 の場合のシミュレーション結果．

表 8.

y → x，x̄ → ȳ ，20–40，df = 32 の場合のシミュレーション結果．

の場合で 89％であった．したがって，1 次抽出単位内における因果の方向性の探索は，これら
の条件であっても高い確率で成功するといえる．これは，1 次抽出単位内の標本数は， C
c=1 Ncであり，C = 25，Nc = 10 − 20 の場合でおよそサンプルサイズが 375 という比較的大きな数になることが理由であると思われる．
一方，1 次抽出単位間については，およそ 90％の正判別確率を得るためには， C = 400 が必要となる．したがって，1 次抽出単位間における因果の方向性の探索のためには， 1 次抽出単位の数を比較的大きくする必要がある．これは，1 次抽出単位内の標本数が C そのものである
ことに理由がある．また，1 次抽出単位内・間で方向が同じ場合には，逆の場合よりも 1 次抽
出単位内・間ともに正判別確率が若干高いこともみてとれる．

2 次と 3 次の積率を用いたマルチレベルモデルの提案

5.

217

展望および今後の課題

シミュレーション研究では 2 変数 x と y について，因果の方向性の検討を行ったが，nnSEM
では擬似相関についても分析可能である．また，導出した 3 次の積率の期待値は 2 変数 x と yに関するものであるが，変数 z 等が加わった場合には E[sxwywzw ] および E[sxbybzb ] を求めれば，3 変数以上の分析を行うことも可能である． x と z ，y と z 間の 3 次の積率の期待値は x とy の間の 3 次の積率の期待値と同じ構造となる．
マルチレベルモデルでは，ランダム切片モデルやランダム傾きモデル（Raudenbush and Bryk,
2002; Kreft and de Leeuw, 1998）が分析に用いられることが多い．学校から生徒を抽出し，生徒を標本として勉強時間（独立変数）と成績（従属変数）の関係を分析する場合を例に挙げると，
ランダム切片とは勉強時間の学校内平均であり，これは学校ごとに異なる．また，ランダム傾
きとは，学校内における成績に対する勉強時間の傾きであり，これも学校ごとに異なる．これ
らのランダム係数を学校レベルにおいて変数として扱うモデルが，ランダム切片モデルやラン
ダム傾きモデルである．Mplus ではランダム係数を潜在変数として扱うことで SEM の枠組み
でモデルを表現している． nnSEM を使った本論文ではランダム係数が扱われておらず，これは今後の課題としたい．
2 段抽出モデルはこれまで本論文で挙げた例以外にも， Patz et al.（2002）など，エッセーな
どの厳密に客観的な採点が困難なテストの採点を複数の評定者が行った場合に，評定者の厳しさや一貫性をモデルに組み込みつつ受験者を評価する方法としても使われている．あるいは，双生児データから表現型（観測変数）
に対する遺伝と環境の影響を推定する行動遺伝学の分野では，Guo and Wang（2002）は双生児データを，1 次抽出単位を家庭，2 次抽出単位を人とした 2
段抽出データとして分析している．以上のように，2 段抽出モデルはさまざまな分野に適用可能な方法であり，（3.1）式から（3.4）式を利用することで 3 次の積率を利用した分析をそれらの分野のデータに対して実行できる．
nnSEM を使った本方法は，上に挙げたように pw = 2，pb = 0 かつランダム係数のない単純なモデルを扱ったのみである．したがって，今後のモデル発展および，実データへ適用した応用研究が望まれる．
付 録A
ここでは，確率変数 X, Y に関して E[(X̄ − µX )3 ]，E[(X̄ − µX )2 (Ȳ − µY )] を求める．なお，
付録 A においては，x1 , x2 , . . . を確率変数 X の実現値とする．また，µX と µY はそれぞれ Xと Y の母平均である．
E[(X̄ − µX )3 ] については，まず，サイズ U の母集団から 2 個の標本 X1 , X2 を抽出したときを考えると，3 
X1 + X2− µX(A.1) E[(X̄ − µX )3 ] = E
2

3 
(X1 − µX ) + (X2 − µX )=E
2
1
= E (X1 − µX )3 + 3(X1 − µX )2 (X2 − µX ) + 3(X1 − µX )(X2 − µX )2
8

+(X2 − µX )3
1= E 2σX 3 + 3(X1 − µX )2 (X2 − µX ) + 3(X1 − µX )(X2 − µX )2
8

統計数理

218

第 58 巻

第 2 号 2010となる．ここで，σX 3 は確率変数 X の 3 次の母集団積率である．(X1 − µX )2 (X2 − µX ) と(X1 − µX )(X2 − µX )2 について考える．母集団サイズが U のとき，(Xi − µX )2 (Xj − µX ) の並
べ方は U (U − 1) 通りある．したがって，x1 , x2 , . . . , xU を確率変数 X の実現値とするとき，

E 3(X1 − µX )2 (X2 − µX )
1
E 3(x1 − µX )2 (x2 − µX ) + · · · + 3(xU − µX )(xU −1 − µX )2=
U (U − 1)1=
E {(x1 − µX ) + · · · + (xU − µ)}3 − (x1 − µX )3 − · · · − (xU − µX )3
U (U − 1)1
E −(x1 − µX )3 − (x2 − µX )3 − · · · − (xU − µX )3=
U (U − 1)−1
N σX 3=
U (U − 1)−1
σ 3=
U −1 Xとなる．E 3(X1 − µX )(X2 − µX )2 についても同様である．よって，これを（A.1）式に代入す
ると，


1
−2
E[(X̄ − µX )3 ] = E 2σX 3 +
σX 3
8
U −1
U − 2 σX3=
U −1 4となる．この結果を使って，サイズ U の母集団から n 個の標本を抽出した場合には，
E[(X̄ − µX )3 ]
3 
X1 + · · · + Xn
− µX
=E
n

3 
(X1 − µX ) + · · · + (Xn − µX )=E
n1
E[(X1 − µX )3 + · · · + (Xn − µX )3n3+ 3 (X1 − µX )2 (X2 − µX ) + · · · + (XN−1 − µX )2 (XN − µX )
+ 3 (X1 − µX )(X2 − µX )2 + · · · + (XU −1 − µX )(XU − µX )2


σX 3
1
= 3 E nσX 3 − n(n − 1)n
U −1
U − n σX 3=
U − 1 n2=となる．ここで，U がある程度大きい場合には(A.2)

E[(X̄ − µX )3 ] 

σX 3
n2

となる．同様の計算により，
(A.3)

となる．

E[(X̄ − µX )2 (Ȳ − µY )] σX 2 Yn22 次と 3 次の積率を用いたマルチレベルモデルの提案

219

参 考 文 献Bentler, P. M.（1983）. Some contributions to eﬃcient statistics in structural models: Speciﬁcationand estimation of moment structures, Psychometrika, 48, 493–517.
. EQS 6 Structural Equations Program Manual, Multivariate Software, Inc.,
Bentler, P. M.（2006）
Encino, California.
Bollen, K. A.（1989）. Structural Equations with Latent Variables, John Wiley & Sons, New York.
. Covariance structures, Topics in Applied Multivariate Analysis（ed. D. M.
Browne, M. W.（1982）
Hawkins）, 72–141, Cambridge University Press, Cambridge.
Browne, M. W.（1984）. Asymptotically distribution-free methods for the analysis of covariance structures, British Journal of Mathematical and Statistical Psychology, 9, 665–672.
Comon, P.（1994）. Independent component analysis, A new concept?, Signal Processing, 36, 287–
314.
. Kendall’s Library of Statistics 9: Multilevel Statistical Models, Edward Arnold,
Goldstein, H.（2002）
London.
. The mixed or multilevel model for behavior genetic analysis, Behavior
Guo, G. and Wang, J.（2002）
Genetics, 32, 37–49.
Hoyer, P. O., Shimizu, S., Kerminen, A. and Palviainen, M.（2008）. Estimation of causal eﬀects using
linear non-Gaussian causal models with hidden variables, International Journal of Approximate Reasoning, 49, 362–378.
. Causal inference using nonnormality, Proceedings of InternationalKano, Y. and Shimizu, S.（2003）
Symposium on Science of Modeling The 30th Anniversary of the Information Criterion（AIC ） , 261–270.
. Introducing Multilevel Modeling, Sage Publications Ltd, California.
Kreft, I. and de Leeuw, J.（1998）
Lacerda, G., Spirtes, P., Ramsey, J. and Hoyer, P. O.（2008）. Discovering cyclic causal models by
independent components analysis, Proceedings of 24th Conference on Uncertainty in Artificial, Helsinki.
Intelligence（UAI2008 ）Mooijaart, A.（1985）. Factor analysis for non-normal variables, Psychometrika, 50, 323–342.
. Multilevel covariance structure analysis, Sociological Methods and Research,
Muthén, B. O.（1994）
22, 376–398.
Muthén, L. K. and Muthén, B. O.（2010）. Mplus User’s Guide（Version 6 ）, Muthén & Muthén, Los
Angeles, California.
. The hierarchical rater model
Patz, R. J., Junker, B. W., Johson, M. S. and Mariano. L. T.（2002）
for rated test items and its application to large-scale educational assessment data, Journal of
Educational and Behavioral Statistics, 27, 341–384.
Rasbash, J., Charlton, C., Browne, W. J., Healy, M. and Cameron, B.（2009）. MLwiN Version 2.1,
Centre for Multilevel Modelling, University of Bristol, Bristol.
Raudenbush, S. W. and Bryk, A. S.（2002）. Hierarchical Linear Models: Applications and Data
Analysis Methods, Sage Publications, London.
. HLM 6 for Windows, Scientiﬁc Software
Raudenbush, S. W., Bryk, A. S and Congdon, R.（2004）
International, Inc., Lincolnwood, Illinois.
Shimizu, S. and Hyvärinen, A.（2008）. Discovery of linear non-Gaussian acyclic models in the presence of latent classes, Proceedings of 14th International Conference on Neural Information, 752–761.
Processing（ICONIP2007 ）. Use of non-normality in structural equation modeling: Application
Shimizu, S. and Kano, Y.（2008）
to direction of causation, Journal of Statistical Planning and Inference, 138, 3483–3491.
Shimizu, S., Hoyer, P. O., Hyvärinen, A. and Kerminen, A.（2006）. A linear non-Gaussian acyclic

220

統計数理

第 58 巻

第 2 号 2010

model for causal discovery, Journal of Machine Learning Research, 7, 2003–2030.
Shimizu, S., Hoyer, P. O. and Hyvärinen, A.（2009）. Estimation of linear non-Gaussian acyclic models
for latent factors, Neurocomputing, 72, 2024–2027.
. Kendall’s Advanced Theory of Statistics, Vol. 1 ;Stuart, A., Ord, J. K. and Kendall, M. G.（1994）
Distribution Theory, 6th ed., Edward Arnold, London.

豊田秀樹（1998）. 『調査法講義』，朝倉書店，東京．
豊田秀樹（2000）. 『共分散構造分析［応用編］』，朝倉書店，東京．
豊田秀樹（2007）. 『共分散構造分析［理論編］』，朝倉書店，東京．
土屋隆裕（2009）. 『概説標本調査』，朝倉書店，東京．

Proceedings of the Institute of Statistical Mathematics Vol. 58, No. 2, 207–221 (2010)

221

A Multilevel Model Using 2nd and 3rd Order Moments
Koken Ozaki1 , Kentaro Nakamura2 and Hiroto Murohashi3
1 The Institute of Statistical Mathematics
2 Faculty of Management, Saitama Gakuen University
3 Ochanomizu Research Center for Human Development and Education, Ochanomizu University

Multi-stage sampling is used to collect data, such as social survey data, educational
survey data, and so on. The multilevel model is appropriate for analyzing this kind of data,and is often used in sociology and pedagogy. Structural Equation Modeling (SEM), which
is widely used in social sciences, can accommodate a multilevel model using multi-group
analysis. In SEM, (means and) covariances are used as information. However, a method
named non-normal Structural Equation Modeling (nnSEM) has recently been developed.
The merits of nnSEM are that it can not only handle non-normally distributed variables,
but can also statistically judge the direction of causation between cross sectional variables.
In this paper, we develope a two-stage model whithin the framework of nnSEM.The model
can judge the direction of causations between cross sectional variables in both sampling
units. Simulation studies were performed to examine the characteristics of the model.

Key words: Multilevel model, two-stage sampling, SEM, nnSEM.