統計数理（2014）
第 62 巻 第 1 号 103–122
c 2014 統計数理研究所


特集「疫学研究のデザインとデータ解析：最近の理論的展開と実践」［総合報告］

欠測データに対するセミパラメトリックな解析法——その理論的背景について——

逸見 昌之†（受付 2013 年 12 月 10 日；改訂 2014 年 3 月 12 日；採択 3 月 12 日）

要旨

疫学におけるコホート研究のような観察研究では，様々な要因により統計解析の結果に偏りが生じる可能性があるが，そのうちの 1 つはデータの欠測によるものである．統計的な欠測デー
タの解析法には，大きく分けて最尤法や多重代入法などのパラメトリックな手法と，逆確率重み付け法や二重頑健推定法などのセミパラメトリックな手法があるが，本稿では後者の方に焦点を当てる．後者の方法は，近年主に医学統計学の分野で急速に発展し，実際の適用例はまだそれほど多くはないが，有力な方法として注目されており，今後その需要はさらに高まってくるものと思われる．セミパラメトリックな統計手法は，特に医学統計学の分野では他にも様々な場面で用いられているが，本稿ではその理論的背景の理解にも資するように，セミパラメト
リック推測の一般論に基づいた解説を行う．
キーワード：MAR，影響関数，推定関数，局外接空間，漸近有効性．
1.

はじめに

疫学におけるコホート研究のような観察研究では，様々な要因により統計解析の結果に偏りが生じる可能性がある． 例えば，ある曝露（喫煙など）がある疾患（肺癌など）の原因であるかどうかを調べるような因果推論における交絡などはその代表的なものであるが，データの欠測もまた，結果にバイアスをもたらし得る要因として考慮すべきものである．コホート研究などの疫学研究では，その参加者から常に（最後まで）データが観測できるとは限らない． このとき，もし欠測の仕方がその研究で対象としている結果変数（ある疾患に罹るかどうかなど）に影響のある
ものであれば，つまり，例えば二群比較の際，曝露群において悪い結果となる参加者ほど欠測しやすい傾向にあれば，欠測を考慮しない解析では曝露効果を過小評価してしまう．欠測データを取り扱うための統計的方法は，大きく分けて 2 種類のものがある．1 つは最尤法や多重代入法などのパラメトリックな手法であり，もう 1 つは
（Rubin, 1987; Carpenter and Kenward, 2013）
逆確率重み付け法（inverse probability weighting; Robins et al., 1994; Seaman and White, 2011）
や二重頑健推定法（doubly robust estimation; Scharfstein et al., 1999; Bang and Robins, 2005）などと呼ばれるセミパラメトリックな手法である． 後者の方法は近年，主に医学統計学の分野で急速に発展し，実問題への適用例はまだそれほど多くはないが，有力な方法として注目されており，今後その需要はさらに高まってくるものと思われる．本稿の主な目的は，これらの方法をセミパラメトリック推測の一般論の観点から解説し，その理論的な背景を明らかにすること†

統計数理研究所：〒 190–8562 東京都立川市緑町 10–3

104

統計数理

第 62 巻

第 1 号 2014

である．数学的な議論の詳細（厳密な正則条件等）にはあまり深く立ち入らず，手法導出の道筋をはっきりさせることに重点を置くが，特に理論的な立場から，二重頑健推定法の仕組みとその役割を理解することを 1 つの目標とする．
本稿の構成は以下の通りである． 第 2 節ではまず，統計的な欠測データ解析の出発点として欠測メカニズムの分類について述べ，単純な方法の問題点に触れた後，パラメトリックな解析法として最尤法について簡単に述べる．そして，それと対比しながらセミパラメトリックな解
析法の必要性について触れ，ごく簡単な場合に，逆確率重み付け法と二重頑健推定法を紹介する． 続いて第 3 節では，理論的な取り扱いの準備として，必要な範囲でセミパラメトリック推測の一般論の概説を行う． そしてそれをもとに第 4 節では，二重頑健推定量の導出を中心に，セミパラメトリックな欠測データ解析法について理論的な立場から解説する．最後に第 5 節として，注意事項や関連事項等について述べる．
2.

欠測データ解析の基本事項

2.1 欠測メカニズムの分類
欠測データの解析を行う際に最も留意すべきことの 1 つは，データの欠測がどのようなメカ
ニズムによって生じているかということである．ここでは，欠測データ解析の出発点としてまによって提唱されて以来，広く受け入れられている欠測メカニズムの 3 つの分
ず，Rubin（1976）
類について述べる． ある（疫学）
研究において，n 人の対象者から m 次元の多変量データ(2.1)

Z 1 = (Z11 , . . . , Z1m ), . . . , Z n = (Zn1 , . . . , Znm )を測定する状況を想定し，これらはある確率分布からのランダムサンプル（互いに独立で同一の確率分布に従う確率ベクトル）であるとする．但し，これらのデータは必ずしも全て観測されるとは限らないとし，例えば，対象者 1 に対しては Z11 から Z1m までが観測されるが，対象者 nに対しては Zn1 と Zn2 しか観測されないということが起こり得るとする． このとき，観測の指
示変数 Rij を，Zij が観測されれば 1, 欠測すれば 0 の値を取る確率変数として定義し，対象者 iに対する観測データ（Z i の各成分のうち，実際に観測される部分）
を Z i と Ri1 , . . . , Rim の関数として O(Z i , Ri1 , . . . , Rim ) と表すことにすると，データの欠測メカニズムは以下のように分類．
される（Rubin, 1976; Little and Rubin, 2002）• 完全にランダムな欠測（Missing Completely At Random, MCAR）∀r1 ∈ {0, 1}, . . . , ∀rm ∈ {0, 1},
P (Ri1 = r1 , . . . , Rim = rm |Z i ) = P (Ri1 = r1 , . . . , Rim = rm )• ランダムな欠測（Missing At Random, MAR)∀r1 ∈ {0, 1}, . . . , ∀rm ∈ {0, 1},
P (Ri1 = r1 , . . . , Rim = rm |Z i ) = P (Ri1 = r1 , . . . , Rim = rm |O(Z i , r1 , . . . , rm ))• ランダムでない欠測（Not Missing At Random, NMAR)∃r1 ∈ {0, 1}, . . . , ∃rm ∈ {0, 1} s.t.
P (Ri1 = r1 , . . . , Rim = rm |Z i ) = P (Ri1 = r1 , . . . , Rim = rm |O(Z i , r1 , . . . , rm ))「完全にランダムな欠測（MCAR）」とは，文字通りデータの欠測が完全にランダムに生じていて，上述のようにデータの各観測（あるいは欠測）パターンの確率がデータそのものの値には全く依らないことを意味する．この場合は，観測されたデータのみ用いて解析を行っても統計的推測の結果にバイアスは生じないが（但しデータの不足による精度の低下は起こる），データの（MAR）」
欠測を意図的に制御しない限り実際の場面では稀である． 2 番目の「ランダムな欠測欠測データに対するセミパラメトリックな解析法

105

は，しばしば「観測（あるいは欠測）
確率が観測データのみに依存する」と言い表されるが，より正確に表すと上述のようになる．最も単純な場合として，m = 2 で Zi1 の方は常に観測されるという状況を考えてみると，この場合，観測の指示変数は Ri2 の方だけ考えればよく，MAR は
P (Ri2 = r|Zi1 , Zi2 ) = P (Ri2 = r|Zi1 ) (r = 0, 1) と表される．これはつまり，Zi2 が観測される
Zi1 の値のみによって決まることを意味している．例えば，
かどうかの傾向は（常に観測される）ある集団に対する健康診断において血圧の値を（日を改めて）2 回に分けて測定する場合，1 回目の血圧 Zi1 は全員に対して測定されるが，2 回目の血圧 Zi2 を測定するかどうかは 1 回目の結果を見て健診対象者自身が決められるものとすると，Zi2 の欠測のメカニズムは MAR であると考」は，文字通り MAR でないということである
えられる． 最後の「ランダムでない欠測（NMAR）
が，これはつまりデータのある観測パターンの確率が，観測データだけでなく欠測データの値にも依存してしまうということを意味する．
統計的な欠測データ解析法の多くは MAR の仮定の下で構築されているが，データの欠測メ
カニズムが MAR であるかどうかは観測されたデータだけからは（原理的に）検証できない．したがって，データ取得のデザインなどから MAR であることが保証されていない限りは，まずMAR を仮定して解析を行ったとしても，その仮定からのずれの影響を調べる感度解析を行うことが望ましいとされている．本稿における手法の解説は基本的に MAR の仮定の下で行っていくが，感度解析などについては，例えば第 5 節で挙げる文献等を参照されたい．
2.2 単純な方法の問題点
近年までの疫学研究では，欠測データの取り扱い方として多くの場合，Complete-Case 解析や単一代入法といった単純な方法が用いられているが（Eekhout et al., 2012），これらは非常に限られた状況でしか正当化されない．例えば，上述の血圧測定の例において（記号を見やすくするために）Xi = Zi1 , Yi = Zi2 , Ri = Ri2 (i = 1, . . . , n) とおき，2 回目の血圧 Yi の平均 μ = E(Yi )に興味があるとすると，Complete-Case 解析による μ の推定量は


