🔖 キーワード索引

本ページは 欠損メカニズム（Missing Data Mechanism）を多角的に解説します。 上のチップは、 検索・関連語の手がかりです。

🔖 索引・このページの読み方

本ページは 欠損メカニズム（Missing Data Mechanism）を 12+ セクションで多角的に解説します。

📂 分野：データ処理 ／ 関連語：MCAR / MAR / MNAR / 多重代入法 / Heckman

🎨 直感で掴む — 具体例で理解する

3 つのメカニズムを具体例で：

例：健康診断データで 「肥満気味の人ほど体重欄を空欄」にすると、 これは MNAR。 単純に「他の人の体重平均」で埋めると、 全体の平均体重を 過小評価してしまいます。

🎨 直感で掴む（詳説）

「データが穴あきになる原因」を 3 つに分類します。

🎭 比喩 — アンケートの茶封筒

1000 通配って 200 通無回答のとき：

• MCAR：郵便事故でランダム紛失 → 残り 800 で不偏推定
• MAR：老人ホーム居住者が無回答（年齢で説明可） → MICE で復元
• MNAR：高所得者が拒否 → 外部情報なしには補正不能

仮定の強さ：MCAR ⊂ MAR ⊂ MNAR。 強い仮定ほど分析は楽だが、 現実は MNAR 寄り。

📐 定義・数式

3 メカニズムの数学的定義（$Y$ = 観測したい値、 $M$ = 欠損指示子、 $X$ = 他の観測変数）：

📐 数式・定義（詳説）

$Y$ = 観測したい変数、 $M$ = 欠損指示子、 $X$ = 補助変数。 各メカニズムは $P(M\mid Y,X)$ の依存性で定義：

MICE が正当化されるのは MAR まで。 MNAR では Heckman 選択モデルが必要。

尤度関数

完全データを $(Y, M)$ とすると、 $p(Y_{obs}, M \mid \theta, \psi) = \int p(Y \mid \theta) p(M \mid Y, \psi) dY_{mis}$。 MAR では $\psi$ を無視できる（ignorable）。

📐 MCAR / MAR / MNAR の数式定義（Rubin の枠組み）

完全データを $Y=(Y_{\text{obs}}, Y_{\text{mis}})$、 欠損インジケータを $R$（$R_i=1$ なら観測、 $0$ なら欠損）、 補助変数を $X$ とおく。 Rubin (1976) は欠損メカニズムを次の 3 段階で分類した。

MCAR (Missing Completely At Random)

欠損確率がデータの値にも観測値にも一切依存しない。 例: アンケート用紙を風で飛ばされた、 サーバが落ちて記録できなかった。 観測されたデータは全データの単純無作為標本になる。

MAR (Missing At Random)

欠損確率は 観測されている変数 には依存するが、 欠損している値自体 には依存しない。 例: 「男性の方が収入を答えないことが多い」が、 答えなかった男性の収入分布は答えた男性と同じ。

MNAR (Missing Not At Random)

欠損確率が 欠損している値そのもの に依存する。 例: 「収入が高い人ほど収入を答えない」。 観測データだけからは厳密に検証できず、 感度分析が必須。

数式を言葉で読み解く: なぜ MAR は「実用上 OK」と言われるのか

$P(R\mid Y_{\text{obs}}, Y_{\text{mis}}, X) = P(R\mid Y_{\text{obs}}, X)$ という条件は、 補助変数 $X$ と観測値 $Y_{\text{obs}}$ で条件づけたときに欠損が「もはやランダム」になることを意味する。 だから、 「条件付け先の変数」さえ十分にモデルへ入れれば、 多重代入や最尤法（FIML）で 不偏推定が可能 になる。 これが MICE/Amelia/mice パッケージが「MAR を仮定する」と書く理由。

📐 Rubin の理論 — 「無視可能性」「適切性」と Bayes 整合

Rubin (1976, Biometrika) が定義した欠損メカニズムの枠組みは、 単なる分類ではなく「いつ欠損プロセスを無視できるか」を厳密に与える 無視可能性 (ignorability) の理論である。 ここでは数式と言葉の対応を最小限に絞って読み解く。

完備データと観測パターン

完備データを $Y = (Y_{obs}, Y_{mis})$、 欠損指示子を $R$（$R_{ij}=1$ なら $Y_{ij}$ が観測）とする。 欠損メカニズムは条件付き分布 $f(R \mid Y, \phi)$ で定義され、 パラメータ $\phi$ を持つ。 ここで重要なのは、 $f$ が $Y_{mis}$ にどう依存するか である。

3 メカニズムの厳密定義

無視可能性 (ignorability) の定義

尤度ベース推論で欠損プロセスを「無視できる」とは、 (A) MAR が成り立ち、 (B) $\theta$（興味のあるパラメータ）と $\phi$（欠損機構パラメータ）が distinct（分離可能）であること。 この 2 条件下で観測尤度は $L(\theta \mid Y_{obs}) = \int f(Y_{obs}, Y_{mis} \mid \theta) dY_{mis}$ で書け、 $\phi$ を明示的に推定しなくてよい。

Bayesian の場合さらに事前 $p(\theta, \phi) = p(\theta) p(\phi)$ の分離（priori independence）が無視可能性の前提に追加される。 実務的にはこれが厳しく、 弱情報事前でも distinct でない場合があるため、 感度分析を併用する。

適切性 (proper) と MI の正当化

Rubin の「適切な多重代入 (proper MI)」は次の 3 条件を要求する: (1) 各代入が posterior 予測分布からの draw、 (2) m → ∞ で pooled 推定が完備データ ML と一致、 (3) Rubin の variance 結合式が完備データ漸近分散の 不偏推定になる。 sklearn の IterativeImputer(sample_posterior=True) は (1) を満たし、 (2)(3) は大標本で近似的に満たす。

Rubin の分散統合式（再掲）

$T = U + (1 + 1/m) B$、 ただし $U$ は within-imputation 分散の平均（各補完での通常標準誤差の二乗の平均）、 $B$ は between-imputation 分散（補完間の点推定のばらつき）、 $m$ は代入回数。 自由度は $\nu = (m-1)(1 + U / ((1+1/m) B))^2$ で、 m が大きく B が小さいと自由度は無限大に近づき正規近似で良い。 m が小さいと t 近似が必要。

💡 覚えるべき相対関係: $B/U$ 比は「欠損による情報損失の指標」、 $\lambda = (1+1/m)B / T$ は「全分散に占める欠損由来の割合」。 $\lambda > 0.5$ なら m を増やす、 ドメイン知識を追加する、 または対象パラメータの推定をあきらめる、 のいずれかが必要。

📝 論文・レポートでの記述テンプレート（欠損メカニズム編）

学術論文・社内分析報告書・行政白書のいずれでも、 欠損処理は「再現可能なレベル」で記述することが求められる。 以下は STROBE / CONSORT / Sterne et al. (2009, BMJ) のガイドラインに準拠したテンプレート。 そのまま雛形として使えるよう、 SSDSE 文脈での記述例を併記する。

テンプレート 1: Methods 節（欠損処理の宣言）

テンプレート 2: Results 節（欠損診断の報告）

テンプレート 3: Discussion 節（限界の明示）

テンプレート 4: 補足表（Supplementary Table）

💡 査読対策: 査読者からの最頻指摘は「MNAR を検討したか」「m が小さすぎる」「Little's test だけで判定していないか」「complete case との比較がない」の 4 点。 上記テンプレートを使うとこの 4 点すべてに先回り回答できる。

🎓 発展トピック — 縦断・パネル・因果推論での欠損

縦断データ（longitudinal）の欠損

繰返し測定では 3 種類の欠損形態を区別する: (a) intermittent（一時的に欠損後復帰）、 (b) dropout / attrition（脱落以降全欠損）、 (c) death（観測機会自体が消失）。 dropout の場合、 「dropout 確率が過去の観測値で説明できる」なら MAR、 「dropout 時点の未観測値で説明される」なら MNAR。 mixed-effects model + FIML は MAR を仮定するため、 MNAR が疑われるときは pattern-mixture / shared-parameter model に切り替える。

パネルデータの attrition と inverse-probability weighting

家計パネル調査などで「離脱確率を観測共変量で予測」し、 残った個体に 1/離脱確率 の重みをかけて推定する IPW (Inverse Probability Weighting) は、 MAR 仮定下で不偏推定を得る古典的手段。 doubly robust estimator（IPW + アウトカム回帰モデルの両方を併用）は片方のモデルが正しければ不偏で、 実務では IPW + MI のハイブリッドが推奨される。

