統計数理（2005）
第 53 巻 第 1 号 133–157
c 2005 統計数理研究所


特集「日本人の国民性調査 50 年」［研究ノート］

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題「日本人の国民性調査」データへの適用
原 宏和

†

（受付 2004 年 7 月 30 日；改訂 2005 年 2 月 4 日）

要

旨

データに正規性を仮定した下での共分散構造モデルの選択問題に関する赤池情報量規準（AIC ）
は，実際に観測値が従う分布である真の分布が正規分布でなく未知の非正規分布であった場合，真の分布の尖度に依存した定数バイアスを持つ．本論文の目的は，部分的に交差確認法（クロス・バリデーション）を適用することでバイアスを補正した，修正 AIC を選択規準に用いる
ことにより真の分布の非正規性の影響が小さくなるような選択法を提案することである．この修正 AIC のバイアスは AIC や T IC のバイアスよりも小さくなることを漸近展開式とシミュ
レーション実験により確かめる．また，実解析への適用例として国民性調査第 11 回全国調査のデータに「日本人の満足感」に関する因子分析モデルをあてはめ，その最適モデルを探索する．
キーワード： 赤池情報量規準，因子分析モデル，交差確認法，非正規性，モデル選択．
1.

はじめに

日本人の国民性調査データに関する解析における重要な目的の一つに，類似した複数の質問
項目の背後にあってそれらの質問項目に対する共通の回答傾向を生ぜしめる要因を探索するということがある．しかしながら，このような要因は直接観測されないことが多い．そのため，意味のある因果関係を，観測変数間の関係としてよりもむしろ，潜在変数（直接観測されない変数）間の関係として表現し，その潜在変数を確率変数とおく．そうすると，潜在変数間の因果関係の探索問題は，観測変数の共分散構造分析などに代表される，共分散行列の構造に関する推定問題に帰着され，目的とする因果関係の推論が可能となる．以上のことから，国民性調査データに関する解析において共分散構造モデルは非常に重要なモデルであると言え，前田（1995）においてもその有用性が論じられている．共分散構造モデルにおいて，仮定した因果関係にはそれに対応した共分散構造があり，それらのどの構造を選べば良いかということは非常に重要な問題となる．
観測値が本当に従う分布である真の分布はほとんどの場合未知であり，一般的な統計的解析
手法では観測値が正規分布に従っているという仮定の下で議論を進めているものが多い．他の分布に較べ最尤推定法が簡便であり，また従来のソフトウエアでは正規性の仮定の下で提案さ†

筑波大学大学院 システム情報工学研究科：〒305–8573 茨城県つくば市天王台 1–1–1；
yanagi@sk.tsukuba.ac.jp

134

統計数理 第 53 巻

第 1 号 2005

れた統計的手法が主流になっていることを考えると，解析の初期段階で取敢えず正規性を仮定することは，理論的には問題はあるかもしれないが，実用面から考えると解析の妥当な第一歩であると考える．正規性の仮定の下で最も一般的なモデル選択法は，Akaike（1973）により提案された赤池情報量規準（AIC ）が最も小さくなるモデルを最適なモデルとする手法である．AIC
は予測に関するカルバック＝ライブラーの距離に基づくリスクの推定量の一つであり，標本対数尤度関数の −2 倍にリスクに対するバイアス補正項である「独立パラメータ数の 2 倍」を足すことにより定義され，その式の簡便性から様々な分野で応用されている．ところが，AIC は
そのバイアス補正項のシンプルさから様々な状況で無視できないバイアスを持つことがある．
そのため，Sugiura（1978），Bozdogan（1987），Hurvich and Tsai（1989），Fujikoshi and Satoh（1997）等の多くの研究者により，そのようなバアイスを補正した情報量規準が多く提案されている．
しかしながら，これらのバイアスを補正した情報量規準は仮定したモデルの分布と真の分布が同じであるという条件の下で提案されたもので，もしそれらの分布間に違いが生じればまた異なったバイアスが発生する．国民性調査データなど社会学，心理学の分野において興味のあるデータはカテゴリカルデータであることが多いため，正規性の仮定を満たしていないものが多い．そのため，仮定した分布と真の分布の違いの影響を受けやすい手法でモデル選択を行え
ば，解析結果にバイアスが生じることも十分考えられる．
竹内（1976）によって提案された T IC は，候補のモデルの分布と真の分布が異なるという仮定の下で導出された情報量規準であり，真の分布が独立に同一な分布であれば AIC が持つような定数バイアスを補正しているということが知られている．この T IC では，最尤推定量の漸近正規性を用いてバイアス補正項を計算している．ところが多変量モデルでは，観測値の次元が標本の大きさに較べ大きくなれば推定量の正規分布への収束が遅くなり，結果的に T ICのバイアスが大きくなることがある．また T IC のバイアス補正項は真の分布の尖度に依存しているため，標本数が巨大でなければうまくバイアスを補正することができないということもFujikoshi et al.（2005）で報告されている．
情報量規準を用いる以外に仮説検定を行うことで共分散行列の構造を選択する方法もある．
観測値が正規分布に従っているという仮定の下で導出された検定統計量の帰無分布は，真の分
布が正規分布であった場合，漸近的に χ2 分布に従うことが知られており，χ2 近似により検定を行うことがほとんどである．しかしながら，その漸近帰無分布は，正規性は仮定したが真の分布は非正規であるといった状況の下では真の分布の尖度に基づく重みを持った χ2 分布の重み付き和の分布になる（Satorra and Bentler（1988），Yuan and Bentler（1999），Yanagihara
et al.（2005）等参照）．そのため，真の分布が正規分布ではないときに χ2 近似を用いて検定を行うと，パーセント点がずれ，解析者が望む有意水準と実際の有意水準が大きく異なってしまい，検定結果に誤りを与えてしまうことがある．特に，正規性の下での尤度比検定では，真の分布の尖度が大きくなるほど検定統計量の帰無仮説の下での漸近平均と漸近分散は χ2 分布の平均と分散よりも大きくなるので，χ2 分布を使った検定は非保守的な検定になる傾向がある
ことが Yanagihara et al.（2005）で述べられている．真の分布の非正規性の影響を小さくするような修正をほどこした検定統計量も Browne（1984），Satorra and Bentler（1988），Kano et al.
（1990），Kano（1992）等により提案されている．しかしながら，それらの多くは尖度の推定量に依存した補正法に基づくものであり，その補正がうまく機能するためには非常に大きな標本数が必要となる（Hu et al.（1992），Yuan and Bentler（1998）等参照）．以上のことから，仮説検定で共分散構造を選択する場合，大きな危険性が存在していると考えられる．
本論文では，Yanagihara（2004a）でのアイデアを基に部分的に交差確認法（クロス・バリデー
ション）を用いることで AIC のバイアスを補正し，その新しい情報量規準を使うことで非正非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

135

規性の影響を受けにくい選択法を提案する．真の分布がどのような分布であってもバイアスを補正しているこの情報量規準を使う選択法では，真の分布の非正規性の影響が小さくなっている．以下，第 2 章では，本論文で用いたモデルと今までに提案されている情報量規準 AIC とT IC を紹介し，新しいバイアス修正情報量規準である修正 AIC（Corrected AIC：AICc ）を提案する．第 3 章では，AICc は AIC や T IC よりも小さいバイアスを持ち，またリスクが最も小さいモデルを選ぶ頻度も改善していることをシミュレーションにより確かめる．第 4 章では，
第 11 回国民性調査データ（坂元 他（2004））に「日本人の満足感」の構造（前田（1995））に関する因子分析モデルをあてはめ，その最適モデルの探索に新しい情報量規準を用いた解析例を紹介する．
2.

共分散構造モデルとモデル選択規準

本章では，共分散構造モデルとその最適モデルの選択のための従来の情報量規準を紹介し，更に従来の情報量規準の問題点を改善した新しい情報量規準を提案する．
2.1 カルバック ライブラーの距離に基づくリスク
第 i 番目の個体 (i = 1, . . . , n) に対する p × 1 観測値ベクトルを y i = (yi1 , . . . , yip ) とする．
共分散構造モデルでは，観測値ベクトル y i に q × 1 未知パラメータベクトル θ = (θ1 , . . . , θq )
で表現できる共分散行列 Σ(θ) を持つ候補のモデル，
M : y i ∼ i.i.d. Np (µ, Σ(θ)),(2.1)(i = 1, . . . , n),をあてはめる．ただし µ = (µ1 , . . . , µp ) は p × 1 未知パラメータベクトルである．本論文では，候補のモデルに正規性を仮定するが，真の分布は正規分布に従っているかどうかわからないという状況を考える．すなわち，真のモデルをM∗ : y i = µ + Σ1/2
∗ εi ,(2.2)(i = 1, . . . , n),とする．ただし，それぞれの誤差ベクトル εi = (εi1 , . . . , εip ) は独立に同一な平均 E[εi ] = 0，
共分散行列 Cov[εi ] = I p をもつ未知の分布に従うとし，Σ∗ は未知の真の共分散行列であるとする．
今，ξ = (µ , θ ) を全ての未知パラメータを並べた (p + q) × 1 ベクトルとし，f (y i |ξ) を平
均 µ 共分散行列 Σ(θ) を持つ p 変量正規分布の密度関数，
(2.3)