n
n

1 
μ̂CC =
Ri Yi
Ri
N=
N i=1
i=1

で与えられるが，これは欠測メカニズムが MCAR でないと（漸近的に）バイアスが生じ得る．実μ̂CC は E(Yi |Ri = 1) に確率収束し，漸近バイアス
際，大数の法則より（n → ∞ のとき）
E(Yi |Ri = 1) − E(Yi ) = {E(Yi |Ri = 1) − E(Yi |Ri = 0)}(1 − p) (p = P (Ri = 1))

は MCAR 以外の場合では 0 とは限らない．一方，単一代入法については，例えば回帰代入法による μ の推定量は以下のように与えられる．
1
{Ri Yi + (1 − Ri )m(Xi ; β̂)}.
n i=1
n

(2.2)

μ̂RI =

ここで β̂ は，ある回帰モデル E(Yi |Xi ) = m(Xi ; β) （例えば線型回帰なら，m(Xi ; β) =に対して，β を「観測データだけ」を用いて推定したときの（最小 2β0 + β1 Xi , β = (β0 , β1 )）乗）推定量である．つまり，推定量 μ̂RI は，Yi が欠測しているところについては回帰モデルに基づく予測値で補完してから，全体で標本平均を取ったものである．このとき，欠測のメカニズム
一致性を持ち漸近が MAR で，かつ回帰モデルが正しく特定されていれば，μ̂RI は（μ に対して）的なバイアスは 0 となる．しかしながら，推定の対象を例えば Yi の 2 次モーメント ν = E(Yi2 )に置き換えると，回帰代入法によるその推定量

統計数理

106

第 62 巻

第 1 号 2014

1[Ri Yi2 + (1 − Ri ){m(Xi ; β̂)}2 ]n i=1nν̂RI =は，MAR の下で回帰モデルが正しく特定されていても，もはや一致性を持たない．回帰代入法
以外にも平均値代入法や Hot deck などの方法があるが，一般にこのような単一代入法は正当性が非常に限られるので，注意が必要である．
2.3 パラメトリックな解析法 — 最尤法 —冒頭でも述べたように，本稿の主な目的は欠測データに対するセミパラメトリックな手法について解説することであるが，その前に，より基本的なパラメトリックな手法として，最尤法について簡単に述べておく．
に対して Z i の分布（確率密度関数）が，あるパラメトリックモデル
多変量データ（2.1）
S = {p(z; θ) | θ ∈ Θ ⊂ Rq }に属していると仮定し，Z i の各成分の（条件付き）
同時観測確率もまた
P (Ri1 = r1 , . . . , Rim = rm |Z i = z) = ω(r1 , . . . , rm , z; ψ) (ψ ∈ Ψ ⊂ Rr )と（パラメトリックに）モデル化されているとすると，全対象者の観測データ
Z oi = O(Z i , Ri1 , . . . , Rim ), Ri1 , . . . , Rim (i = 1, . . . , n)に基づくパラメータ θ, ψ の尤度はL(θ, ψ) =n 
p(Z i ; θ)ω(Ri1 , . . . , Rim , Z i ; ψ)dZ mii=1で与えられる．但し，Z m（つまり，Z mi は Z i における欠測部分を表す
i = O(Z i , 1 − Ri1 , . . . , 1 −．ここでもし，データの欠測メカニズムが MAR であるとすると，ω(Ri1 , . . . , Rim , Z i ; ψ)Rim )）は Z oi にのみ依存し Z mi には依らないのでn
L(θ, ψ) =po (Z oi ; θ)ω(Ri1 , . . . , Rim , Z i ; ψ), po (Z oi ; θ) = p(Z i ; θ)dZ mi .
i=1

したがって，MAR に加えて，2 つのパラメータ θ, ψ が互いに無関係に変化し得るとすると，パラメータ θ の最尤推定量は単に
n
n 


Lo (θ) =po (Z oi ; θ) =p(Z i ; θ)dZ m
i
i=1

i=1

を最大化することによって得られる．これは，各対象者に対して多変量データ Z i のどの部分
尤度であり，通常「観測デーが欠測するかが予め分かっているときの（パラメータ θ に対する）と呼ばれる．このように，（a）欠測メカニズムが MAR でかタ尤度（observed-data likelihood）」（distinct）であれば，最尤法においては（同時）観測確率つ（b）パラメータ θ と ψ が互いに無関係
P (Ri1 = r1 , . . . , Rim = rm |Z i ) をモデル化する必要はなく，この意味で（a）と（b）は欠測が無視可能（ignorable）であるための条件と言われるが，これはパラメータの推測法に依存する概念であることに注意すべきである（実際，これから述べるセミパラメトリック推測の場合は，これらの条件が成り立っても欠測は無視可能でない）．MAR という概念は，もともと Rubin（1976）によって，最尤法やベイズ法などのパラメトリックな方法で欠測データを扱うことを想定して導で見るように，セミパラメトリックな方法を用いる際にも重要入されたが，以降の節（第 4 節）欠測データに対するセミパラメトリックな解析法

107な役割を果たす．
2.4 IPW 法と DR 法 — 単純な場合 —観測データから最尤法によって完全データ Z i (i = 1, . . . , n) の分布に関する推測を行う際には，
上述の p(z; θ) のように Z i の分布（確率密度関数）が完全にパラメトリックにモデル化されている必要があるが，実際に推測の対象となる興味あるパラメータは，θ の一部（あるいは多対一関数）であることが多い．例えば，2.2 節で触れた血圧測定の例で見てみると，興味あるパラメータは 2回目の血圧 Yi の平均 μ = E(Yi ) であるが，MAR（P (Ri = r|Xi , Yi ) = P (Ri = r|Xi ) (r = 0, 1)）の下で最尤法によって μ の推定を行おうとすると，Yi だけでなく Xi と Yi の同時分布をモデル化する必要がある．そこで，仮にその同時分布を平均ベクトル m ，分散共分散行列 Σ の 2 変がパラメータ θ に相当
量正規分布とすれば，m と Σ（の各成分をベクトルとしてまとめたもの）
し，興味あるパラメータ μ はその一部となる．さらにもし，ここで MAR を成立させている Xi
がもっと高次元のものであれば，同時分布をモデル化するためのパラメータ数も増大し，モデルの誤特定の可能性も増してくる．一方，セミパラメトリックな方法では，完全データの同時分布をフルにモデル化する必要はない．例えばこの例の場合，セミパラメトリック推定量の 1つである μ の逆確率重み付け推定量（IPW 推定量）は，以下のように与えられる．
1  Ri Yi.
n i=1 e(Xi )nμ̃IP W =ここで e(Xi ) は，Xi が与えられた下で Yi が観測される条件付き確率を表し（i.e. e(Xi ) =，しばしば傾向スコア（propensity score）と呼ばれる．この推定量は，CompleteP (Ri = 1|Xi )）（Ri = 1）の各対象者に対して，Yi の値をその観測確率（傾向スコア）の逆数倍することによっCaseて欠測値を疑似的に復元し，それらについて標本平均を取るという考え方に基づいているが，これによって，μ̃IP W は MAR の下で μ の不偏推定量となる．しかしながら傾向スコア e(Xi ) は
通常未知なので，その場合は e(Xi ) を観測データから，例えばロジスティック回帰モデル
e(Xi ) =exp(α0 + α1 Xi )(= e(Xi ; α), α = (α0 , α1 ))1 + exp(α0 + α1 Xi )等で推定する必要があり，実際は
1  Ri Yi
n i=1 e(Xi ; α̂)nμ̂IP W =の形で用いることが多い．ここで，α̂ は上のモデルに基づく α の最尤推定量を表すが，このとき μ̂IP W は MAR の下で μ の一致推定量になる．IPW 推定量は直観的に分かりやすく，完全データに対するモデル化も不要であるが，通常は傾向スコアに対するモデル化が必要となるのでそれが誤特定されれば，（漸近的に）
バイアスが生じる．また，そのモデルが正しく特定され（推定量の構成そのものには）用いてていたとしても，Yi が欠測している対象者の Xi の情報はいないので，μ に対する推定効率（漸近分散）は一般にあまり良くない．これらの点を改善するであるが，それは以下のように与えられる．
のが二重頑健推定量（DR 推定量）


