ベイズデータ解析３

概要

標本調査において欠測はつきものである．観測単位が欠測している場合 (unit nonresponse)，call-back や follow-up などの調査を行うか，それができない場合は 荷重校正 (calibration weighting) が可能である．一方で，項目が欠測している場合 (item nonresponse)，代入法 (imputation) が用いられる．

1 有限標本論の概要

1.1 設定

\([N]\) を母集団とする．

\([N]\) の部分集合の全体 \(P([N])\) 上の確率分布を 抽出計画 (sampling design) といい，ある既知の抽出分布に従って得られる標本 \(S\subset[N]\) を 確率標本 (probability sample) という．日本語では 無作為抽出標本 などとも呼ばれる．

抽出計画の例

抽出計画 \(\mathcal{L}[S]\) には，

単純無作為抽出 (SRS: Simple Random Sampling)
系統無作為抽出 (Systematic Random Sampling)
層別抽出 (Stratified Sampling)：母集団を層別し，各層の間では独立な抽出を行う．層ごとに抽出計画は異なっても良い．条件付きランダム化 (conditional randomization) ともいう (Section 2.2 Hernán and Robins, 2020, p. 17)．
クラスター抽出 (Cluster Random Sampling)：クラスターがまずランダム抽出され，そのクラスター内の全構成員が標本に加わる．クラスターのことを PSU (Primary Sampling Unit) ともいう．クラスター内でもランダム抽出が行われた場合，２段階クラスターサンプリング という．

などの方法が存在する．

例えば，日本の国勢調査は２段階の層別抽出である．

確率標本 \(S\) では（１次の）包含確率 \[ \pi_i:=\operatorname{P}[i\in S]=\operatorname{P}[I=1],\qquad I:=1_{i\in S}. \] が定まる．

狭義の「確率標本」

先ほど，\(\mathcal{P}([N])\) 上の確率変数を確率標本と呼ぶとしたが，正確に \(S\) が 確率標本 と呼ばれるためには，\(\pi_i>0\) が母集団 \(i\in[N]\) の全域で成り立つことが必要である (Kim, 2024, p. 12)．

1.2 Horvitz-Thompson 推定量

確率標本 \(S\in\mathcal{L}(\Omega;\mathcal{P}([N]))\) に対しては，ある量 \(y\) についての母集団の総和 \[ Y:=\sum_{i=1}^Ny_i \] が \[ \widehat{Y}_\mathrm{HT}:=\sum_{i\in S}\frac{y_i}{\pi_i} \] により不偏推定できる．

\(\widehat{Y}_\mathrm{HT}\) は (Horvitz and Thompson, 1952) 推定量と呼ばれる．¹

(Sen, 1953)-(Yates and Grundy, 1953)

Horvitz-Thompson 推定量の分散は次で与えられる：

\[ \mathrm{V}[\widehat{Y}_\mathrm{HT}]=\sum_{i,j=1}^N\biggr(\pi_{ij}-\pi_i\pi_j\biggl)\frac{y_iy_j}{\pi_i\pi_j}. \] ただし， \[ \pi_{ij}:=\operatorname{P}[i\in S,j\in S]=\operatorname{P}[I=1=J]. \] を２次の包含確率という．

Horvitz-Thompson 推定量の要点には，「計画した欠損ならば，重みづけによって不偏推定量を得ることができる」という点にある．

そこで抽出計画が不明な場合もこれを推定し，バイアスを補正しようとするアプローチを傾向スコアの方法，または 擬似ランダム化 (pseudo-randomization) の方法という．

1.3 「校正」：効率の改善に向けて

HT 推定量は確率標本 \(S\) の分布，すなわち抽出計画に依らずに不偏性を持つ．

これを計画不偏性 (design-unbiasedness) というが，この性質を持つ線型な推定量は HT に限られる．

しかし，HT 推定量はいつでも分散が最小というわけではない．

計画不偏性は bias-variance trade-off の観点からは欠点でもあり，それゆえ抽出計画に関する情報を用いて分散を低減することも考えられる．

特に，HT 推定量の荷重 \((\pi_i^{-1})\) を，補助変数 \(x_i\) に関する 外部一致性 \[ \sum_{i\in S}w_ix_i=\overline{x} \] を保ちながら新しいもの \((w_i)\) に変更するものが多く考えられた．

「外部一致性」の別名

有限標本論は普遍的な統計推測の基礎であると言える．

実際，この外部一致性の条件は多くの分野で考慮されており，種々の名前が与えられている．

有限標本論：校正条件 (calibration condition / benchmarking property) （第 2.3 節）
欠測データ解析：共変量バランシング (covariate balancing) (Imai and Ratkovic, 2014)
機械学習（継続学習）：共変量シフト (covariate shift) (Shimodaira, 2000)

このアプローチを 荷重校正 (calibration weighting) という．

次章にてこれ以降，種々の荷重校正推定量を紹介する．

2 回帰推定量

2.1 はじめに

前述の通り，補助変数 \(x\) が母集団上で知られている場合に，ここから抽出計画に対する追加情報を抽出して推定量に組み込むことで，計画的に欠測させられたデータ（＝確率標本）に対する不偏推定量 (Horvitz and Thompson, 1952) （以降 HT 推定量という）の効率を改善することを考える．

以降，補助変数 \(x\in\mathbb{R}^p\) は母集団上で既知であるとし，その総和を \[ X:=\sum_{i\in[N]}x_i \] で表す．

2.2 比による校正

補助変数の次元が \(p=1\) のとき，最も安直には \(X\) の HT 推定量から，真の値 \(X\) との「ズレ方」を用いて，\(Y\) の推定量を「校正」することができる．