1
2π

f (y i |ξ) =

p/2





1
|Σ(θ)|−1/2 exp − (yi − µ) Σ(θ)−1 (y i − µ) ,2とする．このとき，候補のモデル M（2.1）の下での対数尤度関数は，(2.4)

(ξ|Y ) =n


log{f (y i |ξ)}

i=1

n
1
np
log(2π) − log |Σ(θ)| −(y − µ) Σ(θ)−1 (yi − µ),222 i=1 in=−となる．ただし Y = (y 1 , . . . , y n ) である．これにより，候補のモデル M の下での µ の最尤
推定量は，
1
y ,
n i=1 i
n

ȳ =

136

統計数理 第 53 巻となり，θ の最尤推定量は，

第 1 号 2005θ̂ = arg min log |Σ(θ)| + tr(SΣ(θ)−1 ) ,

(2.5)

θ

となる．ただし S は標本共分散行列で1
(y − ȳ)(y i − ȳ) ,n i=1 inS=である．
複数の候補のモデルに対してモデルの適合度を比較するためには，何かしらのリスクを定義する必要がある．最終的には，それぞれのモデルにおけるリスクを比較し，リスクが最小になるようなモデルを最適なモデルとみなす．本論文では，そのようなリスクとして予測に基づく
カルバック ライブラーの距離を用いる．今，p × 1 ベクトル ui (i = 1, . . . , n) を互いに独立に，
また y i とは独立に同一の分布に従うような将来の観測値ベクトルとすると，候補のモデル Mの下でのリスクは，
(2.6)





RKL = np log(2π) + nE∗Y log |Σ(θ̂)|
+

n


∗ ∗



EU EY (ui − ȳ) Σ(θ̂)−1 (ui − ȳ)



i=1



= np log(2π) + nE∗Y log |Σ(θ̂)|









+nE∗Y tr(Σ∗ Σ(θ̂)−1 ) + nE∗Y (ȳ − µ) Σ(θ̂)−1 (ȳ − µ) ,と定義される．ただし，E∗ は真のモデル M∗（2.2）の下での期待値を表すとし，U = (u1 , . . . , un )
である．
2.2 AIC 規準
Akaike（1973）により提案された AIC はリスク（2.6）の推定量の 1 つである．その AIC は，
（2.4）式の対数尤度関数に ξ = ξ̂ を代入した標本対数尤度関数 (ξ̂|Y ) をξ̂ = (ȳ  , θ̂ ) とすると，
用いて，
(2.7)

AIC = −2(ξ̂|Y ) + 2(p + q)
= np log(2π) + n log |Σ(θ̂)| + tr(SΣ(θ̂)−1 ) + 2(p + q),と定義される．共分散構造モデルの下位モデルである因子分析モデルでの AIC に関しては，
Akaike（1987）により詳しく議論されているのでそちらを参照して欲しい．式（2.7）からわかる
ように，AIC は標本対数尤度の −2 倍，すなわち，−2(ξ̂|Y ) でリスク RKL（2.6）を推定したときに生じるバイアスを 2(p + q) で補正した推定量となっている．つまり AIC では，
RKL − E∗Y [−2(ξ̂|Y )] = nE∗Y [tr(Σ∗ Σ(θ̂)−1 )]+nE∗Y [(ȳ − µ) Σ(θ̂)−1 (ȳ − µ)] − nE∗Y [tr(SΣ(θ̂)−1 )]≈ 2(p + q),と近似していることになる．しかしながら上記の式からわかるように，バイアスは真の分布に依存しているにも拘らず，AIC のバイアス補正項はどのような分布に対しても一定である．その結果，AIC はリスク RKL に対し真の分布に依存するバイアスを持つことになる．
AIC のバイアスを計算すると以下のようになる．今，候補のモデル M（2.1）の下で標本数 n を−1/2が収束する値を θ0 ，つまり，θ 0 = limn→∞ θ̂ とし，Ω = Σ1/2
大きくしたときに θ̂（2.5）∗ Σ(θ 0 )

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

137と置く．このとき，AIC のリスク RKL に対するバイアス BAIC = RKL − E∗Y [AIC] は，(2.8)

BAIC = 2tr(Ω Ω) + 2tr Ω Ω ⊗ Ω Ω D(θ 0 |Σ∗ )+tr Γκ {Ω ⊗ Ω} D(θ0 |Σ∗ ) Ω ⊗ Ω

− 2(p + q) + O(n−1 ),となる（導出に関しては Appendix A.3 を参照）．ただし，D(·|·) は Appendix A.2 の（A.6）式で定義される p2 × p2 行列である．また，δab を a = b のとき 1 で a = b であれば 0 となるクロネッカーデルタ，e(p)を第 j 番目の成分が 1 で残りの成分が 0 となる p × 1 ベクトル， pa1 ,a2 ,...
jp
p
は a1 =1 a2 =1 · · · を示す表記とすると，Γκ は Yanagihara et al.（2005）で定義される，
(2.9)

p


Γκ =

(p) 

(p) 

κabcd e(p)a eb

⊗ e(p)
c ed,

a,b,c,dであるような真の分布の尖度に依存する p2 × p2 行列である．ここでの κabcd は εi の四次キュムラントであり，ε = (ε1 , . . . , εp ) を εi と同一の分布に従う確率変数とすると，(2.10)

κabcd = E∗Y [εa εb εc εd ] − δab δcd − δac δbd − δad δbc ,である．今 κ(1)を4
(2.11)

(1)

κ4 =

p


κaabb = E∗Y [(ε ε)2 ] − p(p + 2),a,bのように定義される多変量尖度（Mardia（1970），Isogai（1983）等）とすると，（2.9）と（2.11）は(1)
κ4 = tr(Γκ ) という関係式を満たす．このことから，AIC のバイアスは，真の分布の尖度に依存していること，更に候補のモデル M の下での共分散行列の期待値と真のモデルの共分散行列の構造の違いにも依存していることがわかる．ここで更に特別な形を考える．今，P A を行
列 A の列ベクトルで張る空間への射影行列 P A = A(A A)−1 A ，また ∆(·) を Appendix A.1の（A.1）式で定義される Σ(θ) の一次の偏微分に関する p2 × q 行列であるとする．このとき，
Σ∗ = Σ(θ0 ) であれば，Ω = I p ，D(θ0 |Σ∗ ) = P ∆(θ0 ) となることがわかる．このことにより，
Σ∗ = Σ(θ 0 ) であれば，式（2.8）での AIC のバイアスは以下のようになる．
BAIC = tr(Γκ P ∆(θ0 ) ) + O(n−1 ).

もし εi が正規分布に従えば，κabcd = 0 である．よって，式（2.7）で与えられる AIC は，真のモデルと候補のモデルが一致する，すなわち，εi が正規分布でかつ Σ(θ0 ) = Σ∗ であれば定数バイアスを持たないことがわかる．しかしながら，真のモデルと候補のモデルが異なれば定数バイアスを持つ．さらに真の分布の尖度が大きくなれば AIC のバイアスも大きくなる傾向があることもわかる．
2.3 T IC 規準
竹内（1976）により提案された T IC は，AIC が持つようなモデルの不一致によるバイアス√が修正された情報量規準である．この T IC では， n(θ̂ − θ0 ) が漸近的に (p + q) 変量正規分布に従うことを利用して，AIC のバイアス補正項である 2(p + q) を再評価している．今，(p + q) × (p + q) 行列 J (ξ̂) と I(ξ̂) を（2.3）式の p 変量正規分布の密度関数 f (y i |ξ) と（2.4）式の対数尤度関数 (ξ|Y ) を用いて以下のように定義する．
(2.12)1
J (ξ̂) = −
n



 

∂2

(ξ|Y )

∂ξ∂ξ

,

ξ=ξ̂

138

統計数理 第 53 巻

第 1 号 2005



1 ∂
∂

I(ξ̂) =
log{f (y i |ξ)}  log{f (y i |ξ)}