n
Ri Yi
1
Ri
μ̂DR =+ 1−
m(Xi ; β̂) .
n i=1 e(Xi ; α̂)e(Xi ; α̂)
ここで β̂ は，ある回帰モデル E(Yi |Xi ) = m(Xi ; β) に対して，β を観測データから（最小二乗）
推定したものである．推定量 μ̂DR は，形式的には回帰代入法による推定量（2.2）において，Ri
を Ri /ê(Xi ) に置き換えたものになっているが，MAR の下で以下の性質を持つ．

108

統計数理

第 62 巻

第 1 号 2014

（1） 傾向スコア e(Xi ) に対するモデル e(Xi ; α) と条件付き期待値 E(Yi |Xi ) に対するモデル
m(Xi ; β) のうち，少なくともどちらか一方が正しく特定されていれば，DR（回帰モデル）
推定量 μ̂DR は μ に対して一致性を持つ．
（2） どちらのモデルも正しく特定されていれば，DR 推定量 μ̂DR の漸近分散は IPW 推定量μ̂IP W よりも小さくなり，さらにセミパラメトリック漸近有効となる．
（1）の性質は μ̂DR の二重頑健性（double robustness）と呼ばれ，二重頑健推定量という名称はこにおいて「セミパラメトリック漸近有効」というのは，MAR（と傾の性質に由来する．また，（2）
向スコアに対するパラメトリックモデル）
によって規定されるセミパラメトリックモデルの下推定量のあるクラスを考えたとき，μ̂DR の漸近分
で，μ に対する（一致性と漸近正規性を持つ）
散が，そのクラスに属する推定量の漸近分散の下限（semiparametric eﬃciency bound）に達するという意味であるが，これは次節で明らかになる．
疫学や医学研究で用いられる統計的手法においては，経時データ解析の一般化推定方程式（GEE）や生存時間解析の Cox 回帰モデルなど，データの欠測が全くなかったとしても，そもそもセミパラメトリックな統計モデルを用いることが多い．このような場合は，欠測に対処する際にそもそも最尤法などのパラメトリックな方法を用いることはできず，その意味でも欠測デー
タに対するセミパラメトリックな解析法は重要な役割を果たす．
一般論を述べる前にごく簡単に歴史的なことについて触れておくと，IPW 推定量の基になる考え方自体は古くからあり，例えば標本調査の分野における Horvitz-Thompson 推定量（Horvitz
and Thompson, 1952）なども同様の考え方に基づいているが，欠測データ解析に（セミパラメトリック推測の観点を伴って）導入されたのは比較的最近であり，Robins et al.（1994）が初期のものである．また二重頑健推定量は，その後，セミパラメトリック推測理論に基づく欠測データ
解析法の研究の進展に伴って見い出されたものであるが，Scharfstein et al.（1999）によって最初に導入されて以来，現在でも多くの研究が行われている．
3.

セミパラメトリック推測の一般論

次節において，欠測データに対するセミパラメトリックな解析法，特に二重頑健推定量についてその導出を中心に解説するが，ここではまずその準備として，必要な範囲でセミパラメト
リック推測の一般事項について述べる．
（なお，本節および次節の議論は主に Tsiatis, 2006 に基）づく．ここでは概略のみ述べるが，より詳しくは Tsiatis, 2006 を参照されたい．
に対して Z i の分布（確率密度関数）が，あるセミパラメトリックモデル
多変量データ（2.1）(3.1)

S = {p(z; θ, η) | θ ∈ Θ ⊂ Rr , η ∈ H}に属しているとする（但しここでは，Z 1 , . . . , Z n は全て観測されるとする）．ここで，θ は有限（未次元（r 次元）の興味あるパラメータ，η は無限次元の局外パラメータで，応用上はしばしば
知の）密度関数や回帰関数などとして現れる．例えば，Y を（連続的な）結果変数，X を（連続的な）説明変数ベクトル，θ を有限次元の回帰パラメータとして，回帰モデル
E(Y |X ) = m(X ; θ)

(3.2)

を考えたとき，Z = (Y, X ) の確率密度関数はp(z; θ, η1 , η2 ) = η1 (y − m(x; θ), x)η2 (x) (z = (y, x))と書ける．但しここで，η1 , η2 は

欠測データに対するセミパラメトリックな解析法
η1 (, x)d = 1,109
η1 (, x)d = 0 (∀x),

η2 (x)dx = 1

を満たす
（未知の）
非負値関数である（η1 は X = x の下での  = Y − E(Y |X ) の条件付き密度関．これはつまり，（3.2）のモデルが，θ を興味あるパラ
数，η2 は X の周辺密度関数に相当する）
局外パラメータとするセミパラメトリックモデルを規定していメータ，η1 , η2 を（無限次元の）ることを意味する．
セミパラメトリック推測論の 1 つの主題は，与えられたセミパラメトリックモデル（3.1）の下で，興味ある（有限次元）パラメータ θ に対する「最良の」推定量を見い出すことである．そのためにはまず，どのような推定量のクラスで考えるかということを定めなくてはならないが，ここでは次のような関係式を満たす推定量 θ̂ を考える．
(3.3)

√

1 
n(θ̂ − θ) = √φ(Z i , θ, η) + op (1) (∀θ ∈ Θ, ∀η ∈ H).
n i=1n

ここで，φ(z, θ, η) は
E{φ(Z i , θ, η)} = 0, E{ φ(Z i , θ, η) 2 } < ∞ (∀θ ∈ Θ, ∀η ∈ H)

を満たすある r 次元のベクトル値関数であり，推定量 θ̂ の影響関数（inﬂuence function）と呼ば．式（3.3）を満たす推定量 θ̂ はれる（但し，op (1) は n → ∞ のときに 0 に確率収束する項を表す）
通常，漸近線形推定量（asymptotic linear estimator）と呼ばれるが，影響関数 φ(z, θ, η) は θ̂ に対して一意に定まる（ことが示せる）ので，漸近線形推定量は影響関数によって特徴付けられる．
また，大数の法則と中心極限定理により，漸近線形推定量 θ̂ に対して一致性と漸近正規性が成り立つ．すなわち，n → ∞ のとき√
一致性：θ̂ → θ（確率収束）, 漸近正規性： n(θ̂ − θ) → N (0, Avar(θ̂))（分布収束）.

ここで，Avar(θ̂) は θ̂ の漸近分散共分散行列と呼ばれ，
Avar(θ̂) = E{φ(Z i , θ, η)φ(Z i , θ, η)T }で与えられる．右辺はちょうど影響関数の分散共分散行列に相当するが，漸近線形推定量についてはこの量が（対称行列の順序関係の意味で）小さいほど，漸近的に精度の良い推定量である．
漸近線形推定量の典型的な例は，推定関数（あるいは推定方程式）から得られる推定量（M - 推定量）である． 推定関数とは，例えば，以下のような条件
∂u
E{u(Z i , θ)} = 0, E{ u(Z i , θ) 2 } < ∞, det E(Z i , θ) = 0 (∀θ ∈ Θ, ∀η ∈ H)
∂θを満たす r 次元ベクトル値関数 u(z, θ) のことであるが，ランダムサンプル Z 1 , . . . , Z n に対し，
推定方程式
n


u(Z i , θ) = 0

i=1の解として得られる θ の推定量を θ̂ とすると，ある適当な正則条件の下で θ̂ は一致性を持ち，また以下の関係式が成り立つ．
n  
−1
√
1 
∂u
n(θ̂ − θ) = √u(Z i , θ) + o（
(Z i , θ)Ep 1) (∀θ ∈ Θ, ∀η ∈ H).
∂θ
n i=1

統計数理

110

つまり，M - 推定量 θ̂ は E∂u(Z i , θ)∂θ

