観察データの推測の限界−揺らぎモデルアプローチ

観察データの推測の限界−揺らぎモデルアプローチ−
江口真透 （統計数理研究所, 総合大学院大学統計科学専攻）
1. はじめに
観察研究から得られたデータに対して標準的な統計的方法を単純に適用することは不適
切な結論を導く恐れがある．このことは標準的な統計的方法が正当化される暗黙の仮定を
観察データが満たさないことから生じる問題である．このようにデータ適用性を巡って統
計学の基本的なパラダイムにこの重大な問題が議論された文献がある [1, 3, 6, 7, 11]．ミッ
シング，打ち切り，未観測交絡，競合リスクなど不完全観測を伴うデータにはランダム・サ
ンプルの仮定への乖離が生じる可能性があり，しかもほとんどの場合，その乖離の有無は
データからは検証できない [9]．
バイアスに対する研究は膨大な文献に著され，様々なアプローチが展開されている．こ
の発表において紹介するアプローチは，仮定された統計モデルのミス・スペシケイション
を表現するために，統計モデル包含する管状の近傍を考察した公式 [4] を利用している．モ
デルからの乖離を記述するパラメータを導き，そのパラメータによる感度解析をすること
から始まった [2]．例題としてヒロシマ被爆者の生存解析において遅延登録によるバイアス
の問題を考察した [10]．この研究を発展させて，全ての選択バイアスの影響を許容した信
頼領域の構成法を提案した [3]．その補正公式は簡明で最尤推定量の漸近分散行列を 2 倍に
することによって得られる． 間接喫煙から肺がんの危険性を評価する問題が動機を与えて
いる例として考察される．
2. モデルの不確定性
必ずしも全てが観測できない仮想的な確率変数 Z を考えよう．その統計モデル
M = {f = f (z, θ) : θ ∈ Θ}
が与えられたとする．統計モデル M から仮想的な観測ベクトル z 1 , · · · , z n が得られると
考える．実際には観測ベクトルの分布が正確に統計モデル M にあると考えるよりも，管
状近傍
N (M ) = {g : min D(g , f ) < }
f ∈M
にあると考えた方が多くの場合で正確である．ここで D は分布の隔たりを表す尺度を表す．
（ここで，また以下でも混乱が生じる恐れがないときは分布と密度を混同して使う．
）また，
これからの議論では θ の推測には尤度法を考えることから D はクルバック・ライブラーの
ダイバージェンスを考えることが自然である [12]．確率変数 Z が連続分布を持つならば，
N (M ) は関数空間となるが，十分小さい を選べば，
L (f ) = {g ∈ N (M ) : D(g , f ) = ∗min D(g , f ∗ )}
f ∈M
は分布 f を貫くリーフとなり，統計モデル M 上の全ての分布 f をとってリーフ L (f )
を重ね合わせると Z の可能な分布全体の中で葉層を成すことが分かる．このようなリーフ
L (f ) は
g = g (z, , u ) = f (z, θ) exp{u (z, θ) − 12 2 }
(1)
と書ける [2]. ここで u は分布 f の下で平均 0, 分散 1 を持つ確率変数とする．以下の議
論では, u が更にモデル M のスコアーベクトル s と無相関であると仮定する．このよ
うな仮定からデータ分布 (1) を持つ完全データ z 1 , · · · , z n の最尤推定量 θ̂
はバイアスを
またないことが判る：
Ef (θ̂ ) − θ = O(n−1 ) + O(2 ).
(2)
ここで Ef はデータが f からのランダムサンプルとしたとき期待値を表す．この関係が
意味あるためには
1
= O(n− 2 )
(3)
でなくてはならない．これを概パラメトリックな仮定と呼び [5]，以下の議論を通してこの
仮定をおく．次節では現実には不完全データしか得られない状況を考える．このとき最尤
推定量は，観測の不完全性によって生じたモデルの不確定性によって (2) のような漸近不
偏性が保たれなくなり，選択バイアスが現われることを考察しよう．
3. 不完全な観測
前節では完全データのモデルの不確定性は理想的に最尤推定量のバイアスは漸近的に
無視可能となる枠組みが示された． この節ではミッシング，打ち切りなどの不完全観測を
伴う観察によるデータによって，この枠組みがどのように誘導されるか考えよう．データ
の不完全性によって，もたらされる推定量の分散とバイアスの影響を表す漸近的な公式を
与える．
不完全データ確率ベクトル Y が完全データ確率ベクトル Z に対して Y = h(Z) と書か
れたとする．ここで h は多対１の関数とする．このとき Y のモデル分布は, Z のモデル分
布 f (z, θ) から誘導されて，形式的に
f (y, θ) =
−1 ()
f (z, θ)dz
(4)
と書かれる．ここで h−1 (y) = {z : h(z) = y} を表す．このようにして Z のフィッシャー
情報行列 I は不完全変換 h によって Y のフィッシャー情報行列 I は情報の損失を受け
る．実際，Y のスコアーベクトル s は s の Y = y の条件付期待値で表され，情報損失
は I − I = var(s − s ) となる．もし不完全データ y 1 , · · · , y n がモデル分布 f に従っ
ているならば，θ の最尤推定量 θ̂ は θ の漸近的不偏推定量で漸近正規性
√
D
n(θ̂ − θ) −→ N(0, I−1)
(5)
−1
を持つ．ここで N(0, I−1
) は平均ベクトル 0, 分散行列 I を持つ正規分布を表す．この性
質より，点推定 θ̂ を中心とする水準 α を持つ信頼領域は
C0 (rα ) = {δ : n(θ̂ − θ)T I (θ̂ − θ) ≤ rα2 }
(6)
として与えられる．ここで rα2 はカイ 2 乗分布の α 分位点とする．これらの統計的方法の漸
近有効性も標準的な考察から得られる．しかしながら，これらの性質は，不完全データが
モデル分布 f に従っているという仮定よって支持されるものである．観察データがミッシ
ングなどの不完全観測を伴って得られたときに，この仮定が成立している保証はない．
これから，モデルの不確定性を考慮して完全データの分布を (1) で定義した g である
と仮定すると，上の式 (4) の考察と同様に，Y のデータ分布は
g (y, θ, , u ) = f (y, θ) exp{u (y, θ) − 12 2 ρ }
(7)
となる．ここで u は u の条件 Y = y の条件付期待値で，ρ は u の分散を表す．θ の
最尤推定量 θ̂ はモデル分布 f の仮定の下では，上の議論のように標準的な漸近性 (5) を
持つが，データ分布 g であれば，どうなるだろうか？実は，完全データの場合に仮定され
た漸近不偏性 (2) は不完全関数 h によって崩される．一方で漸近分散はオーダーの仮定 (3)
の下では影響を受けないことが示される．実際には，その影響は高次のオーダーに吸収さ
れる．
不完全データによる最尤推定量 θ̂ のバイアスは
b = Ef (θ̂ − θ̂ ) = 2 Ef {u I−1 s }
(8)
と表現される．バイアスベクター b の長さは次の不等式を持つことが主張される．
定理 1. 完全データのモデル分布 (1) において，モデル不確定性を与える u が集合
Uf = {u = u (z, θ) : Ef (u ) = 0, Ef (u2 ) = 1, Ef (u s ) = 0}
にあると仮定する．このとき，任意の u ∈ Uf に対して, (8) で定められる最尤推定量 θ̂ のバイアスベクトル b は
2
−1
−1
bT
I b ≤ I (I − I )
(9)
を満たす．ここで は作用素ノルムを表す．等号は，ある定数ベクトル ω が存在して
u (y, θ) = ω T s (y, θ),
ω T (I−1 − I −1 )−1 ω = 1
(10)
のときに限り成立する．
注意 1. 定理 1 の不等式 (9) の等号を満たす (10) で与えられる Y のデータ分布 g は
g (y, , u ) = f (y, θ + ω) + O(3 )
(11)
となる．この最悪の選択バイアスをもたらす場合はモデルの識別性が不成立であり，ω が
推定不可能であることが結論される．このように選択バイアスの問題は 2 重の困難が横た
わる．一つは観察データを標準的な方法で解析すると無視できないバイアスが生じて標準
的な推測が壊れてしまう困難さであり，もう一つはその最悪となるケース (10) に近づくと，
そのバイアスが推定不能となるという困難さに会う．すなわち，観察データの解析者は，い
つも，大きな誤りを導くかもしれな標準的な統計方法を妥当性を疑いながら使わなければ
ならない困難がある．
この問題の実際的な対処法として，ω を推定する代わりに可能な値を与えてデータ分
布 (7) の下で θ を推定する方法が提案されている [2]．このようなアプローチを一般に感度
分析と呼び，ω を感度パラメータと呼ばれる [9]. しかしながら感度分析で分かることは，
感度パラメータ ω を適当に与えて統計方法がどの位の影響を受けるかを計ることしかでき
ない．
次の節では，この感度分析をある反事実仮想を行うことによって，自動的に選択バイア
スの影響の受け入れるアプローチを紹介する．これによって，選択バイアスがあったとし
ても信頼水準や P 値を保証するように統計量を補正する簡単な方法が導かれる．
3. 信頼領域の最悪評価
前節で議論されたように分散の影響は漸近的に無視可能であり，古典的な水準 α の信
頼領域の形状は保たれるがバイアスの影響によって中心が揺らいでしまう．この選択バイ
アスを許容して水準 α を保つような信頼領域の補正を提案する．この補正された信頼領域
は検出不可能な選択バイアスによって揺らいだ古典的な信頼領域を常に包含することが示
される．
最尤推定量 θ̂ の廻りの古典的な水準 α の信頼領域 C0 (rα ) は (6) のように与えられる．
前の節の枠組みに従うと観測できるのは不完全データベクトル Y であって完全データベク
トル Z は観測できない．にもかかわらず，仮想的に完全観測値ベクトル z 1 , · · · , z n が観測
されて，それに基づく最尤推定量 θ̂ が得られたとしよう．このように θ の 2 つの最尤推
定量 θ̂ と θ̂ の同時分布は (1) で与えられるデータ分布 g の下で以下のような漸近正規
性を持つ．
√
n
θ̂ − θ
θ̂ − θ̂
D
−→ N
√
n
b
b
,
I−1 ,
I−1 − I −1
−1
−1
I − I , I−1 − I −1
ここで，2 節と同様に (3) の仮定する．これより，次の条件付き漸近分布と漸近周辺分布
√
√
D
n(θ̂ − θ)|( n(θ̂ − θ̂ ) = u) −→ N(u, I −1 ))
√
D
n(θ̂ − θ̂ ) −→ N(u, I−1 − I −1 )
が得られる．このように 2 つの水準 α の信頼領域
C(u, rα ) = {θ : (θ − θ̂ − u)T I (θ − θ̂ − u) ≤ rα2 }
B(rα ) = {u : uT (I−1 − I −1 )−1 u ≤ rα2 }
が考えられる．この議論は観測不能な Z からの最尤推定 θ̂ を用いたものであるから現実
のデータ解析には意味を成さない．実際には条件付けをする u はデータから決めることが
できないがここでは u ∈ B(rα ) という仮定をする．このことは最尤推定量 θ̂ を導いた仮
想的な実現値ベクトル z 1 , · · · , z n は十分モデル分布 f にフィットしていることを意味す
る．もしこの仮定が満たされてないのなら，理想的なモデル分布から乖離したデータ分布
g に従う Z が不完全変換 Y = h(Z) に伴って Y のデータ分布 g が誘導されるという枠
組みの出発点の漸近不偏性 (2) に矛盾する．しかし，漸近不偏性を満たすようにデータ分
布 g が (1) とパラメトライズできるのでこの仮定は自然である．漸近的な選択バイアスが
u = 0 であるという仮説が水準 α で採択されるという条件付き信頼領域は
C1 (rα ) =
C(u, rα )
∈B(rα )
(12)
で表される．以上の議論に対して次の被覆定理が主張される．
定理 2. 上式 (12) で定義された C1 (rα ) に対して
√
C1 (rα ) ⊆ C0 ( 2rα )
が成立する．
注意 2. 定理 2 は古典的な信頼領域の半径を
√
2 倍すると，起こりうる全ての選択バイア
スのケースに対して信頼水準 α を保証することが主張される．これは，モデル分布の下で
最尤推定量の漸近分散を 2 倍して信頼楕円を構成したことと等価になる．これより「分散
2 倍ルール」と呼ぶ．
注意 3. もし変換 h : X → Y による情報損失がないならば, データ分布 g もまた漸近不
偏性を満たすので選択バイアスは生じない. このことは C1 (rα ) も自然に C0 (rα ) となるこ
とに対応する．このように分散 2 倍ルールは最悪評価を与えている．一般に情報損失行列
√
I −1 I の最小固有値が 12 以下ならば C1 (rα ) が C0 ( 2rα ) の境界に接する．すなわち最小の
最悪評価を与えていることになる．下の図 1 では 4 点で境界に交わっていることが見れる．
1
0.5
-1
-0.5
0.5
1
-0.5
-1
√
Figure 1: C1 (rα ) ⊆ C0 ( 2rα )
4. 終わりに
この発表では完全な観測から不完全性を伴ったときに (3) のオーダーでモデルからの弱
い乖離を与える揺らぎモデルアプローチを考察した．可能な全ての方向のモデルの乖離に
対して漸近的に信頼水準 α を保証する簡便な信頼領域の構成を分散 2 倍ルールによって与
えた．しかしながら，意味のある可能な全てのモデルへの乖離を大域的に特徴付けること
は困難で一般論は未完成である．メタ解析の出版バイアスの内容では大域的な最悪評価が
与えられている [8] が，このアプローチとの融合が今後の課題として挙げられる．
社会が情報化時代の深化する局面を向え，様々な重要なデータは観察研究から従来の規
模や範囲を超えて生産されるようになってきている．今後，この問題は更に，重要になる
と思われる．観察研究の推測に対して数理統計学が長い歴史を通して培った多くの定理が
適用できない．ここに数理統計の今後の方向を導く潜在的な可能性があるが，未だその見
通しの良い考えは提出されてないように思える．このように多くの未解決問題が山積され
ている現状で観察研究のデータ解析の実際に即応した数理統計学の研究者の理論的貢献が
切望されている．
References
[1] Cochran, W. G. The Planning of Observational Studies of Human Populations. J. Royal
Statist. Soc. A, 128 (1965) 234-255
[2] Copas, J. and Eguchi, S. Local sensitivity approximation for selectivity bias. J. Royal
Statist. Soc. B 63 (2001) 871-895.
[3] Copas, J. and Eguchi, S. Local model uncertainty and incomplete data bias (with
discussion). J. Royal Statist. Soc. B, 67 (2005) 459-512 .
[4] Eguchi, S. and Copas, J. A class of logistic-type discriminant functions, Biometrika,
89, 1-22 (2002).
[5] Eguchi, S. and Copas, J. A class of local likelihood methods and near-parametric asymptotics. J. Royal Statist. Soc. B 60, (1998) 709-724.
[6] Greenland, S. (2005) Multiple-bias modelling for analysis of observational data (with
discussion). J. Royal Statist. Soc. A, 168, 267-306.
[7] Heckman, J.J. Sample bias as a specification error. Econometrica 47, 1 (1979) 153-162.
[8] Henmi, M., Copas, J. and Eguchi, S. Confidence intervals and P-values for meta analysis
with publication bias. Biometrics 63 (2007) 475-482.
[9] Little, R. J. A. and Rubin, D. B. (2002) Statistical Analysis with Missing Data. 2nd
edn. New York: Wiley.
[10] Matsuura, M. and Eguchi, S. Modeling late entry bias in survival analysis. Biometrics,
61, 559-566 (2005).
[11] Rosenbaum, P. R. and Rubin, D. B. Reducing bias in observational studies using
subclassification on the propensity score. J. Amer. Statist. Assoc. 79, 387 (1984) 516524.
[12] White, H. (1982) Maximum likelihood estimation of misspecified models. Econometrica, 50, 1-25.