統計学 学習ノート

はじめに

高年収の3人が答えないと、平均は134万円も低くなる

「アンケートの空欄、その人ごと消してしまっていいの？」

例えば、10人に年収を聞いたとしましょう。本当の年収は、200万円から1,000万円まで、平均520万円です。ところが、年収の高い3人が答えてくれませんでした。残った7人で平均をとると、約386万円です。全員の平均より134万円も低くなりました。

計算はどこも間違っていません。問題は、「空欄がなぜ空欄なのか」を考えなかったことです。実は、データの空欄の扱いは、空欄の 理由 で正解が変わります。

そこで登場するのが、今回の主役 欠測（missing data）です。欠測とは、本来は値があるはずなのに、そのデータが手に入っていないことです。そして、そのようなデータを不完全データと呼びます。

前回の記事では、すべてのセルが埋まった分割表をモデルで表しました。今回の記事では、空欄が生まれる仕組みである欠測メカニズム、空欄を捨てる削除法、空欄を埋める補完法、そして欠測があっても最尤推定を行うEMアルゴリズムについて解説します。

欠測メカニズム：空欄の理由を3つに分ける

空席の理由が、くじか、見える名簿か、見えない点そのものか

はじめに、欠測が起きる理由の分類について解説します。

欠測メカニズム とは、どの値が欠測するかを決めている確率の仕組みのことです。データ全体を、観測された部分 $Y_{\text{obs}}$ と欠測した部分 $Y_{\text{mis}}$ に分けます。さらに、値が観測されたら1、欠測したら0をとる変数 $R$ を用意します。欠測メカニズムは、$R$ の条件つき分布です。

これは3種類に分けられます。

【欠測メカニズムの3分類】

1. MCAR（完全にランダムな欠測）: $P(R \mid Y_{\text{obs}}, Y_{\text{mis}}) = P(R)$
2. MAR（ランダムな欠測）: $P(R \mid Y_{\text{obs}}, Y_{\text{mis}}) = P(R \mid Y_{\text{obs}})$
3. MNAR（ランダムでない欠測）: $R$ の確率が $Y_{\text{mis}}$ に依存する

言葉で言うと、MCARはデータと無関係に抜ける場合です。MARは、見えている値を見れば抜け方が説明できる場合です。MNARは、見えていない値そのものが、抜ける理由になっている場合です。

教室の空席に例えてみましょう。くじで休みの人を決めたならMCARです。部活の大会がある人が休んだなら、部活の名簿を見れば説明できるのでMARです。テストの点が悪そうな人だけが休んだなら、休んだ理由は受けなかったテストの点そのものなので、MNARです。

具体的に見てみましょう。健康診断の問診票で、体重の欄に空欄があるとします。年齢はすべて記入されています。

【体重の空欄の例】

• 集計機が故障して、無作為に一部のデータが消えた： MCAR
• 若い人ほど空欄が多いが、同じ年齢の中では体重と無関係： MAR
• 同じ年齢の中でも、体重が重い人ほど書かなかった： MNAR

冒頭の年収の例は、年収の高い人が年収を答えなかったので、MNARです。

ここがポイントです。MARなら、欠測の仕組みをモデルに入れなくても、観測データの尤度だけで正しく最尤推定ができます。欠測の仕組みの母数を $\psi$、データの母数を $\theta$ とし、2つは互いに無関係だとすると、観測できるものの同時分布は次のように書けます。

$$f(Y_{\text{obs}}, R \mid \theta, \psi) = \int f(Y_{\text{obs}}, Y_{\text{mis}} \mid \theta)\, f(R \mid Y_{\text{obs}}, \psi)\, dY_{\text{mis}} = f(Y_{\text{obs}} \mid \theta)\, f(R \mid Y_{\text{obs}}, \psi)$$

MARの仮定により、$R$ の確率が積分の外へ出ます。$\theta$ を含むのは前の因数だけなので、後ろの因数は無視できます。これを、欠測メカニズムが 無視可能 であると言います。

一方で、MARとMNARの区別は、データだけからは確かめられません。見えていない値の影響を、見えている値で調べることはできないからです。そこで実務では、欠測の原因になりそうな変数も分析に入れてMARに近づけ、MNARのときに結果がどれだけ変わるかを感度分析で調べます。

試験のポイント：「欠測の起きやすさが、観測された別の変数で説明できる」ならMAR、「欠測している値そのものに依存する」ならMNARです。

以上が、欠測メカニズムについての解説です。

削除法：空欄のある行を捨てる