n i=1 ∂ξ∂ξn(2.13).

ξ=ξ̂

T IC の一般公式（竹内（1976））から，候補のモデル M（2.1）における T IC は，
T IC = −2(ξ̂|Y ) + 2tr(I(ξ̂)J (ξ̂)−1 ),(2.14)と定義される．上記のトレース部分を更に計算すると，(2.15)T IC = np log(2π) + n log |Σ(θ̂)| + (n + 2)tr(SΣ(θ̂)−1 )1 (ε̂ ⊗ ε̂i )D(θ̂|S)(ε̂i ⊗ ε̂i )n i=1 in+−vec(Σ(θ̂)−1/2 SΣ(θ̂)−1/2 ) D(θ̂|S)vec(Σ(θ̂)−1/2 SΣ(θ̂)−1/2 ),となる（導出に関する詳細は Appendix A.2 を参照）．ただし ε̂i は候補のモデル M の下での標準化された残差ベクトルの推定量，(2.16)

ε̂i = Σ(θ̂)−1/2 (y i − ȳ),である．式（2.15）から，T IC は，AIC とは異なり，多変量尖度の推定量に依存し真の分布と仮定した分布の違いを考慮に入れた情報量規準となっていることがわかる．実際に εi が独立に同一な分布に従えば，T IC は AIC が持つような定数バイアスを修正しているということが知られている．つまり，BT IC を T IC のバイアスとすると，
BT IC = RKL − E∗Y [T IC] = O(n−1 ),となる．
ところが，理論上ではバイアスが補正されていると雖も，実際の状況において T IC のバイ
アス補正項がうまく機能しないことがある．そのような T IC の欠点は，多変量線形モデルにおける変数選択の場合ではあるが，Yanagihara（2004a），Fujikoshi et al.（2005）にも述べられている．一般的に，多変量分布の尖度は次元 p が大きくなればなるほど大きくなる．しかしながら，Mardia（1970）に代表される従来の最小二乗残差に基づく推定量により大きい尖度の推定を行えば，巨大な標本数が無ければ極端な過小推定をしていまい，その結果非常に大きなバイアスが生じるということが知られている（Yanagihara（2004b））．以上のようなことが（2.15）
式の T IC でも起き，適度な標本数であっても非正規性に対する補正がうまく行われないこと√がある．また T IC は，標本数 n に関してリスク RKL（2.6）を n(θ̂ − θ0 ) でテーラー展開し，その第一項目による近似を基にした推定量である．もし p が大きくなれば，無視した第二項目
以降が大きくなり，このテーラー展開による近似の精度は悪くなる．そのため，次元が大きい
ときに T IC のバイアスが大きくなるという現象が起こる．一般的に，共分散構造モデルの次元は大きくならざるをえない．以上のことから，通常の状況であれば T IC と AIC は結果的にほとんど変わらない状況が往々に起こりうる．
2.4 クロス・バリデーションを用いた情報量規準
2.2 と 2.3 章からわかるように，AIC では真の分布によりバイアスが変動し，T IC では多変量尖度の推定がうまくいかなかったり次元の増大により θ̂ の正規近似が悪くなったりして，バ
イアスがうまく補正されないという状況が起こることがあり問題がある．そのため，候補のモデルにおける分布と真の分布との相違に関する影響がより小さくなるような情報量規準が必要となる．本論文では，Yanagihara（2004a）でのアイデアを用い，リスク RKL（2.6）における未来の観測値ベクトルに関する部分だけをクロス・バリデーションによる予測残差に置き換える

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

139

ことで新しい情報量規準を構成する．クロス・バリデーションを用いた情報量規準の構成は多くの研究者により行われている．特に Stone（1974）が有名であり，その規準量は漸近的に AICと同値であることが Stone（1977）により示されている．また，正規多変量線形モデルと正規一
般化多変量分散分析（GMANOVA）モデルにおいて，Fujikoshi et al（．2003）によりクロス・バリ
デーションを用いた情報量規準のバイアスを補正した新しい規準量も提案されている．
新しい情報量規準である修正 AIC（Corrected AIC：AICc ）を以下のように定義する．まず，ȳ [−i] と S [−i] を第 i 番目の個体に関する観測値ベクトル y i を取り除いたデータにより得た以下のような推定量とする．
1 
n
1
y,
yj =
ȳ −
n−1
n−1
n−1 i
n

(2.17)ȳ [−i] =j=i(2.18)

S [−i] =n

1(y − ȳ [−i] )(y j − ȳ [−i] )
n − 1 j=i j

n
=
n−1



1
S−
(y − ȳ)(yi − ȳ)
n−1 i



.

上記の式変形については Fujikoshi et al（．2003）等を参照してほしい．今，θ̂ [−i] を第 i 番目の個体に関する観測値ベクトル y i を取り除いたデータにより得た θ の最尤推定量とする．つまり，(2.19)θ̂[−i] = arg min log |Σ(θ)| + tr(S [−i] Σ(θ)−1 ) .
θ

このとき AICc を，(2.20)

AICc = np log(2π) + n log |Σ(θ̂)|
+

n 


y i − ȳ [−i]Σ(θ̂ [−i] )−1 y i − ȳ [−i]i=1

= np log(2π) + n log |Σ(θ̂)|
+

n
n−1

2 
n

(y i − ȳ) Σ(θ̂[−i] )−1 (y i − ȳ)

i=1と定義する．AICc のバイアスは T IC と同じ
BAICc = RKL − E∗Y [AICc ] = O(n−1 ),である（証明は Appendix A.5 を参照）．しかし，T IC と異なって個別に尖度の推定をする必要がないため，AICc では尖度の推定量による過小推定の問題は発生しない．更に，テーラー展
開を使っていないため，次元の増大に対するバイアス項への影響は T IC に較べ小さくなる．
共分散構造の選択問題においては，Cudeck and Browne（1983），Browne and Cudeck（1989），
De Gooijer（1995）等でクロス・バリデーションを用いたモデル選択の指標が提案されている．
これらの指標では，データを二つ（もしくは複数）の calibration 標本と validation 標本に分けそれらの標本から指標値を計算している．しかしながら，AICc ではそのような手順を踏む必
要はない．
式（2.20）からもわかるように，AICc は第 i 番目の個体に関する観測値ベクトル y i を取り除いたデータにより得た θ の最尤推定量である θ̂ [−i]（2.19）に依存している．そのため，AICcを得るためには n 回の繰り返しが必要となる．もともと θ̂ を得るためには，いくらかの繰り返し計算が必要であり，それだけで多くの計算時間を必要とする．そのため，標本数 n が大きい時には，何らかの手法を用いて計算時間の短縮をはかったほうが良い．標本数 n が十分大き140

統計数理 第 53 巻

第 1 号 2005

い時には，Fujisawa and Eguchi（2003）のような近似公式を使い，θ̂ [−i] の計算時間を短縮することもできる．その近似公式は，
1
H(θ̂|S)−1 ∆(θ̂) vec(Σ(θ̂)−1/2 W i Σ(θ̂)−1/2 ),(2.21)
θ̂ [−i] ≈ θ̂ +n−1である（導出に関する詳細は Appendix A.4 を参照）．ただし H(·|·) は Appendix A.1 での（A.5）
式で定義される尤度関数の二次偏微分に関する行列であり，W i は，
n(2.22)
Wi = S −(y − ȳ)(yi − ȳ) ,n−1 iである．しかしながら，近似公式（2.21）において，観測値の次元 p が標本数 n に較べ大きいときは，第二項目以降の次元の影響が大きくなり，近似の精度が悪くなることがある．またこれは T IC でも同じことが言えるのだが，H(θ̂|S) が正則でない場合に近似公式（2.21）は使えないものとなる．更に，正則の場合でも H(θ̂|S)−1 が不安定な値になることもある．以上のことから，この近似公式（2.21）を使用するときは十分に注意しなければならない．
3.

数値実験

第 2 章では，共分散構造モデルの選択問題に関する従来の情報量規準のバイアスに関する問
題点を考え，それらを改善するような新しい情報量規準を提案した．本章では，新しく提案した AICc のバイアスが実際に AIC や T IC のバイアスと較べ小さくなっていることをシミュ
レーション実験により確かめる．
3.1 実験に用いたモデル
本論文の数値実験では，真のモデル M∗（2.2）の共分散行列として Hu et al.（1992）で使用している以下のような 15 次元の因子分析モデルを使用した．
Σ∗ = Λ∗ Λ∗ + Ψ∗ .

(3.1)