\[ \widehat{Y}_{\mathrm{R}}:=\widehat{Y}_\mathrm{HT}\frac{X}{\widehat{X}_\mathrm{HT}} \] とできるだろう．

この推定量は ratio estimator などと呼ばれ，性能の代わりにバイアスが生じてしまう．

一般に，\(X,Y\) が正の相関を持つとき大きな分散低減が得られる (Deng and Wu, 1987), (Kim, 2024, p. 92)．

\(x_i=1\) と取った場合を Hajék 推定量ともいう．Hajék 推定量が HT 推定量よりも推奨される状況が (Särndal et al., 1992, p. 182) にリストされている．

この推定量は昔は計算の簡単さから使われていたが，一般の次の回帰推定量の方が MSE が小さいことが知られている (Deng and Wu, 1987)．

2.3 回帰推定量

超母集団模型 \[ Y=X^\top\beta+\epsilon,\qquad\epsilon\overset{\text{iid}}{\sim}(0,\sigma^2) \] を想定し，得られている標本のみから \(\widehat{\beta}\) を推定する．こうして得られる \[ \widehat{y}_i:=x_i^\top\widehat{\beta},\qquad\widehat{\beta}:=\left(\sum_{i\in S}\pi_i^{-1}x_ix_i^\top\right)^{-1}\sum_{i\in S}\pi_i^{-1}x_iy_i \] の総和が，\(Y\) に対する 回帰推定量 (regression estimator) と呼ばれる．²

これは \((y_i)\in\mathbb{R}^n\) に関する線型推定量になっている．加えて，外部一致性 \[ \sum_{i\in S}w_ix_i=\overline{x} \tag{1}\] を満たす荷重 \[ w_i:=\overline{X}^\top\left(\sum_{i\in S}\pi_i^{-1}x_ix_i^\top\right)^{-1}\pi_i^{-1}x_i \] に関して， \[ \widehat{Y}_{\mathrm{reg}}=\sum_{i\in S}w_iy_i \] という形の線型推定量になっている．

式 (1) を 外部一致性 (external consistency)，または 校正条件 (calibration / benchmarking property) (Deville and Särndal, 1992) という．

回帰推定量は \(X,Y\) の関係に依らず一致性を持ち，\(X,Y\) の間の相関の絶対値が大きいほど分散低減効果が高くなる (Kim, 2024, p. 95)．³

2.4 事後層別化

事後層別化 (post-stratification / stratification after selection) は標本抽出の結果を見て標本を層別化する手法であるが，回帰推定量の特別な場合と見れる．

母集団が \(G\) 個の層に分けられるとする：\(N=N_1+\cdots+N_G\)．

このとき，\(i\in[N]\) 番目の単位が層 \(g\in[G]\) に属するかどうかの指示変数 \(x_{ig}\in2\) のベクトル \(x_i:=(x_{i1},\cdots,x_{iG})^\top\in2^G\) に関する回帰推定量 \[\begin{align*} \widehat{Y}_{\mathrm{post}}&:=\sum_{i=1}^Nx_i^\top\left(\sum_{i\in S}\pi_i^{-1}x_ix_i^\top\right)^{-1}\sum_{i\in S}\pi_i^{-1}x_iy_i\\ &=\sum_{g=1}^G\sum_{i\in S_g}\pi_i^{-1}\frac{N_g}{\widehat{N}_g}y_i,\qquad\widehat{N}_g:=\sum_{i\in S}\pi_i^{-1}x_{ig}. \end{align*}\] を事後層別化推定量という．

MRP (Multilevel Regression and Post-stratification) (Gelman and Little, 1997), (Gelman, 2014) は事後層別化の階層モデル・縮小推定版である．

2.5 ランキング法／繰り返し比例的フィッティング法

(Deming and Stephan, 1940) では 1940 年の国勢調査の結果の分析を考えていた．

特に，基本的な情報は全数調査されるが，詳細な情報は標本調査でしか得られない状況下で，母集団の \(I\times J\) 分割表の各セル \(U_{ij}\) の値 \(N_{ij}\) の推定を考えていた．

ただし，周辺和 \(N_{i-},N_{-j}\) は全数調査で得られているとする．

このとき，\(N_{ij}\) の推定量の候補として \[ \frac{n_{ij}}{n_{i-}}N_i,\quad\frac{n_{ij}}{n_{-j}}N_{-j},\quad\frac{n_{ij}}{n}N \] の３つが考えられる．３番目が良いと考えるかもしれないが，その結果得られる分割表は周辺和を保存しない．

この問題は次のような形でも現れる：指示変数 \[ x_k=(x_{1-k},\cdots,x_{I-k},x_{-1k},\cdots,x_{-Jk}),\qquad x_{ijk}:=1_{U_{ij}}(k), \] に基づく事後層別化推定量 \[ \widehat{Y}_{\mathrm{post}}=\sum_{i\in S}\pi_i^{-1}g_i(S)y_i,\qquad g_i(S):=\left(\sum_{k=1}^Nx_k\right)^\top\left(\sum_{k\in S}\pi_k^{-1}x_kx_k^\top\right)^{-1}x_i \] を考えたいが，これが \(\operatorname{rank}\left(\sum_{k\in S}\pi_k^{-1}x_kx_k^\top\right)=I+J-1\) であるため，一意な表示を持たない．

