欠測データ：平均で埋めると信頼区間が嘘をつく【第31回】

はじめに

第29章は不完全データの統計処理です。要するに、データに空白があるときどうするかという話です。

この章に入る前、私はここを軽く見ていました。空白があるなら、その行を捨てるか、平均で埋めるか、どちらかでいいだろうと思っていたのです。特に平均で埋めるのは無害に見えました。全体の平均を動かさないのだから、大きな害はないはずだ、と。

手を動かしてみたら、この予想は半分当たっていて、半分は完全に外れていました。

平均代入で平均は当たります。真の平均 452.0 万円に対して、推定値 451.67。ここは予想通りでした。

外れていたのはその先です。同じ推定に付ける95%信頼区間の被覆率が 81.5% まで落ちました。95%と名乗る区間が、5回に1回近く外すということです。しかもこれは、欠測行を捨てた場合（94.1%）より悪い。埋めたことで、データを捨てるより悪い状態になっていました。

原因を追いかけて分かったのは、ひとつの単純な事実でした。

表に数字を1つ書き込むという行為が、その人の分散をゼロにしてしまう。

この一文がこの章の背骨で、平均代入がダメな理由も、EM アルゴリズムが何をしているのかも、多重代入法で「なぜ複数回埋めるのか」も、全部ここから出てきます。

もうひとつ収穫がありました。MCAR・MAR・MNAR という3つの略語です。名前が似ていて区別がつかず、この章の最初の関門だと思っていたのですが、同じデータで欠測の仕組みだけを変えた3枚の散布図を並べたら一気に整理できました。しかも「MAR かどうかは検定できない」という事実が、なぜできないのかまで含めて数値で示せました。

いつものように、出てくる数値はすべて自分で計算し、理論値と突き合わせています。

この回で扱う用語

何に使うのか：空白のある表を前にして

数式の前に、この章が生まれた事情から入ります。

アンケートを 2000 人に配ったとします。聞いたのは年齢と年収です。年齢は名簿から分かるので全員埋まりました。年収は 600 人分が空白でした。

平均年収を出したい。素朴な手は2つです。

1. 空白の行を捨てて 1400 人で平均を出す（リストワイズ削除）
2. 空白を「答えてくれた人の平均」で埋める（平均代入）

どちらも、初めて聞くと十分まともに思えます。この章の全体は「この2つがいつどう壊れるのか」を明らかにして、壊れない道具を用意する話です。

そして壊れ方は、欠測がどういう仕組みで起きたかによって変わります。だから最初に欠測の分類が来ます。分類が先に来るのは、この章の構成として自然なのです。

欠測の仕組みは3種類ある

同じ 2000 人のデータで、回答率をどれも約70%にそろえて、「誰が答えなかったか」の決まり方だけを変えました。データの中身はまったく同じです。

灰色が無回答、青が回答ありです。黒破線が全員の真の平均、赤実線が回答者だけの平均。左のパネルだけ赤線と黒線が重なっています。

3つの略語はこう整理できます。

実測値です。真の平均年収は 452.8 万円でした。

MCAR はほぼ当たります。MAR は若い人が抜けたので +38.3 万円高く出て、MNAR は高収入が抜けたので −45.5 万円低く出ました。

MAR と MNAR の分かれ目は「偏りの原因が手元にあるか」

MCAR とそれ以外の差はすぐ分かります。難しいのは MAR と MNAR の区別で、ここが一番間違えられるところです。

MAR の行をもう一度見てください。回答者の年齢平均が 46.1 歳です（全体は 41.4 歳）。あからさまに高齢に偏っています。

そしてこの偏りは、データを見れば分かります。年齢は全員について観測できているからです。だから「年齢で重み付けし直す」「年齢を説明変数に入れる」という対処ができます。

MNAR はどうでしょうか。回答者の年齢平均は 38.0 歳で、これも偏ってはいます。しかし偏りの本当の原因は「年収が高い人が答えなかった」ことで、その年収は欠測しています。手元のデータをどう眺めても「高収入層が何人抜けたか」は分かりません。直す手がかりがデータの中に存在しないのです。

言い換えるとこうなります。

定義を式で書くと

$R$ を「回答したか」を表す 0/1 の変数、$Y_{\text{obs}}$ を観測できた値、$Y_{\text{mis}}$ を欠測した値とします。

MAR の式が言っているのは「$Y_{\text{obs}}$（年齢）を知っている状態なら、$Y_{\text{mis}}$（年収）を追加で知っても欠測確率の予想は変わらない」ということです。これは第2回でやった条件付き独立と同じ形です。

判定基準は変数の名前ではなく「何を測ったか」

MAR の原因は「年齢である」に限りません。観測できている変数なら何でもよいのです。性別でも居住地でも、それらの組み合わせでも。「見えている情報のどれかで説明がつく」なら MAR です。

そして MNAR も「原因が年収そのもの」に限りません。知りたい変数と関係している未観測の何かが原因なら、すべて MNAR になります。たとえば「慎重な性格の人は答えない」で、性格が年収と相関しており、性格を測っていないなら MNAR です。

ここで注意すべき点があります。未観測の原因があること自体は MNAR の条件ではありません。 その原因が知りたい変数とつながっているかどうかが分かれ目です。

具体例で確かめます。健康診断のデータで、体重は全員測り、血圧は一部が欠測しました。欠測の原因は「その日の血圧計の故障」で、検査室の記録は持っていません。

これは MCAR です。故障は、その人の血圧の値とも、年齢や体重とも無関係に起きています。だから欠測した人たちは全体からランダムに抜き取られた集団になります。検査室の記録を持っていなくても MCAR のままで、MCAR は他の変数の助けを必要としない唯一のケースだからです。機器の故障・記録用紙の紛失・データ入力の取りこぼしは、教科書が MCAR の例として挙げる定番です。

ではこの例を MNAR にひねってみます。

混んでいる時間帯は昼休みに来る現役の会社員が多く、空いている時間帯は平日昼間に来られる高齢者が多い。そして血圧は年齢と強く相関する。

こうすると、故障の起きた時間帯が「若くて血圧が低い人」を選択的に落とすことになり、MCAR が崩れます。このとき年齢を測っていれば MAR、年齢も測っていなければ MNAR です。