ただし，因子負荷量行列 Λ∗ と 独自因子分散行列 Ψ∗ は，
Λ∗ = (I 3 ⊗ l∗ )Q1/2∗ ,Ψ∗ = diag(13 ⊗ w∗ ),である．ここでの 1k は全ての成分が 1 である k × 1 ベクトルであり，l∗ ，Q∗ ，w∗ はそれぞれ

 
0.70



0.70 l∗ = 
0.75 
,
 0.80 
1.0 0.3 0.4
Q =  0.3 1.0 0.5  ,
∗

0.4

0.5

1.0

0.80



0.5100
0.5100 w∗ = 
0.4375 , 0.3600 
0.3600である．非正規誤差は，Yuan and Bentler（1997）で提案されている，以下のようなモデルにより発生させた．
データモデル．x1 , . . . , xk (k ≥ p) はそれぞれ独立に同一な平均 0 分散 1 を持つ分布に従う
確率変数とし，x = (x1 , . . . , xk ) とする．また r は，x と独立で，分散 1 を持つ分布に従う確率変数とする．このとき，データを(3.2)

ε = rA x,

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

141で発生させる．ただし A = (a1 , . . . , ak ) は rank(A) = p である k × p 行列で A A = I p を満たすものとする．このとき，式（2.11）で定義される ε の多変量尖度は，(1)

κ4 = β(η − 3)k(ai ai )2 + (β − 1)p(p + 2),i=1となる．ただし E[x4i ] = η であり E[r 4 ] = β である．
χ2f を自由度 f の χ2 分布に従う確率変数とし，16 × 15 行列 B を



B=

I 15
115



(I 15 + 115 115 )−1/2 ,とする．本論文では，（3.2）式のデータモデルに対して，以下のような x，A，r により誤差分布を発生させた．この設定は Yanagihara（2004b）での数値実験と同じものである．
誤差分布 1 （正規分布）: xi ∼ i.i.d. N(0, 1)，r = 1，A = I 15（κ(1)．
4 = 0）
 2(1)
誤差分布 2 （t-分布）: xi ∼ i.i.d. N(0, 1)，r = 6/χ8 ，A = I 15（κ4 = 127.5）．
√
誤差分布 3 （一様分布）: xi = wi /(5/ 3)，wi ∼ i.i.d. U(−5, 5)，r = 1，A = B（κ(1)=
4−16.875）．
√
誤差分布 4 （χ2 分布）: xi = (wi − 4)/2 2，wi = χ24 ，r = 6/χ28 ，A = B（κ(1)．
4 = 190.781）
√
1/8
1/2
1/4
誤差分布 5 （対数正規分布）: xi = (wi − e )/(ee− 1)，log wi ∼ i.i.d. N(0, 1/4)，
(1)
r = 6/χ28 ，A = B（κ4 = 251.920）．
3.2 実験結果
1/2
3.1 章での式（3.1）と（3.2）を用いて発生させた擬似データ y i = 115 + Σ∗ εi (i = 1, . . . , 200)に対して，因子分析モデル，
Σ(θ) = ΛΛ + Ψ,(Λ = (λ1 , . . . , λm ), Ψ = diag(ψ1 , . . . , ψp )),を当てはめ，その因子数 m を AIC ，T IC ，AICc を用いて決定する実験を 500 回繰り返し行った．このモデルでは m = 3 のときが真のモデルである．Λ と Ψ の実際の推定アルゴリズムは，
柳井 他（（1990），pp. 49–92）等を参考にしていただきたい．因子数 m の探索範囲は 1 ≤ m ≤ 6とした．図 1 には，それぞれの誤差モデルで得たリスクと AIC ，T IC ，AICc の平均を示している．図中の横軸の目盛りは候補のモデルでの因子数に対応している．更に，リスクは実線で，AIC の平均は〇で，T IC の平均は□で，AICc の平均は△で表している．また，表 1 はそれ
ぞれの誤差モデルにおけるリスクとバイアス，モデルの選択確率（度数（％））を示している．
図と表から，AIC のバイアスは真の分布の尖度が大きくなれば大きくなっていることがわ
かる．またそれに伴い，尖度が大きくなればリスクが最小になるモデル（ここでは真のモデル（m = 3））を選択する確率が低くなっている．T IC は AIC に較べバイアスが小さくなっていることが見てとれるが，その大きさはまだまだ大きく，十分にバイアスを補正しているとは言い難い．このような現象は，誤差分布 2，4，5 のように尖度が大きい分布において顕著に現れている．その原因は尖度の過小推定にあると考えられる．また，T IC のバイアスは AIC のバイ
アスよりも小さくなっていると雖も，モデル選択確率は逆に AIC よりも悪くなっている．一方，AICc のバイアスはほとんどの場合に最も小さくなっている．それは特に誤差分布 5 で顕著な結果として現れている．更に，モデルの選択確率は，尖度が小さい分布である誤差分布 1と 3 でこそ AIC よりも低い確率になっているがその差は小さいものであり，尖度が大きい分
布では AIC よりもはるかに高い確率となっている．また，すべての分布において真のモデル

142

統計数理 第 53 巻

図 1.

第 1 号 2005

リスク RKL と AIC ，T IC ，AICc の平均．

を選ぶ確率が 90％を超えているため，AICc を用いたモデル選択法では，AIC や T IC を用いた手法と異なり選択確率における分布間の違いはあまり無い．
以上のことから，AIC を使った選択法は，真の分布が正規分布，または真の分布の尖度が小さい場合にはほとんど問題ないが，真の分布の尖度が大きくなればそのバイアスも大きくな
り，リスクが最小になるモデルを選ぶ確率が下がる．一方，AICc はすべて誤差分布においてバイアスを最も良く補正しており，それを用いた選択法は，リスクが最小になるモデルを選ぶ
確率もどの分布においてもあまり変わらない．そのため，AICc を用いた選択手法は，候補のモデルと真のモデルの分布間の相違に関する影響を最も受けにくい手法であると言える．
4.

国民性データへの適用

本章では，国民性調査第 11 回全国調査 M 型調査票のデータ（坂元 他（2004））に「日本人の満足感」（前田（1995））に関する因子分析モデルを当てはめ，2.4 章で提案した新しい情報量規準
AICc を用いて潜在因子の最適構造の探索を行う．
4.1 データと候補のモデル
本論文の解析では， 前田（1995）と同様の表 2 で示した 14 項目を変数として使用した．具体的な質問項目文については，坂元 他（2004）を参照して欲しい．それぞれの項目において その他 と ‘D.K.’ 回答が一つでもある個体を取り除いた 835 名のデータを用いて解析を行った．

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題
表 1.

143

数値実験の結果（リスク，バイアスの平均，モデル選択確率）．

また， 健康状態 （#2.80）と記したものは，前田（1995）と同様に，五つの下位項目，すなわ
ち， 頭痛・偏頭痛 （#2.80a）， 背中の痛み （#2.80b）， いらいら （#2.80c）， うつ状態（#2.80d）， 不眠症 （#2.80e），に対してここ一ヶ月の間に かかったことがない と回答した
場合に 1 ポイントを与えて，その合計点を新たな変数と定義したものである．
分析の目的は，上記のデータに対して因子分析モデル，(4.1)

Σ(θ) = ΛΦΛ + Ψ,をあてはめ，「日本人の満足感」に関した潜在因子の最適構造を探索することである．ただし，
Λ は成分に部分的に 0 を持つ 14 × m 因子負荷量行列，Φ は対角成分が 1 であるような m × m

144

統計数理 第 53 巻

第 1 号 2005

表 2. 分析に使用した項目と関連する構成概念の内容．