\(g_i(S)\) の候補のうち，次を満たす \(g_i\) を選ぶことが目標である： \[ \sum_{k\in S}\frac{g_k}{\pi_k}x_{i-k}=\sum_{k=1}^Nx_{i-k}=N_{i-}, \tag{2}\] \[ \sum_{k\in S}\frac{g_k}{\pi_k}x_{-jk}=\sum_{k=1}^Nx_{-jk}=N_{-j}. \tag{3}\]

(Iterative Proportional Fitting / Ranking algorithm Deming and Stephan, 1940)

\(g^{(0)}_k\gets1\) と初期化する．
\(x_{i-k}=1\) すなわち \(k\in U_{i-}\) であるとき， \[ g^{(t+1)}_k\gets g_k^{(t)}\frac{\sum_{k=1}^Nx_{i-k}}{\sum_{k\in S}\frac{g^{(t)}_k}{\pi_k}x_{i-k}}. \] これにより条件 (2) が満たされる．
\(z_{-jk}=1\) すなわち \(k\in U_{-j}\) であるとき， \[ g^{(t+2)}_k\gets g_k^{(t+1)}\frac{\sum_{k=1}^Nx_{-jk}}{\sum_{k\in S}\frac{g^{(t+1)}_k}{\pi_k}x_{-jk}}. \] これにより条件 (3) が満たされる．
収束するまで繰り返す．

これは特定の目的関数を最小化することに等しい．(Deming and Stephan, 1940, p. 428), (Zieschang, 1990), (Jean-Claude Deville and Sautory, 1993) も参照．

3 荷重校正推定量

3.1 はじめに

回帰推定量は \(X\) から \(Y\) に関する情報を抽出することで，HT 推定量の効率を改善することができる方法である．

しかし，HT のもう一つの魅力的な性質であった 計画一致性 (design consistency) が失われている．

回帰推定量の性質である 外部一致性 (external consistency) を保ちながら，別の解を見つけることで，回帰推定量を一般化する形で計画一致性を持つ効率的な推定量を構成することを考える．

実はこの方法は，モデリングの観点からは \(X,Y\) の間のモデルを，標本レベルから母集団レベルに一般化することに相当する．こうして考えられる超母集団モデルを 一般化回帰モデル (GREG: Generalized Regression) という．

このような方法で HT 推定量を改善した計画一致性を持つ推定量を model-assisted estimator，特に特定の制約下最適化問題の解として与えられるものを 校正推定量 (calibrated estimator) という．

校正推定量は計画一致性を持つために，傾向スコアの推定に成功していれば不偏性が保証される．この性質は二重頑健推定量の構成の基礎となる．

3.2 差分推定量

補助的な量 \(y_i^{(0)}\) が母集団全体で観測されている場合， \[ \widehat{Y}_{\mathrm{diff}}:=\sum_{i=1}^Ny_i^{(0)}+\sum_{i\in S}\pi_i^{-1}\left(y_i-y_i^{(0)}\right) \] は 差分推定量 (difference estimator) と呼ばれる．

HT 推定量同様不偏であるが，分散の値は変化し，特に \(y_i^{(0)}\) が \(y_i\) の良い近似であるほど分散が小さくなる (Kim, 2024, p. 99)．

この \(y_i\) の proxy とも言える量 \(y_i^{(0)}\) を，他の共変量 \(x_i\) から回帰により構成することで，回帰推定量（第 2.3 節）よりも複雑な \(x_i,y_i\) 関係もうまく取り込んだ分散低減が可能になる．

このように（暗黙裡にでも）モデルを用いており，加えてモデルの特定が成功しているかに依らず HT 推定量を改善できる方法を model-assisted estimation といい，校正推定量の基本的な考え方である．

3.3 一般化回帰モデルと射影推定量

まず母集団 \([N]\) に応用 \(Y\) のモデルを当てはめる： \[ y_i=x_i^\top\beta+\epsilon_i,\qquad\epsilon_i\overset{\text{iid}}{\sim}(0,c_i(x_i)\sigma^2). \tag{4}\] このように母集団に置かれるモデルを 超母集団モデル (superpopulation model) (Isaki and Fuller, 1982) という．

特に式 (4) の Gauss-Markov 型の超母集団モデルを 一般化回帰モデル (GREG: Generalized Regression) ともいう．

これを解いて得る推定量 \(\widehat{y}_i=x_i^\top\widehat{\beta}_c\) の総和として得られる推定量 \[ \widehat{Y}_{\mathrm{P}}:=\sum_{i=1}^N\widehat{y}_i \] を（モデルベースの） 射影推定量 (projection estimator) という．

射影推定量は計画一致性を持つとは限らない．

仮に GREG モデルで \[ \frac{c_i}{\pi_i}\parallel x_i \] が成り立つならば，内部バイアス校正 (IBC: Internally Biased Calibration) (Firth and Bennett, 1998) 条件 \[ \sum_{i\in S}\frac{1}{\pi_i}(y_i-\widehat{y}_i)=0 \] が成り立つ．

この IBC が，射影推定量が抽出計画に依らずに一致性を持つための十分条件である (補題9.1 Kim, 2024, p. 100)．

3.4 一般化最小二乗法 (GLS)

当然 GREG モデルが IBC 条件を満たすとは限らない．

そのような場合でも計画一致性を持つような推定量を考えたい．実は， \[ \widehat{Y}_{\mathrm{GREG}}:=\widehat{Y}_\mathrm{HT}+\biggr(X-\widehat{X}_\mathrm{HT}\biggl)^\top\widehat{\beta}_c \] は計画一致性を持つ．

これは 一般化回帰推定量 (GREG: Generalized Regression Estimator) または計量経済学において GLS (Generalized Least Squares) (Aitken, 1936) と呼ばれる．⁴