第 62 巻−1

第 1 号 2014

u(z, θ) を影響関数とする漸近線形推定量であり，それゆえ漸近正規性が成り立ち，漸近分散共分散行列は(3.4)

Avar(θ̂) ∂u= E(Z i , θ)
∂θ

−1

 
∂u
E{u(Z i , θ)u(Z i , θ) } E(Z i , θ)∂θ
T−T

で与えられる． パラメータ θ の（漸近的な）信頼領域などを構成するためには，（3.4）をデータから推定する必要があるが，通常，（その形から）サンドイッチ推定量と呼ばれる以下のような推定量が用いられる．
(3.5)

ˆ θ̂)
Avar(
 
∂u
= Ê
(Z i , θ̂)
∂θ

−1

 
∂u
Ê{u(Z i , θ̂)u(Z i , θ̂)T } Ê(Z i , θ̂)
∂θ

−T

.

但し，Ê は経験分布による期待値を表す．すなわち Z i の関数 f (Z i ) に対して Ê{f (Z i )} =n1
（推定方程式）は，セミパラメトリックモデルの下でのパラメー
i=1 f (Z i ) である． 推定関数nタ推定法としてよく用いられるものであり，次節で述べる欠測データに対するセミパラメトリッ
クな解析法においても重要な役割を果たすが，推定関数について詳しくは例えば，van der Vaart（2000）の第 5 章，Boos and Stefanski（2013）の第 7 章等を参照されたい．
さて，漸近線形推定量は一致性と漸近正規性を持ち，また M - 推定量がそのクラスに含まれるということから自然な要請であると考えられるが，これだけでは例えば，超有効性と呼ばれる性質を持つ極めて不自然な挙動を示す推定量まで含まれてしまうことになる．そこで，そのようなものを除外するためにある正則条件を付加するのが通常であり，その条件を満たす漸近（RAL とは Regular Asymptotic Linear の略．ここで線形推定量のことを RAL 推定量と呼ぶ．
）
は，正則条件の内容については述べないが，詳しくは Tsiatis, 2006 を参照のこと．
が与えられたとき，興味あるパラメータ θ の RAL 推定量（あセミパラメトリックモデル（3.1）
るいはそれを特徴付ける影響関数）
はどのようなものになるであろうか？ 次にそれを見い出す
ための手がかりになる事実（性質）について述べるが，まず，そのために必要となるいくつかの基礎概念を定義しておく．
セミパラメトリックモデル S の各点 p(z; θ, η) に対し
定義 1.（パラメトリックサブモデル）
p(z; θ, η) ∈ Ssub ⊂ S(3.6)

を満たすパラメトリックモデル
Ssub = {p(z; θ, γ) | θ ∈ Θ ⊂ Rr , γ ∈ Γ ⊂ Rs } (s はある自然数)

を，p(z; θ, η) における S のパラメトリックサブモデルという．
セミパラメトリックモデル S の各点 p(z; θ, η) に対し，上記のパラメト
定義 2.（局外接空間）
リックサブモデル Ssub の局外接空間をTθN,γ（Ssub ) = {λT sγ (z, θ, γ) | λ ∈ Rs }として定義する． ここで，γ は p(z; θ, η) に対応するもの（つまり p(z; θ, η) = p(z; θ, γ)）であり，sγ (z, θ, γ) は γ に関するスコア関数，すなわちsγ (z, θ, γ) =∂log p(z; θ, γ)∂γ

欠測データに対するセミパラメトリックな解析法

111である．このときTθN,η（S) = ∪Ssub TθN,γ（Ssub )(3.7)

を，p(z; θ, η) における S の局外接空間（nuisance tangent space）という．但し，ここでパラメト
リックサブモデル Ssub の局外接空間 TθN,γ（Ssub ) は< h1 , h2 >θ,η = E{h1 (Z)h2 (Z)} (∀h1 , h2 ∈ Hθ,η )を内積とするヒルベルト空間(3.8)

Hθ,η = {h | E{h(Z)} = 0, E[{h(Z)}2 ] < ∞}（Z は p(z; θ, η) を密度関数とする確率ベクトル，h は Z の値域で定義された実数値可測関数）に含まれるとし，∪Ssub は全てのパラメトリックサブモデル Ssub について和（集合）を取ること（つまりその集合の要素から成る無限点列の極限を，上付きの横線は Hθ,η において閉包を取る（3.7）の右辺は，一般であるような Hθ,η の点を全て含める）ことを意味するものとする． また，には Hθ,η の線形部分空間になるとは限らないが，ここではそのようになる場合についてのみ考える．
（実際，応用上現れる多くの例で，そのようになっている．
）
さて，セミパラメトリックモデルの局外接空間と RAL 推定量の影響関数の間には，次のような関係がある．
命題 1.（影響関数の満たすべき条件）セミパラメトリックモデル（3.1）の下で，パラメータ θの RAL 推定量 θ̂ の影響関数 φ(z, θ, η) は，必ず以下の条件を満たす．
(3.9)

φ(z, θ, η) ⊥ TθN,η (S),(3.10)

E{φ(Z, θ, η)sθ (Z, θ, η)T } = Ir (r 次の単位行列).

ここで，（3.9）は影響関数 φ(z, θ, η) の各成分が（ヒルベルト空間 Hθ,η の内積に関して）
局外接空N（3.10）における sθ (z, θ, η) は θ に関するスコア関数を
間 Tθ,η (S) に直交していることを意味し，表す．すなわちsθ (z, θ, η) =∂log p(z; θ, η).
∂θ（3.9）の条件は別の言い方をすれば，影響関数の各成分は常に局外接空間の（Hθ,η における）
直交補空間に属するということであるが，このことから逆に，RAL 推定量の影響関数がどのようなものであるかを知るにはまず，その直交補空間を求めればよいということが示唆される．また，（多次元の）
比例関係にある推定関数から得られる M - 推定量については影響関数と推定関数がので，局外接空間の直交補空間を求めることは推定関数を求めるための手がかりにもなる．
（実
際，推定関数の各成分もその直交補空間に属している．
）
さて，RAL 推定量の漸近分散共分散行列は影響関数の分散共分散行列で与えられたので，条（3.10）
を満たす φ(z, θ, η) のうち，その分散共分散行列が（対称行列の順序関係の意味
件（3.9）とで）最小となるものが見い出せれば，それが漸近分散共分散行列の「下限」を与えることになる．
これについて以下のことが成り立つ．
（3.1）の下で，パ
命題 2.（セミパラメトリック Cramer-Rao 不等式）セミパラメトリックモデル
ラメータ θ の任意の RAL 推定量 θ̂ の漸近分散共分散行列 Avar(θ̂) に対し，以下の不等式が成り立つ．
(3.11)E
T −1
Avar(θ̂) ≥ [E{sE.
θ (Z, θ, η)sθ (Z, θ, η) }]

統計数理

112

第 62 巻

第 1 号 2014

ここで，sE
（の各成分）
を局外接空間 TθN,η (S) の直交
θ (z, θ, η) は θ に関するスコア関数 sθ (z, θ, η)
補空間に直交射影したものとして定義され，θ に関する有効スコア関数（eﬃcient score function）と呼ばれる．また，等号は θ̂ の影響関数がE
T −1 EφE (z, θ, η) = [E{sEsθ (z, θ, η)θ (Z, θ, η)sθ (Z, θ, η) }]で与えられるときに成立し，この φE (z, θ, η) を θ に関する有効影響関数（eﬃcient inﬂuence funcと呼ぶ．
（有効影響関数は，条件（3.9）と（3.10）を満たす φ(z, θ, η) のうち，最小の分散共分tion）
散行列を持つものである．但しこれは，影響関数として対応する RAL 推定量が存在するか否の右辺はしばしば，セミパラメトリック効率限界（semiかに依らないことに注意．
）
不等式（3.11）と呼ばれるが，これはパラメトリックモデルにおける（漸近分散にparametric eﬃciency bound）
対する）Cramer-Rao の下限のセミパラメトリック版に相当する．
有効スコア関数は一般に未知の（無限次元）
局外パラメータにも依存するが，例えば生存時間のような特別なセミパラメトリックモデル
解析における比例ハザードモデル（Cox 回帰モデル）に対しては，興味ある（有限次元の）パラメータのみに依存し，そのような場合には有効スコア関数を基にした推定関数によって，常に漸近有効な RAL 推定量（大域的（漸近）有効推定量）が得られる． しかしながら，例えば経時データ解析における一般化推定方程式（GEE）のセミパラメ
トリックモデル（制限モーメントモデル）
では，有効スコア関数は未知の局外パラメータにも依存する部分（結果変数の分散共分散行列）を含み，有効スコア関数を基に推定方程式を構成する
際にはその部分を推定する必要がある．ただこの場合は，それを推定するためのモデルが誤特は得ら定されていたとしても，（漸近有効性は失われるが）少なくとも一致推定量（RAL 推定量）れるという性質を有しており，このような推定量は，そのモデルが正しく特定されているときのみ漸近有効性を持つという意味で，局所（漸近）有効推定量と呼ばれる．これらについて詳しの第 4–5 章を参照されたい．
くは，Tsiatis（2006）
4.