因子間相関行列，また，Ψ は Ψ = diag(ψ1 , . . . , ψ14 ) の 14 × 14 独自因子分散行列である．前田（1995）では，「日本人の満足感」というものを「個人生活への満足感」と「社会への満足感」という二つの構成概念によって検討している．解析で取り扱ったデータでは，「個人生活への満足感」が潜在因子になっていると考えられる変数は 1，8，9，10，11，12，であり，「社会への満足感」が潜在因子となっていると考えられる変数は 2，3，4，5，6，7，13，14 である．本論
文では，「日本人の満足感」は「個人の満足感」と「社会への満足感」の二つの構成概念によって表現されていると捉え，表 2 における変数の構成概念から「個人の満足感」と「社会への満足感」を複数に分割した潜在因子によって「日本人の満足感」を説明することを試みる．実際のあてはめに用いたモデルは以下の六つのモデルである．
モデル 1：潜在因子は， 個人の満足感 と 社会への満足感 （m = 2）．
モデル 2：潜在因子は， 個人の満足感 と 帰属意識 と 日本の物質的豊かさへの評価（m = 3）．
モデル 3：潜在因子は， 個人の満足感 と 豊かさへの評価 と 社会への満足 （m = 3）．
モデル 4：潜在因子は， 個人の満足感 と 生活水準 と 社会への評価 （m = 3）．
モデル 5：潜在因子は， 個人の満足感 と 生活水準 と 日本の物質的豊かさへの評価と 社会への満足 （m = 4）．
モデル 6：潜在因子は， 健康状態 と 生活水準 と 日本の物質的豊かさへの評価 と個人生活への満足 と 社会への満足 （m = 5）．
それぞれのモデルにおけるパス図を図 2 に，因子負荷量行列 Λ（4.1）のパターンを Appendix
A.6 に示している．
4.2 解析結果
表 3 は，候補のモデル 1 ∼ 6 を当てはめたときのそれぞれの情報量規準 AIC ，T IC ，AICcとモデルの適合度 GFI の結果である．また，図 3 はそれぞれのモデルにおける情報量規準の変動を示している．図 3 において，AIC は〇，T IC は□，AICc は△で表されている．また，図中の横軸の目盛りは使用したモデルの番号に対応している．表 3 と図 3 から，候補のモデル

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

図 2.

それぞれのモデルでのパス図．

145

146

統計数理 第 53 巻

第 1 号 2005

表 3. それぞれの候補のモデルでの情報量規準の比較．

図 3.

それぞれのモデルでの情報量規準．

の中で最適なモデルは，すべてのモデルの中で最も GFI が大きいモデル 5 であることがわかる．
表 3 と図 3 から，それぞれのモデルにおいて，情報量規準 AIC ，T IC ，AICc の差があまり無いことに気づく．この原因には，データの尖度があまり大きくないということが関係していると予想される．使用したデータの変数別のヒストグラムを図 4 に示す．この図から，変数
1 は大きく正規分布から逸脱しているが，それ以外の変数ではそれほど正規分布から逸脱していないことが見てとれる．実際に真の分布が正規分布のときに不偏推定量になるような多変量
尖度の推定量（Mardia（1970），Browne（1984）等参照），
(4.2)

(1)

κ̂4 =

(n + 1)   2ε̃i ε̃i − p(p + 2),n(n − 1) i=1n−1/2を計算してみると，κ̂(1)(yi − ȳ) である．この数値からは
4 = 26.67 となる．ただし ε̃i = S
データが正規分布から大きく逸脱しているとは言えない．しかしながら，上記の推定量のような従来の最小二乗残差に基づく推定量では，データの非正規性に対してそれほど敏感に反応しない（Yanagihara（2004b））ため，データの非正規性を過小に評価していることが考えられる．
そのため，Yanagihara（2004b）で提案されている（4.2）よりも非正規性に対して敏感に反応する
多変量尖度の推定量でもデータの非正規性を評価してみる．Yanagihara（2004b）ではチューニ(1)
ングパラメータ α によって残差ベクトルの長さを調節した多変量尖度の推定量 κ̃4 (α) を提案している．この推定量は α = 0 のとき 式（4.2）での推定量と一致し，α = 0.5, 1 のときはそれ非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

図 4.

ぞれ，
(1)

κ̃4 (0.5) =(4.4)

κ̃4 (1) =(1)

それぞれの変数のヒストグラム．

(ε̃i ε̃i )(n − p − 3) − p(p + 2),n(n − 1) i=1 1 − ε̃i ε̃i /(n − 1)
2

n

(4.3)

147

(n − p − 3)(n − p − 5) n(n − 1)2i=1nε̃i ε̃i1 − ε̃i ε̃i /(n − 1)

2

− p(p + 2),(1)となる．それぞれ計算すると，κ̃(1)
（4.2）
（
，4.3）
（
，4.4）に4 (0.5) = 28.12， κ̃4 (1) = 29.61 となる．
おいて，p(p + 2) = 224 であることに注意すると，データの尖度は大きくなく，ほとんど正規
分布とみなせることがわかる（厳密性を求めるならここで仮説検定を行うべきであろうが，本論文ではそれを省略する）．データが正規分布からそれほど逸脱していないので，それぞれの情報量規準の結果がさほど変わらないものになったと推測できる．
最後に，情報量規準により最適なモデルであると選択されたモデル 5 での解析結果を図 5 に示す．この結果から，本論文で取り扱ったデータでは，「日本人の満足感」は，「個人の満足感」と「社会への満足感」を三つに分割した，「生活水準」，「日本の物質的豊かさへの評価」，「社会への満足」という因子で説明されていることになる．また，それぞれの潜在変数間の相関はそれほど高くないという結果も得た．このモデルでの GFI は 0.964 であり，現在求められている， GFI が 0.95 付近 という規準（狩野 · 三浦（2002），p. 141）を満たしている．しかしながら，それぞれの因子負荷には 0.26 等の値が低いものも存在し，さらに説明力が高いモデルを探索する必要があることを指摘したい．

148

統計数理 第 53 巻

図 5.

5.

第 1 号 2005

最適なモデルにおける解析結果．

おわりに従来の統計手法の多くは正規性の仮定の下で提案されたものである．統計解析のソフトウエ
アでは正規性の下で提案された手法が一般的であることや真の分布は未知であることを考えると，解析の初期段階においてこのような手法を用いることはある意味妥当ではある．しかしながら，国民性全国調査によるデータはカテゴリカルデータであり，明らかに正規分布に従っていない．そのため，従来の正規性の仮定に基づく手法をそのままデータに適用したときに，仮定した正規分布とデータが従う真の分布の相違が解析結果にどのような影響を与えるかを調べ，もし必要であればその非正規性の影響を小さくするような修正を施す必要がある．本論文
では，分布に正規性を仮定したときの共分散構造モデルの選択問題に関する AIC のバイアスを漸近展開式から評価し，非正規性の影響が小さくなるような新しいバイアス補正情報量規準である AICc を提案した．実際のシミュレーション実験により，この AICc のバイアスは AIC
や T IC のバイアスよりも小さくなり，更に AICc を用いたモデル選択ではリスクを最小にするモデルを選ぶ確率がどのような分布に関してもほぼ同じくらいになることが確かめられた．
国民性調査データを使った解析例においては，正規分布からさほど外れていないデータであったため，AIC ，T IC ，AICc のどの情報量規準もほぼ同じ値になり，すべての情報量規準が同じモデルを最適なモデルと選んだ．この結果は，本論文で提案した情報量規準をモデルの最適化に用いた場合，データが正規分布に 近い 場合でも，分布に正規性を仮定して導出された AIC を用いた結果と同じような結果となるということを示している．これは非常に重要な結果である．もしデータが正規分布に 近い 場合，非正規性の下でバイアスを補正した情報量規準が正規性の仮定の下での情報量規準とまったく違う結果を与えてしまうのであれば，バイアスを補正した情報量規準を安心して使用することができない．その点を考えれば，4 章
での解析結果は非常に重要な結果である．
しかしながら，前田（1995）で述べているように，4 章で選ばれたモデルは候補のモデルの中で最適なモデルであり，本当の意味での真のモデルを示しているのではない．また，4 章の最後にも述べたが，因子負荷の問題も残っており，さらに説明力が高いモデルを探索する必要がある．また，カテゴリカルデータに正規分布の仮定をしても良いか ? という問題もある．カテゴ
リー数が大きければさほど問題はないが，厳密な意味で言えば，因子分析モデルもカテゴリカ

149

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

ルデータにおける因子分析モデル（柳井 他（（2000），pp. 129–145）等）を用いる必要がある．カテ
ゴリカルデータでは，観測値ベクトル y i に E[y i ] = µ，Cov[y i ] = Σ(θ) である離散分布を仮定する方法，yi は平均 µ，共分散行列 Σ(θ) を持つ連続分布からの潜在変数をカテゴライズした変数であると仮定する方法等が考えられる．前者の場合は問題無いが，後者の場合 E[y i ] = µ，
Cov[y i ] = Σ(θ) であるため，本論文で用いた推定法では推定結果にバイアスが生じ，AICc でも正確にバイアスを補正できない．このようなモデルで推定を行った結果と，正規性の仮定の下で推定を行った結果との違いも調べる必要がある．また，本論文での数値実験では，モデル選択法での非正規性の影響を調べるということを中心に行ったが，その他にも標本数，因子負荷量の大きさ等も影響を与えていると考えられる（Bandalos（1993））．更に，仮定した分布と真の分布が異なった場合に情報量規準のバイアスを補正するその他の手法として，ブートスト
ラップ法を用いた方法があり，Ishiguro et al.（1997）ではブートストラップ法でバイアスを補正した情報量規準である EIC を提案している．本論文では EIC との比較を行っていないが，多変量モデルに関していえば，標本数が次元に対して十分大きくなければブートストラップ法
は良い近似を与えないということが Wakaki et al.（2002）や Fujikoshi et al.（2005）等で報告されている．これらが検証できる数値実験も行う必要がある．以上のような点がこれからの課題
になる．
謝