一般化回帰推定量は次の最適化による特徴付けがある： \[ \widehat{Y}_{\mathrm{GREG}}=\sum_{i\in S}\widehat{\omega}_iy_i,\qquad\widehat{\omega}_i:=\pi_i^{-1}+\left(X-\widehat{X}_\mathrm{HT}\right)^\top\left(\sum_{i\in S}\frac{1}{c_i}x_ix_i^\top\right)^{-1}\frac{x_i}{c_i}. \] この荷重 \(\widehat{\omega}_i\) は，校正条件 (calibration constraint) （式 (1) との違いに注意）を満たすものの中で \[ Q(\omega):=\sum_{i\in S}(\omega_i-d_i)^2c_i,\qquad d_i:=\pi_i^{-1},\quad\operatorname{subject to}\sum_{i\in S}\omega_ix_i=\sum_{i=1}^Nx_i. \] を最小にするものとも特徴付けられる (Kim, 2024, p. 102)．

特に，\(\widehat{w}_i\xrightarrow[n\to\infty]{\mathrm{p}}d_i\)．

3.5 校正推定量

一般に，校正条件制約を満たす \((\omega_i)\) のうち，凸関数 \(G\) が定める目的関数 \[ Q(\omega):=\sum_{i\in S}d_iG\left(\frac{\omega_i}{d_i}\right)c_i \] を最小にするものを 校正荷重 (calibration weight)，校正荷重に関する線型推定量を 校正推定量 (calibration estimator) という (Deville and Särndal, 1992), (Kim, 2024, p. 103)．

ほとんどの校正推定量は漸近的に GREG 推定量に一致する．

一般に，有限母集団に対する確率標本からの一様最小分散不偏推定量 (UMVUE) は存在しない (Godambe and Joshi, 1965) が，GREG 推定量は「期待漸近分散」の下界を達成する (Isaki and Fuller, 1982)．

3.6 最適校正推定量

特に， \[ Q(\omega)=\sum_{i\in S}\omega_i^2c_i \] を最小化するものは 最適校正推定量 (optimal calibrated estimator) と呼ばれる (Kim, 2024, p. 110)．

これはモデルの視点からは \(x\) を拡張して人工的に IBC 条件を満たすようにした射影推定量（第 3.3 節）とも見れる．

最適校正推定量は超母集団モデル (4) が誤特定されている場合に GREG 推定量より良い性能を示す (Kim, 2024, p. 112)．

GREG モデルより一般的な超母集団モデルに対しての同様の手続きは モデル校正 (model calibration) (Wu and Sitter, 2001) と呼ばれている．この方法では \(X,Y\) の関係を推定し，\(Y\) の線型推定量を \(m(X)\) の形で構成してから，最適構成推定量の議論に還元する．

3.7 一般化エントロピー法

最適構成推定量の構成に倣い， \[ Q(\omega):=\sum_{i\in S}G(\omega_i)c_i\qquad\operatorname{subject to}\sum_{i\in S}\omega_ig(d_i)c_i=\sum_{i=1}^Ng(d_i)c_i \] の最小化により校正荷重を構成する方法を 一般化エントロピー法 (generalized entropy method) (Kwon et al., 2024) という．

これは目的関数には計画荷重 \(d_i=\pi_i^{-1}\) が入っていないが，制約条件に入っていることで計画一致性を達成している．

超母集団モデルである GREG モデルが正しく特定されているならば (Godambe and Joshi, 1965) の下界を達成するが，そうでなくとも一致性は保たれる上に，一般の校正推定量（第 3.5 節）よりも分散は小さいである (Kwon et al., 2024)．⁵

4 欠測データの扱い

4.1 はじめに

観測単位が欠測している場合 (unit nonresponse)，call-back / follow-up 調査を行うか，それができない場合は次の２つの対処が可能である：

単位欠測の扱い

欠測メカニズムを抑える共変量は見えている場合（MAR 条件），傾向スコア推定量が利用可能（第 4.2 節）．これは欠測メカニズムのモデリングに基づく．
一般の校正推定量に対しても，

単位欠測の場合は，２段階の標本抽出と状況が似ているのである．さらには，非確率標本（調査観察データ，ビッグデータなど）の扱いとも似通う．これについては次稿も参照．

一方で，項目が欠測している場合 (item nonresponse)，代入法 (imputation) が用いられる．⁶

現状は多重代入法（第 5.2 節）が主流であると言える (Buuren, 2018)．

4.2 傾向スコア推定量

標本の観測 \(Y_i\) は，\(\delta_i=0\) のとき欠損しているとする．

4.2.1 MAR 条件：欠測のメカニズムを抑える共変量が観測できている

加えて，標本全体についてある変数 \(X\) が観測できており，これについて次の条件が成り立つとする：

(MAR condition Rubin, 1976)⁷

欠測の指示変数 \(\delta\) について， \[ \operatorname{P}[\delta=1|X,Y]=\operatorname{P}[\delta=1|X]=:p(X) \] が成り立つ．

これは条件付き独立性 \(\delta\perp\!\!\!\perp Y\mid X\) よりも弱い条件で，MAR (Missing At Random) の条件と呼ばれる．⁸

4.2.2 欠測メカニズムの推定

欠測確率 \(p(x):=\operatorname{P}[\delta=1|X=x]\) にノンパラメトリックなモデル \(p_\phi(x)\) を課したとする．

このとき，パラメータ \(\phi\) は擬似最尤推定量 \(\widehat{\phi}\) により一致推定をすることができる．

4.2.3 傾向スコア推定量

仮に母平均 \[ Y:=\sum_{i=1}^Ny_i \] が推定対象であったとしよう．