因果推論との接続: 欠損 = 観測されない潜在アウトカム

ランダム化比較試験 (RCT) では「処置群の潜在アウトカム」と「対照群の潜在アウトカム」のうち、 各個体で一方が常に欠損する。 ATE（平均処置効果）の推定はこの欠損を「ランダム化により MCAR」と見做して可能になる。 観察研究での confounder 調整は MAR 化のプロセスと等価で、 propensity score matching は欠損アウトカムの補完を間接的に行っているとも解釈できる。 この視点は Rubin の「potential outcome framework」と直結する。

機械学習での欠損

XGBoost / LightGBM などツリーモデルは欠損を「default direction」として学習でき、 単純な mean imputation より高精度になる場合が多い。 しかしこれは 予測タスク での話で、 因果推定や統計推論では Rubin の枠組みに従った MI が依然として必要。 「予測精度が良いから補完は不要」という主張は 推論目的では誤り。

深層学習での欠損: GAIN, MIWAE

GAIN (Yoon et al., 2018) は GAN を使った補完で MAR 下で MICE を上回ることが報告されている。 MIWAE (Mattei & Frellsen, 2019) は変分オートエンコーダによる補完で MNAR にも一定の頑健性を示す。 ただし両者とも「不確実性の定量化」は工夫が必要で、 Rubin の規則に従った標準誤差結合は別途実装する必要がある。

💡 発展まとめ: 欠損は「データの問題」ではなく 研究設計の一部。 設計段階で欠損メカニズムを予測し、 補助変数を用意し、 縦断追跡で attrition 対策を埋め込むのが「欠損に強い研究」の本質。 補完手法の選択はその後の話に過ぎない。

❓ よくある誤解と追加 FAQ

🪜 欠損メカニズム判定 — 12 ステップの実務マニュアル

新しいデータセットを受け取り、 「欠損があります」と気づいた瞬間から「補完済みデータで本解析を始める」までの 12 ステップを並べる。 各ステップで「何を見て」「何を記録するか」を具体化したので、 そのままチェックリストとして使える。

1. 欠損率の表化: df.isna().mean() で列ごと、 df.isna().sum(axis=1) で行ごとの欠損数を確認。 欠損率 0% の列と 30% 超の列を分類。
2. 欠損パターン可視化: missingno.matrix(df) で「どの行のどの列が欠損しているか」を画素的に確認。 縦縞 → 全変数同時欠損、 散在 → 個別欠損。
3. 共変量との関連検定: 各欠損列について「欠損 vs 観測共変量」のロジスティック回帰。 係数が有意 → MCAR ではない。
4. 欠損あり/なし群の平均比較: t 検定または Mann-Whitney で 2 群の他変数の差を確認。 有意差 → MAR 以上。
5. Little's MCAR test: statsmodels 公式実装は無いが、 R の misty::na.test() や Python の impyute で代替。 棄却 → MCAR ではない。 棄却なし ≠ MCAR。
6. ドメイン知識の聞き取り: データ取得担当者・調査企画者に「欠損の典型理由」を 30 分インタビュー。 これが MAR/MNAR 判定の最大根拠になる。
7. 補助変数の探索: 「欠損を説明しそうな変数」が他にないか検討。 ある場合は補完モデルに必ず含める。
8. MAR 仮定の宣言: ここまでの証拠を文書化し「MAR を仮定する」を明示。 後の感度分析で MNAR を検討する。
9. MI 実行: m = max(20, 100 × γ) で IterativeImputer(sample_posterior=True) を回す。 各補完で本解析を実行。
10. Rubin 規則で統合: 点推定、 標準誤差、 自由度を統合。 B/U 比と λ を報告。
11. complete case との対照: 同じ解析を listwise でも実行。 結果が大きく違うなら MI の根拠が強化される。
12. MNAR 感度分析: δ-adjustment や pattern-mixture で MNAR シナリオを検討。 結論の頑健性を報告。

💡 運用 Tips: 上記 12 ステップを Jupyter Notebook の章立てとして並べておくと、 新規データセットに対して「セル順に走らせるだけで再現可能な欠損診断レポートが出る」状態を作れる。 SSDSE 系の公的統計、 医療コホート、 マーケティング調査いずれにも転用可能。

🚫 アンチパターン集 — やってはいけない 8 つのこと

1. 欠損行を黙って削除: 欠損率と削除後の標本サイズを必ず記録。 削除しただけで論文に書かない、 はインテグリティ違反。
2. 平均代入で点推定だけ報告: 標準誤差が過小で偽の有意性が出る。 mean imputation を使うなら理由と注意を明記。
3. カテゴリ変数に「Unknown」を新カテゴリ化: 一見便利だが、 元のカテゴリ分布が歪む。 統計推論では特に危険。
4. m=1 の単一代入: MI のメリット（不確実性の定量化）が消える。 m ≥ 5、 推奨は m ≥ 20。
5. MAR 仮定を根拠なく宣言: 「MAR とした」とだけ書く論文が多いが、 ロジスティック診断と感度分析で裏付けを示さないと査読落ち。
6. MNAR 感度分析を省略: 「MAR を仮定したので感度分析は不要」は誤り。 MAR が間違っていたときの影響評価が必須。
7. 補完後の値を「実測値」として扱う: 補完値の信頼区間は実測より広いことを忘れる。 「補完済み」と明記する。
8. 機械学習の予測精度だけで補完法を選ぶ: 推論目的では Rubin 規則に従えるかどうかが第一基準。 RMSE 最小化ではない。

⚠️ 業務で見かけるパターン: 「欠損は KNNImputer で埋めたので大丈夫」だけの記述。 これは 補完手法の選択理由、 メカニズムの根拠、 感度分析 がすべて欠けている。 査読・監査で必ず指摘される。

📜 欠損メカニズム理論の歴史 — Rubin 以前と以後

欠損データ問題は古くから認識されていたが、 1970 年代の Rubin 論文以前は「経験則ベースの補完」が主流だった。 ここでは主要な転換点を年表で示す。

Rubin の 1976 論文以前は、 「欠損は厄介だが何とか埋めれば良い」程度の認識が一般的で、 mean imputation や last-observation-carried-forward (LOCF) が標準だった。 Rubin の貢献は (1) 欠損プロセス自体を確率モデルとして扱う視点を導入したこと、 (2) 「いつ無視できるか」の数学的条件を与えたこと、 (3) 多重代入で不確実性を定量化する道筋を開いたこと、 にある。 これは現代統計学・データサイエンス全体に影響を与え、 観察研究の因果推論や Bayesian inference の発展とも交差している。

💡 歴史から学ぶ教訓: 「最先端の補完手法（GAIN, MIWAE）を使えば良い」のではなく、 Rubin の理論枠組みの上で 不確実性を正しく定量化 することこそが本質。 1976 年の枠組みは 50 年経った今も全く色褪せず、 最新の深層生成モデルも結局この理論で評価される。

⚖️ 主要補完手法の徹底比較 — メカニズム別の妥当性マトリクス

補完手法は 10 種類以上あり、 メカニズム別に妥当性が大きく異なる。 「とりあえず KNN」ではなく、 「データのメカニズム」と「解析目的（推論 vs 予測）」の組合せで選ぶ。

○ = 不偏で適切な不確実性定量化が可能、 △ = 条件付きで可、 × = 偏り or 不確実性過小評価。 MNAR の列は「感度分析と組み合わせるか、 欠損機構を明示モデル化するか」が必要なため、 完全な ○ は無い。

💡 選択フローチャート: (1) 推論目的か予測目的か → 推論なら MI / FIML、 予測なら ML 系も可。 (2) 欠損機構は MAR で十分か → MNAR 疑いなら pattern-mixture を追加。 (3) データ規模 → 小規模 (~1000 行) なら MI、 大規模なら MICE + 並列化または深層補完。 (4) カテゴリ変数の有無 → 多いなら R mice、 連続中心なら sklearn。

🗂 用語整理 — 欠損関連の必須概念ミニ辞書

論文や教科書で頻出する欠損関連の専門用語をまとめる。 各用語の英語表記・略語・日本語訳を併記し、 さらに「いつ使う言葉か」を補足した。

💡 用語の使い分けポイント: 統計系論文では Rubin の MCAR/MAR/MNAR、 計量経済では「missing-not-at-random」「selection on observables」、 機械学習では「missingness as a feature」「imputation」などの用語がよく見られる。 同じ概念でも分野ごとに表現が異なるため、 査読時は用語統一に注意。