各列1割の欠測でも、行ごと捨てると4割の個体を失う

次に、最も手軽な方法である削除法について説明します。

削除法 は、欠測のあるデータを分析から取り除く方法です。虫食いのあるページを、本から破り捨てて読むようなものです。代表的なものは2つあります。

【2つの削除法】

• リストワイズ削除（完全ケース分析）： 1つでも欠測のある個体を丸ごと除く
• ペアワイズ削除： 計算に使う2変数がそろっている個体を、計算ごとに使う

リストワイズ削除は単純で、多くの統計ソフトの標準の動作です。ただし、弱点が2つあります。

1つ目は、データが大きく減ることです。例えば、5つの変数がそれぞれ独立に10%ずつ欠測すると、すべてがそろう個体の割合は $0.9^5 \approx 0.59$ です。1変数あたりの欠測は1割でも、4割の個体を失います。

2つ目は、MCARでないと偏りが生まれることです。MCARなら、残った個体は全体の無作為な部分標本なので、推定量に偏りは出ません。標準誤差が大きくなるだけです。しかし、MARやMNARでは、残った個体の顔ぶれが全体と違います。

ペアワイズ削除はデータを多く残せますが、相関係数ごとに使う個体が変わります。そのため、相関行列の全体で矛盾が起き、分散が負になるような行列が出ることもあります。

例題：削除法の偏りを計算する

100人の生徒がいて、中間テストの成績で上位50人と下位50人に分かれます。期末テストの平均は、上位が80点、下位が60点でした。全員の平均は70点です。ところが、期末テストの欠席者が上位で10%、下位で50%出ました。中間テストの成績は全員分があるとして、次を求めましょう。

(1) 期末テストを受けた人だけの平均点　(2) 群ごとの平均から全体を推定した値

【解き方】

1. 受けた人数を求める： 上位45人、下位25人
2. 受けた人の平均： $(45 \times 80 + 25 \times 60) \div 70 \approx 72.9$
3. 群の人数比で平均する： $0.5 \times 80 + 0.5 \times 60 = 70$

答えは (1) 約72.9点、(2) 70点です。欠席のしやすさが中間テストの成績で説明できるMARなので、リストワイズ削除は偏りますが、観測された中間テストの情報を使えば偏りを直せます。「MARなら観測値を手がかりにできる」とは、こういうことです。

以上が、削除法についての説明です。

補完法：空欄を埋めて計算する

平均値で埋めると、平均は同じでも分散が縮む

続いて、空欄に値を埋める補完法について解説します。

補完法（代入法、imputation）は、欠測した値をもっともらしい値で埋めてから、ふつうの分析を行う方法です。ジグソーパズルの抜けたピースを、周りの絵柄から作って埋めるイメージです。1つの値で埋める単一代入と、複数の値で埋める多重代入があります。1つずつ見ていきましょう。

1. 単一代入

代表的な方法は次のとおりです。

【単一代入の方法】

• 平均値代入： 観測値の平均で埋める
• 回帰代入： 他の変数からの回帰の予測値で埋める
• ホットデック法： 似た個体の実際の値を借りて埋める
• 確率的回帰代入： 回帰の予測値に、残差の分だけ乱数を足して埋める

一見よさそうな平均値代入には、大きな弱点があります。ばらつきを小さく見せてしまうのです。

例えば、10人のテストのうち8人の点数が 52, 58, 60, 62, 68, 70, 74, 76 で、2人が欠測だとします。8人の平均は65点、不偏分散は約69.7です。欠測の2人を65点で埋めると、平均は65点のままですが、偏差の平方和が増えないまま分母が増えます。

$$s^2_{\text{代入後}} = \frac{\sum_{\text{観測}} (y_i - \bar{y})^2}{n - 1} = \frac{n_{\text{obs}} - 1}{n - 1}\, s^2_{\text{obs}}$$

ここでは $7/9$ 倍の約54.2に縮みました。分散が小さく見えるので、標準誤差も信頼区間も小さくなりすぎます。さらに、他の変数との相関も弱まります。

回帰代入は、平均値代入より賢いですが、埋めた値が回帰直線の上にぴったり乗るので、今度は相関を強く見せます。どちらも、「埋めた値には誤差がある」という事実を無視していることが問題です。

2. 多重代入

そこで登場するのが 多重代入法（multiple imputation）です。ここでは、乱数を使って補完したデータセットを $m$ 個作り、それぞれを分析して、結果をまとめます。1人の予想屋に頼らず、$m$ 人の予想屋に予想を出してもらい、意見のばらつきも結果に含めるイメージです。