このとき，推定された \(\widehat{\phi}\) を元に，次の推定量が構成できる：

\[ \widehat{Y}_\mathrm{PS}:=\sum_{i\in\delta^{-1}(1)}\frac{1}{\pi_i}\frac{y_i}{p_{\widehat{\phi}}(x_i)}. \]

命題（傾向スコア推定量の一致性）⁹

欠測確率 \(p\) のモデル \(p_\phi(x)\) の特定に成功しているとき，ある正則性に関する条件が満たされる限り，傾向スコア推定量 \(\widehat{Y}_\mathrm{PS}\) は一致推定量に \(n^{-1}\) のオーダーで漸近する．

4.3 校正推定量

ある校正荷重 \((d_i)\) に関して，計画一致性を持つ推定量 \[ \widehat{Y}=\sum_{i\in S}d_iy_i \] を考えているが，単位欠測により特定の \(y_i\) が得られず，計算できないものとする．

この場合でも，応答があった部分標本 \[ S_R:=\delta^{-1}(1) \] 上の校正推定量 \[ \widehat{Y}_\omega:=\sum_{i\in S_R}\omega_iy_i \] であって，欠測メカニズム \(p(x)\) の特定か，または超母集団モデル \[ y_i=x_i^\top\beta+\epsilon_i,\qquad\epsilon_i\overset{\text{iid}}{\sim}(0,c_i\sigma^2) \] の特定に成功すれば一致性を持つ，二重頑健なものを構成できる (Kim and Haziza, 2014)．

4.4 代入法とその不偏性条件

項目非反応がある場合，代入値を \(y_i^*\) として \[ \widehat{Y}_{\mathrm{I}}:=\sum_{i\in S}\frac{1}{\pi_i}\biggr(\delta_iy_i+(1-\delta_i)y_i^*\biggl) \] による推定が試みられる．

代入 \(y_i^*\) を行うことでリストワイズの削除をするよりも推定の効率を上げることができる．

(代入推定量の不偏性 Kim, 2024, p. 162)

\[ \operatorname{E}[Y^*|\delta=0]=\operatorname{E}[Y|\delta=1] \] が成り立つならば，\(\widehat{Y}_\mathrm{I}\) は不偏推定量である．

この条件は，標本内で MAR 条件（第 4.2.1 節）が成り立つとき： \[ Y|(X,\delta=1)=Y|(X,\delta=0), \tag{5}\] \(Y^*\) を \(Y|(X,\delta)\) からのサンプリングで代入すれば達成される．

さらに強い条件 \[ \delta\perp\!\!\!\perp Y\mid X \] が成り立つとき，標本内の MAR 条件が成り立つ．

換言すれば代入法において，欠測の原因 \(X\) を突き止め，欠測したグループにおける \(Y\) の値 \(Y|(X,\delta=1)\) にモデル (outcome model) を立て，そこからサンプリングをすることを目指す．

4.5 回帰による代入

仮に共変量 \(X\) が \(Y\) と強い相関を持つとする．このように線型回帰模型を背後に想定することが適切な場合は，よく次のような手続きで代入がされる．

まず共変量により母集団を \([N]=N_1+\cdots+N_G\) 個に層別化し，それぞれの層で \[ Y_i=X_i^\top\beta+\epsilon_i,\qquad\epsilon_i\overset{\text{iid}}{\sim}(0,\sigma^2) \tag{6}\] というセミパラメトリック回帰モデルを考える．

次に推定されたモデルを用いて，\(\epsilon_i^*\sim(0,\sigma^2)\) を残差 \[ \widehat{\epsilon}_i:=y_i-x_i^\top\widehat{\beta} \] の分布から（リ）サンプリングし， \[ y_i^*\gets x_i^\top\widehat{\beta}+\epsilon_i^* \] を代入値とする．

以上の手続きは 確率的回帰代入法 (stochastic regression imputation) と呼ばれる．平均を代入する場合は単に回帰代入法または条件付き平均代入法 (conditional mean imputation) という．

\(Y\) と強い相関を持つ補助変数 \(X\) がいつでも見つかるとは限らない．

その場合は Gauss-Markov モデル (6) を一般の統計モデルに一般化すれば良い．

4.6 マッチングによる代入

層の中の他のセルをランダムに選んでその値を代入する hot deck imputation や，セルの加重平均を代入する fractional hot deck も同様の考え方に基づく (Fuller and Kim, 2005)．

このような手法は マッチング と呼ばれ，カーネル法と関連が深い (Cheng, 1994)．加重平均は対象のセルとの関連度を「距離」によって測り，距離を計算するのに使われる変数はキーともいう (高井啓二 et al., 2016, p. 110)．傾向スコアマッチングでは傾向スコアがキーである．

最も単純には同一データセット内の最も似ている単位を持ち出してその値を代入するのがマッチングであるが，最も洗練された方法としては類似度に依存して関連度を自動的に重みづけて，データセット全体で加重平均をとっても良いわけである．

他の標本の値を参考にする場合は cold deck imputation という．

などの Least squares method も同様の考え方に基づく (Little, 1992)．

4.7 母集団モデルによる代入法

一方で，母集団上での \(Y,X\) の関係についてモデルを立てて \(Y|X\) からサンプリングをすることも考えられる．

注（無情報サンプリング条件）

母集団の分布と標本の分布が一致するとき，無情報サンプリング (noninformative sampling) が実施されたという．そうでない場合は informative sampling という．

サンプリングが無情報であるための十分条件には \[ \operatorname{P}[I=1|X,Y]=\operatorname{P}[I=1|X] \] が挙げられる．(Sugden and Smith, 1984) はこれを無情報サンプリング条件という．