補足: 関連する近接概念

truncation（切断）とcensoring（打ち切り）は欠損とは区別される。 truncation は「特定の条件を満たす個体しか観測されない」（例: 所得 0 円以下は標本に入らない）状況で、 標本選択バイアスを引き起こす。 censoring は「観測はあるが値が一部しか分からない」（例: 生存時間で追跡終了時点まで生存していた個体）。 これらは別個の統計理論（survival analysis, Tobit model）で扱われ、 単純に MI を当てるのではなく専用手法を選ぶこと。

measurement error（測定誤差）と欠損も近縁概念で、 「観測されているが真値とずれている」のが測定誤差、 「観測されていない」のが欠損。 ただし測定誤差の極端版（誤差が無限大 = 観測不能）として欠損を扱う統一的枠組みもある。

🔬 記号・要素の読み解き

$Y$ — 観測したい完全データ

$Y = (Y_{obs}, Y_{mis})$ と分解。

$M$ — 欠損指示子

$M_{ij}=1$ で欠損。

$X$ — 補助変数

MAR では「条件付ければ $Y$ と独立」。

$\theta$ — 推定対象

母平均、 回帰係数など。

$\psi$ — 欠損モデル係数

MNAR で同時推定が必要。

ignorable

MCAR・MAR は無視可、 MNAR は不可。

識別性

観測データのみで MAR/MNAR を統計検定不能。

Rubin's rules

多重代入を統合する公式。

🔬 多重代入 (MI) と MNAR 感度分析 — Rubin の規則で結論を頑健化

MAR の妥当性が確認できたら多重代入 (Multiple Imputation, MI) を実行する。 MI は欠損値を 1 つの値ではなく m 個の確率的代入で補い、 各補完データセットで解析し、 Rubin の規則で統合する。 これにより「欠損による不確実性」が標準誤差に正しく反映される。

ステップ 3: scikit-learn IterativeImputer で m=20 の MI

このコードでやること: sklearn.experimental.IterativeImputer を sample_posterior=True で起動し、 20 個の補完版を生成。 各版で「人口 → 貯蓄」の回帰を当て、 Rubin の規則で統合する。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

import numpy as np
import pandas as pd

# df_mnar はこのあとのブロックで作っているので、ここでも用意しておく
# MNAR: 値そのものの大きさが欠測しやすさを決める（消費支出が低い県ほど欠けやすい）
_rng = np.random.default_rng(0)
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)
df_mnar = df.copy()
_p = 1 - (df_mnar['L3221'].rank(pct=True))     # 低いほど欠測確率が高い
df_mnar.loc[_rng.random(len(df_mnar)) < _p * 0.4, 'L3221'] = np.nan

from sklearn.experimental import enable_iterative_imputer  # noqa
from sklearn.impute import IterativeImputer
import statsmodels.api as sm
import numpy as np
import pandas as pd

cols = ['A1101','A1301','A4103','L3221']
data = df_mnar[cols].values

m = 20
betas, ses = [], []
for seed in range(m):
    imp = IterativeImputer(sample_posterior=True, random_state=seed, max_iter=20)
    completed = imp.fit_transform(data)
    cdf = pd.DataFrame(completed, columns=cols)

    X = sm.add_constant(cdf[['A1101']])
    y = cdf['L3221']
    res = sm.OLS(y, X).fit()
    betas.append(res.params['A1101'])
    ses.append(res.bse['A1101'])

beta_bar = np.mean(betas)
U = np.mean(np.array(ses)**2)             # within-imputation 分散
B = np.var(betas, ddof=1)                 # between-imputation 分散
T = U + (1 + 1/m) * B                     # 総分散
se_total = np.sqrt(T)
print(f'pooled beta = {beta_bar:.4f}')
print(f'pooled SE   = {se_total:.4f}')
print(f'between/within 比 = {B/U:.3f}')
print(f'欠損による情報損失率 λ = {(1 + 1/m) * B / T:.3f}')

📤 実行例:

💬 pooled β = 0.0026 は「総人口 A1101 が 1 人増えると消費支出 L3221 が 0.0026 円/月 増える」＝ 人口 100 万人あたり月 2,600 円という関係を示す（どちらも生の単位のままなので係数が小さく出る）。 pooled SE = 0.0013 なので t = 0.0026/0.0013 ≒ 2.0 で、 かろうじて有意水準 5% に届く程度。 between/within 比 0.155 は欠損由来の分散が within 分散の 15.5% を上乗せしていることを意味し、 これを無視（= 単一代入）すると標準誤差を 7% 程度過小評価する。 λ=0.140 は「総分散の 14% が欠損由来」という Rubin の指標で、 0.5 を大きく下回るので m=20 の代入回数で十分。

ステップ 4: MNAR 感度分析 — δ-adjustment

このコードでやること: 各補完値に「欠損値は MAR 補完よりも δ だけ大きい」というシフトを加え、 δ を [-1, +2] の範囲で動かして結論（β の符号と有意性）の頑健性を確認する。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

import numpy as np
import pandas as pd

# 各 δ に対して MI 後 OLS を再計算
deltas = np.arange(-1.0, 2.01, 0.5)
result = []
for delta in deltas:
    bs = []
    for seed in range(m):
        imp = IterativeImputer(sample_posterior=True, random_state=seed, max_iter=20)
        completed = imp.fit_transform(data)
        cdf = pd.DataFrame(completed, columns=cols)
        # MNAR シフト: 元欠損行の貯蓄に δ × σ を加える
        sigma = cdf['L3221'].std()
        mask = df_mnar['L3221'].isna().values
        cdf.loc[mask, 'L3221'] += delta * sigma
        X = sm.add_constant(cdf[['A1101']])
        res = sm.OLS(cdf['L3221'], X).fit()
        bs.append(res.params['A1101'])
    result.append((delta, np.mean(bs)))

for d, b in result:
    print(f'δ = {d:+.1f} → pooled beta = {b:+.4f}')

📤 実行例:

💬 δ を -1σ から +2σ まで動かしても β はすべて正の符号を保ち、 値も 0.0024-0.0031 の範囲（δ=0 の 0.0026 に対して ±20% 程度の振れ幅）。 結論「人口が大きい県ほど消費支出が多い」は MNAR 仮定の幅広いシナリオで頑健。 報告には「δ ∈ [-1, +2] σ で β > 0 が保たれた」と書けば査読でも説明可能。

💡 感度分析の報告ルール: 「δ ∈ [a, b] で結論が変わらない」と明示し、 結論が反転する tipping point（符号が逆転する δ）を併記すること。 査読者・規制当局が安心する。

🌍 分野別ケースで掴む — メカニズムは分野で典型形が決まる

欠損メカニズムは「データの性格」と密接に関連する。 同じ「MAR」でも、 医療と金融、 教育と Web 行動ログでは 典型的な発生機構 が異なる。 分野別パターンを掴むことで、 新しいデータセットに出会ったとき「最初にどこを疑うか」が即座にわかる。

事例 1: 公的統計（SSDSE-B-2026 風）

SSDSE-B-2026 自体には欠損はほぼ無いが、 同種の「県別経済指標」では 小規模県のサンプル数不足が原因で一部の細目（例: 個別業種の出荷額）が「-」記号で欠損化されることがある。 これは「人口」という観測変数で説明できるため MAR と扱える。 公的統計では「秘匿」と「未集計」の区別が重要で、 メタデータ（注記）を必ず読むこと。

事例 2: 医療コホート研究

5 年追跡の高血圧コホートで「血圧の値」が脱落者で欠損したとする。 「血圧が高い被験者ほど脱落しやすい」と判断したい場合、 直接は確認できないが、 (a) 脱落直前の血圧と非脱落者の血圧分布の比較、 (b) 「脱落 vs 観測変数」のロジスティック回帰、 (c) 主治医インタビューでの「来院理由」確認、 を組み合わせて MAR/MNAR を判定する。 MNAR が疑われる場合は δ-adjustment や pattern-mixture を必ず提示する。

事例 3: 金融・年収アンケート

家計調査で「年収」項目の無回答率が約 15-25% に達することは珍しくない。 高所得層と無所得層の 両側 で欠損が増える「U 字パターン」が観測される。 この場合、 線形ロジスティックでは説明できないため、 二次項を入れたモデルや有限混合モデルで欠損機構を推定し、 selection model 系の MNAR フレームで再解析する。

事例 4: IoT センサー