$j$ 番目のデータセットの推定値を $\hat{Q}_j$、その分散を $U_j$ とします。全体の推定値と分散は、ルービンのルール でまとめます。

$$\bar{Q} = \frac{1}{m}\sum_{j=1}^{m} \hat{Q}_j, \qquad T = \bar{U} + \left(1 + \frac{1}{m}\right) B$$

$\bar{U}$ は $U_j$ の平均で、データセットの中の分散です。$B = \frac{1}{m-1}\sum_j (\hat{Q}_j - \bar{Q})^2$ は、データセットの間の分散です。言葉で言うと、$B$ が「埋め方が分からないことによる不確かさ」を表し、それを足すことで標準誤差が小さくなりすぎるのを防いでいます。

例えば、$m = 5$ で平均点の推定値が 70.2, 71.0, 69.8, 70.6, 70.4、それぞれの分散が 1.2, 1.1, 1.3, 1.2, 1.2 だったとします。$\bar{Q} = 70.4$、$\bar{U} = 1.2$、$B = 0.2$ なので、$T = 1.2 + 1.2 \times 0.2 = 1.44$ です。標準誤差は $\sqrt{1.44} = 1.2$ となります。

試験のポイント：平均値代入は平均を保っても、分散と相関を小さくします。多重代入は補完の不確かさを $B$ として反映する方法です。

以上が、補完法についての解説です。

EMアルゴリズム：埋めて、解いて、繰り返す

下書きと清書を繰り返すと、50から始めた値が65に落ち着く

最後に、欠測があるときの最尤推定を見ていきましょう。

最尤推定をしたいのに、欠測があると尤度の式が複雑になります。ここには、ニワトリと卵のような関係があります。母数 $\theta$ が分かれば欠測値を予想できますし、欠測値が分かれば $\theta$ を最尤推定できます。どちらも分からないなら、仮の値から始めて、交互に更新すればよいのです。

これが EMアルゴリズム（Expectation-Maximization algorithm）です。完全データの対数尤度を $\ell_c(\theta) = \log f(Y_{\text{obs}}, Y_{\text{mis}} \mid \theta)$ として、次の2つのステップを繰り返します。

【EMアルゴリズムの手順】

1. Eステップ： 現在の値 $\theta^{(t)}$ のもとで、$\ell_c$ の条件つき期待値を求める
2. Mステップ： その期待値を最大にする $\theta$ を $\theta^{(t+1)}$ とする

$$Q(\theta \mid \theta^{(t)}) = E\left[\ell_c(\theta) \mid Y_{\text{obs}}, \theta^{(t)}\right], \qquad \theta^{(t+1)} = \arg\max_{\theta} Q(\theta \mid \theta^{(t)})$$

下書きと清書の繰り返しに例えられます。今の案で空欄を下書きし、下書きを前提に清書して案を更新します。すると、案はだんだん落ち着いていきます。

1. 正規分布の平均

先ほどの10人のテストで、$\mu$ を推定してみましょう。8人の合計は520点です。Eステップでは、欠測2人の期待値を $\mu^{(t)}$ とします。Mステップでは、10人分の平均をとります。

$$\mu^{(t+1)} = \frac{520 + 2\mu^{(t)}}{10}$$

$\mu^{(0)} = 50$ から始めると、62.0、64.4、64.88、64.976と進み、8人の平均の65に収束します。1回ごとに、目標との差が0.2倍になっています。0.2は欠測の割合です。一般に、EMアルゴリズムの収束の速さは、欠測した情報の割合で決まります。欠測が多いほど、ゆっくり進みます。

ただし、EMアルゴリズムは「欠測値を埋める」方法そのものではありません。埋めるのは、完全データの十分統計量の条件つき期待値です。例えば、分散を推定するときは $y^2$ の期待値が必要で、これは $\{E[y]\}^2$ に条件つき分散を足した値です。平均値代入のようにばらつきが縮まないのは、このためです。

例題：打ち切りのある寿命をEMアルゴリズムで推定する

寿命が平均 $\theta$ の指数分布に従う電球を10個調べました。7個は途中で切れ、その寿命の合計は350時間でした。残りの3個は100時間の時点で切れておらず、そこで観測をやめました。生存時間の記事で学んだ打ち切りです。$\theta$ を求めましょう。

【解き方】

1. Eステップ： 指数分布の無記憶性から、$E[X \mid X > 100] = 100 + \theta^{(t)}$
2. Mステップ： 完全データの指数分布の最尤推定量は標本平均なので、次の更新式になる