この下では母集団のモデルと標本のモデルとは一致するが，一般にはこの２つは厳密に峻別しなければ混乱の源である．

標本内の MAR 条件 (5) だけでなく，母集団上で MAR 条件が成り立つ場合は，\(Y|X\) の尤度を \(f_\theta(y|x)\) としてモデリングをし，これを \[ \ell(\theta):=\sum_{i\in S}w_i\delta_i\log f_\theta(y_i|x_i) \] の最大化によって \(M\)-推定することが考えられる．¹⁰

ただし，\(w_i\) は \(Y\) の計画一致性を持つ校正推定量を定める校正荷重であるとする．\(w_i\) の存在は標本と母集団のズレに起因する．

最終的に学習されたモデル \(f_\theta(y|x_i)\) からのサンプリングによって代入値 \(y_i^*\) を生成する．

このモデル \(f_\theta(y|x_i)\) を当てはまりの度合いを見ながらベイズ推論によって得る方法もよく取られるようになっている (C. K. Enders et al., 2020)．

母集団上の MAR 条件が成り立たない場合は \(Y|(X,\delta=0)\) のモデリングを考える必要がある．

5 多重代入法

5.1 はじめに

ベイズの観点からは，欠測データとパラメータとは違いがない (Chapter 18 Gelman et al., 2014, p. 449)．

ベイズ事後分布は欠測データとパラメータの上に同時に定まり，欠測データに関して積分をすることで最終的な推論が実行される．

これを模倣する形で提案されたのが 多重代入法 (MI: Multiple Imputation) (Rubin, 1978), (Rubin, 1987) である．

多重代入法ではベイズ事後分布から補完値を複数生成し，複数の擬似完全データに関して同じ解析を実行し，最後に結果を平均する．

擬似完全データに対する解析が一貫したベイズ推論であった場合，この一連の手続きによって（近似的な）ベイズ推論が実行されることになる．

しかしデータの補完とその後の擬似完全データ解析は 融和性 (congeniality) を保つ限り別の方法を用いても良いように拡張された (Meng, 1994)．

このことにより多重代入法は広く使われるようになっている．

5.2 多重代入法

多重代入法では，モデルベースの代入法（第 4.7 節）をさらに推し進める．

本来の推定量 \[ \widehat{Y}=\sum_{i\in S}w_iy_i \] を代入推定量 \[ \widehat{Y}_\mathrm{I}=\sum_{i\in S}w_i\biggr(\delta_iy_i+(1-\delta_i)y_i^*\biggl),\qquad y_i^*\sim f_\theta(y_i|x_i) \] で模倣する際，ベイズ事後予測分布で \[ y_i^*\sim f(y_i|y_{\text{obs}}) \] によって補間することが理想的である．

(Multiple Imputation Rubin, 1978)

事後予測分布から補間値を \(M\) 個生成する： \[ y_i^{(j)}\sim f(y_i|y_{\text{obs}}),\qquad j\in[M]. \]
それぞれの補間値について推定量 \(\widehat{Y}^{(j)}\) を計算し，その平均を最終的な推定値とする： \[\newcommand{\MI}{\mathrm{MI}} \widehat{Y}_\MI:=\frac{1}{M}\sum_{j=1}^M\widehat{Y}^{(j)}. \]

(Royston and White, 2011) は \(M\approx10^3\) を推奨している．

5.3 連鎖方程式による多重代入

多重代入法において事後予測分布から補間値を生成することは，\(Y\) に関してモデルを立てる必要があるためネックになりがちである．

相互条件付き識別性 (FCS: Fully Conditional Specification) (Stef Van Buuren and Rubin, 2006) が成り立つモデルについては，モデルの具体的な形に依らない Gibbs サンプラーによるサンプリングが可能になる．

これを 連鎖方程式による多重代入 (MICE: Multiple Imputation by Chained Equations) (Buuren and Groothuis-Oudshoorn, 2011) といい，R 言語 mice パッケージで実装されている．

その実用性も相まってか，近年の Lancet 誌，New England Journal of Medicine 誌のレビューでは，欠測データの取り扱いに最も多く用いられている手法は MICE であるという報告もある(Hayati Rezvan et al., 2015)．
(久史, 2017, p. 75)

5.4 その他の代入法

ランダムな欠損ではなく，計画された大規模な欠損がある場合は，two-phase sampling の考え方を応用することができる (Kim, 2024, p. 173)．

なお，全ての代入法はモデル \(Y|(X,\delta)\) の特定を間違えると，\(\widehat{Y}\) の不偏性が失われることに注意 (Hayati Rezvan et al., 2015)．

5.5 代入をしない

代入をせず，欠測しているなら欠測したままで最尤推定を実行することも考えられる．

このアプローチは 完全情報最尤推定 (FIML: Full Information Maximum Likelihood)，より最近では　pairwise likelihood estimation とも呼ばれる．¹¹

欠測が \(Y\) に依存しない場合，この「最尤推定量」は MAR の下で一致性と漸近正規性を持つ．¹²

ただし，推定されたモデルから，欠測値を代入してから結果を出してももちろん良い．ベイズの観点からは，モデルの平均を取ってから予測することに当たる．¹³

(1.6節 Buuren, 2018) も参照．

5.6 欠測値をどう扱うべきか？

いつでも多重代入法を使えば良いというものではない．

例えば \((X,Y)\) の関数関係が知りたい回帰分析の状況下で被説明変数 \(Y\) の欠損は，これを無視してリストワイズ消去をした complete-case analysis が代入法と等価になる．