同じ「血圧計の故障」でも、故障のタイミングが患者の種類と相関しているかどうかで分類が変わる。 そしてこれは少し不思議な性質を含んでいます。

MAR か MNAR かは、データの性質ではなく「何を一緒に集めたか」で決まる。

同じ現象でも、年齢を聞いていれば MAR、聞いていなければ MNAR です。これは後で実務の話をするときに効いてきます。

名前についての警告

MAR（ランダムな欠測）は、ランダムではありません。 上の図の真ん中は、若い人が明らかに狙って抜けています。

「他の観測変数で条件付けたあとはランダム」という意味なのですが、この命名は統計学で最も分かりにくいものの一つと言われています。私は MAR = 条件付きランダム と読み替えることにしました。

MCAR は検定できる。MAR と MNAR は原理的に区別できない

分類が分かったところで、当然の疑問が来ます。実際のデータで、どれなのか判定できるのでしょうか。

答えは「MCAR は判定できる。MAR と MNAR は判定できない」です。しかも後半は「まだ良い方法が発明されていない」のではなく、原理的に不可能でした。

MCAR の検定

やることは単純です。回答した人としなかった人に分けて、両方について観測できている変数（年齢）を比べる。MCAR なら無回答は完全にランダムなので、2群の年齢は一致するはずです。第13回のウェルチのt検定で比べました。

MCAR かどうかは判定できました。ただしMAR と MNAR はどちらも同じ「棄却」に落ちます。この検定は両者を区別しません。

なお、複数の変数をまとめて検定する方法としてリトルの MCAR 検定（Little's MCAR test）があり、実務ではこちらが使われます。ただし区別できるのは MCAR かどうかまでで、上の限界は変わりません。

区別できないことの実演

手元のデータを1つ固定します。1400 人は年齢と年収の両方、600 人は年齢だけ。

この同じデータを生む、2つの世界を考えます。

答えは 35.9 万円も違うのに、手元のデータ（上の表の3行）はAとBで完全に同一です。

検定とは、観測データを入力にして判定を出す関数です。入力が同じなら出力も同じになります。

だから MAR と MNAR を見分ける検定は作れない。

これは検定の設計が下手だからではなく、情報が存在しないからです。判別に必要な情報が、まさに欠測している部分に入っています。

実務では、MAR か MNAR かは分野の知識で仮定するしかありません。「この調査で年収を答えなかった人は、年齢と職種を揃えれば答えた人と同じと考えて差し支えないか」を、統計ではなく調査設計の知識で判断することになります。

リストワイズ削除はいつ許されるのか

「欠測がある行を丸ごと捨てる」という最も単純な方法が、どこまで許されるのかを測ります。4000 回のシミュレーションで、真の平均年収 452.0 万円、真の傾き 8.0 です。

図は1回のシミュレーションで、表の値は 4000 回の平均です。MAR のパネルで赤と黒がほぼ重なっているのが要点です。

図に出ている傾きは1標本の値なので、表の 4000 回平均とは少しずれます。MAR は図で 7.70、4000 回平均は 7.9996（この 7.70 は真値 8.00 から −1.34 標準偏差の位置にある通常のばらつきの範囲内です）。MNAR は図で 7.18、4000 回平均は 7.0867。MNAR のずれは何度繰り返しても消えませんが、MAR のずれは繰り返せば消えます。この違いが偏りとばらつきの区別です。

MAR で平均は壊れるのに傾きは壊れない

MAR の行を見てください。平均は +39.68 も外すのに、傾きは −0.0004 しか外していません。

同じデータ、同じ削除操作なのに、壊れ方が違います。つまり

「MAR なら削除はダメ」は、推定したい量とモデルの正しさを指定しないと判定できない。

理由はこうです。削除によって年齢の分布は歪みます（若い人が減る）。しかし各年齢での年収の条件付き分布は歪んでいません。欠測確率が年齢だけで決まっているので、同じ年齢の中では回答者と無回答者が同じ分布を持つからです。

そして回帰の傾きは、この条件付き分布だけで決まる量です。だから無傷。一方、平均は「年齢の分布」の重みで加重平均する量なので、そこが歪めば直撃を受けます。

ただしこれは MAR 一般の性質ではない

ここは危ないので、条件を崩した反例を自分で測りました。真の傾きを 2.0 にした別のデータで、4000回のシミュレーションです。

③が一番怖いところです。真の関係が曲線なのに直線を当てると、その傾きは「$x$ の分布のどこを重く見るか」で変わってしまいます。全データに直線を当てたときの傾きは 1.9957 なのに、削除すると 3.4487。$x$ の分布が歪むことが、ここで直撃します。

②については前提をひとつ書き添えておきます。上の 1.8471 は $z = y + \text{誤差}$ とした設定で、この $z$ は $x$ とも相関しています。$z$ を $x$ と完全に無相関にする（$y$ の誤差項にだけ相関させる）と、傾きは 2.0000 に戻り、ずれるのは切片だけになりました。「$y$ と相関する観測変数に依存すれば必ず傾きが偏る」わけではないという点は、私も試して初めて分かりました。

逆に、誤差が不等分散でも傾きは無傷でした（実測 1.9983）。効率は落ちますが偏りは出ません。

用語についても正確に書いておきます。「欠測確率が、回帰モデルに入っている共変量だけで決まる」は MAR の特殊ケースで、共変量依存の欠測（covariate-dependent missingness）と呼ばれます。だから正しい言い方は「MAR だから傾きは大丈夫」ではなく、

共変量依存の欠測で、かつモデルが正しく特定されているから大丈夫。

になります。私は当初これを「MAR の定義そのもの」と書いていて、そこが誤りでした。

この区別は実務で意味を持ちます。「平均を知りたい」のか「関係を知りたい」のかで、許される手抜きが違うのです。①②を確認でき、モデルの形にも根拠があるなら、回帰係数だけが目的なら完全ケース分析が使える場面があります（③のモデルの正しさはデータから確認できる性質ではないので、ここは分野の知識に頼ることになります）。

MCAR でも代償はある

偏りが消えても、無料ではありません。

比は 1.1912 で、理論値 $\sqrt{2000/1400} = 1.1950$ とよく合います。精度はルートで悪化するわけです。第10回の「精度はルートでしか改善しない」の裏返しです。