温度センサーが 40℃ を超えると誤動作で NaN を返す仕様があったとする。 これは「観測温度（直前の値）」で説明できる MAR。 隣接時刻の値で Kalman smoother 補完するのが標準。 一方、 機械振動でセンサー自体が故障し全期間欠損になる場合は MNAR に近いため、 別個のセンサー（冗長系）からの補完か当該機械の除外を検討する。

事例 5: パネル調査の attrition

JHPS（日本家計パネル調査）のような複数年追跡では、 ベースライン特性（年齢・初期所得・地域）が attrition の予測因子となることが多く、 MAR 仮定が概ね妥当。 ただし「政治意識調査」など態度変化と離脱が連動するケースでは MNAR を疑い、 inverse-probability weighting (IPW) や joint model を組み合わせる。

💡 分野横断の鉄則: 欠損メカニズムの仮定は「データだけ」では決められない。 必ず ドメイン知識（業務担当者・調査企画者・センサー仕様）と組み合わせて判定する。 純粋に統計的検定だけで MAR を主張するのは危険。

🧮 数値例・実値計算

所得アンケート（真平均 650万、 高所得 MNAR）：

Rubin's rules で 3 個の補完統合

$\hat\theta_1=540, \hat\theta_2=555, \hat\theta_3=548$、 各分散 $U_k=900,950,920$。

$\bar\theta = 547.67$、 $\bar U = 923.3$、 $B = 56.3$、 $T = \bar U + (1+1/m)B = 998.4$、 $SE \approx 31.6$、 95% CI = $(485.7, 609.6)$

🧮 SSDSE で 3 メカニズム比較（実値シミュレーション）

🎯 このコードでやること: SSDSE-B-2026 の「総人口」を完全データとみなし、 MCAR / MAR / MNAR の 3 種類の欠損メカニズムを意図的に注入し、 mean imputation の平均がどれだけバイアスを受けるか比較する。

📥 入力: SSDSE-B-2026、 総人口(A1101)・15歳未満人口(A1301) の 2 列（総人口を欠損させ、 15歳未満人口を欠損確率のドライバに使う）。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

import pandas as pd
import numpy as np
from scipy.stats import bernoulli

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1]); df = df[df['SSDSE-B-2026'] == 2023].dropna()
pop  = df['A1101'].values
gdp  = df['A1301'].values
true_mean = pop.mean()

rng = np.random.default_rng(0)

# MCAR: 50% を完全無作為に欠損
mcar = pop.copy().astype(float)
mcar[rng.random(len(pop)) < 0.5] = np.nan

# MAR: 観測変数 A1301(15歳未満人口) が中央値以上の県ほど欠損しやすい
mar  = pop.copy().astype(float)
prob = 0.2 + 0.6 * (gdp > np.median(gdp))
mar[rng.random(len(pop)) < prob] = np.nan

# MNAR: 人口が大きいほど欠損しやすい（値自体に依存）
mnar = pop.copy().astype(float)
prob = 0.2 + 0.6 * (pop > np.median(pop))
mnar[rng.random(len(pop)) < prob] = np.nan

for name, x in [('MCAR', mcar), ('MAR', mar), ('MNAR', mnar)]:
    obs   = x[~np.isnan(x)]
    bias  = obs.mean() - true_mean
    print(f'{name:5s}: 残存={len(obs):2d}  観測平均={obs.mean():,.0f}  真値={true_mean:,.0f}  bias={bias:+,.0f}')

📤 実行例:

💬 結果の読み方: MCAR は 50% を無作為に欠損させても観測平均が真値の周りに散らばるだけで、 バイアスは相対的に小さい（ここでは +255,451。 n=47 の標本変動によるもので、 n が増えれば 0 に近づく）。 一方 MAR・MNAR では「値の大きい県」が体系的に消えるため観測平均が大きくバイアスする。 今回は観測変数（15歳未満人口 A1301）が総人口と強く相関するため MAR のバイアスが最大になったが、 本質は MAR は補助変数（15歳未満人口）で補正できるのに対し、 MNAR は欠損値自身に依存するため補助変数では打ち消せない点にある。

🧮 SSDSE-B-2026 で欠損メカニズム判定の実値ワークフロー

ここでは SSDSE-B-2026（47 都道府県 × 統計指標）を素材に、 「欠損列が来たら次に何をするか」を 6 ステップで実行する。 全ステップ実値ベースで、 教科書的に「MAR を仮定して MI」と決めつけない手順を示す。

ステップ 0: 入力と前提

このコードでやること: SSDSE-B-2026 を読み込み、 まず欠損率と欠損パターンを表化する。 仮想的に「消費支出(L3221)の上位 5 県だけ欠損」を作って MNAR っぽい状況を再現する（実データ自体は完備しているため、 教材として人工的に欠損を入れる）。

📥 入力データ（SSDSE-B-2026 抜粋、 単位は元データに準拠）:

ステップ 1: 欠損率と欠損パターン表

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

import pandas as pd
import numpy as np

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', skiprows=[1], encoding='cp932')
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 2023 年の 47 都道府県に絞る

# 教材用: 消費支出(L3221)の上位 5 県を欠損化（MNAR シナリオ）
top5 = df.nlargest(5, 'L3221').index
df_mnar = df.copy()
df_mnar.loc[top5, 'L3221'] = np.nan

# 欠損率の確認
missing_rate = df_mnar.isna().mean().sort_values(ascending=False)
print('--- 欠損率 ---')
print(missing_rate[missing_rate > 0])
print()
print('--- 欠損行の他変数の要約 ---')
print(df_mnar[df_mnar['L3221'].isna()][['A1101','A1301','A4103']].describe())

📤 実行例:

💬 欠損率は約 10.6%（5/47）。 欠損 5 県の総人口(A1101)・15歳未満人口(A1301) は 全 47 県の平均より明らかに大きい（例: 欠損群の総人口平均 ≈ 520 万人）。 つまり「欠損が観測変数（人口規模）と相関している」サインがあり、 単純な MCAR ではない。 まず MAR を疑い、 さらに「消費支出そのものが大きい県で欠損」という前提を観測データからは確認できないため MNAR の可能性も保持する。

ステップ 2: ロジスティック回帰で MAR シグナルを定量化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

from sklearn.linear_model import LogisticRegression
import numpy as np

# 欠損ダミー
y = df_mnar['L3221'].isna().astype(int)
# 観測変数のみで予測（標準化）
X = df_mnar[['A1101','A1301','A4103']].copy()
X = (X - X.mean()) / X.std(ddof=0)

clf = LogisticRegression(max_iter=1000, penalty='l2', C=1.0)
clf.fit(X, y)

print('係数 (z-score 単位):')
for name, coef in zip(X.columns, clf.coef_[0]):
    print(f'  {name:12s}: {coef:+.3f}')
print(f'切片        : {clf.intercept_[0]:+.3f}')
print(f'pseudo R^2  : {clf.score(X, y):.3f}')

📤 実行例:

💬 総人口(A1101) の係数が正（+0.284）で、合計特殊出生率(A4103) が負（−0.448）→ 観測変数で欠損傾向がある程度説明でき、 MAR と矛盾しない。 15歳未満人口(A1301) は +0.160 と小さく、 総人口とほぼ同じ情報を持つため単独では効きにくい。 ここでの pseudo R² は clf.score（= 正解率）で 0.915 だが、 47 県中 42 県が「欠損なし」なので全部「欠損なし」と答えるだけで 0.894 になる。 つまり 0.915 という値はほとんど下駄で、 説明力が高く見えること自体は MAR の証明にはならない点に注意。 さらに「消費支出そのもの」を説明変数に入れていない（欠損しているため観測できない）ので、 MNAR を完全には排除できない。 次のステップで多重代入を実行し、 後に MNAR 感度分析を回す。

🧮 数式に値を入れて手で計算する: MCAR/MAR/MNAR の判定

合成データで欠損パターン別に平均推定値を計算する。

Step 1: シナリオ別の欠損率

Step 2: バイアス

🐍 Python で再現

1
2
3
4

truth = 50
obs = {'MCAR':49.8, 'MAR':49.2, 'MNAR':40.0}
for k, v in obs.items():
    print(f"{k}: バイアス={v-truth}")

📤 実行結果

💬 手計算 (Step 2) と Python 出力が完全一致。 MNAR が最も深刻。

🧮 SSDSE-B-2026 実値で Step 1〜5 — 東北 6 県 食料費で 3 メカニズムを再現