他にも complete-case analysis や代入をしない方がむしろ適切な場合は多い (2.7節 Buuren, 2018)．

6 終わりに

6.1 多重代入法について

パッケージに実装される都合上，多重代入法はパラメトリックな手法であるという言説があるが，必ずしもそうである必要はない．この場合，傾向スコア推定量や校正推定量がセミパラメトリックな手法と呼ばれる．

また多重代入法が代入に使われたのちに，後続の解析は全く違うモデルが使われることもあり，このような場合は 融和性 (congeniality) (Meng, 1994) の議論が必要になる．

特に公的統計においては，後続のタスクが One Number Principle に従うようにするために欠測のあるデータは代入してしまい，擬似完全データを作成することもあり得る (Section 13.1 Kim, 2024, p. 161)．

この観点から見れば，多重代入法とは擬似完全データを複数作ることで後続推定精度営みというようにも思える．

いずれの場合も，多重代入法の「代入法」としての側面が強調されるあまり理論的背景が捨象され，また多重代入法の実際の使われ方が使用可能なパッケージでの実装方式に強く依存され，元来の手法の数理的本体が見失われている状況と言うことができるだろう．

Bayesian inference draws no distinction between missing data and parameters; both are uncertain, and they have a joint posterior distribution, conditional on observed data. (Chapter 18 Gelman et al., 2014, p. 449)

6.2 MAR について

MAR は現行で最も緩い条件である．

そして，MAR が成立しているかは確認したがく，感度分析などが推奨される (逸見昌之, 2014)．

欠測するかどうか \(\delta\) が，欠測する所の値 \(Y\) に依存している場合，これを MNAR という．この場合のセミパラメトリック最適な推定法は (Morikawa and Kim, 2021) などが提案されている．

References

Aitken, A. C. (1936). IV.—on least squares and linear combination of observations. Proceedings of the Royal Society of Edinburgh, 55, 42–48.

Berg, E., Kim, J.-K., and Skinner, C. (2016). Imputation Under Informative Sampling. Journal of Survey Statistics and Methodology, 4(4), 436–462.

Buuren, S. van. (2018). Flexible Imputation of Missing Data. Boca Raton, FL.: CRC Press.

Buuren, S. van, and Groothuis-Oudshoorn, K. (2011). Mice: Multivariate imputation by chained equations in r. Journal of Statistical Software, 45(3), 1–67.

Cheng, P. E. (1994). Nonparametric estimation of mean functionals with data missing at random. Journal of the American Statistical Association, 89(425), 81–87.

Deming, W. E., and Stephan, F. F. (1940). On a least squares adjustment of a sampled frequency table when the expected marginal totals are known. The Annals of Mathematical Statistics, 11(4), 427–444.

Deng, L.-Y., and Wu, C. F. J. (1987). Estimation of variance of the regression estimator. Journal of the American Statistical Association, 82(398), 568–576.

Deville, J.-C., and Särndal, C.-E. (1992). Calibration estimators in survey sampling. Journal of the American Statistical Association, 87(418), 376–382.

Enders, Craig K., and Bandalos, D. L. (2001). The relative performance of full information maximum likelihood estimation for missing data in structural equation models. Structural Equation Modeling: A Multidisciplinary Journal, 8(3), 430–457.

Enders, C. K., Du, H., and Keller, B. T. (2020). A Model-based Imputation Procedure for Multilevel Regression Models with Random Coefficients, Interaction Effects, and Nonlinear Terms. Psychological Methods, 25(1), 88–112.

Firth, D., and Bennett, K. E. (1998). Robust models in probability sampling. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1), 3–21.

Fuller, W. A., and Kim, J.-K. (2005). Hot deck imputation for the response model. Survey Methodology, 31(2), 139–149.

Gelman, A. (2014). How Bayesian Analysis Cracked the Red-State, Blue-State Problem. Statistical Science, 29(1), 26–35.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2014). Bayesian data analysis. Boca Raton : CRC Press.

Gelman, A., and Little, T. (1997). Poststratification into many categories using hierarchical logistic regression. Survey Methodology, (1997002).

Godambe, V. P., and Joshi, V. M. (1965). Admissibility and bayes estimation in sampling finite populations. i. The Annals of Mathematical Statistics, 36(6), 1707–1722.

Hayashi, F. (2000). Econometrics. Princeton University Press.

Hayati Rezvan, P., Lee, K. J., and Simpson, J. A. (2015). The rise of multiple imputation: A review of the reporting and implementation of the method in medical research. BMC Medical Research Methodology, 15(1), 30.

Hernán, M. A., and Robins, J. M. (2020). Causal Inference: What If. Boca Raton: Chapman & Hall/CRC.

Horvitz, D. G., and Thompson, D. J. (1952). A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association, 47(260), 663–685.

Imai, K., and Ratkovic, M. (2014). Covariate balancing propensity score. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 76(1), 243–263.

Isaki, C. T., and Fuller, W. A. (1982). Survey design under the regression superpopulation model. Journal of the American Statistical Association, 77(377), 89–96.

Jean-Claude Deville, C.-E. S., and Sautory, O. (1993). Generalized raking procedures in survey sampling. Journal of the American Statistical Association, 88(423), 1013–1020.

Kim, J. K. (2024). Statistics in survey sampling.

Kim, J. K., and Haziza, D. (2014). DOUBLY ROBUST INFERENCE WITH MISSING DATA IN SURVEY SAMPLING. Statistica Sinica, 24(1), 375–394.

Kwon, Y., Kim, J. K., and Qiu, Y. (2024). Debiased calibration estimation using generalized entropy in survey sampling.