そして多変量ではこれが深刻になります。10 変数がそれぞれ独立に 10% 欠測すると、全部揃っている行は $0.9^{10} \approx 35\%$ しか残りません。1変数あたりの欠測が少なくても、変数が増えると全滅に近づくのがリストワイズ削除の第2の問題です。

なお、この問題を避けるためにペアワイズ削除（相関を計算するときはその2変数が揃っている行だけ使う）という方法もありますが、こちらは組み合わせごとに使う人が違うため、分散共分散行列が正定値でなくなるという別の困りごとを生みます。第24回で扱った固有値分解が使えなくなる場合があり、多変量解析の前処理としては危険です。

平均代入がダメな理由

ここが今回いちばん学びのあった部分です。設定は MCAR で 30% 欠測、つまり最も条件の良いケースにしました。偏りの心配がない状況で試して、それでも壊れるかを見ます。n=200、4000 回のシミュレーションです。

赤い×が埋めた値です。図の中央（平均代入）では1本の水平線に並び、右（回帰代入）では1本の斜め線に並びます（図のパネル番号①②③は上の表の行番号とは対応していません。図は真のデータ・平均代入・回帰代入の3枚です）。どちらも「そこにばらつきがない」ことが問題です（図は1回のシミュレーションなので、値は上の表の4000回平均とは少しずれます）。

まず、平均は当たっている

③の平均は 451.67。当たっています。

これは当然で、$n$ 個の平均に、その平均自身を何個足しても平均は変わりません。だから「平均代入は無害に見える」という直感は錯覚ではなく、平均だけを見るなら本当に無害です。

壊れるのはそれ以外の全部でした。

壊れたもの1：ばらつきが16.7%過小

SD 100.38 は真値 120.50 の 83.3% です。理論値がぴったり出ます。欠測割合を 30%（残り 70%）とすると、

$120.50 \times 0.8367 = 100.82$ で、実測 100.38 と合います。

理由は図②の赤い×が全部同じ高さに並んでいることです。平方和 $\sum (y_i - \bar y)^2$ を作るとき、埋めた60人の寄与はちょうど 0 になります（埋めた値が平均そのものなので差がゼロ）。分子は観測された 140 人分のまま、分母だけ 200 になる。だから $140/200 = 0.70$ 倍です。

壊れたもの2：信頼区間が嘘をつく

ここが最も深刻でした。

平均代入が報告する標準誤差は 7.098 です。ところが推定値の本当のばらつきは 10.399（リストワイズ削除と同じ値になります。情報量は増えていないので当然です）。1.47 倍も自分を過大評価しているわけです。

その結果、95%信頼区間の被覆率が 81.5% に落ちます。

そして表を縦に見てください。欠測なし（94.8%）より悪い。削除（94.1%）より悪い。 つまり

平均代入は、データを捨てるより悪い。

捨てれば少なくとも嘘はつきません。埋めると 200 人分のデータがあるかのように振る舞ってしまう。実際には 140 人分の情報しかないのに、そう見えなくなるのです。

第11回で信頼区間を扱ったとき「信頼区間は作り方に付いている数字」という話をしましたが、これはその作り方が壊れている例です。推定値が正しくても、精度の申告が嘘なら結論は信用できません。

壊れたもの3：相関と回帰係数が逆方向に壊れる

2つの代入法が逆方向に壊れているのが面白いところです。

数字の関係を先に確認しておきます。平均代入では共分散が 0.70 倍になります（埋めた60人の $y_i - \bar y$ がゼロなので $\sum(x_i-\bar x)(y_i-\bar y)$ への寄与が消え、分子は140人分のまま分母が200になる）。$s_x$ は不変、$s_y$ は $\sqrt{0.70}$ 倍。だから

同じ現象なのに傾きと相関で係数が違うのは、相関だけが $s_y$ の縮小に助けられて半分だけ戻るからです。

平均代入では、埋めた点が水平に並びます。つまり年齢と年収の関係がない点を60個追加したことになる。だから相関も傾きも薄まります。

回帰代入では、埋めた点が回帰直線上にぴったり乗ります。つまり誤差ゼロの完璧な点を60個追加したことになる。だから相関が真値より高く出ます。

なお回帰代入で傾きだけが当たるのは、追加した点が既存の回帰直線の上に乗っているので直線の向きを変えないからです。当てはまりの良さ（相関）は水増しされますが、直線そのものは動きません。

回帰代入も不十分だが、惜しいところまで来ている

④の行を見ると、平均 451.73 も傾き 7.99 も合っており、平均代入より明らかにましです。しかし SD 111.91（真値 120.50）と標準誤差 7.913（本当は 9.387）はまだ過小で、被覆率も 90.1%。

足りないものは1つだけでした。

埋めた値に「予測の誤差」が入っていない。

48 歳の人の年収は「たぶん 510 万円」ですが、本当は $510 \pm 80$ くらいのはずです。回帰代入は予測の中心だけを置いて、$\pm 80$ を捨てている。だから点が直線に張り付き、ばらつきが足りなくなります。

ここから解決策が2方向に分かれます。

1. 誤差を足して埋める。 ただし1回では乱数任せなので、乱数を変えて何回も埋める → 多重代入法
2. そもそも埋めるのをやめる。 埋めずに、欠測を持ったまま尤度を最大化する → EM アルゴリズム / 完全情報最尤法

MCAR でも代入した方が得だという話

ここでひとつ、私が誤解していた点を書いておきます。

「MCAR なら削除しても偏らないのだから、それで十分」と思っていました。偏りについては正しいのですが、精度については損をしています。

MCAR・欠測30%で、平均年収を推定した結果（6000回、n=2000）です。

どちらも偏りはゼロですが、分散が18.5%減っています。標準誤差は「欠測なし」比で 1.20 倍から 1.08 倍まで戻りました。

なぜかというと、捨てた600人にも情報が残っているからです。年収は分かりませんが、年齢は分かっている。そして年齢と年収の相関は $r = 0.747$ です。「48歳の人の年収は不明だが、たぶん510万円くらい」という予測は、まったくの無情報より確実にましです。捨てるとこの情報を自分から放棄することになります。

つまり MCAR での正しい理解は、こうです。

削除しても答えは合う。ただしデータを無駄にしている。