セミパラメトリックな解析法 — 二重頑健推定量の導出 —2.4 節において，ごく簡単な場合に，欠測データに対するセミパラメトリック推定量として逆
確率重み付け推定量（IPW 推定量）と二重頑健推定量（DR 推定量）を紹介したが，本節ではセミパラメトリック推測の一般論に基づき，さらにより一般の設定の下でこれらの推定量について述べる．
4.1 観測データセミパラメトリックモデル
2.3 節では多変量データ（2.1）
に対し，Z i の確率密度関数がパラメトリックモデルに属しているとしたが，ここでは以下のようなセミパラメトリックモデル
S F = {p(z; β, η) | β ∈ B ⊂ Rp , η ∈ H}に属しているとする．ここで，β は有限次元（p 次元）の興味あるパラメータ，η は無限次元の局（つまり Zij が観測されれば 1, 欠測すれば 0 の外パラメータである．Rij を Zij の観測指示変数
値を取る確率変数）とし，データの欠測メカニズムが MAR であるとすると，2.3 節と同様の議論により，対象者 i に対する観測データ
Z oi = O(Z i , Ri ), Ri = (Ri1 , . . . , Rim )Tの（同時）確率密度関数は欠測データに対するセミパラメトリックな解析法

113

q(z o , r) = po (z o ; β, η)ω(z o , r)と書ける．但しここで



(4.1)

po (z o ; β, η) =p(z; β, η)dz m (z = (z o , z m )),(4.2)

ω(z o , r) = P (Ri = r|O(Z i , r) = z o )である． すなわち，観測データの分布もまた，これらによって規定されるセミパラメトリック
モデルに属する． 興味あるパラメータ β の二重頑健推定量は，このセミパラメトリックモデルに対して前節の一般論を適用することにより導かれるが，以下ではそのことについて順を追って説明する．
4.2 観測データ局外接空間とその直交補空間 1 — データ観測確率が既知の場合 —まず，データの（条件付き）観測確率（4.2）
が既知の場合について考える．前節の命題 1 によれ
ば，興味あるパラメータに対する RAL 推定量の（影響関数の）クラスを知るには，まずセミパラ
メトリックモデルの局外接空間とその直交補空間を求めることが重要であったが，観測データ
セミパラメトリックモデル
SO = {q(z o , r; β, η) | β ∈ B ⊂ Rp , η ∈ H}の局外接空間およびその直交補空間は，以下のように与えられる．
命題 3.（観測データ局外接空間とその直交補空間— データ観測確率が既知の場合 —）完全デー
局外接空間を T N F ，完全データヒルベルト空間
タモデル S F の（点 p(z; β, η) における）
HF = {h | E{h(Z)} = 0, E[{h(Z)}2 ] < ∞}における T N F の直交補空間を ΛF とし，また，観測データモデル SO の（点 q(z o , r; β, η) における）局外接空間を TηN O ，観測データヒルベルト空間
HO = {g | E{g(Z o , R)} = 0, E[{g(Z o , R)}2 ] < ∞}F
Oにおける TηN O の直交補空間を ΛO（線形写像）
η とする．このとき，H から H への線形作用素
O
F∗
L と H から H への線形作用素 L をそれぞれL(h)(Z o , R) = E{h(Z)|Z o , R} (∀h ∈ HF )L∗ (g)(Z) = E{g(Z o , R)|Z} (∀g ∈ HO )によって定義すると，(4.3)

TηN O = L(T N F ) = {L(h) | h ∈ T N F }

(4.4)

∗ −1
ΛO
(ΛF ) = {g ∈ HO | L∗ (g) ∈ ΛF }
η =（L )

が成り立つ．
∗（h は ΛF の任意の要素）の解集合を（h ∈ ΛF について）
式（4.4）は，ΛOη が線形方程式 L (g) = h
全て併合したものであることを示しているが，I(R = r) を R = r であれば 1，そうでなければ0 の値を取る 2 値の確率変数とし，1 = (1, . . . , 1)T ∈ Rm , また P (R = 1|Z) > 0 と仮定すると

g(Z o , R) =I(R = 1)h(Z)P (R = 1|Z)

114

統計数理

第 62 巻

第 1 号 2014

がその線形方程式の 1 つの解（特殊解）であることから
(4.5)

ΛO
η =

I(R = 1)ΛF
+ Ker（L∗ )P (R = 1|Z)が成り立つ．ここで，Ker（L∗ ) は線形作用素 L∗ の核，つまり Ker（L∗ ) = {h ∈ HO |L∗ (h) = 0}であるが，L∗ の定義と観測指示変数ベクトル R が離散的である（有限個の値しか取らない）
ことから，Ker（L∗ ) の要素は一般に次のように書ける．

I(R = 1)P (R = r|Z)(4.6)
k(Z o , R) =hr (Z or ).
I(R = r) −P (R = 1|Z)r =1


但し， r=1 は 1 以外の全ての r の値に関して和を取ることを意味し，hr (Z or ) は Z or = O(Z, r)のある実数値関数を表す．式（4.6）の表現は一般にはこれ以上簡単にならないが，単純な場合として，完全データ Z が Z = (V T , W T )T のように予め 2 つの部分に分かれ，V は常に観測され（全ての成分が）欠測するかのどちらかであるという場合（2 水準の欠測）
るが W は観測されるかは以下のように表される．
について考えてみると，
（4.6）

(4.7)

k(Z o , R) =R − e(V )h(V ).
e(V )但し，R は Z（の全ての成分）が観測されれば 1, そうでなければ（つまり W が欠測すれば）0の値を取る 2 値の観測指示変数であり，e(V ) = P (R = 1|V ), また h(V ) は V のある実数値関数
である．
4.3 観測データ局外接空間とその直交補空間 2 — データ観測確率が未知の場合 —さて，次にデータの観測確率（4.2）
が未知の場合について考える．このときは，（4.2）を観測データから推定する必要があるが，そのためにそれをパラメトリックにモデル化し，一般にP (Ri = r|O(Z i , r) = z o ) = ω(z o , r; α) (α ∈ A ⊂ Rq )と表すことにすると，観測データモデル(4.8)

S O = {q(z o , r; β, α, η) | β ∈ B ⊂ Rp , α ∈ A ⊂ Rq , η ∈ H}

(4.9)

但し，q(z o , r; β, α, η) = po (z o ; β, η)ω(z o , r; α)の局外接空間とその直交補空間は以下のように与えられる．
命題 4.（観測データ局外接空間とその直交補空間— データ観測確率が未知の場合 —）観測デー
局外接空間を T N O , 観測データヒルベルト空間 HO
タモデル S O の（点 q(z o , r; β, α, η) における）
NO
O
における Tの直交補空間を Λ とすると，(4.10)T N O = TηN O ⊕ TαN O , TηN O ⊥ TαN O(4.11)

ΛO = {h − Π(h|TαN O ) | h ∈ ΛO
η }が成り立つ．ここで，TαN O は S O の（有限次元の）
局外パラメータ α に関する接空間，すなわちTαN O = {λT sα (z o , r, α) | λ ∈ Rq }但し，sα (z o , r, α) =∂log ω(z o , r; α)∂αであり，また Π h|TαN O は，ヒルベルト空間 HO における，h の TαN O への直交射影を表す．

欠測データに対するセミパラメトリックな解析法

115