合成データではなく、 SSDSE-B-2026 の東北 6 県 食料費 (千円/月) で MCAR/MAR/MNAR の観測平均バイアスを 手で求めて Python と一致確認する。

Step 1: 真の値 (全 6 県観測できた場合の平均)

Step 2: MCAR — ランダムに偶数番目 (秋田・岩手・山形) を欠損

Step 3: MAR — 「観測される県名」に依存。 北寄り 2 県 (青森・秋田) のみ欠損

Step 4: MNAR — 「食料費そのもの」が小さいと欠損 (x < 80)

Step 5: バイアスの大小比較

🐍 Python で再現 — Step 1〜5 と一致確認

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

import numpy as np
# SSDSE-B-2026 東北 6 県 食料費 (千円/月)
prefs = ['青森', '秋田', '福島', '岩手', '宮城', '山形']
x = np.array([77.90, 78.12, 78.98, 82.00, 83.83, 84.11])
mu_true = x.mean()                                  # Step 1

mcar_mask = np.array([False, True, False, True, False, True])  # Step 2 (秋田/岩手/山形 欠損)
mar_mask  = np.array([True, True, False, False, False, False])  # Step 3 (北 2 県 欠損)
mnar_mask = x < 80                                       # Step 4 (x<80 欠損)

for name, mask in [('MCAR', mcar_mask), ('MAR ', mar_mask), ('MNAR', mnar_mask)]:
    obs_mean = x[~mask].mean()
    print(f"{name} 観測平均={obs_mean:7.4f}  バイアス={obs_mean-mu_true:+7.4f}")
print(f"真の平均 = {mu_true:.4f}")

📤 実行結果

💬 Step 2-4 の手計算 (80.2367 / 82.2300 / 83.3133) と Python 出力が小数 4 桁まで一致。 SSDSE-B の実値でも MNAR が最大の +2.49 バイアスを発生させ、 「値そのもので脱落する」メカニズムが最も補正困難であることが 実データの数値として確認できる。

📌 このページのコードで使う範囲について： 補完の比較は 最新年度の 47 行 × 7 列に絞っています。 SSDSE-B-2026 は 564 行 × 110 列あり、そのまま多重代入（MICE）にかけると 1 回の補完に数分かかってブラウザでは実行できません。 列を絞っても、MCAR / MAR / MNAR で補完のふるまいが変わるという 結論そのものは変わりません。 自分のデータで試すときは、まず小さく絞って挙動を確かめ、 そのうえで列数を増やしていくのが安全です。

🐍 Python 実装例

最小コードで動かしてみる例：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

import pandas as pd
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]        # 最新年度の 47 行だけを使う

# ── 練習用の欠損をここで作ります（架空）─────────────────────
# SSDSE-B-2026 には欠損値が 1 件もありません。そのままでは補完しても
# 何も変わらないので、このページの主題である 3 つの仕組みを人工的に作ります。
import numpy as np

_rng = np.random.default_rng(0)
_n = len(df)
# MCAR: 消費支出を、他の変数と無関係に 15% 落とす
df.loc[_rng.random(_n) < 0.15, 'L3221'] = np.nan
# MAR: 教育費を、15 歳未満人口（A1301）が中央値未満の県ほど落ちやすくする
_p = np.where(df['A1301'] < df['A1301'].median(), 0.30, 0.05)
df.loc[_rng.random(_n) < _p, 'L322108'] = np.nan
# MNAR: 住宅地の標準価格を、値そのものが高い県ほど落ちやすくする
_p = np.where(df['C5401'] > df['C5401'].quantile(0.75), 0.40, 0.05)
df.loc[_rng.random(_n) < _p, 'C5401'] = np.nan
# ────────────────────────────────────────────────────

_rate = df.isnull().mean()
print('欠損率（0 より大きい列だけ）:')
print(_rate[_rate > 0].round(3))

# 多重代入法（MICE） — MAR を仮定
# 110 列すべてを説明変数にすると 1 回の補完に数分かかる。
# 欠損のある 3 列と、それを説明できる少数の列だけに絞る。
_use = ['L3221', 'L322108', 'C5401', 'A1101', 'A1301', 'A1303', 'A4101']
imputer = IterativeImputer(max_iter=10, random_state=0)
df_filled = imputer.fit_transform(df[_use])
print('補完後の形:', df_filled.shape)

🐍 Python 実装（4+ サンプル）

① 欠損パターン可視化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

import pandas as pd
import numpy as np
import os
import missingno as msno
import matplotlib.pyplot as plt

os.makedirs('figures', exist_ok=True)   # 保存先が無いと savefig は失敗する

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', header=1)

# ── 練習用の欠損をここで作ります（架空）─────────────────────
# SSDSE-B-2026 には欠損値が 1 件もないので、可視化する対象を人工的に作ります。
import numpy as np

_rng = np.random.default_rng(0)
_n = len(df)
df.loc[_rng.random(_n) < 0.15, '消費支出（二人以上の世帯）'] = np.nan          # MCAR
_p = np.where(df['15歳未満人口'] < df['15歳未満人口'].median(), 0.30, 0.05)
df.loc[_rng.random(_n) < _p, '教育費（二人以上の世帯）'] = np.nan              # MAR
_p = np.where(df['標準価格（平均価格）（住宅地）']
              > df['標準価格（平均価格）（住宅地）'].quantile(0.75), 0.40, 0.05)
df.loc[_rng.random(_n) < _p, '標準価格（平均価格）（住宅地）'] = np.nan         # MNAR
# ────────────────────────────────────────────────────

print('欠損率:', df.isnull().mean().sort_values(ascending=False).head(10))
msno.matrix(df); plt.savefig('figures/missing_matrix.png', dpi=120)
msno.heatmap(df); plt.savefig('figures/missing_heatmap.png', dpi=120)

② 簡易 MCAR 検定

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

# ── この抜粋で使うデータを用意します（練習用の欠損は架空）──
# 110 列すべてを説明変数にすると 1 列あたり数十秒かかる。
# 最新年度の 47 行 × 欠損に関係しそうな 7 列だけで見る。
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler

_USE = ['L3221', 'L322108', 'C5401', 'A1101', 'A1301', 'A1303', 'A4101']

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023][_USE].copy()
_rng = np.random.default_rng(0)
_n = len(df)
df.loc[_rng.random(_n) < 0.15, 'L3221'] = np.nan                       # MCAR
_p = np.where(df['A1301'] < df['A1301'].median(), 0.30, 0.05)
df.loc[_rng.random(_n) < _p, 'L322108'] = np.nan                       # MAR
_p = np.where(df['C5401'] > df['C5401'].quantile(0.75), 0.40, 0.05)
df.loc[_rng.random(_n) < _p, 'C5401'] = np.nan                         # MNAR

num = df.select_dtypes(include='number')
miss_flag = num.isnull().astype(int)
for col in num.columns[num.isnull().any()][:5]:
    others = num.drop(columns=col).fillna(num.mean())
    others = StandardScaler().fit_transform(others)   # 収束を速くする
    lr = LogisticRegression(max_iter=500).fit(others, miss_flag[col])
    print(col, '説明力:', round(lr.score(others, miss_flag[col]), 3))

③ 補完 3 種比較

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

# ── この抜粋で使うデータを用意します（練習用の欠損は架空）──
# 110 列すべてを MICE に渡すと 1 回の補完に数分かかるので、
# 最新年度の 47 行 × 補完に使う 7 列だけに絞る。
import numpy as np
import pandas as pd

_USE = ['L3221', 'L322108', 'C5401', 'A1101', 'A1301', 'A1303', 'A4101']

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023][_USE].copy()
_rng = np.random.default_rng(0)
_n = len(df)
df.loc[_rng.random(_n) < 0.15, 'L3221'] = np.nan                       # MCAR
_p = np.where(df['A1301'] < df['A1301'].median(), 0.30, 0.05)
df.loc[_rng.random(_n) < _p, 'L322108'] = np.nan                       # MAR
_p = np.where(df['C5401'] > df['C5401'].quantile(0.75), 0.40, 0.05)
df.loc[_rng.random(_n) < _p, 'C5401'] = np.nan                         # MNAR
num = df.select_dtypes(include='number')

from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import SimpleImputer, KNNImputer, IterativeImputer
imp_mean = SimpleImputer(strategy='mean').fit_transform(num)
imp_knn  = KNNImputer(n_neighbors=5).fit_transform(num)
imp_mice = IterativeImputer(max_iter=10, random_state=0).fit_transform(num)
# 欠損を作る前の「本当の平均」も並べて、どの補完がそれに近いかを見る
_true = (pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
         .query('`SSDSE-B-2026` == 2023')[_USE].mean())