そして無駄にした分を追加調査で買い戻すには費用がかかります。この章の道具は、追加調査なしで精度を回収する技術でもあるわけです。

そして MAR なら、偏りそのものが消える

もうひとつ大事な点です。ここからは平均の話です（傾きについては前節のとおり、条件がそろえば削除でも無傷でした）。3000 回のシミュレーションで、削除と「年齢を使った回帰代入」を比べました（真の平均 452.0 万円）。

MAR の行を見てください。削除だと +39.60 外すのに、年齢で埋めると +0.02。ほぼ完全に直りました。

これは当然のことで、MAR の定義は「欠測の原因が観測されている変数で説明できる」でした。原因が手元にあるなら、それを使って埋めればよい。 削除が MAR で壊れるのは、MAR が直せないからではなく、削除という手段がその情報を使っていないからです。

MNAR でも −46.48 が −22.63 に半減しています。完全には直らないが、何もしないよりはましという位置づけです。

ここで用語を1つ区別しておきます。「代入すると良くなる」の「良くなる」が2つの別物を指しています。

前者はおまけ、後者は必須です。MAR が「正しい結果を出せる境界線」になります。

EM アルゴリズム

「E ステップと M ステップを繰り返す」という説明を何度も読んだのですが、それでは何も見えませんでした。そこで n=10 の例で反復を1行ずつ数値で追いました。

使うデータ

x は全員わかっていて、y は後半4人が欠測です。しかも欠測しているのは x が大きい側だけという、いかにも危なそうな配置にしました。

2変量正規分布を仮定して、$\mu_y$（y の平均）・$\sigma_{yy}$（y の分散）・$\sigma_{xy}$（x と y の共分散）を推定します。x 側は全員観測なので $\mu_x = 4.2$、$\sigma_{xx} = 1.32$ で確定です。

2つのステップが何をしているか

この2つは互いの出口が互いの入口になっています。 パラメータが分かれば埋められる、埋まればパラメータが出る。しかし最初はどちらも手元にありません。

だから適当な値から始めて往復させる。

これが EM の全体像です。「鶏と卵」の状況を、片方を仮に置くことで回し始める、という発想です。

反復を数値で追う

初期値をわざと外して、$\mu_y = 3.0$、$\sigma_{yy} = 2.0$、$\sigma_{xy} = 0$（x と y は無関係だと思い込んでいる状態）から始めます。

上段で動きが見えます。反復1では傾きが 0 なので、4人を全部 3.0 という同じ高さに置いています。これは平均代入と同じ状態です。反復4で少し斜めになり、反復80では直線に沿って階段状に並びます。

そして赤いひげが縮んでいくのも見てください。ひげの長さは $\pm\sqrt{\sigma^2_{y \mid x}}$ なので、条件付き分散が $2.00 \to 0.054$ と小さくなるのに対応して、ひげは $\pm 1.414 \to \pm 0.232$ と短くなります（パネルのタイトルに出ている 2.0000 や 0.0540 は分散の値で、ひげの長さはその平方根です）。EM は「埋める値」と「その値に対する自信」を同時に育てています。

下段では、$\sigma_{yy}$ が一度 1.29 から 0.46 付近まで下がってから 1.14 まで戻る非単調な動きをしています（表に載せていない反復7が谷底で、最小値は反復6の 0.4569 と反復7の 0.4768 のあいだ、図では約 0.45 に見える位置です）。単調に改善するわけではないのが面白いところです（後述しますが、単調非減少が保証されているのは尤度であってパラメータではありません）。

収束先は最尤推定値そのもの

この欠測パターン（x が全観測で y が単調に欠測）では、運良く最尤推定値が閉じた式で書けます。完全ケース6人で回帰を推定すると $\hat b = 0.907143$、切片 $0.899048$、残差分散 $0.054032$ で、ここから直接計算できます。

12桁一致しました。 ここが私の最大の疑問への答えでした。

EM は新しい推定量ではない。第10回でやった最尤推定値を求めるための計算手順である。

答えは最尤法と同じもので、そこに至る道が反復になっているだけです。「EM 推定量」という別の推定量があるわけではありません。

他の方法と並べると差が見えます。

削除だと $\mu_y$ が 3.98、EM は 4.71 で、0.73＝73万円分の差です。欠測しているのは x が大きい（年齢が高い）側で、x と y には強い正の関係がある。だから「捨てた4人は高い側にいたはずだ」という情報を、x を通じて回収しているわけです。

肝は E ステップの「+ 条件付き分散」

ここが平均代入・回帰代入と決定的に違う点です。M ステップで分散を計算するには $E[y^2]$ が必要ですが、

と、右辺に第2項が付きます。これは新しい式ではなく、第3回の分散の定義を移項しただけです。

に、全部「$\mid x$」を付けただけ。$\sigma^2_{y \mid x}$ は $V[y \mid x]$ の別名です。

数値で確かめます。7人目（$x = 4.8$）は、収束時には「5.253 を中心に、分散 0.054 でばらつく」という分布だと見なされています。仮にこの人の年収が等確率で2通り、$5.253 \pm 0.2324$（$0.2324 = \sqrt{0.054}$）だとしましょう。つまり 5.0206 か 5.4854 です。

一方、中心の値を2乗しただけなら $5.253^2 = 27.594$ です。差は $27.648 - 27.594 = 0.054$ で、ぴったり条件付き分散になります。

y の平均は 5.253 だが、$y^2$ の平均は $5.253^2$ ではない。 第10回で扱った $E[g(y)] \ne g(E[y])$ の話です。そして $g(y) = y^2$ は凸関数なので、イェンセンの不等式から $E[y^2] \ge (E[y])^2$ と必ず上に外れる。その外れ幅がちょうど分散です。

表に数字を書くと分散がゼロになる

M ステップで分散を計算するとき、平方和が必要で、そこに $y^2$ が入ります。

7人目に何を足すか。

• 回帰代入：表に 5.253 と書いたので、$5.253^2 = 27.594$ を足す
• EM：$27.648$ を足す（$= 27.594 + 0.054$）

回帰代入は1人あたり 0.054 ずつ足し忘れている。 これが積み上がって分散が過小になっていました。前節の「SD 111.91 vs 真値 120.50」の正体はこれです。