ちなみに，局外接空間 T N O が（4.10）のように直交分解されることは，観測データの密度関数が（4.9）のように分解されることからの帰結であり，（4.11）もその直交分解から自然に導かれる．
4.4 観測データ影響関数とその最適性
一般に興味あるパラメータの RAL 推定量の影響関数は，局外接空間の直交補空間に属する（3.10）
も満たす必要があるが，そのことも考慮すると（4.5），（4.11）から，ことに加えて命題 1 の（観測データ影響関数）
は，一般に観測データに基づくパラメータ β の RAL 推定量の影響関数
次のように表されることが分かる．
(4.12)oOoNO
φO (Z o , R, β, α, η) = φO (Z , R, β, α, η) − Π{φ  (Z , R, β, α, η)|Tα }

ここで，
o
φO
 (Z , R, β, α, η) =I(R = 1)φF (Z, β, η)+ k(Z o , R, α)ω(Z, 1; α)であり，φF (Z, β, η) は完全データに基づく β の RAL 推定量の影響関数（但し，命題 1 の（3.9），k(Z o , R, α) はその各成分が Ker（L∗ ) に属する Rp 値と（3.10）を満たすという意味において）の右辺第 2 項は，成分ごとに直交射影を取るものとする．
関数である．また（4.12）
さて，RAL 推定量の漸近分散共分散行列はその影響関数の分散共分散行列によって与えられたので，（4.12）において任意性のある完全データ影響関数 φF (Z, β, η) と関数 k(Z o , R, α) を変化させたときに，観測データ影響関数 φO (Z o , R, β, α, η) の分散共分散行列が最も小さくなる場
合が最適な場合となる（つまりそのときに，（4.12）は有効影響関数となる）．しかしながら，その意味で最適な場合は一般に複雑となるので，まず，完全データ影響関数を 1 つ固定したときの最適性について考える．これについて，以下のことが成り立つ．
命題 5.（最適な観測データ影響関数 — 完全データ影響関数を固定した場合 —）完全データ影（4.12）によって与えられる観測データ影響関数の分散響関数 φF (Z, β, η) を 1 つ固定したとき，共分散行列が最小となるのは
I(R = 1)φF (Z, β, η) ∗k(Z o , R, α) = −Π
 Ker（L )ω(Z, 1; α)の場合である．このとき，その観測データ影響関数は(4.13)φO (Z o , R, β, α, η) =I(R = 1)φF (Z, β, η)−Πω(Z, 1; α)I(R = 1)φF (Z, β, η) ∗ Ker（L )ω(Z, 1; α)となる．
式（4.13）の右辺第 2 項は一般には複雑で明示的には表せないが，特別な場合，例えば 4.2 節での場合には，Ker（L∗ ) の要素が（4.7）のように簡単な形になることから，触れた「2 水準の欠測」（4.13）は以下のように表される．
(4.14)φO (Z o , R, β, α, η) =RφF (Z, β, η)R − e(V ; α)−E{φF (Z, β, η)|V }.
e(V ; α)e(V ; α)

ここで，Z = (V T , W T )T であり，V は常に観測され，W は観測されるか（全ての要素が）欠測するかのどちらかである．また，R は Z（の全ての成分）が観測されるかどうかを表す 2 値の指（ロジスティック
示変数で，e(V ; α) は傾向スコア P（R = 1|V ) に対するパラメトリックモデル
回帰モデル等）
である．

116

統計数理

第 62 巻

第 1 号 2014

4.5 二重頑健推定量 — 2 水準の欠測の場合 —完全データ影響関数を固定するということは，実際には例えば，完全データに基づく推定関
数（完全データ推定関数）
を 1 つ与えることに対応する．すなわち，興味あるパラメータ β に対する完全データ推定関数 uF (Z, β) が与えられると，（推定方程式を通して得られる RAL 推定
量の影響関数として）完全データ影響関数
  F
−1
∂u
(4.15)φFuF (Z, β)(Z, β)u (Z, β, η) = E∂βが定まるが，このとき（4.14）の代わりに(4.16)

uO (Z o , R, β) =RuF (Z, β)R − e(V ; α)−E{uF (Z, β)|V }e(V ; α)e(V ; α)によって β に対する観測データ推定関数を与えると，これから得られる β の RAL 推定量の影
響関数は(4.17)oφOu (Z , R, β, α, η) =RφF
R − e(V ; α)u (Z, β, η)−E{φFu (Z, β, η)|V }e(V ; α)e(V ; α)となる． 但しこれは（4.16）
において，（E{uF (Z, β)|V } の部分の計算に必要な）V が与えられた下での W の条件付き分布とパラメータ α が既知であると仮定した場合のものである．しかしながら実際にはこれらは未知なので，推定関数（4.16）を用いるためにはこれらを観測データか（密度関数）にら推定する必要があるが，そのために V が与えられた下での W の条件付き分布
対しても，何らかのパラメトリックモデル p(w|v; ξ) (ξ ∈ Ξ ⊂ Rd ) を想定する． ここで，W は常に観測されるとは限らないが，その欠測メカニズムは MAR であると仮定しているので，Wと R は V が与えられた下で条件付き独立であり，その結果，V が与えられた下での W の条
件付き分布は V と R (R = 1) が与えられた下での条件付き分布と等しくなる．すなわち，パラメータ ξ は V と W がともに観測された場合（complete-case）のデータのみから推定可能である． このとき，観測データによる 2 つのパラメータ α, ξ の最尤推定量をそれぞれ α̂, ξ̂ とすると，推定方程式n Ri − e(V i ; α̂)Ri uF (Z i , β)(4.18)−
E{uF (Z i , β)|V i ; ξ̂} = 0e(Ve(Vi ; α̂)i ; α̂)i=1
(但し E{uF (Z i , β)|V i ; ξ̂} = uF (V i , w, β)p(w|V i ; ξ̂)dw)の解 β̂ DR について，以下のことが成り立つ．
命題 6.（二重頑健推定量 — 2 水準の欠測の場合 —）（1） 傾向スコア P（R = 1|V ) に対するモデル e(V ; α) と，V が与えられた下での W の条件付き分布に対するモデル p(w|v; ξ) のうち，少なくともどちらか一方が正しく特定されていれば，β̂ DR は β に対する RAL 推定量となる．
（2） どちらのモデルも正しく特定されていれば，β に対する RAL 推定量 β̂ DR の影響関数は（4.17）で与えられる．
の左辺第 1 項だけに着目すると，（完全データ推定関数 uF (Z, β) に基
さて，推定方程式（4.18）づく）β の IPW 推定方程式(4.19)n
Ri uF (Z i , β)i=1e(V i ; α̂)=0欠測データに対するセミパラメトリックな解析法

117が得られるが，この解として定まる IPW 推定量 β̂ IP W の影響関数は，（傾向スコアに対するモデルが正しく特定されていれば）観測データ影響関数の一般形（4.12）においてoφF (Z, β, η) = φFu (Z, β, η), k(Z , R, α) = 0としたものになる． 但し，φF（4.15）で与えられる，完全データ推定関数 uF (Z, β)u (Z, β, η) はに対応する影響関数である． よって，命題 5 と命題 6 の（2）より，傾向スコアに対するモデルと V が与えられた下での W の条件付き分布に対するモデルのどちらも正しく特定されていれ（少なくとも）IPW 推定量 β̂ IP W の漸近分散共分
ば，DR 推定量 β̂ DR の漸近分散共分散行列は散行列よりも小さくなる．
第 3 節でも述べたように，パラメータ β に対する（漸近的な）信頼領域などを得るためには，
DR 推定量 β̂ DR の漸近分散共分散行列とその推定量が必要となるが，推定方程式（4.18）において，パラメータ α と ξ の推定の影響を考慮しながら DR 推定量 β̂ DR の漸近分散共分散行列
Avar(β̂ DR ) を計算し，（それに対して）サンドイッチ推定量（3.5）を適用すれば，傾向スコアに対するモデルと V が与えられた下での W の条件付き分布に対するモデルのうち，少なくともど
ちらか一方が正しく特定されているという状況の下で，Avar(β̂ DR ) の一致推定量が得られる．
の第 10 章を参照されたい．
それについて詳しくは，Tsiatis（2006）
2.4 節において，最も単純な 2 水準の欠測問題の設定で IPW 推定量と DR 推定量を紹介したが，それらはそれぞれ完全データ推定関数 uF (Y, μ) = Y − μ に対する IPW 推定方程式，DR 推
定方程式の解として得られる推定量となっており，そこで述べた DR 推定量の（1）二重頑健性と（2）漸近有効性は命題 6 の特別な場合に当たる． 但し，その問題設定では，MAR と傾向スコア（つまモデルのみが仮定され，完全データモデル S F に対しては何も特別な制約を加えていない
り完全データに対してはノンパラメトリックモデルが想定されている）ので，興味あるパラメー
タ μ に対する完全データに基づく RAL 推定量の影響関数は唯一つしか存在しない．したがっは有効影響関数（漸近分散が最小となる RAL 推定量て，この場合においては，影響関数（4.17）（2）
は，DR 推定量が（2 つのモデルが正しく特定されていれば）
セの影響関数）であり，命題 6 のミパラメトリック漸近有効な推定量であることを意味する． しかしながら一般には，影響関数（4.17）は有効影響関数であるとは限らず，DR 推定量も（2 つのモデルが正しく特定されていたとしても）セミパラメトリック漸近有効な推定量になるとは限らない．次節の例は，そのような場合についてのものである．
4.6 ロジスティック回帰モデルにおける共変量欠測
2 水準の欠測の例として，ここではロジスティック回帰モデルにおける共変量の欠測問題について考え，それを通して二重頑健推定量の構成法をもう少し具体的に見てみることにする．
Y を（1 または 0 の値を取る）2 値の結果変数とし，共変量（説明変数）については簡単のために実数値を取る 2 つの連続変数 X1 , X2 を考え，これらについてロジスティック回帰モデル(4.20)