辞

因子分析に関して有益なコメントを頂いた 筑波大学 上市秀雄先生，ハワイ大学 林健太郎先生，本論文の内容に対して適切なコメントを頂いた 統計数理研究所 坂元慶行先生に厚く御礼申し上げます．またこの特集号の編集委員，さらに二人の査読者の方々には大変丁寧に論文を読んで戴き，改訂のための助言を戴きました．ここに深く感謝の意を表します．
Appendix
A.1 Σ(θ) の偏微分に関する表記
共分散行列 Σ(θ) に関する標準化された j 次偏微分行列を
Σ̇a1 ···aj (θ ∗ ) = Σ(θ)−1/2








∂j
Σ(θ) Σ(θ)−1/2 ,∂θa1 · · · ∂θaj
θ=θ∗とする．この行列を用い，以下のような p2 × q の一次偏微分行列を定義する．
(A.1)∆(θ∗ ) = (vec(Σ̇1 (θ∗ )), . . . , vec(Σ̇q (θ∗ )))

∂
 .
= Σ(θ)−1/2 ⊗ Σ(θ)−1/2
vec(Σ(θ))∂θθ=θ∗また θ に関する対数尤度関数をnnnplog(2π) − log |Σ(θ)| − tr(SΣ(θ)−1 ),(A.2)
ln (θ|S) = −222とし，この対数尤度に対し，一次と二次の偏微分に基づく q × 1 ベクトル，q × q 行列を，2g(θ∗ |S) = −
n






∂
ln (θ|S)
∂θ
θ=θ∗



,

2
H(θ∗ |S) = −
n






∂2
ln (θ|S)
∂θ∂θ 
θ=θ ∗,と定義する．今，候補のモデル M（2.1）の下での共分散行列で標準化した残差行列を(A.3)

R(θ|S) = Σ(θ)−1/2 {Σ(θ) − S} Σ(θ)−1/2 ,150

統計数理 第 53 巻

第 1 号 2005とおくと，共分散構造モデルにおける g(θ|S) と H(θ|S) は以下のようになる．
(A.4)

g(θ|S) = ∆(θ) vec(R(θ|S)),(A.5)

H(θ|S) = ∆(θ) ∆(θ) − 2∆(θ) {I p ⊗ R(θ|S)} ∆(θ)
+

q





(q) 

tr R(θ|S)Σ̇ab (θ) e(q)
a eb

.

a,b

ただし，e(q)は第 j 番目の成分が 1 で残りの成分が 0 となる q × 1 ベクトルである．
j
A.2 T IC の導出
竹内（1976）によって提案された T IC の一般公式は，式（2.14）で定義される．共分散構造モ
デルにおける J (ξ̂)（2.12）と I(ξ̂)（2.13）は，


Σ(θ̂)−1J (ξ̂) = O1
H(θ̂|S)
2

O



,


Î 11
I(ξ̂) =  
Î 12


,

Î 12
Î 22

である．ただし O は全ての成分が 0 であるような行列であり，Î 11 = Σ(θ̂)−1 SΣ(θ̂)−1 ,1Î 12 = Σ(θ̂)−1/2
2
1
Î 22 = ∆(θ̂)
4



 
n

1
ε̂i ⊗ ε̂i ε̂in i=1∆(θ̂),

vec(R(θ̂|S))vec(R(θ)|S) +1
ε̂i ε̂i ⊗ ε̂i ε̂i
n i=1
n



−vec(I p − R(θ̂|S))vec(I p − R(θ̂|S))∆(θ̂),

である．ここでの ε̂i は式（2.16）で定義される候補のモデル M（2.1）の下での残差ベクトルの推定量であり，R(·|·) は式（A.3）で定義した候補のモデル M（2.1）の下での共分散行列で標準化した残差行列とする．今，θ̂ は最尤推定量であるので，g(θ̂|S) = ∆(θ̂) vec(R(θ̂|S)) = 0 を満たすことを利用すると Î 22 は，Î 22 =1∆(θ̂)
4 n1ε̂i ε̂i ⊗ ε̂i ε̂in i=1−vec(Σ(θ̂)−1/2SΣ(θ̂)

−1/2)vec(Σ(θ̂)−1/2SΣ(θ̂)−1/2 )∆(θ̂),となる．ここで，D(θ|S) = ∆(θ)H(θ|S)−1 ∆(θ) ,(A.6)という p2 × p2 行列関数を定義すると，求める式（2.14）中のトレースは，
2tr(I(ξ̂)J (ξ̂)−1 )= 2tr(SΣ(θ̂)−1 ) +1 (ε̂ ⊗ ε̂i )D(θ̂|S)(ε̂i ⊗ ε̂i )n i=1 in−vec(Σ(θ̂)−1/2 SΣ(θ̂)−1/2 ) D(θ̂|S)vec(Σ(θ̂)−1/2 SΣ(θ̂)−1/2 ),となる．よって，式（2.15）を得ることができる．

151

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

A.3 AIC のバイアス
観測値ベクトル y i が独立同一分布に従う場合，T IC のリスク RKL（2.6）に対するバイアス
は RKL − E∗Y [T IC] = O(n−1 ) となる（竹内（1976））．この関係式を使うと AIC のバイアスは，
BAIC = RKL − E∗Y [AIC](A.7)= E∗Y [T IC] − E∗Y [AIC] + O(n−1 )= 2E∗Y [tr(I(ξ̂)J (ξ̂)−1 )] − 2(p + q) + O(n−1 )= B1 − B2 + B3 − 2(p + q) + O(n−1 ),となる．ただし，

B1 = 2E∗Y tr(SΣ(θ̂)−1 ) ,

B2 = E∗Y vec(Σ(θ̂)−1/2 SΣ(θ̂)−1/2 ) D(θ̂|S)vec(Σ(θ̂)−1/2 SΣ(θ̂)−1/2 ) ,
1 ∗  EY (ε̂i ⊗ ε̂i )D(θ̂|S)(ε̂i ⊗ ε̂i ) ,
n i=1
n

B3 =

である．今，E∗Y [θ̂] → θ0 ，E∗Y [S] → Σ∗ (n → ∞) であることを利用すると，



∗



(A.8)

E∗Y tr(SΣ(θ̂)−1 ) = tr(Ω Ω) + O(n−1 ),(A.9)

EY vec(Σ(θ̂)−1/2 SΣ(θ̂)−1/2 ) D(θ̂|S)vec(Σ(θ̂)−1/2 SΣ(θ̂)−1/2 )= vec(Ω Ω) D(θ0 |Σ∗ )vec(Ω Ω) + O(n−1 ),1/2となることがわかる．ただし Ω = Σ∗ Σ(θ0 )−1/2 である．また同様に，ε̂i = Ω εi + Op (n−1/2 )を利用すると，1 ∗  EY (ε̂i ⊗ ε̂i )D(θ̂|S)(ε̂i ⊗ ε̂i )n i=1n(A.10)

= tr  {Ω ⊗ Ω} D(θ0 |Σ∗ ) Ω ⊗ Ω
 p




(p)
(p)
∗
(p)
(p)
 + O(n−1 ),×EY [εa εb εc εd ] ea eb ⊗ ec ed
a,b,c,dを得ることができる．ここで，K p を Magnus and Neudecker（（1999），pp. 46–47）で定義されている，任意の p × p 行列 C に対して K p vec(C) = vec(C  ) を満たす p2 × p2 の交換行列（commutation matrix）とする．Yanagihara et al.（2005）より，この交換行列 K p は，
Kp =

p


(p) 

e(p)
a eb

(p)

⊗ eb e(p)
a



,

a,b

となる．これと四次キュムラント（2.10）の関係式より，
p


(p) 

E∗Y [εa εb εc εd ] e(p)
a eb

(p) 

⊗ e(p)c ed= Γκ + I p2 + K p + vec(I p )vec(I p ) ,a,b,c,dを得る．ただし Γκ は式（2.9）で定義される真の分布の尖度に依存する p2 × p2 行列である．上記の期待値を式（A.10）に代入すると以下のようになる．
tr

 





Ω ⊗ Ω D(θ0 |Σ∗ ) {Ω ⊗ Ω} Γκ + I p2 + K p + vec(I p )vec(I p )