$$\theta^{(t+1)} = \frac{350 + 3(100 + \theta^{(t)})}{10}$$

3. $\theta^{(0)} = 50$ から計算する： 80、89、91.7、92.51、……
4. 収束先を求める： $\theta = (650 + 3\theta)/10$ を解いて、$\theta = 650/7 \approx 92.9$

答えは約92.9時間です。これは、総観察時間650時間を、切れた個数7で割った値です。打ち切りのある指数分布の最尤推定量と一致しています。

2. なぜ尤度が増えるのか

EMアルゴリズムの大切な性質は、更新のたびに、観測データの対数尤度が減らない ことです。導出の要点を確認しましょう。

条件つき確率の定義から $f(Y_{\text{obs}} \mid \theta) = f(Y_{\text{obs}}, Y_{\text{mis}} \mid \theta) / f(Y_{\text{mis}} \mid Y_{\text{obs}}, \theta)$ です。対数をとり、$\theta^{(t)}$ のもとで $Y_{\text{mis}}$ について期待値をとると、次の式が得られます。

$$\log f(Y_{\text{obs}} \mid \theta) = Q(\theta \mid \theta^{(t)}) - H(\theta \mid \theta^{(t)})$$

ここで $H(\theta \mid \theta^{(t)}) = E[\log f(Y_{\text{mis}} \mid Y_{\text{obs}}, \theta) \mid Y_{\text{obs}}, \theta^{(t)}]$ です。イェンセンの不等式から、$H$ は $\theta = \theta^{(t)}$ で最大になります。つまり、Mステップで $Q$ が増えれば、$H$ は増えないので、差である対数尤度は必ず増えるか同じままです。

一方で、弱点もあります。収束先は局所的な最大値や鞍点のこともあるので、初期値を変えて試します。収束が遅い場合もあります。標準誤差は自動では出ないので、別の計算が必要です。

補足（ここは読み飛ばしてOK）： EMアルゴリズムは、欠測がない場面でも使えます。例えば、2つの正規分布の混合モデルでは、「どちらの分布から来たか」というラベルを欠測とみなします。Eステップで所属の確率を計算し、Mステップで重みつきの平均と分散を計算します。

試験のポイント：Eステップは完全データの対数尤度の条件つき期待値の計算、Mステップはその最大化です。対数尤度は単調に増加しますが、大域的な最大値に届く保証はありません。

以上が、EMアルゴリズムについての解説です。

細かい式を忘れてしまっても問題ありません。ここでは、「まず空欄の理由を考える」「平均値で埋めるとばらつきが縮む」「EMは期待して最大化を繰り返す」という点だけ押さえておきましょう。

今回のまとめ

空欄の理由を考えてから、消す・埋める・EMを選ぶ

それでは、今回の内容を振り返ります。

今回は、不完全データの扱いについて紹介しました。空欄を消すことも埋めることも、裏では「空欄の理由」について何かを仮定しています。手法を選ぶ前に、その仮定が現場に合っているかを考えることが大切です。

今回のポイントを振り返りましょう。

振り返り

• 欠測メカニズムはMCAR、MAR、MNARの3つ。見えない値そのものが理由ならMNAR。
• MARなら欠測の仕組みを無視して、観測データの尤度で最尤推定ができる。
• 削除法は手軽だが、データが減り、MCARでなければ偏る。
• 補完法のうち平均値代入は、分散と相関を小さくする。多重代入は不確かさを $B$ として足す。
• EMアルゴリズムは下書きと清書の繰り返しで、対数尤度を単調に増やす。

キーワードのおさらい

• 欠測: 本来あるはずの値が手に入っていないこと
• 欠測メカニズム: どの値が欠測するかを決める確率の仕組み（MCAR、MAR、MNAR）
• 削除法: 欠測のある個体を分析から除く方法（リストワイズ、ペアワイズ）
• 補完法（代入法）: 欠測をもっともらしい値で埋めてから分析する方法
• 多重代入法: 補完を $m$ 回行い、結果をルービンのルールでまとめる方法
• EMアルゴリズム: Eステップ（条件つき期待値）とMステップ（最大化）を繰り返す最尤推定の方法

さて、今回まででさまざまなモデルの当てはめ方を学びました。でも、当てはめられるモデルはいくつもあります。変数を増やせば、データにはいくらでも合わせられそうです。どのモデルを選べばよいのでしょうか？

次回は、「AIC」「BIC」「クロスバリデーション」によるモデル選択について紹介します。

それでは、次の記事へ進んでいきましょう。