Little, R. J. A. (1992). Regression with missing x’s: A review. Journal of the American Statistical Association, 87(420), 1227–1237.

Meng, X.-L. (1994). Multiple-Imputation Inferences with Uncongenial Sources of Input. Statistical Science, 9(4), 538–558.

Morikawa, K., and Kim, J. K. (2021). Semiparametric optimal estimation with nonignorable nonresponse data. The Annals of Statistics, 49(5), 2991–3014.

Royston, P., and White, I. R. (2011). Multiple imputation by chained equations (MICE): Implementation in stata. Journal of Statistical Software, 45(4), 1–20.

Rubin, D. B. (1976). Inference and missing data. Biometrika, 63(3), 581–592.

Rubin, D. B. (1978). Multiple imputations in sample surveys - a phenomenological bayesian approach to nonresponse. Proceedings of the Survey Research Methods Section, ASA, 20–28.

Rubin, D. B. (1987). Multiple Imputation for Nonresponse in Surveys. John Wiley & Sons.

Särndal, C.-E., Swensson, B., and Wretman, J. (1992). Model assisted survey sampling. Springer New York.

Sen, A. R. (1953). On the estimate of the variance in sampling with varying probabilities. Journal of the Indian Society of Agricultural Statistics, 5, 119–127.

Shimodaira, H. (2000). Improving predictive inference under covariate shift by weighting the log-likelihood function. Journal of Statistical Planning and Inference, 90(2), 227–244.

Stef Van Buuren, C. G. M. G.-O., J. P. L. Brand, and Rubin, D. B. (2006). Fully conditional specification in multivariate imputation. Journal of Statistical Computation and Simulation, 76(12), 1049–1064.

Sugden, R. A., and Smith, T. M. F. (1984). Ignorable and informative designs in survey sampling inference. Biometrika, 71(3), 495–506.

Wu, C., and Sitter, R. R. (2001). A model-calibration approach to using complete auxiliary information from survey data. Journal of the American Statistical Association, 96(453), 185–193.

Yates, F., and Grundy, P. M. (1953). Selection without replacement from within strata with probability proportional to size. Journal of the Royal Statistical Society: Series B (Methodological), 15(2), 253–261.

Zieschang, K. D. (1990). Sample weighting methods and estimation of totals in the consumer expenditure survey. Journal of the American Statistical Association, 85(412), 986–1001.

久史. (2017). 連鎖方程式による多重代入法. 応用統計学, 46(2), 67–86.

狩野裕. (2019). 欠測データ解析のmissとmyth. 日本統計学会誌, 48(2), 199–214.

逸見昌之. (2014). 欠測データに対するセミパラメトリックな解析法――その理論的背景について――. 統計数理, 62(1), 103–122.

高井啓二, 星野崇宏, and 野間久史. (2016). 欠測データの統計科学：医学と社会科学への応用,Vol. 1. 岩波書店.

Footnotes

Inverse probability weighting estimator ともいう (Hernán and Robins, 2020, p. 22)．↩︎
結果的に，Weighted Least Squares と同じ形になっている．WLS は誤差分散が既知の形 \(W^{-1}:=\mathrm{diag}(\sigma_1^2,\cdots,\sigma_n^2)\) をしている場合の最良線型不偏推定量 (BLUE) である．一般に最小二乗法は広い設定で BLUE を与え続け，一般の既知の分散 \(V(X)\) を持つ場合は GLS (Generalized Least Squares) と呼ばれる．\(V(X)\) が既知である場合などなく，一般にはこれの推定から始める必要があり，これは Feasible GLS と呼ばれる (Hayashi, 2000, p. 59)．↩︎
この抽出計画に依らない性質を以て，(Särndal et al., 1992) は model-assisted 推定量と呼んでいる．model-dependent 推定量とは対照的である．↩︎
この２つの類似性は (Zieschang, 1990) が指摘している．一般の回帰分析の設定下では “GLS is more efficient than OLS under heteroscedasticity (also spelled heteroskedasticity) or autocorrelation” などと説明される．↩︎
ただし，余分な項があるために，正しく特定されている下では校正推定量よりもやや分散が大きい．↩︎
総務省統計局では，Imputation の訳語として「補定」を用いる．↩︎
最も古典的な形のものであり，母集団上の条件であることから，population MAR とも呼ばれる．母集団上の MAR と抽出計画の無視可能性 (Sugden and Smith, 1984) との２条件が成り立つとき，標本の MAR が成り立つ (Berg et al., 2016)．↩︎
\(Y\to X\to\delta\) が Markov 連鎖をなす，とも換言できる．↩︎
(Kim, 2024, p. 154) 定理12.1も参照．↩︎
一方で，重み付き推定方程式の解として定まる \(Z\)-推定量として構成することもできる．(5.2節高井啓二 et al., 2016, p. 163)．↩︎
完全情報最尤推定の言葉は初期の構造方程式モデリングプログラム AMOS に組み込まれて有名になっていた (Craig K. Enders and Bandalos, 2001)．直接尤度 (direct likelihood) または観測尤度 (observed likelihood) の方法ともいう (狩野裕, 2019)．完全尤度 (full likelihood) の用語は (高井啓二 et al., 2016) など．↩︎
(狩野裕, 2019) に素晴らしい解説がある．日本語の文献としては (高井啓二 et al., 2016) もあり，第５章で推定方程式の観点から解説されている．↩︎
そういえば Bayes 的な integral out に関して doubly robust という考え方はないのか？doubly robust の Bayesian counterpart はなんだろう？↩︎