logit{P (Y = 1|X1 , X2 )} = β0 + β1 X1 + β2 X2が成り立っているものとする． 但し，Y と X1 は常に観測されるが，X2 については観測されるとは限らないという場合を想定し，X2 の観測指示変数を R として（つまり X2 が観測されれば R = 1, 欠測すれば R = 0），その欠測メカニズムは MAR であると仮定する． 完全データ, Yn , X1n , X2n ) とすると，これらがもし全（Y, X1 , X2 ）のランダムサンプルを（Y1 , X11 , X21 ), . . . （における回帰パラメータ β = (β0 , β1 , β2 )T の最尤推定量は推定方て観測されたならば，（4.20）
程式（スコア方程式）

118(4.21)

統計数理n


第 62 巻

第 1 号 2014

{Yi − eY (X i ; β)}X i = 0

i=1の解として与えられる． 但し，X i = (X1i , X2i )T , X i = (1, X1i , X2i )T , eY (X i ; β) = 1/{1 +exp（ − β T X i )} (i = 1, . . . , n) である． しかしながら，データ (Yi , X1i , X2i ) のうち X2i は n 個
体全てで観測されるとは限らないので，この（完全データ）推定方程式をベースとして，観測デー（4.18）を構成すると，それは以下のように与えられる．
タに基づく DR 推定方程式
n Ri uF (Yi , X i ; β)Ri − eR (Yi , X1i ; α̂)(4.22)−
E{uF (Yi , X i ; β)|Yi , X1i ; ξ̂} = 0ee
R (Yi , X1i ; α̂)R (Yi , X1i ; α̂)i=1
ここで，uF (Y, X ; β) = {Y − eY (X ; β)}X , eR (Y, X1 ; α) は傾向スコア P (R = 1|Y, X1 ) に対するパラメトリックモデル（ロジスティック回帰モデル等）
で，α̂ はそのモデルにおけるパラメータα の最尤推定量である．また ξ̂ は，Y と X1 が与えられた下での X2 の条件付き分布（密度関数）
に対するあるパラメトリックモデル p(x2 |y, x1 ; ξ) におけるパラメータ ξ の最尤推定量であり，
(4.23)
E{uF (Yi , X i ; β)|Yi , X1i ; ξ̂} = uF (Yi , X1i , x2 ; β)p(x2 |Yi , X1i ; ξ̂)dx2
である．但しここで重要なのは，パラメトリックモデル p(x2 |y, x1 ; ξ) は，完全データに対して（元々）仮定されているロジスティック回帰モデル（4.20）と整合するように設定されなければならないということである． 例えば，Y = 0, Y = 1 の下での X の条件付き分布が，それぞれ平均ベクトルは異なるが，分散共分散行列は共通の二変量正規分布であると想定すると，ロジス
ティック回帰モデル（4.20）
が成立し，またこのとき，Y と X1 が与えられた下での X2 の条件付
き分布（密度関数）は以下のように与えられる．
(4.24)p(x2 |y, x1 ; ξ) = n(x2 ; ξ0 + ξ1 x1 + ξ2 y, ξσ2 ) (ξ = (ξ0 , ξ1 , ξ2 , ξσ2 )T ).

但し，右辺は平均 ξ0 + ξ1 x1 + ξ2 y, 分散 ξσ2 の正規分布の密度関数を表すものとする． つまり，（密度関数）
に対する
二重頑健推定のための，Y と X1 が与えられた下での X2 の条件付き分布のように設定しておけば，特に「Y = 0, Y = 1 の下での X「作業モデル」
として，例えば（4.24）の条件付き分布が，それぞれ平均ベクトルは異なるが，分散共分散行列は共通の二変量正規分
布である」という場合には，元々のロジスティック回帰モデル（4.20）と両立し，矛盾がない． モ
デル（4.24）のパラメータ ξ は，MAR の下では Complete-Case（Yi , X1i , X2i の全てが観測されている個体 i）のデータのみを用いて最尤法で推定されるが，今の場合は正規分布を仮定しているので，標準的な線型回帰の方法（最小二乗法）によって推定できる．一方，（4.23）の積分計算は簡単ではないが，数値積分やモンテカルロ法などにより実行可能である．
は最尤推定量を導くスコア方程式であり，完全データに対してはモ
さて，推定方程式（4.21）
デル（4.20）の下で最適なものであるので，それをベースとして観測データに対する DR 推定方（推定方程式の構
程式を構成するのは自然に思われるが，その解として得られる DR 推定量は，成に必要な 2 つのモデルが正しく特定されていたとしても）最適，つまりセミパラメトリック漸（4.22）はあくまで，完全データ推定関数を
近有効になるとは一般には限らない．DR 推定方程式
uF (Y, X ; β) に固定したときに最適な（漸近分散共分散行列が最小の）推定量を与えるものであり，それがセミパラメトリック漸近有効になる保証はなく，実際，そのようにはなっていない．
このことはつまり，漸近有効な推定量を与えるような完全データ推定関数（あるいは影響関数）が別に存在する可能性を示唆しているが，（この場合も含めて）一般にそれは複雑なものとなる．
この問題については，ここではこれ以上述べないが，詳しくは Tsiatis（2006）の第 11 章以降を参照されたい．

欠測データに対するセミパラメトリックな解析法

5.

119

おわりに