言い方を変えると、これが今回の核心になります。

表に1つの数字を書き込むという行為が、そもそも「その人の分散はゼロ」を意味してしまう。

だから EM は表に数字を書かず、「この人は 5.253 ± 0.2324 の分布」という幅ごと持ち回ります。E ステップの正式な説明が「値を代入するのではなく対数尤度の期待値を計算する」なのは、$y$ と $y^2$ の両方の期待値を扱っているという意味でした。

第18回の IRLS と似ているのか

どちらも「解析的に解けないから反復する」のは同じですが、解けない理由が違います。

表の「保証なし」については補足が必要でした。実際に測ってみると、失敗の仕方が2種類ありました。

ロジスティック回帰の対数尤度は正準リンクなので凹であり、普通のデータを普通の初期値で回せばまず単調に増加します。「保証なし」が効いてくるのは、初期値を大きく外してフルステップで行き過ぎる場合で、実務ではステップ幅を縮めるだけで防げます。

もうひとつの完全分離は性質が違って、尤度は下がらないのに推定値が発散するという失敗です。初期値の問題ではなく最尤推定値が有限の値として存在しないためで、これは EM に替えても解決しません（第18回で扱った Hauck-Donner 効果の背景です）。

ロジスティック回帰の場合、尤度は $\prod \pi_i^{y_i}(1-\pi_i)^{1-y_i}$ と完全に書けます。データは全部揃っている。困るのは微分して0にした連立方程式が $\beta$ について解けないという、純粋に代数の問題です。

欠測の場合はもっと手前で困ります。尤度を書こうとすると 7 人目の $y_7$ が式の中に必要なのに、その値がない。式が書けないので微分もできない。

そこで EM は $y_7$ を1つの数字ではなく分布として持ち、その分布で尤度の期待値を取ります。「$y_7$ が 5.02 なら尤度はこう、5.49 ならこう、それを重み付けして平均」という形にすれば、$y_7$ が式から消えてパラメータだけの関数になる。そうなればいつもの最尤推定が使える。 M ステップが「いつもの平均・分散の式を当てるだけ」で済んだ理由です。

遅いが頑丈というトレードオフ

今回 80 回かかりました。ニュートン法なら5回程度で済む場面です。ではなぜ EM を使うのでしょうか。

1. 微分がいらない。 2階微分（ヘッセ行列）を導出・実装するのは面倒で、欠測パターンが複雑だと現実的でありません
2. 絶対に下がらない。 ニュートン法は初期値が悪いとステップが行き過ぎて発散します（上の実測で対数尤度が $-768 \to -1.7\times10^{28}$）。EM は初期値をどれだけ外しても（今回は $\sigma_{xy} = 0$ から始めました）尤度が単調非減少で、必ずどこかに着きます
3. 1回の更新が「いつもの推定式」で書ける。 実装が驚くほど簡単です

安定性を速度で買っているわけです。

ただし2に関して弱点もあります。尤度が単調非減少であることは保証されますが、たどり着くのが最大値だとは保証されていません。 一般に保証されるのは停留点に着くことまでで、局所最大だけでなく鞍点に捕まることもあります。混合正規分布のあてはめでは実際に起きます。第24回で扱った「停留点と最大値は別物」という話と同じ構図です。対策は初期値を複数試すことです。

なお EM の反復は $\theta^{(1)} \to \theta^{(2)} \to \theta^{(3)} \to \cdots$ と、現在の値だけから次の値が決まる写像の繰り返しです。第22回のマルコフ連鎖で扱った「$\pi P = \pi$ を満たす定常分布に収束する」のと同じ構図で、EM の収束先も $\theta = M(E(\theta))$ という不動点です。$\sigma_{yy}$ が非単調に動いたのは、この写像が不動点に向かう途中の挙動でした。

EM は欠測専用の道具ではない

応用範囲はこの章より広いです。「観測できない何かがある」という形をしていれば全部 EM に乗ります。

私は研究で隠れマルコフモデルを使ったことがあるのですが、そのライブラリの中身のパラメータ推定は EM でした（この文脈ではバウム・ウェルチ法と呼ばれます）。「各時刻の隠れ状態」が欠測データで、E ステップで「時刻 $t$ に状態 $k$ にいた確率」を計算し、M ステップで推移確率と出力確率を更新する。今日追った反復と構造が同じです。

第26回のクラスター分析で扱った k-means 法も、EM の親戚です。「各点がどのクラスターに属するか」を割り当てる段階と、「クラスターの中心を計算し直す」段階の往復になっており、E と M の構造そのものです。違いは、k-means が所属を 0/1 で断定する（ハード割り当て）のに対し、混合正規分布の EM は「70% はクラスター1、30% はクラスター2」という確率で持つ（ソフト割り当て）点です。ここでも「1つに決めてしまうか、分布として持つか」という今回の主題が顔を出します。

打ち切りは欠測なのか

第19回で扱った打ち切りデータとの関係です。答えは「欠測に見えるが、実際には欠測していない（区間として観測されている）」でした。この違いが、打ち切りだけが直せる理由になります。

契約期間（解約までの月数）が平均24ヶ月の指数分布で、観測期間は12ヶ月とします。12ヶ月時点でまだ解約していない人は打ち切りです（このケースでは 60.6% が打ち切りになります）。4000 回、n=500 のシミュレーションです。

②が第19回で見た「打ち切り捨てで77%短くなる」現象です（$5.50/24.0 = 23\%$）。打ち切りが 60.6% なので、リストワイズ削除すると6割を捨てることになります。しかも捨てられるのは長く続いた人ばかりなので、残った4割は短命な人の集まりです。

一見 MNAR に見えるのに直せる理由

「解約月」という変数だけを見ると、長く続いた人が欠測している。つまり欠測が値そのものに依存しているので、これは MNAR の形です。前半でやった「年収が高い人が答えない」と同じ構造です。

ところが打ち切りには、通常の欠測にない情報が2つあります。

打ち切られた人の尤度への貢献を、密度 $f(t)$ ではなく $S(12) = P(T > 12)$ にする。「値は知らないが12以上だ」という情報を、そのまま尤度に書けるわけです。

つまり打ち切りは欠測に見える状況の中で最も恵まれたケースです。第19回で打ち切りをちゃんと扱えたのは、欠測の仕組みが既知だったからでした。