+ O(n−1 ).

152

統計数理 第 53 巻

第 1 号 2005

ここで，Magnus and Neudecker（（1999），p. 47）より，任意の p × p 行列 C に対して，K p (C ⊗C) = (C ⊗ C)K p であり，また，もし C が対称行列であれば K p vec(C) = vec(C) を満たすということを利用すると，∆(θ) の定義式（A.1）より，Σ̇α (θ) (α = 1, . . . , q) は対称行列である
ので，
K p (Ω ⊗ Ω)∆(θ0 ) = (Ω ⊗ Ω)∆(θ0 ),となることがわかる．D(θ 0 |Σ∗ ) = ∆(θ 0 )H(θ0 |Σ∗ )−1 ∆(θ0 ) であるので，上記の式より，
K p (Ω ⊗ Ω)D(θ0 |Σ∗ ) = (Ω ⊗ Ω)D(θ0 |Σ∗ ),である．よって，
1 ∗  EY (ε̂i ⊗ ε̂i )D(θ̂|S)(ε̂i ⊗ ε̂i )n i=1n(A.11)

= tr Γκ {Ω ⊗ Ω} D(θ0 |Σ∗ ) Ω ⊗ Ω
+2tr

 





Ω Ω ⊗ Ω Ω D(θ0 |Σ∗ )+vec(Ω Ω) D(θ0 |Σ∗ )vec(Ω Ω) + O(n−1 ),である．式（A.8），（A.9），（A.11）を（A.7）に代入すれば，AIC のバイアスは（2.8）式のようになることがわかる．
A.4 クロス・バリデーションに関する近似公式
第 i 番目の個体の観測値ベクトル y i (i = 1, . . . , n) を除いたデータで定義される θ の対数尤度関数は，(A.12)

ln−1 (θ|S [−i] ) = −n−1 log(2π) + log |Σ(θ)| + tr(S [−i] Σ(θ)−1 ) ,2である．ここで，式（2.18）と W i（2.22）から S [−i] = S + W i /(n − 1) であることがわかるので，式（A.12）は式（A.2）でのすべてのデータを使った対数尤度関数 ln (θ|S) と W i を使って，ln−1 (θ|S [−i] ) =n−1
1ln (θ|S) − tr(W i Σ(θ)−1 ),n2と書き換えることができる．ここで，θ̂[−i] は ln−1 (θ|S [−i] ) を最大にする解であることを利用
すると，




∂
ln−1 (θ|S [−i] )
∂θ
θ=θ̂

[−i]

n−1
1
=−
g(θ̂[−i] |S) + ∆(θ̂[−i] ) vec(Σ(θ̂[−i] )−1/2 W i Σ(θ̂[−i] )−1/2 ) = 0,22である．ただし g(·|·) は式（A.4）で与えられるすべてのデータに基づく対数尤度関数 ln (·|·) の一次微分に関するベクトルであり，0 はすべての成分が 0 であるようなベクトルである．今，
z i = (zi1 , . . . , ziq ) = θ̂[−i] − θ̂ とすると，テーラー展開により上記の式は，−

n−11g(θ̂|S) + H(θ̂|S)z i + ∆(θ̂) vec(Σ(θ̂)−1/2 W i Σ(θ̂)−1/2 ) + O( z i ) = 0,22

( z i → 0) と展開できる．ここで，θ̂ は最尤推定量であるので g(θ̂|S) = 0 となることに注意
すれば，
zi =

1
H(θ̂|S)−1 ∆(θ̂) vec(Σ(θ̂)−1/2 W i Σ(θ̂)−1/2 ) + O( z i ) .
n−1

153

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

W i = Op (1)，S = Σ∗ +Op (n−1/2 )，θ̂ = θ0 +Op (n−1/2 ) であるので，上記の式から z i = Op (n−1 )であることがわかる．よって，(A.13)

θ̂[−i] = θ̂ +1
H(θ̂|S)−1 ∆(θ̂) vec(Σ(θ̂)−1/2 W i Σ(θ̂)−1/2 ) + Op (n−2 ),n−1となり，この式の第二項目まで考慮に入れたものがクロス・バリデーションの近似公式（2.21）となる．
A.5 AICc のバイアス式（2.17）から ȳ [−i] − µ は，ȳ [−i] − µ =n1(ȳ − µ) −(y − µ),n−1n−1 iと変形できる．この式より AICc の期待値は，



n
 

∗

E∗Y [AICc ] = np log(2π) + nE∗Y log |Σ(θ̂)| +
+

2 
n
n
n−1


∗

EY tr(Σ∗ Σ(θ̂[−i] )−1 )i=1EY (ȳ − µ) Σ(θ̂[−i] )−1 (ȳ − µ)





i=1
n


 ∗
2n

−1
−
E
(ȳ
−
µ)
Σ(
θ̂
)
(y
−
µ)
+ O(n−1 ),[−i]Yi(n − 1)2 i=1となることがわかる．今，式（A.13）を，θ̂ [−i] − θ̂ =1
di + Op (n−2 ),nと変形する．ただし，
di = H(θ̂|S)−1 ∆(θ̂) vec(Σ(θ̂)−1/2 {S − (y i − ȳ)(y i − ȳ) }Σ(θ̂)−1/2 ),である．ここで，nS = ni=1 (y i − ȳ)(yi − ȳ) であることに気がつけば，これを利用すると，テーラー展開により任意の p × p 行列 C に関して，ni=1 di = 0 となる．

tr(CΣ(θ̂ [−i] )−1 ) = tr(CΣ(θ̂)−1 )1+ vec(Σ(θ̂)−1/2 CΣ(θ̂)−1/2 ) ∆(θ̂) di + Op (n−2 ),nであるので，ntr(CΣ(θ̂[−i] )−1 ) = ntr(CΣ(θ̂)−1 ) + Op (n−1 ),i=1となることがわかる．今，(ȳ − µ) Σ(θ̂[−i] )−1 (ȳ − µ) = tr((ȳ − µ) (ȳ − µ) Σ(θ̂[−i] )−1 ) である
ので，
n

i=1
n



tr(Σ∗ Σ(θ̂[−i] )−1 ) = ntr(Σ∗ Σ(θ̂)−1 ) + Op (n−1 ),(ȳ − µ) Σ(θ̂[−i] )−1 (ȳ − µ) = n (ȳ − µ) Σ(θ̂)−1 (ȳ − µ) + Op (n−1 ),

i=1

n

i=1

(ȳ − µ) Σ(θ̂[−i] )−1 (y i − µ) = n (ȳ − µ) Σ(θ̂)−1 (ȳ − µ) + Op (n−1 ),154

統計数理 第 53 巻

第 1 号 2005

である．以上より，n2 /(n − 1)2 = 1 + 2/n + O(n−2 )，−2n/(n − 1)2 = −2/n + O(n−2 ) である
ことを考慮すると，



E∗Y [AICc ] = np log(2π) + nE∗Y log |Σ(θ̂)|









+nE∗Y tr(Σ∗ Σ(θ̂)−1 ) + nE∗Y (ȳ − µ) Σ(θ̂)−1 (ȳ − µ) + O(n−1 ).

よって式（2.6）のリスク RKL と比較すれば，AICc のバイアスのオーダーは O(n−1 ) となることがわかる．
A.6 実解析での因子負荷量行列のパターン
4 章での実解析で用いた候補のモデル（4.1）における因子負荷量行列のパターンは以下のようなものである．ただし因子負荷量行列 Λj の添え字は，それぞれのモデルの番号に対応するも
のである．
λ1
0 
 0 λ2  0 λ3  00 λλ4  0 λ5  0 λ6 
Λ1 = 
 λ8 07  , λ9 0  λ10 0  λ11 0  λ12 0 
0
λ13






Λ4 = 






0

λ14

λ1
0
0
0
0
0
0
λ8
λ9
λ10
λ11
λ12
0
0

0
λ2
λ3
λ4
0
0
0
0
0
0
0
0
0
0

0
0
0
0
λ5
λ6
λ7
0
0
0
0
0
λ13
λ14

λ1
0
0 
 0 λ2 0 
 0 λ3 0 
 0 λ4 0 
 00 00 λλ5  0 0 λ6 
Λ2 = 
 λ8 0 07  , λ9 0 0  λ10 0 0  λ11 0 0  λ12 0 0 
0
λ13 0
0

λ14

0

0

0

0

λ1
0
0 
 0 λ2 0 
 0 λ3 0 
 0 λ4 00  00 λλ5 0 
 0 λ6 0 
