想定読者：データ分析・機械学習を始めたばかり／dropna()や平均補完でとりあえず進めたことがある人
この記事のゴール：
- 欠損（missing） と MCAR/MAR/MNAR の意味を「用語から」理解する
- ダミーデータで “同じ欠損率でも、欠損の仕組み次第で結論が逆転する” を体験する
- 何から疑い、どう対処するかの 実務チェックリスト を持ち帰る
TL;DR（結論だけ先に）
- 欠損には「偶然抜けた」だけでなく、何らかの理由で偏って抜けるケースがあります。
- 同じ欠損率（例：35%欠損）でも、欠損が
- MCAR（完全にランダム）なら、単純に欠損行を捨てても大きくはズレにくい
- MAR（観測できる変数に依存）なら、単純比較は逆転しうるが、適切に条件（共変量）を入れると回復しうる
- MNAR（欠損している値そのものに依存）だと、条件を入れても結論が逆転したままになりうる
- この記事の再現コード（seed固定）では、真の効果が +0.4（有効） なのに：
- MAR：単純比較は -0.62（有害に見える）／回帰で補正すると +0.44（回復）
- MNAR：単純比較は -0.79（有害）／回帰で補正しても -0.24（まだ有害）
1. まず「欠損」と「欠損メカニズム」を用語から
欠損（Missing data）
データの一部が NaN などで抜けている状態です。
- センサー不調
- アンケート未回答
- 医療ならフォローアップ脱落（dropout）
- ログ欠落、通信途絶 など
欠損メカニズム（Missingness mechanism）
欠損が「どういう仕組みで起きているか」を分類したのが MCAR / MAR / MNAR です。
欠損しているかどうかを表す変数を R（観測できたら1、欠損なら0）とすると、ざっくり次のイメージです。
MCAR（Missing Completely At Random：完全にランダム）
欠損が、データの値と無関係に起きる（たまたま抜けた）
例：ランダムにログが落ちた、ランダムに測定が失敗した
MAR（Missing At Random：観測できるものには依存）
欠損が、観測できる変数（年齢、重症度、群など）には依存するが、
欠損している値そのものには直接依存しない（※条件づけると）
例：重症度が高いほど通院が途切れやすい（重症度は観測できる）
MNAR（Missing Not At Random：欠損している値そのものに依存）
欠損が、欠損している値（本当の値）に依存する
→ 欠損しているからこそ、その値が見えない（厄介）
例：アウトカムが悪い（/良い）人ほどフォローアップに来ない
→ 「来なかった人のアウトカム」は見えないので、データだけでは判定しにくい
重要：
実務では MCARかMARかMNARかを“データだけで証明する”のは原理的に難しい です。
だからこそ「欠損理由の記録」「ドメイン知識」「感度分析」が大事になります。
2. 最小実験：同じ欠損率でも結論が逆転するのを再現する
この記事の設定
治療の効果（t=1が効くかどうか）を推定するトイモデルで再現します。
-
t：治療（0=なし、1=あり）※ここではランダムに割り付け -
x：観測できる共変量（例：重症度） -
y_true：本当のアウトカム（欠損前） -
y：観測されたアウトカム（欠損後はNaN）
真の世界はこう作ります：
-
真の治療効果：
+0.4（本当は「効く」） -
xもアウトカムに影響する（現実っぽく）
そして 欠損率はどれも約35% に揃えたうえで、
- MCAR：ランダム欠損
- MAR：観測できる
xとtに依存して欠損 - MNAR：観測できない
y_trueに依存して欠損（=現実では見えない）
を作り、推定結果がどう変わるか見ます。
3. Google Colab 実行コード（コピペでOK）
3.1 import & 補助関数
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
def sigmoid(z):
return 1 / (1 + np.exp(-z))
3.2 ベースデータ（欠損なし）を作る
def make_base_data(n=5000, seed=0, beta_t=0.4, beta_x=1.0, noise=1.0):
"""
t: treatment (0/1) ... randomized
x: observed covariate
y_true: true outcome (before missingness)
"""
rng = np.random.default_rng(seed)
x = rng.normal(0, 1, n)
t = rng.integers(0, 2, n) # randomized treatment
y = beta_t * t + beta_x * x + rng.normal(0, noise, n)
return pd.DataFrame({"x": x, "t": t, "y_true": y})
base = make_base_data()
display(base.head())
print("True effect (beta_t) = 0.4")
print("n =", len(base))
print("t=1 ratio =", base["t"].mean())
3.3 欠損を作る（MCAR / MAR / MNAR）
ポイント：欠損率が揃うように、ロジスティック関数の切片を 二分探索で調整しています。
（「仕組みだけ変えて欠損率は同じ」にしたいからです）
def apply_missing_MCAR(df, miss_rate=0.35, seed=2):
rng = np.random.default_rng(seed)
m = rng.random(len(df)) < miss_rate
out = df.copy()
out["y"] = out["y_true"].copy()
out.loc[m, "y"] = np.nan
out["missing"] = m.astype(int)
out["p_miss"] = miss_rate
return out
def apply_missing_MAR(df, miss_rate=0.35, strength=3.0, seed=2):
"""
MAR: missing depends on observed x and t (NOT on y_true directly).
Here we make it extreme on purpose:
- when t=1, larger x -> more likely missing
- when t=0, smaller x -> more likely missing
"""
rng = np.random.default_rng(seed)
x = df["x"].values
t = df["t"].values
z0 = strength * (2*t - 1) * x # depends on observed variables only
# tune intercept b so average missing rate is ~ miss_rate
lo, hi = -10, 10
for _ in range(60):
mid = (lo + hi) / 2
if sigmoid(mid + z0).mean() > miss_rate:
hi = mid
else:
lo = mid
b = (lo + hi) / 2
p = sigmoid(b + z0)
m = rng.random(len(df)) < p
out = df.copy()
out["y"] = out["y_true"].copy()
out.loc[m, "y"] = np.nan
out["missing"] = m.astype(int)
out["p_miss"] = p
return out
def apply_missing_MNAR(df, miss_rate=0.35, strength=1.2, seed=2):
"""
MNAR: missing depends on y_true itself (unobserved when missing).
Again, we exaggerate on purpose:
- when t=1, larger y_true -> more likely missing
- when t=0, smaller y_true -> more likely missing
"""
rng = np.random.default_rng(seed)
y = df["y_true"].values
t = df["t"].values
z0 = strength * (2*t - 1) * y # depends on y_true => MNAR
lo, hi = -10, 10
for _ in range(60):
mid = (lo + hi) / 2
if sigmoid(mid + z0).mean() > miss_rate:
hi = mid
else:
lo = mid
b = (lo + hi) / 2
p = sigmoid(b + z0)
m = rng.random(len(df)) < p
out = df.copy()
out["y"] = out["y_true"].copy()
out.loc[m, "y"] = np.nan
out["missing"] = m.astype(int)
out["p_miss"] = p
return out
df_mcar = apply_missing_MCAR(base)
df_mar = apply_missing_MAR(base)
df_mnar = apply_missing_MNAR(base)
for name, d in [("MCAR", df_mcar), ("MAR", df_mar), ("MNAR", df_mnar)]:
print(name, "missing rate =", d["y"].isna().mean())
3.4 「結論」を2通りで出す（単純比較 vs 回帰で補正）
ここでは初心者がやりがちな2つを比較します：
-
単純比較：観測できた
yだけでmean(t=1)-mean(t=0) -
回帰で補正：
y ~ t + xの回帰で、tの係数を見る（xで条件づけ）
def estimate_effects(df, label):
if label == "Full":
obs = df.copy()
obs["y"] = obs["y_true"]
miss_rate = 0.0
else:
obs = df.dropna(subset=["y"]).copy()
miss_rate = df["y"].isna().mean()
# 1) naive difference in means
naive = obs.loc[obs.t==1, "y"].mean() - obs.loc[obs.t==0, "y"].mean()
# 2) regression adjustment
X = obs[["t", "x"]].values
y = obs["y"].values
lr = LinearRegression().fit(X, y)
beta_t = lr.coef_[0]
beta_x = lr.coef_[1]
return {
"dataset": label,
"n_obs": len(obs),
"miss_rate": miss_rate,
"naive_diff": naive,
"reg_beta_t": beta_t,
"reg_beta_x": beta_x,
}
results = []
results.append(estimate_effects(base, "Full"))
results.append(estimate_effects(df_mcar, "MCAR"))
results.append(estimate_effects(df_mar, "MAR"))
results.append(estimate_effects(df_mnar, "MNAR"))
res_df = pd.DataFrame(results)
display(res_df)
3.5 欠損の仕組み＆推定結果の比較
(A) MCAR：xに依存しない（だいたいフラット）
def plot_missing_vs_x_quantile(df, title, filename, q=10):
tmp = df.copy()
tmp["x_bin"] = pd.qcut(tmp["x"], q=q, duplicates="drop")
g = tmp.groupby(["t", "x_bin"], observed=True)["missing"].mean().reset_index()
g["x_mid"] = g["x_bin"].apply(lambda b: (b.left + b.right)/2).astype(float)
plt.figure(figsize=(6,4))
for tval in [0, 1]:
sub = g[g["t"]==tval].sort_values("x_mid")
plt.plot(sub["x_mid"], sub["missing"], marker="o", label=f"t={tval}")
plt.ylim(0, 1)
plt.xlabel("x (quantile bin midpoint)")
plt.ylabel("Missing rate of y")
plt.title(title)
plt.legend()
plt.tight_layout()
plt.savefig(filename, dpi=200)
plt.show()
plot_missing_vs_x_quantile(df_mcar, "MCAR: missing does NOT depend on x", "mcar_missing_vs_x.png")
print("saved: mcar_missing_vs_x.png")
(B) MAR：観測できるxに依存（t=0とt=1で逆方向に偏らせる）
plot_missing_vs_x_quantile(df_mar, "MAR: missing depends on observed x (and t)", "mar_missing_vs_x.png")
print("saved: mar_missing_vs_x.png")
(C) MNAR：本当はy_trueに依存（※現実では“欠損したy”は見えないのがミソ）
def plot_missing_vs_ytrue_quantile(df, title, filename, q=10):
tmp = df.copy()
tmp["y_bin"] = pd.qcut(tmp["y_true"], q=q, duplicates="drop")
g = tmp.groupby(["t", "y_bin"], observed=True)["missing"].mean().reset_index()
g["y_mid"] = g["y_bin"].apply(lambda b: (b.left + b.right)/2).astype(float)
plt.figure(figsize=(6,4))
for tval in [0, 1]:
sub = g[g["t"]==tval].sort_values("y_mid")
plt.plot(sub["y_mid"], sub["missing"], marker="o", label=f"t={tval}")
plt.ylim(0, 1)
plt.xlabel("y_true (quantile bin midpoint)")
plt.ylabel("Missing rate of y")
plt.title(title)
plt.legend()
plt.tight_layout()
plt.savefig(filename, dpi=200)
plt.show()
plot_missing_vs_ytrue_quantile(df_mnar, "MNAR: missing depends on unobserved y (here y_true)", "mnar_missing_vs_ytrue.png")
print("saved: mnar_missing_vs_ytrue.png")
(D) 推定結果の比較（ここが“結論逆転”の本丸）
beta_t_true = 0.4
labels = res_df["dataset"].tolist()
xpos = np.arange(len(labels))
width = 0.35
plt.figure(figsize=(8,4.5))
plt.bar(xpos - width/2, res_df["naive_diff"], width, label="Naive: mean(t=1)-mean(t=0)")
plt.bar(xpos + width/2, res_df["reg_beta_t"], width, label="Regression: coef of t (y ~ t + x)")
plt.axhline(beta_t_true, linestyle="--", linewidth=2, label=f"True effect = {beta_t_true}")
plt.xticks(xpos, labels)
plt.ylabel("Estimated treatment effect")
plt.title("Missing-data mechanism can flip the conclusion")
plt.legend()
plt.tight_layout()
plt.savefig("effect_comparison.png", dpi=200)
plt.show()
print("saved: effect_comparison.png")
3.6 （Colab用）図をまとめてダウンロード
from google.colab import files
for fn in ["mcar_missing_vs_x.png", "mar_missing_vs_x.png", "mnar_missing_vs_ytrue.png", "effect_comparison.png"]:
files.download(fn)
4. 結果の図
欠損の仕組み（MCAR / MAR / MNAR）
推定結果（結論が逆転する様子）
5. 結果の解釈
- Full（欠損なし）：治療効果は +0.4 付近（真の効果どおり）
- MCAR：欠損行を捨てても、単純比較も回帰もだいたい +0.4（ズレにくい）
-
MAR：
- 単純比較は 負（「治療は有害」に見える）←結論が逆転
- でも
xを入れて条件づける回帰だと 正（ほぼ回復）
-
MNAR：
- 単純比較は 負（有害に見える）
-
xを入れても回帰が 負のまま（結論が戻らない）
なぜMARで回復して、MNARで回復しないのか？
-
MAR：欠損は観測できる
xで説明できる
→y ~ t + xのように 欠損に関係する変数を入れて条件づけすれば、推定が回復しうる -
MNAR：欠損が
y自体に依存している
→ 欠損しているyは見えないので、データからは補正しきれない（仮定が必要）
6. チェックリスト（まず何をすべき？）
6.1 まず「欠損理由」をログ/設計で押さえる（最重要）
- 欠損は「通信・計測エラー」か「人が来ない（脱落）」か？
- 欠損は特定の条件（年齢、重症度、時間帯、装置、施設、群）で増えていないか？
- 欠損の定義は統一されているか？（空文字/0/NaN混在してないか）
6.2 MCARっぽいとき
- 欠損行削除（complete-case）でも大きくは歪みにくい
ただし サンプル数が減って不確実性が増える（検出力低下）
6.3 MARっぽいとき
- 欠損に関係する変数（今回なら
x）を モデルに含める - 予測なら：
SimpleImputer+MissingIndicator、あるいは（状況次第で）IterativeImputer/ 多重代入 - 効果推定なら：回帰で条件づけ、あるいは重み付け（IPW）や多重代入など
6.4 MNARが疑われるとき（要注意）
- 「MNARをデータだけで否定できない」 を前提にする
- 可能なら：
- 欠損理由を追加で取る（来院しない理由、測定できない理由）
- 感度分析（仮定を振って結果がどれだけ動くか）
- 欠損メカニズム自体をモデル化（選択モデル、パターン混合モデル等）
- あるいは「欠損そのものが情報（予兆）」として扱う（ビジネス/医療ではよくある）
7. まとめ
- 欠損は「NaNを埋める/捨てる」前に、欠損の仕組み（MCAR/MAR/MNAR）を疑う
- 同じ欠損率でも、欠損の仕組みが違うと 結論が逆転しうる
- MARなら観測できる変数で条件づけて回復できる場合がある
- MNARなら追加情報・仮定・感度分析なしに“正解”を主張しにくい