用語の整理（ここは私が混同していました）

当初、私はこれを「値そのものが欠測を決めるが仕組みが既知なので、無視できる（ignorable）欠測に分類される」と理解していました。これは誤りでした。

ルービンの枠組みでの ignorability の定義は「MAR ＋ パラメータの分離性」です。つまり MNAR かつ ignorable という組み合わせは存在しません。では打ち切りは何なのでしょうか。

答えは「打ち切りは欠測ではなく粗視化（coarsened）データ」です。値が消えているのではなく、「$T > 12$ という区間」に粗くなっているだけ。この枠組みは CAR（Coarsening At Random＝ランダムな粗視化） と呼ばれます。条件は「粗視化の仕組みが値に依存しない」ではありません（12ヶ月打ち切りは $T>12$ かどうかという値に依存して起きています）。正しくは粗視化の確率が、観測された粗いデータ（「$T>12$」という情報）だけで決まり、その区間内のどこに真値があるかには依存しないことです。MAR の粗視化版に当たる条件で、これが成り立つと尤度が分離して $S(12)$ を掛ける形が正当化されます。

「値そのものが欠測を決めるのに直せる」の正体は、そもそも欠測していないことでした。 区間として観測されているのです。

「仮定なしで直る」も言い過ぎでした

必要な仮定が2つあります。

ひとつは無情報打ち切り（non-informative / independent censoring）で、これが欠測データ論の MAR に対応します。全員が12ヶ月で観測終了する管理打ち切りなら自動的に満たされますが、「早く解約しそうな不満のある人が、先に連絡が取れなくなって脱落する」ような打ち切りだと崩れます。実際に試すと 24ヶ月が 46.01ヶ月になりました。

もうひとつは分布の特定です。上の 24.09 は「指数分布」を正しく仮定できたから出た数字でした。

ワイブルで 63.88 は 2.7 倍の過大です。60.6% が打ち切られている状況では、観測できているのは分布の左端だけで、右側の形はほぼ全部を仮定が決めている。カプラン・マイヤー法に逃げても、最後の観測点より先の平均は識別できません。

正しい言い方はこうなります。

打ち切りは、区間情報が残るので仮定の量が少なくて済む。ただしゼロではない。

逆向きの応用

第19回では打ち切りを生存時間解析の道具として扱いましたが、通常の欠測を打ち切り風に扱うという発想もあります。

アンケートで年収を数値で聞く代わりに「400万円未満／400〜600万／600万以上」と区間で聞けば、答えてもらいやすくなり、しかも打ち切りと同じ形の情報が手に入ります。値は不明でも区間が分かるので、尤度に正しく書ける。

「答えてくれない」を「区間だけ答えてくれる」に設計変更すると、欠測問題が打ち切り問題に格下げされる。 判別不能な MAR / MNAR の問いを、検証しやすい仮定に置き換えられるわけです。実務ではかなり有効な手だと思います。

多重代入法：なぜ複数回埋めるのか

私の疑問は「1回でよいのではないか」でした。乱数で誤差を足すなら、1回埋めれば十分な気がします。

結論から書くと、推定値については1回で十分でした。複数回必要なのは別の理由です。

まず確率的回帰代入

回帰代入に足りないのは「予測の誤差」でした。そこで埋める値を

埋める値 = 予測の中心 $\hat a + \hat b x$ + 乱数 $\varepsilon \sim N(0, \hat s^2)$

にします。これが確率的回帰代入です。点が直線に張り付かず、本物のデータと同じように散らばります。

これを $m$ 回繰り返して結果をまとめるのが多重代入法です。$m$ を変えて測りました（3000回、n=200、MCAR 30%欠測、真の平均 452.0）。

推定値の列を見てください。1回でも m=100 でも 452 前後で変わりません。

変わるのは報告するSEです。1回だけだと 8.502 と報告するのに本当は 9.957。過小評価なので被覆率が 90.9% に落ちます。m=2 にした瞬間 9.278 になり、被覆率が 94.2% に回復します。

複数回埋める理由は、標準誤差を正しく報告するためである。

そして m=2 でほぼ効果が出て、m=5 で十分でした。m=100 にしてもほとんど改善しません。教科書に「$m = 5$ 程度でよい」と書かれている理由が、この表で納得できました。

1回だと足りない理由

1回で埋めると、埋めた表は「本物のデータ」に見えます。そこから普通に標準誤差を計算すると、200人分のデータがある前提の値が出ます。しかし 57 人分の値は自分が乱数で作ったものなので、そこに情報はありません。

欠測の不確実性には2種類あります。

②を見るには、答えを複数持つ必要があります。だから乱数を変えて何通りも埋めて、答えがどれくらいばらけるかを実測する。第11回のブートストラップと同じ「ばらつきを知りたければ何回もやってみる」という発想です。

観測 143 人（青）は5枚とも同一で、赤い×だけが毎回違います。その結果、5枚の平均が 443.43 / 443.91 / 444.31 / 439.12 / 441.41 とばらけます。このばらつきそのものが②の測定値です。

ルービンの公式は何を足しているのか

上の5枚から出た数値です（この標本では観測143人・欠測57人・残差SD 74.81）。

推定値は単純平均です。

分散は2つを足します。これがルービンの公式です。

ルービンの公式が足しているのは、①と②である。

「もしデータが揃っていたら残っていた誤差」＋「揃っていないことによる追加の誤差」です。

1回の代入では $\bar U$ しか計算できません。$B$ は複数の答えがないと定義できない量です。これが「複数回埋める」の答えでした。

$(1 + 1/m)$ という補正がついているのは、$m$ 回という有限回で $B$ を推定していること自体の誤差を織り込むためです。$m$ が小さいほど大きくなり（$m=5$ なら1.2倍）、$m \to \infty$ で1に近づきます。

情報の損失割合

$B$ と $T$ の比が、欠測でどれだけ損したかを教えてくれます。

57人分（28.5%）が欠測しているのに、全分散のうち欠測に由来する分は 7.35%（同じ設定を繰り返したときの平均は 13.0%）でした。欠測割合の 28.5% よりずっと小さく収まるのは、年齢という強い予測変数（相関 0.75）があるおかげで、欠測した人についてもかなりのことが分かるからです。