」という特別な欠測パターンの下で二重頑健推
前節では，「2 水準の欠測（two-level missingness）「単調な欠測（monotone missing定量を導いたが，応用上もう 1 つ重要な欠測のパターンとして，」と呼ばれるものがある． これは例えば経時データ解析において，一度脱落（dropout）した
ness）
対象者は二度と対象とならない（つまり，脱落した時点の直前までしかデータが観測されず，以後の時点では全て欠測する）というような欠測パターンのことであるが，この場合も，観測デーの右辺第 2 項（Ker（L∗ ) への直交射影の部分）を明示的に表すことが可能であタ影響関数（4.13）
り，その表現に基づいて二重頑健推定量が導出される．
（完全データ推定方程式として一般化推を用いる場合が，その典型例である．
）しかしながら，欠測のパターンが非単調
定方程式
（GEE）
な場合
（一般の場合）には，原理的には二重頑健推定量を考えることはできるが，既に述べたようの右辺第 2 項は一般には明示的には表せず計算も単純ではなく，に，観測データ影響関数（4.13）
またデータの全ての要素が観測される確率（Complete-Case の確率）ω(Z, 1; α) = P (R = 1|Z) の推定も一筋縄には行かないため，応用上はまだ困難な問題が残っている．
「二重頑健推定量」という名称は，その推定量を導く推定方程式を構成するのに必要な 2 つのパラメトリックモデルのうち，どちらか一方が正しく特定されていればその推定量は一致性を持つという性質（二重頑健性）に由来するが，本稿で扱ったセミパラメトリック推測理論の立場から見た場合は，その性質はどちらかと言えば副次的なものであり，理論ではむしろ推定量の漸近的な推定精度の向上（漸近分散を可能な限り小さくするということ）の方に焦点が置かれている． 欠測データ解析法におけるこの「二重頑健性」
については他の立場からも考察がなされて，Kang and Schafer（2007）等を参照さいるが，それについては例えば，Bang and Robins（2005）
れたい． また，二重頑健性という性質自身は一般的なものであり，欠測データ解析法に限らず他の場面にも現れ得るものである．Robins and Rotnitzky（2001）では，そのための条件について，セミパラメトリック推測理論の立場から一般的な考察がなされている．
二重頑健推定量は，二重頑健性と漸近有効性（あるいは漸近的な推定精度の向上）の双方の面
よりは優れていると言えるが，その構成に必要な 2から，逆確率重み付け推定量（IPW 推定量）
つのパラメトリックモデルのどちらも誤特定している場合には，必ずしも良い挙動を示すとは言えない．そこで近年，それらのモデルの誤特定に対してより頑健性の強い「多重頑健推定量」
がいくつか提案されている（Han and Wang, 2013; Chan, 2013; Hattori and Henmi, 2014）． これは，推定対象となる 2 つの要素は変わらないが，その推定のためのモデルとしてそれぞれ複数の候補を許容し，それらのうちのどれか 1 つが正しく特定されていれば，興味あるパラメー
タの一致推定量が得られるという方法である． また，モデルの誤特定の問題以外にも実用上いくつかの問題点が指摘されており，それらに対処するための改良版二重頑健推定量もいくつか提案されている（Kang and Schafer, 2007; Cao et al., 2009; Tan, 2010; Rotnitzky et al., 2012）．
本稿では，欠測のメカニズムが MAR であるという仮定の下で議論を行ったが，2.1 節でも述べたように，その仮定は原理的に観測データからは検証できないものなので，データ取得のためのデザインなどから予め MAR であることが保証されているような場合を除いては，MAR からの乖離の影響を調べる感度解析などを行うことが望ましい．IPW 推定や DR 推定（二重頑健推，Li et al.（2011），
定）に対するこの種の感度解析法については，例えば Scharfstein et al.（1999）
Shen et al.（2011）
等がある．
IPW 推定量や DR 推定量は，（直接的な）
欠測の問題だけでなくしばしば因果推論の問題でも，あるいは反事実結扱われるが，これは，因果推論の問題が潜在的結果変数（potential outcome）と呼ばれるデータの欠測問題として見られることによる．この
果変数（counterfactual variable）
点を通して両者の問題は類似した構造を有しているが，例えば本稿で用いた「傾向スコア」
とい

120

統計数理

第 62 巻

第 1 号 2014

う用語は，もともとは因果推論の文脈で提案されたものである（Rosenbaum and Rubin, 1983）．
（特に後者）については，Tsiatis（2006）の第 13 章や因果推論における IPW 推定量や DR 推定量
Bang and Robins（2005）
等を参照されたい．
欠測データに対する二重頑健推定法は，まだ応用上の課題は残っているものの，完全データの分布全体に対するパラメトリックなモデル化を必要としないセミパラメトリックな解析法としていくつかの有用な性質を持ち，今後も有力な方法として発展していくと思われる．本稿ではセミパラメトリック推測の一般論の立場から，その理論的な側面を中心に解説を行ったが，応用上の観点からは例えば，ロジスティック回帰モデルの共変量欠測については Tchetgen
Tchetgen（2009），一般化推定方程式（GEE）
における経時データの欠測については Seaman and（2009）等が参考になる．
Copas

参

考

文

献

Bang, H. and Robins J. M.（2005）. Doubly robust estimation in missing data and causal inference
models, Biometrics, 61, 962–972.
Boos, D. D. and Stefanski, L. A.（2013）. Essential Statistical Inference, Springer, New York.
Cao, W., Tsiatis, A. A. and Davidian, M.（2009）. Improving eﬃciency and robustness of the doubly
robust estimator for a population mean with incomplete data, Biometrika, 96, 723–734.
Carpenter, J. and Kenward, M.（2013）. Multiple Imputation and Its Application, Wiley, Chichester.
Chan, K. C. G.（2013）. A simple multiply robust estimator for missing response problem, Stat, 2,
143–149.
Eekhout, I., de Boer, M. R., Twisk, J. W. R., de Vet, H. C. W. and Heymans, M. W.（2012）. Missingdata: A systematic review of how they are reported and handled, Epidemiology, 23, 729–732.
. Estimation with missing data: beyond double robustness, Biometrika,
Han, P. and Wang, L.（2013）
100, 417–430.
Hattori, S. and Henmi, M.（2014）. Stratiﬁed doubly robust estimators for the average causal eﬀect,
Biometrics（to appear）.
Horvitz, D. G. and Thompson, D. J.（1952）. A generalization of sampling without replacement from aﬁnite universe, Journal of the American Statistical Association, 47, 663–685.
Kang, J. D. Y. and Schafer, J. L.（2007）. Demystifying double robustness: A comparison of alternative strategies for estimating a population mean from incomplete data, Statistical Science, 22,
523–539.
Li, L., Shen, C., Wu, A. C. and Li, X.（2011）. Propensity score-based sensitivity analysis method for
uncontrolled confounding, American Journal of Epidemiology, 174, 345–353.
. Statistical Analysis with Missing Data, 2nd ed., Wiley, New
Little, R. J. A. and Rubin, D. B.（2002）
Jersey.
Robins, J. M. and Rotnitzky, A.（2001）. Comment on the Bickel and Kwon article, “On double robustness,” Statistica Sinica, 11, 920–936.
Robins, J. M., Rotnitzky, A. and Zhao, L. P.（1994）. Estimation of regression coeﬃcients when some
regressors are not always observed, Journal of the American Statistical Association, 89, 846–
866.
Rosenbaum, P. R. and Rubin, D. B.（1983）. The central role of the propensity score in observational
studies for causal eﬀects, Biometrika, 70, 41–55.
Rotnitzky, A., Lei, Q., Sued, M. and Robins, J. M.（2012）. Improved double-robust estimation in missing
data and causal inference models, Biometrika, 99, 439–456.
Rubin, D. B.（1976）. Inference and missing data, Biometrika, 63, 581–592.

欠測データに対するセミパラメトリックな解析法

121

Rubin, D. B.（1987）. Multiple Imputation for Nonresponse in Surveys, Wiley, New York.
Scharfstein, D. O., Rotnitzky, A. and Robins, J. M.（1999）. Adjusting nonignorable drop-out using
semiparametric nonresponse models, Journal of the American Statistical Association, 94, 1096–
1120.
Seaman, S. and Copas, A.（2009）. Doubly robust generalized estimating equations for longitudinal data,
Statistics in Medicine, 28, 937–955.
Seaman, S. R. and White I. R.（2011）. Review of inverse probability weighting for dealing with missing
data, Statistical Methods in Medical Research, 22, 278–295.
Shen, C., Li, X., Li, L. and Were, M. C.（2011）. Sensitivity analysis for causal inference using inverse
probability weighting, Biometrical Journal, 53, 822–837.
Tan, Z.（2010）. Bounded, eﬃcient and doubly robust estimation with inverse weighting, Biometrika,
97, 661–682.
Tchetgen Tchetgen, E. J.（2009）. A simple implementation of doubly robust estimation in logistic
regression with covariates missing at random, Epidemiology, 20, 391–394.
Tsiatis, A. A.（2006）. Semiparametric Theory and Missing Data, Springer, New York.
van der Vaart, A. W.（2000）. Asymptotic Statistics, Cambridge University Press, Cambridge.

122

Proceedings of the Institute of Statistical Mathematics Vol. 62, No. 1, 103–122 (2014)

Semiparametric Statistical Methods for Missing Data
Masayuki Henmi
The Institute of Statistical Mathematics

In observational studies such as cohort studies in epidemiology, results of statisticalanalysis can be biased for various reasons such as missing data. There are two kinds ofstatistical methods for missing data. One is a parametric method such as the maximum
likelihood method and multiple imputation. The other is a semiparametric method such as
inverse probability weighting and doubly robust estimation. This paper focuses on the latter, which has been rapidly developed in resent years especially in biostatistics. Although
there have been few applications of the semiparametric method in practice, it is expected
as a powerful method for missing data and its demand should increase in the future. Semiparametric statistical methods are widely used for various problems in biostatistics. The
argument in this paper is based on general theory of semiparametric inference in order tohelp our theoretical understanding of other cases as well as missing data problems.

Key words: MAR, influence functions, estimating functions, nuisance tangent spaces, asymptotic efficiency.