result = pd.DataFrame({
    '真値': _true,
    '欠損後': num.mean(),
    '平均補完': pd.DataFrame(imp_mean, columns=num.columns).mean(),
    'KNN': pd.DataFrame(imp_knn, columns=num.columns).mean(),
    'MICE': pd.DataFrame(imp_mice, columns=num.columns).mean()})
# 欠損を作った 3 列だけを見る
print(result.loc[['L3221', 'L322108', 'C5401']].round(1))

📤 実行すると次の出力が得られる：

💬 結果の読み方：3 行がそのまま 3 つの仕組みの教科書的な違いになっている。 MCAR の L3221 はどの補完でも真値 295,856 とほぼ一致（欠損があっても平均は歪まない）。 MAR の L322108 は平均補完だと 10,077 と真値 9,577 より高く出るが、 他の列を使う KNN・MICE は 9,881・9,777 まで戻せる。 MNAR の C5401 は、 値が高い県ほど欠けているので観測平均が 51,045 まで下振れし、 平均補完はそれをそのまま引きずる。 MICE は 52,846 と真値 54,638 の方向へ戻すが、 これは「他の列に価格と相関する情報が残っていた」から成功したのであり、 MNAR なら常にこう救えるわけではない点に注意。

④ 多重代入 m=20

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

# ── この抜粋で使うデータを用意します（練習用の欠損は架空）──
import numpy as np
import statsmodels.api as sm
from statsmodels.imputation.mice import MICE, MICEData

import pandas as pd

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]          # 2023 年の 47 都道府県
df = df.copy()
df['death_rate']   = df['A4200'] / df['A1101'] * 1000       # 人口千人あたり死亡数
df['aging_rate']   = df['A1303'] / df['A1101'] * 100        # 高齢化率 %
df['medical_cost'] = df['L322106'].astype(float)            # 保健医療費
num = df[['death_rate', 'aging_rate', 'medical_cost']].copy()
_rng = np.random.default_rng(0)
num.loc[_rng.random(len(num)) < 0.15, 'medical_cost'] = np.nan   # 練習用の欠損

mice_data = MICEData(num)
mice = MICE('death_rate ~ aging_rate + medical_cost', sm.OLS, mice_data)
print(mice.fit(n_burnin=5, n_imputations=20).summary())

⑤ 欠損フラグを特徴量化

1
2
3

imp = SimpleImputer(strategy='median', add_indicator=True)
X = imp.fit_transform(num)
print('元の列数:', num.shape[1], '/ 補完後:', X.shape[1])

🐍 多重代入（MICE）の実装

🎯 このコードでやること: 上記で作成した MAR 欠損データに対し、 sklearn の IterativeImputer（MICE の単純版）で多重代入を実行し、 単純な mean imputation との平均推定精度を比較する。

📥 入力: SSDSE-B-2026 で MAR 欠損を注入した df。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

import pandas as pd
import numpy as np
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer, SimpleImputer

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1]); df = df[df['SSDSE-B-2026'] == 2023].dropna()
true_mean = df['A1101'].mean()

rng = np.random.default_rng(0)
data = df[['A1101', 'A1301', 'A1303', 'A4101', 'B4101']].copy()
prob = 0.2 + 0.6 * (data['A1301'] > data['A1301'].median())
mask = rng.random(len(data)) < prob
data.loc[mask, 'A1101'] = np.nan

# (1) mean imputation
mean_imp = SimpleImputer(strategy='mean').fit_transform(data)
mean_pop = mean_imp[:, 0].mean()

# (2) iterative imputation (MICE-like)
mice_imp = IterativeImputer(max_iter=20, random_state=0).fit_transform(data)
mice_pop = mice_imp[:, 0].mean()

print(f'真値          : {true_mean:,.0f}')
print(f'mean imputation: {mean_pop:,.0f}  bias={mean_pop-true_mean:+,.0f}')
print(f'MICE (反復補完): {mice_pop:,.0f}  bias={mice_pop-true_mean:+,.0f}')

📤 実行例:

💬 結果の読み方: mean imputation は欠損していない値の平均で穴埋めするので、 MAR で大きく偏った観測平均をそのまま継承する → bias 大。 MICE は 15歳未満人口(A1301) など補助変数を使って欠損値を回帰予測するので、 MAR 仮定が満たされていれば bias が大幅に縮小する。 これが「補助変数を MICE モデルに入れる」ことの効果。

🐍 statsmodels の MICE 実装

🎯 このコードでやること: statsmodels の MICEData を使って、 m=5 個の補完済みデータセットを生成し、 各々で平均を推定して Rubin の規則でプールする「正しい多重代入」を実行する。

📥 入力: 上記と同じ MAR 欠損データ。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import pandas as pd
import numpy as np
from statsmodels.imputation.mice import MICEData

np.random.seed(0)   # MICEData は内部で np.random を使う。 ここで固定すると結果が再現する

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1]); df = df[df['SSDSE-B-2026'] == 2023].dropna()
true_mean = df['A1101'].mean()
rng = np.random.default_rng(0)

data = df[['A1101', 'A1301', 'A1303', 'B4101']].copy()
mask = rng.random(len(data)) < (0.2 + 0.6 * (data['A1301'] > data['A1301'].median()))
data.loc[mask, 'A1101'] = np.nan

mice = MICEData(data)
means, vars_ = [], []
for k in range(5):
    mice.update_all()                 # 1 サイクル分の Gibbs サンプリング
    means.append(mice.data['A1101'].mean())
    vars_.append(mice.data['A1101'].var(ddof=1) / len(data))

m = len(means)
qbar = np.mean(means)
ubar = np.mean(vars_)                 # 平均的な within-imputation variance
B    = np.var(means, ddof=1)          # between-imputation variance
T    = ubar + (1 + 1/m) * B           # Rubin の総合分散
se   = np.sqrt(T)
print(f'pooled mean = {qbar:,.0f}  真値 = {true_mean:,.0f}')
print(f'pooled SE   = {se:,.0f}  95% CI ≈ [{qbar-1.96*se:,.0f}, {qbar+1.96*se:,.0f}]')

📤 実行例:

💬 結果の読み方：MICEData は内部で np.random を使うため、 種を固定しないと実行のたびに pooled mean が 175〜197 万人の範囲で動きます。 上のコードでは冒頭に np.random.seed(0) を置いてあるので、 ここに載せた値は何度実行しても再現します。 中身を見ると pooled mean 1,908,013 は真値 2,645,809 を 大きく下回っていますが、 95% CI の上端 2,642,851 がちょうど真値に届いており、 区間としてはかろうじて真値を含みます。 これは「人口が多い県ほど欠測しやすい」MAR 構造のため点推定が下振れし、 その不確実性を Rubin の総合分散 T が幅として吸収している状態です。 ばらつきを標準誤差に反映させることが多重代入の目的であって、 ばらつきを言い訳にして値を固定しないことではない点に注意してください。

💬 結果の読み方: m=5 個の補完で得た 5 つの推定値を平均し、 内 within 分散 $\bar U$ と外 between 分散 $B$ を Rubin の公式 $T=\bar U+(1+1/m)B$ で合成。 これにより「補完自体が不確実だ」という情報が信頼区間に反映され、 単一代入よりも広い区間になる。 📝 より正確な分析（教材補足）: 今回のように MAR の依存が強く n=47 と小さいと、 プール平均（実測で 175〜197 万人）にも残差バイアスが残り、 真値 2,645,809 が 95% CI の外に出ることがある。 これは「多重代入は万能ではない」ことを示すが、 補完の不確実性を区間幅に誠実に反映するという MI の本質的価値は変わらない。 実務では補助変数の追加や感度分析でこの残差バイアスに対処する。

🐍 欠損パターンの可視化（missingno）

🎯 このコードでやること: 欠損は「どこに、 どんなパターンで」入っているかを目視で診断する。 missingno を使って matrix / heatmap / dendrogram を描画し、 MAR/MNAR の手掛かりを得る。

📥 入力: SSDSE-B-2026 に MAR 欠損を注入したデータ。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

import pandas as pd
import numpy as np
import missingno as msno
import matplotlib.pyplot as plt

# 英字の項目コード（A1101 など）を使うので skiprows=[1] で読む
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].dropna()
rng = np.random.default_rng(0)
data = df[['A1101', 'A1301', 'A1303', 'B4101', 'A4101']].copy()
data.loc[rng.random(len(data)) < (0.2 + 0.6 * (data['A1301'] > data['A1301'].median())),
         'A1101'] = np.nan