ただし数字の読み方に注意が2つ必要でした。

ひとつは用語です。$\lambda$ は正確には「全分散のうち欠測に起因する割合」で、ルービンの FMI（Fraction of Missing Information＝欠測情報割合）は自由度補正を含む $\gamma = (r + 2/(\nu+3))/(r+1)$（ただし $r = (1+1/m)B/\bar U$）という別の式です。$\lambda$ は $\nu$ が大きいときの近似で、この標本では $\lambda = 7.35\%$ に対し FMI は 7.60% でした。近いのですが同じものではありません。

もうひとつはこの 7.35% が1標本の値だという点です。$m=5$ で $B$ を推定しているので、$\lambda$ 自体が激しくばらつきます。同じ設定を800回繰り返すと $\lambda$ の平均は 13.0%、5〜95パーセンタイルは 3.0%〜28.1% でした。7.35% は下振れした1回の引きです。 $m$ が小さいときの $\lambda$ は方向性の目安であって、報告する数字として信用しすぎない方がよさそうです。

この標本でのSEを並べると、こうなります。

多重代入は完全データにかなり近いところまで戻せています。そして削除より良い数字を出すのに、削除より正直です。平均代入は「削除より良い数字を報告するが嘘」でしたから、そこが決定的に違います。

3段階の手順

段階2が「いつもの分析をそのまま」でよいのが、多重代入の実務上の強みです。欠測用に特別な分析手法を覚え直す必要がありません。

EM だと分析ごとに欠測を組み込んだ尤度を書き直す必要がありますが、多重代入は「データを作る」段階で問題を解決してしまうので、後段は普通のツールが使えます。

なお実際の多重代入では、埋める値の生成にベイズ的な手続きを使います。回帰係数 $\hat b$ 自身にも推定誤差があるので、その事後分布からも $\hat b$ を引き直してから予測値を作る、という二段構えです。今回のシミュレーションは $\hat b$ を固定した簡易版なので、正式な多重代入よりわずかに $B$ を小さく見積もっています。多変量の欠測に対しては連鎖方程式による多重代入（MICE）が標準的で、変数ごとに順番に回帰して埋めるループを回します。

完全情報最尤法（FIML）

FIML（Full Information Maximum Likelihood＝完全情報最尤法）は、欠測を埋めず、各人が持っている情報だけで尤度を書いて足し上げる方法です。

（$A$ は両方ある人、$B$ は x だけの人）

「観測されていない変数は、尤度から積分して消す」だけです。 欠測している人を捨てるのでもなく、値を作るのでもなく、その人が知らせてくれた分だけを使う。

FIML は EM とまったく同じ最尤推定値を返します。違いは解き方だけで、EM は反復（E→M→E→M）で解き、FIML はこの $\log L$ を数値最適化で直接最大化します。今回の n=10 の例で EM が出した 4.709048 は、FIML の答えでもあります。

EM の弱点として挙げた「標準誤差が出ない」は実務では重要です。EM は点推定値を返しますが、その精度は教えてくれません（尤度の2階微分を別に計算するか、ブートストラップを回す必要があります）。

FIML と多重代入は標準誤差まで面倒を見てくれるので、実務ではこの2つが主流です。EM は「何が起きているか」を理解するのに最適な形をしていますが、道具としては裏方に回っている印象です。

感度分析：MNAR に対する唯一の誠実な対応

MAR か MNAR かは検定できませんでした。だからと言って諦めるのでもなく、MAR だと決めつけるのでもない第三の道が感度分析です。

やり方は単純で、MNAR の度合いをパラメータ $\delta$ にして振ってみます。

無回答者の年収 = （同じ年齢の回答者からの予測値） + $\delta$ + $\varepsilon$

$\delta = 0$ が MAR の仮定です（無回答者は同じ年齢の回答者と同じ）。$\delta > 0$ なら「無回答者は隠れて年収が高い」という MNAR を表します。

$\delta = 0$ での推定値は 426.5（オレンジの点）です。ところがこのデータの真の $\delta$ は 74 で、真の平均は 450.9（赤破線）でした。

つまり MAR を仮定した多重代入・EM も、実際が MNAR なら外します。ただし灰色の点線（リストワイズ削除の 403.8）はさらに下にあり、描いた範囲のどの $\delta$ でも点推定の線より下にあります（信頼区間の帯まで含めると、$\delta$ が負の側では帯が 403.8 に届きます）。MAR を仮定した推定は、少なくとも削除よりは真値に近い位置から出発しています。

そしてグラフのどの点が正しいかはデータから決められません。だから報告の仕方が変わります。

「無回答者が回答者と同程度（$\delta = 0$）なら平均は 426 万円。無回答者が 100 万円高いなら 459 万円。$\delta = 150$ を超えないと『平均は 475 万円以上』という結論にはならないので、そこまで極端な無回答バイアスが考えにくいなら、平均は 427〜476 万円の範囲と考えられる。」

1つの数字ではなく、条件つきの幅で答える。 この形にすると、読み手が自分の分野の知識で $\delta$ を選べます。

MNAR を検定で解決できない代わりに、判断を「$\delta$ の妥当な範囲はどこか」という答えられる問いに置き換えるのが感度分析の狙いです。

グラフの傾き自体も情報です。緩やかなら「$\delta$ が多少ずれても結論は変わらない＝頑健」、急なら「結論が仮定に大きく依存する＝慎重に扱う」と読みます。

なお MNAR に正面から取り組むモデルとして、ヘックマンの2段階推定（選択モデル）やパターン混合モデルがあります。前者は「回答するかどうか」を決める式と「年収の値」を決める式を同時に推定する方法で、経済学でよく使われます。ただし識別のために「回答するかどうかには影響するが年収には影響しない変数」を見つける必要があり、その仮定が満たされているかは検証できません。MNAR ではどこかに検証不能な仮定が必ず残るという構造は変わりません。

実務：アンケートの無回答をどう扱うか

優先順位をつけると、解析より前の設計段階でできることが圧倒的に効きます。

5について補足すると、代入モデルに入れる変数は最終的な分析に使わないものでも入れてよい（むしろ入れるべき）です。欠測を予測できる変数が多いほど MAR の仮定が現実に近づくからで、これは「モデルは単純な方がよい」という第32回のモデル選択の原則とは逆向きに働きます。代入は予測の問題であって、説明の問題ではないという区別です。