Λ3 = 
 λ8 07 0  , λ9 0 0  λ10 0 0  λ11 0 0  λ12 0 0 
0
0 λ13
0

0

λ14

0

0

0


λ1
0
0
0 
λ1 0
0
0
0 
0
λ
0
λ
0
0
0
0
0 




2
2

 0 λ3 0 0 
 0 λ3 0 0 0  0 λ4 0 0 
 0 λ4 0 0 0  00 00 λλ5 00  00 00 λλ5 00 00 6
 0 0 λ6 0  , Λ5 =  λ 0 07 0  , Λ6 =  00 00 λ07 λ0 00  .
 λ89 0 0 0 
 0 0 0 λ89 0  λ10 0 0 0 
 0 0 0 λ10 0  λ11 0 0 0 
 0 0 0 λ11 0  λ12 0 0 0 
 0 0 0 λ12 0 
 0 0 0 λ13  0 0 0 0 λ13 
λ14

0

λ14

参 考 文 献Akaike, H.（1973）. Information theory and an extension of the maximum likelihood principle, 2nd, 267–281,
International Symposium on Information Theory（eds. B. N. Petrov and F. Csáki）
Akadémiai Kiado, Budapest.
Akaike, H.（1987）. Factor analysis and AIC, Psychometrika, 52, 317–332.
Bandalos, D. L.（1993）. Factors inﬂuencing cross-validation of conﬁrmatory factor analysis models,Multivariate Behavioral Research, 28, 351–342.

非正規データにおける情報量規準を用いた共分散構造モデルの選択問題

155

Bozdogan, H.（1987）. Model selection and Akaike’s information criterion（AIC）: The general theory
and its analytical extensions, Psychometrika, 52, 345–370.
Browne, M. W.（1984）. Asymptotically distribution-free methods for the analysis of covariance structures, British Journal of Mathematical and Statistical Psychology, 37, 62–83.
Browne, M. W. and Cudeck, R.（1989）. Single sample cross-validation indices for covariance structures, Multivariate Behavioral Research, 24, 445–455.
Cudeck, R. and Browne, M. W.（1983）. Cross-validation of covariance structures, Multivariate Behavioral Research, 18, 147–167.
De Gooijer, J. G.（1995）. Cross-validation criteria for covariance structures, Communications in
Statistics, Simulation and Computation, 24, 111–147.
Fujikoshi, Y. and Satoh, K.（1997）. Modiﬁed AIC and Cp in multivariate linear regression, Biometrika,
84, 707–716.
Fujikoshi, Y., Noguchi, T., Ohtaki, M. and Yanagihara, H.（2003）. Corrected versions of crossvalidation criteria for selecting multivariate regression and growth curve models, Annals of
the Institute of Statistical Mathematics, 55, 537–553.
Fujikoshi, Y., Yanagihara, H. and Wakaki, H.（2005）. Bias corrections of some criteria for selection
multivariate linear regression models in a general case, American Journal of Mathematical andManagement Sciences（in press）.
Fujisawa, H. and Eguchi, S.（2003）. Robust estimation in the normal mixture model, Research Memo.,
No. 867, The Institute of Statistical Mathematics, Tokyo.
Hu, L., Bentler, P. M. and Kano, Y.（1992）. Can test statistics in covariance structure analysis betrusted ?, Psychological Bulletin, 112, 351–362.
Hurvich, C. M. and Tsai, C. L.（1989）. Regression and time series model selection in small samples,
Biometrika, 76, 297–307.
Ishiguro, M., Sakamoto, Y. and Kitagawa, G.（1997）. Bootstrapping log likelihood and EIC, an
extension of AIC, Annals of the Institute of Statistical Mathematics, 49, 411–434.
Isogai, T.（1983）. On measures of multivariate skewness and kurtosis, Mathematica Japonica, 28,
251–261.
Kano, Y.（1992）. Robust statistics for test-of-independence and related structural models, Statistics& Probability Letters, 15, 21–26.

狩野 裕，三浦麻子（2002）．『AMOS，EQS，CALIS によるグラフィカル多変量解析 目でみる共分散構造分析 （増補版）』，現代数学社，京都．
Kano, Y., Berkane, M. and Bentler, P. M.（1990）. Covariance structure analysis with heterogeneous
kurtosis parameters, Biometrika, 77, 575–585.
前田忠彦（1995）． 日本人の満足感の構造とその規定因に関する因果モデル 共分散構造分析の「日本人の国民性調査」への適用 ，統計数理，43，141–160．
Magnus, J. R. and Neudecker, H.（1999）. Matrix Diﬀerential Calculus with Applications in Statistics
and Econometrics, Revised Edition, John Wiley & Sons, New York.
Mardia, K. V.（1970）. Measures of multivariate skewness and kurtosis with applications, Biometrika,
57, 519–530.
坂元慶行，中村 隆，前田忠彦，土屋隆裕（2004）． 国民性調査第 11 回全国調査 2003 年全国調査 ，
統計数理研究所研究リポート，No. 92．
Satorra, A. and Bentler, P. M.（1988）. Scaling corrections for chi-square statistics in covariancestructure analysis, Proceedings of the Business and Economic Statistics Section, 308–313,
American Statistical Association, Alexandria.
Stone, M.（1974）. Cross-validatory choice and assessment of Statistical predictions, Journal of theRoyal Statistical Society, Series B, 36, 111–147.

156

統計数理 第 53 巻

第 1 号 2005

Stone, M.（1977）. An asymptotic equivalence of choice of model by cross-validation and Akaike’s
criterion, Journal of the Royal Statistical Society, Series B, 39, 44–47.
Sugiura, N.（1978）. Further analysis of the data by Akaike’s information criterion and the ﬁnite
corrections, Communications in Statistics, Theory and Methods, A7, 13–26.
竹内 啓（1976）． 情報統計量の分布とモデルの適切さの規準，数理科学，153，12–18．
Wakaki, H., Yanagihara, H. and Fujikoshi, Y.（2002）. Asymptotic expansion of the null distributionsof test statistics for multivariate linear hypothesis under nonnormality, Hiroshima Mathematical Journal, 32, 17–50.
Yanagihara, H.（2004a）. Corrected version of AIC for selecting multivariate normal linear regression
models in a general nonnormal case, Tech. Report, No. 04–04, Statistical Research Group,
Hiroshima University, Hiroshima.
Yanagihara, H.（2004b）. A family of estimators for multivariate kurtosis in a nonnormal linear regression model, Tech. Report, No. 04–08, Statistical Research Group, Hiroshima University,
Hiroshima.
Yanagihara, H., Tonda, T. and Matsumoto, C.（2005）. The eﬀects of nonnormality on asymptotic
distributions of some likelihood ratio criteria for testing covariance structures under normal.
assumption, Journal of Multivariate Analysis（in press）
柳井晴夫，繁桝算男，前川眞一，市川雅数（1990）．『因子分析 その理論と方法 』，朝倉出版，東京．
Yuan, K.-H. and Bentler, P. M.（1997）. Generating multivariate distributions with speciﬁed marginal
skewness and kurtosis, SoftStat’ 97
Advances in Statistical Software 6 （eds. W. Bandilla
and F. Faulbaum）, 385–391, Lucius & Lucius, Stuttgart.
Yuan, K.-H. and Bentler, P. M.（1998）. Normal theory based test statistics in structural equation
modelling, British Journal of Mathematical and Statistical Psychology, 51, 289–309.
Yuan, K.-H. and Bentler, P. M.（1999）. On normal theory and associated test statistics in covariancestructure analysis under two classes of nonnormal distributions, Statistica Sinica, 9, 831–853.

Proceedings of the Institute of Statistical Mathematics Vol. 53, No. 1, 133–157 (2005)

157

Selection of Covariance Structure Models in Nonnormal Data
by Using Information Criterion: An Application to Data
from the Survey on the Japanese National Character
Hirokazu Yanagihara
University of Tsukuba

This paper deals with a bias correction of Akaike’s Information Criterion(AIC) fora selection of covariance structure models under a general nonnormal case, when we
assume that an observation is distributed according to the normal distribution. The AIC
under the normal assumption has a constant bias that depends on a kurtosis of the true
distribution. We correct this bias by partly using the cross-validation method. We verify
that our criterion is better than the AIC and T IC by calculating a theoretical bias
and conducting a numerical experiment. Further, an illustrative analysis applying factor
models with respect to Japanese people’s level of satisfaction to an actual data set fromthe eleventh nationwide survey on the Japanese national character is also demonstrated.
Using our new criterion, we search the best model from the candidate ones.

Key words: Akaike’s information criterion, cross-validation, factor analysis model, model selection,
nonnormality.