data.loc[rng.random(len(data)) < 0.3, 'A4101'] = np.nan

# 各列の欠損率
print(data.isna().mean().round(3))
# 欠損相関（同時に欠ける傾向）
print(data.isna().corr().round(2))

# matrix プロット
msno.matrix(data, figsize=(8, 4))
plt.tight_layout()
plt.savefig('missing_matrix.png', dpi=100)

📤 実行例:

💬 結果の読み方: 「総人口(A1101)」が約 43% 欠損、 「出生数(A4101)」が約 28% 欠損。 両者の欠損相関は 0.05 ≈ 0、 つまり「同時に欠ける」傾向はない（独立した欠損パターン）。 もし欠損相関が高ければ、 同一の原因（例: 同じ調査票の同じ問）で欠損している可能性を疑う必要がある。 matrix プロットはこれを 視覚的に 一望できるので、 欠損診断の第一歩として強く推奨される。

数式を言葉で読み解く: 欠損相関と Little の MCAR 検定

Little (1988) の MCAR 検定は、 観測サブセットごとの平均ベクトル $\bar y^{(g)}$ が全体の平均ベクトル $\bar y$ と等しいかを Wald 型統計量で評価する:

$d^2$ が大きく p<0.05 なら「MCAR を棄却」= 少なくとも MAR か MNAR。 ただし MCAR を 棄却できなかった としても MCAR が成立する証拠にはならない（仮説検定の非対称性）。 欠損相関の可視化と組み合わせて総合判断するのが実務的。

論文での記述テンプレ（コピペ可）

このテンプレを基に、 自分のデータの数値に置き換えるだけで論文の Methods 節になる。 「補助変数として何を入れたか」を明記することが、 査読者から最も問われるポイント。

⚠️ よくある落とし穴

⚠️ よくある落とし穴（5+）

⚠️ 欠損処理の 8 つの落とし穴

特に 単一代入と多重代入の違い を理解していない研究者は多い。 「補完したから OK」と片付ける前に、 「補完の不確実性」を信頼区間に反映できているか必ず確認すること。

🌐 関連手法・派生

• Listwise Deletion (CCA)：MCAR でのみ不偏
• Pairwise Deletion：利用可能ペアのみ
• 単一補完（mean/median/mode）：分散縮小
• 回帰補完：他変数から予測
• 確率的回帰補完：予測+ノイズ
• KNN 補完：類似サンプルから借用
• MICE：m 個を Rubin's rules で統合
• EM：欠損データの尤度最大化
• Heckman：MNAR への 2 段階推定
• Pattern-Mixture：欠損パターン別モデル
• FIML：SEM での自然な扱い
• missForest：非線形補完
• GAIN, MIWAE：深層生成モデル

🌐 欠損処理ライブラリ早見表

最初は SimpleImputer でベースラインを出し、 そこから IterativeImputer、 MICEData へとステップアップしていくのが現実的。

🔗 関連用語

🔗 関連用語（前提・並列・発展、 各 4+ 件）

前提として知っておきたい用語（4+）

並列・同カテゴリの用語（4+）

発展・応用先（4+）

📚 関連グループ教材

• 📚 データエンジニアリング
• 📚 データ前処理
• 📚 クレンジング処理
• 📚 データ品質

📜 歴史・背景

Donald B. Rubin (1976) "Inference and Missing Data" で体系化。 それ以前は「削除か平均補完」が標準。 Rubin は MCAR/MAR/MNAR の 3 分類と ignorable 概念を導入。 1987 年の著書で多重代入法が確立。

2000 年代に Schafer・van Buuren が MICE を実装し、 R の mice パッケージで普及。 scikit-learn の IterativeImputer もこの系譜。 GAIN, MIWAE, MissForest など深層補完も活発。

主要文献：Rubin (1976) Biometrika; Little & Rubin (2019) 第 3 版; van Buuren (2018) Flexible Imputation; Heckman (1979) Econometrica; Schafer & Graham (2002) Psychological Methods。

🌍 ケーススタディ

SSDSE-B-2026 では一部指標（小規模県の専門医療、 離島部の通信指標）が欠損。 単純削除では 8 県が消え n=39 に。 MAR と見て多重代入すれば 47 県全部を活用可能。

MIMIC-III 集中治療データベース例：患者症状が重いほど検査値未記録 → 典型的 MNAR。 単純補完では推定が楽観的すぎる。 indicator 変数の追加や重症度モデルとの同時推定が必要。

論文記述例（コピペ用）

❓ FAQ

🎓 深掘り — 高度な論点

サンプル選択との関係

「欠損」の一般化が サンプル選択バイアス。 「観測される事例自体が母集団から偏る」現象。

教育的演習

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

import numpy as np
# 学習用に分布パラメータで生成（本番では公的データを使用）
N = 200
np.random.seed(0)
income = np.random.lognormal(6.0, 0.7, N)
age    = np.random.normal(45, 12, N)

# (1) MCAR
y1 = income.copy(); y1[np.random.rand(N) < 0.3] = np.nan
# (2) MAR
prob = 1 / (1 + np.exp((age - 50)/5))
y2 = income.copy(); y2[np.random.rand(N) < prob] = np.nan
# (3) MNAR
prob3 = (income - income.min()) / (income.max() - income.min())
y3 = income.copy(); y3[np.random.rand(N) < prob3] = np.nan

print('真の平均:', income.mean())
print('MCAR後 :', np.nanmean(y1), '(不偏)')
print('MAR後  :', np.nanmean(y2), '(やや偏る)')
print('MNAR後 :', np.nanmean(y3), '(過小評価)')

補完前後で回帰係数を比較

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

from sklearn.linear_model import LinearRegression
from sklearn.impute import SimpleImputer

X = np.column_stack([age, np.ones(N)])
y_full = income; y_obs = y2  # MAR シナリオ
mask = ~np.isnan(y_obs)

print('真の係数 :', LinearRegression().fit(X, y_full).coef_)
print('Listwise:', LinearRegression().fit(X[mask], y_obs[mask]).coef_)
y_mean = SimpleImputer(strategy='mean').fit_transform(y_obs.reshape(-1,1)).ravel()
print('平均補完 :', LinearRegression().fit(X, y_mean).coef_)
# 観察: Listwise は MAR でもバイアスあり / 平均補完は係数を 0 寄りに

理論的注意：observed-data likelihood

完全データの尤度は $p(Y, M\mid\theta,\psi) = p(Y\mid\theta) p(M\mid Y,\psi)$。 観測データの尤度は $\int p(Y\mid\theta) p(M\mid Y,\psi) dY_{mis}$。 MAR では $p(M\mid Y,\psi) = p(M\mid Y_{obs},\psi)$ なので $\psi$ が $\theta$ と分離可能 → 欠損モデルを無視できる。 MNAR ではこれが成立せず、 同時推定が必要。

主要 R パッケージ早見表

✅ 実務チェックリスト

• □ 「欠損メカニズム」が指す概念を 1 文で説明できるか
• □ 適用条件（前提・制約）を 3 つ以上挙げられるか
• □ 代表的な落とし穴を 5 つ以上避けられるか
• □ Python での実装イディオムを暗記しているか
• □ 関連手法との使い分け基準を持っているか
• □ レポートでの記述例を 1 段落書けるか
• □ 公的データ（SSDSE / e-Stat）で動作確認したか
• □ 同カテゴリの他用語との違いを説明できるか

🔄 標準ワークフロー — 「欠損メカニズム」を使う実務 5 ステップ

推奨ディレクトリ構成

project/
├── data/
│   ├── raw/        # 取得した生データ（git にコミットして再現性確保）
│   └── processed/  # 前処理後の中間ファイル
├── notebooks/      # 探索的分析の Jupyter
├── src/            # 本番コード（モジュール化）
│   ├── preprocess.py
│   ├── analyze.py  # ← 欠損メカニズム の実装
│   └── visualize.py
├── results/        # 出力（CSV, PNG, モデルファイル）
├── docs/           # レポート、 図表説明
├── requirements.txt
└── README.md       # 「目的・データ・手順・結論」を 1 枚で

⚖️ 同種ツール・手法の比較

「欠損メカニズム」と並んで現場で使われる代替・関連手段を整理。 「どれが正解」ではなく「状況に応じて使い分ける」のが実務。

⚡ パフォーマンス・スケーラビリティの目安

高速化テクニック 5 つ