そして、やってはいけないことが1つあります。

平均代入。 これだけは明確に避けます。分散を縮め、相関を薄め、信頼区間で嘘をつく（被覆率 81.5%）。リストワイズ削除の方がまだ正直です。

第21回で扱った標本調査法の無回答バイアスとの関係も整理しておきます。あちらでは重み付け（無回答が多い層の回答者に重みを増やす）という対処を扱いました。これは今回の枠組みで言えば「層という観測変数で条件付ける」ことなので、MAR を仮定した対処です。同じ仮定の上に立っている2つの技術で、重み付けは層単位・多重代入は個人単位という粒度の違いになります。

自分が間違えていたこと

この章で自分の理解が間違っていた点を残しておきます。

1. 平均代入を「無害」だと思っていた

平均が動かないから大きな害はないだろうと考えていました。実際には平均だけが無事で、分散・標準誤差・相関・回帰係数がすべて壊れます。しかも被覆率が 81.5% で、データを捨てた場合（94.1%）より悪い。「何もしないより悪い処置」がありうるという点が、いちばんの発見でした。

2. MCAR を「誰が答えないか分からない状態」と表現した

これは MNAR の説明になってしまいます。MCAR は「どういう人が答えないという偏りが、そもそも存在しない」状態です。「分からない」のは MNAR の困りごとで、MCAR は偏りの有無という一段前の話で脱落しています。

3. 「未観測の原因があれば MNAR」と考えた

血圧計の故障のように、知りたい変数と無関係な未観測の原因なら、記録がなくても MCAR のままです。分かれ目は「未観測かどうか」ではなく「その原因が知りたい変数とつながっているか」でした。

4. 「MAR か MNAR か区別できないから正しい結果は出せない」と考えた

区別できないのは事実ですが、結論は間違いでした。MAR なら削除は偏るが、正しく埋めれば偏らない（+39.60 → +0.02）。そして実務では「MAR を仮定して推定し、MNAR だったらどれだけずれるかを感度分析で示す」という2段構えを取ります。

5. 「MCAR なら削除で十分」と考えた

偏りについては正しいのですが、精度で損をしています。捨てた人にも他の変数（年齢）が残っており、それを使えば分散を 18.5% 減らせました。追加調査なしで精度を回収できる分を、自分から放棄していることになります。

6. 「代入すると良くなる」の中身を区別していなかった

「良くなる」が2つの別物を指していました。MCAR で代入が効くのは精度が上がるからで、MAR で代入が効くのは偏りが消えるからです。前者はおまけ、後者は必須です。同じ「良くなる」で括ると、どちらが致命的かの判断を誤ります。

7. EM を「最尤法とは別の推定法」だと思っていた

EM の収束値は閉じた式の最尤推定値と12桁一致しました。EM は推定量ではなく、最尤推定値に到達するための計算手順です。「EM 推定量」という別物はありません。

8. 多重代入を「複数回埋めれば推定値が良くなる」と誤解した

推定値は 1 回でも m=100 でもほぼ同じ（452 前後）でした。複数回必要なのは標準誤差を正しく報告するためで、$B$（代入間分散）が複数の答えがないと定義できないからです。目的は推定値ではなく精度の申告でした。

9. 「$m$ は大きいほどよい」と思っていた

m=2 でほぼ効果が出て、m=5 で十分でした。m=100 にしても被覆率は改善しません（94.1% → 93.8%）。$(1+1/m)$ の補正が有限回の誤差を織り込んでいるので、$m$ を増やす見返りが早く飽和します。

10. 「MAR なら削除は必ず偏る」と思っていた

推定したい量によります。 MAR＋削除で平均は +39.68 も外しますが、回帰の傾きは −0.0004 しか外しません。削除で歪むのは説明変数の分布で、条件付き分布は歪まないからです。

ただし査読で分かったのは、この性質が MAR 一般のものではないことでした。成立するのは「$y$ だけが欠測」「欠測確率が、回帰モデルに入っている共変量だけで決まる（＝共変量依存の欠測）」「モデルが正しく特定されている」の3条件がそろったときだけです。モデルを誤特定すると傾きも 2.00 → 3.45 と壊れます。 私は当初これを「MAR の定義そのもの」と書いていて、そこが誤りでした。

11. 打ち切りを普通の欠測と同じ枠組みで考えていた

打ち切りは最も恵まれたケースでした。ただし当初の理解には2つ誤りがありました。

第一に、これを「MNAR かつ無視できる（ignorable）欠測」と分類したのが誤りです。ignorability の定義は「MAR ＋ パラメータの分離性」なので、MNAR かつ ignorable は存在しません。打ち切りは欠測ではなく粗視化（coarsened）データで、「$T>12$ という区間」として観測されています。直せる理由は「そもそも欠測していない」ことでした。

第二に「仮定なしで直る」も言い過ぎでした。無情報打ち切りと分布の特定という2つの仮定が必要で、真の分布がワイブルなら 24ヶ月を 63.88ヶ月と答えます。正しくは「仮定の量が少なくて済む。ただしゼロではない」です。

要点まとめ

次回

次回は第30章のモデル選択です。今回の話とは逆方向の問題を扱います。

今回は「情報が足りない」ときにどうするかでした。次回は「使える変数がたくさんある」ときに、どれを使うかを決める話です。第16回で「変数を増やすと決定係数は必ず上がる」という現象を見ましたが、その先にある赤池情報量規準（AIC）やベイズ情報量規準（BIC）、交差検証法を扱います。

そして意外な接続があります。今回の代入モデルでは「予測変数を惜しみなく入れるべき」でしたが、次回のモデル選択では「変数は少ない方がよい」という逆の原則が出てきます。同じデータ、同じ変数でも、予測が目的か説明が目的かで正反対の判断になるという話で、今回の最後に触れた区別がそこで正面から扱われます。

また今回の対数尤度が、次回は AIC の計算式の中に $-2\log L$ として現れます。EM で最大化していたものが、そのままモデル比較の材料になるわけです。

この連載の全体像とこれまでの回は統計検定準1級・独学連載のまとめにあります。