• ベクトル化：for ループを numpy / pandas 関数に置き換え（10〜100 倍速）
• 適切な dtype：int64 → int32、 object → category で半減〜数倍速
• 必要列のみ読み込み：pd.read_csv(..., usecols=[...]) でメモリ削減
• チャンク処理：chunksize=10000 で逐次処理
• 並列化：joblib, dask, multiprocessing で CPU を使い切る

🧪 品質保証 — 再現性・正確性チェック

テストすべき項目

最小限のテストコード例

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

import pytest
import pandas as pd
from src.analyze import missing_mechanism_function  # 仮の関数名

def test_basic_output():
    df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', header=1)
    result = missing_mechanism_function(df)
    assert result is not None
    assert len(result) > 0

def test_handles_missing():
    df = pd.DataFrame({'x': [1, 2, None, 4]})
    result = missing_mechanism_function(df)
    assert result.isnull().sum() == 0  # 仕様: NaN を除いて返す

def test_reproducible():
    df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', header=1)
    r1 = missing_mechanism_function(df)
    r2 = missing_mechanism_function(df)
    pd.testing.assert_frame_equal(r1, r2)

📚 関連論文・教科書（学習ロードマップ）

公的データソース（欠損メカニズム の練習用）

• e-Stat（政府統計の総合窓口） — 各省庁の統計を一元アクセス。 SSDSE-B-2026 の出典でもある
• RESAS（地域経済分析システム） — 観光・人口・産業の地域別データ
• OECD.Stat — 国際比較データ。 GDP・教育・健康など
• World Bank Open Data — 開発指標・経済データ
• e-Gov 法令検索 — 統計法、 個人情報保護法など、 「欠損メカニズム」を扱う際の法的根拠
• 政府 CIO 「データカタログサイト」 — オープンデータの目録

📝 レポート・論文での記述テンプレート

Methods 節（コピペ可）

Results 節（コピペ可）

Discussion 節（コピペ可）

📚 欠損処理を深く学ぶための文献

• Rubin, D.B. (1976) Inference and missing data. Biometrika 63(3): 581-592. — MCAR/MAR/MNAR の原典
• Little, R.J.A. & Rubin, D.B. (2019) Statistical Analysis with Missing Data, 3rd ed. — 教科書の決定版
• van Buuren, S. (2018) Flexible Imputation of Missing Data, 2nd ed. — MICE の作者による実践書（無料 PDF あり）
• Honaker, J., King, G. (2010) What to do about missing values in time-series cross-section data. Amelia の理論的背景
• Sterne, J.A.C. et al. (2009) Multiple imputation for missing data in epidemiological and clinical research. BMJ 338: b2393.

数式を言葉で読み解く: Rubin の規則の心臓部

多重代入で得た $m$ 個の推定量 $\hat Q^{(k)}$ と その分散 $U^{(k)}$ について

$\bar U$ は「もし正解の補完値が分かったときの分散」、 $B$ は「補完値が一意でないことから来る追加の分散」。 両者を合算した $T$ がプールされた分散。 $1/m$ の補正項は「$m$ 有限による補完平均のサンプリング誤差」を反映する。 これが Rubin が 1987 年に確立した 正しい多重代入の分散合算。

🗺 欠損処理ワークフロー（決定木）

💡 大原則: 「欠損は 情報 である」。 欠損のパターン自体が研究対象の挙動を語ることがある（例: 「収入が高い人ほど答えない」という事実が、 すでに収入分布の歪みを示唆している）。 単に補完して進むのではなく、 まず欠損の意味を考えること。

🧪 理解度チェック — 自分の理解を 7 問で診断する

以下の 7 問に答えてから先へ進む。 解答は各設問の 解説 を畳んで確認すること。 「3 メカニズムを言葉だけで説明できるが、 補完手法の選択に詰まる」という人を主な対象としている。

💡 採点ガイド: 5 問以上正答できれば「実務適用フェーズ」に進める。 3-4 問なら本ページの 🌍 ケーススタディ + 📐 数式定義に戻る。 2 問以下なら 🎨 直感セクションから読み直すこと。

🖼 図で掴む — 欠損メカニズムの分布的サイン

3 つのメカニズム（MCAR / MAR / MNAR）は、 単に「ランダムか否か」では区別しきれない。 観測値と欠損値の 分布のずれ を視覚的に確認することで、 仮定の妥当性が見える。 ここでは SSDSE-B-2026 を素材にした典型図を 3 枚並べ、 「何を見れば次のステップが決まるか」を示す。

📌 視覚的診断のチェックリスト: (1) 欠損あり/なし群のヒストグラム重ね合わせ、 (2) 散布図上での欠損行マーカー、 (3) missingno の matrix / heatmap、 (4) ロジスティック回帰で「欠損 vs 観測変数」を回帰し有意性を見る。 4 つすべて見て初めて MAR の妥当性を主張できる。

🔗 隣接手法への橋渡し

欠損機構 (MCAR / MAR / MNAR) は単独で扱う概念ではなく、 上流の欠損パターン可視化 (missingno)、 並列の補完手法選択 (mean / 多重代入 / モデルベース)、 下流の感度分析と組み合わせて初めて妥当な分析になる。 機構を誤判定すると補完がバイアスを生む。

• 上流 (前段処理): 欠損パターンの可視化 — missingno マトリクスで構造を確認
• 並列 (代替/比較): 検出力 / バイアス分析 — 欠損が結果に与える歪みの評価
• 下流 (後段活用): 多重代入 / 重み付き推定 — メカニズムに応じた補完手法選択

SSDSE-B-2026 にあえて欠損を入れて検証する場合、 例えば「高所得県ほど一部列が欠損」(MNAR) を仮定して、 上流で欠損パターン行列を可視化、 中段で MICE 多重代入、 下流で補完ありと完全ケース分析の結果を比較、 という三段構成で頑健性を確認する。

🌳 手法選択フロー

欠損機構の判定と補完方針は、 (1) 観測変数間に欠損依存があるか (Little MCAR 検定)、 (2) 欠損が結果変数と関係するか、 (3) 補完前後で結論が変わるか、 で決まる。 MCAR なら listwise OK、 MAR なら MICE、 MNAR は感度分析必須。

1. 欠損は他変数と無関係か? はい → MCAR、 削除可能
2. 観測変数で説明できるか? はい → MAR、 MICE で多重代入
3. 欠損自体が値に依存するか? はい → MNAR、 感度分析・選択モデルで対応

実務では「完全ケース / 平均補完 / MICE / Pattern-mixture」の 4 法を必ず並べ、 結論が変わらないことを確認するのが標準で、 結論が手法依存なら欠損機構自体を再検討する。

🎮 触って理解する — 欠損メカニズムで平均はどう歪むか

下は 架空データ（実在の統計ではありません）です。 200 人分の x=年齢 と y=所得型スコア（年齢とともに増える相関あり）を生成しています。 欠損メカニズムを 3 つのボタンで切り替え、 欠損率スライダーを動かすと、 観測できる点（濃い点）と欠損した点（グレーの中空◯）、 そして 完全平均・観測平均・回帰代入後平均 の 3 本の水平線がリアルタイムで動きます。

ねらい：MCAR では観測平均が完全平均とほぼ一致（不偏）、 MAR / MNAR では観測平均が下にズレる（バイアス）ことを体感してください。 グラフ上を左右にドラッグしても欠損率を変えられます（タッチ対応）。

なぜ偏るのか（このメカニズムの説明）

補完法の簡易比較 — リストワイズ削除 / 平均代入 / 回帰代入

同じ欠損データに 3 つの対処法を当てた結果です（数値は上のスライダー設定でリアルタイム更新）。 平均代入 は平均は観測平均のまま動かず、 さらに ばらつき（標準偏差）を過小評価 する点に注目。

真値（完全平均）＝ ／ 完全データの標準偏差＝ 。 リストワイズ削除は観測平均と同じ値になります（欠損行を捨てるだけのため）。

💡 ここから持ち帰る 3 点

• MCAR は無料の昼食：観測平均がそのまま不偏。 削除も平均代入も平均は歪めない（ただし効率は落ちる）。
• MAR は x を使えば救える：欠損が観測変数 x に依存するだけなら、 回帰代入（x で y を予測） がバイアスをほぼ消せる。 平均代入は救えない。
• MNAR は観測データだけでは救えない：欠損が y 自身に依存するため、 x を使う回帰代入でも残差バイアスが残る。 感度分析・選択モデルなど外部仮定が要る。