はじめに
こんにちは、事業会社で働いているデータサイエンティストです。
事業会社に限らず、どのような実証分析の分野でも、「欲しい特徴量が欠損している」という問題は避けて通れません。
政治学方法論の例でいえば、アンケートの特定の設問に回答してもらえない場合、その項目は欠損します。さらに言うと、民主主義体制では投票の自由が保障されているため、個人の投票結果データを誰も持っていませんし、そもそも持ってはいけません。
マクロ計量経済学の分野でも、政府機関が何らかの事情で一時的にデータ収集を停止したり、収集するデータの種類を変更したりすれば、同様に欠損が発生します。
ビジネスの現場でもこれは同じです。特に、人が入力するタイプのデータでは、KPIやインセンティブ計算と直接関係のないデータほど欠損しやすい傾向があります。
たとえば営業担当者が商談内容を入力し忘れても、受注さえ取れれば月次目標には計上され、評価上は問題になりません。しかしその結果、データサイエンティストや営業企画担当にとっては貴重な商談内容データが欠損してしまいます。
欠損データの扱い方は非常に奥が深いテーマですが、本記事では専門的な議論には立ち入りません。ここでは、あえて専門用語を避けながら、最もシンプルな状況に対応できるベイズモデルを紹介し、その性能を実際に確認してみます。
ただし注意点として、本記事で紹介するモデルは、いわゆる戦略的あるいは意図的な欠損には対応していません。たとえば、営業の商談データで「失注した場合、上司に怒られるのを避けるために商談内容を記録しない」といったケースでは、欠損の有無が目的変数（この場合は受注の成否）と相関してしまいます。このような状況では、本記事で紹介するモデルはうまく機能しません。
とはいえ、すべての欠損がこのような戦略的な理由で起きるわけではありません。単に「うっかり入力を忘れた」「忙しくて後回しにした」といったケースも多いはずです。社会人として、会社に指示されたことを忘れずにやるのは基本中の基本ですが、それでも人間なので、忘れるときはありますよね💦
そんなふうに、特に戦略的な理由もなくデータが欠損している場合には、ぜひ本記事で紹介するベイズモデルを試してみてください！
データシミュレーション
まずは、シミュレーション用のデータを10万件作成します。
ここでは、3つの異なるセグメント（潜在クラス）が存在する状況を想定しています。それぞれのセグメントには異なる効果（segment_effect）があり、これが目的変数 y に影響を与えます。また、セグメントではない説明変数として x を1つだけ用意し、誤差 e を加えています。
セグメントの一部は観測されず、segment_observed において 確率 0.9 で観測され、確率 0.1 で欠損（0として扱う） という設定にしています。このようにして、「一部のセグメントラベルが欠損している」という現実的な状況を人工的に再現しています。
具体的なコードはこちらです：
set.seed(12345)
my_df <- 100000 |>
rnorm() |>
tibble::tibble(
e = _
) |>
dplyr::mutate(
x = rnorm(dplyr::n()),
segment_true = sample(c(1,2,3), dplyr::n(), replace = TRUE),
segment_effect = dplyr::case_when(
segment_true == 1 ~ -4.66,
segment_true == 2 ~ -0.66,
TRUE ~ 5.33
),
y = 1.2 - 0.5 * x + segment_effect + e,
segment_observed = segment_true |>
purrr::map_int(
\(x){
flag <- rbinom(1, 1, 0.9)
if (flag == 1){
return(x)
} else {
return(0)
}
}
)
)
val_my_df <- my_df |>
dplyr::filter(segment_observed == 0)
このデータを使い、Stanで潜在セグメントを同時に推定するベイズモデルを適用し、欠損したセグメント情報をどこまで正確に復元できるかを検証します。
モデル構造
欠損するセグメントデータに対応するため、本記事ではセグメントの生成過程自体をモデルの一部として推定するという戦略を取ります。つまり、セグメントが観測されている場合も、観測されていない場合も、どちらも確率的な生成過程として一貫して扱い、その不確実性をベイズ的に取り込むことで、欠損した情報を自然に補完します。
具体的には、まずセグメントの生成過程を次のように定義します。
$$
\alpha \sim Gamma(1, 1)
$$
$$
\theta \sim Dirichlet(\alpha)
$$
$$
セグメント_{i} \sim Categorical(\theta)
$$
ここで、$\theta$ は各セグメント（潜在クラス）の生起確率を表し、$\alpha$ はそのディリクレ分布のハイパーパラメータです。この部分がどのセグメントがどれくらいの割合で存在するのかをモデル化しています。
次に、セグメントごとに異なる回帰係数（固定効果） $\beta_{セグメント}$ を持たせます。
$$
y_{i} \sim Normal(切片 + \beta_{x} x_{i} + \beta_{セグメント_{i}}, \sigma_{y})
$$
もしセグメントが観測されている場合（segment_observed != 0）、そのラベルをそのまま使います。 一方、欠損している場合（segment_observed == 0）は、その観測値のセグメントを一つの確率変数として扱います。
ただし、Stanをはじめとする多くの機械学習・統計モデリング用パッケージは、基本的に連続型パラメータの推定に最適化されており、離散型パラメータ（たとえばセグメントのラベル）の直接的な推定には対応していません。そのため、欠損したセグメントに対しては、すべてのセグメントの可能性を考慮し、それぞれの尤度を総当たりで評価したうえで周辺化を行う必要があります。
Stanでは、この周辺化を効率的に実装するために、以下のように log_sum_exp 関数を利用します。これにより、全てのセグメント候補の尤度を対数空間で安定的に合計し、数値的な安定性を保ちながら確率的な混合を表現できます。
Stanでの実装
Stanでの実装コードは以下のとおりです。各種係数には事前分布を設定していますが、本記事の主眼はモデル構造の説明にあるため、細部についてはコード内のコメントをご参照ください。
また、セグメントはカテゴリ変数であるため、そのままでは回帰モデルに直接投入できません。通常、このような場合は ダミー変数化を行い、どれか 1 つのセグメントを基準値として除外するか、あるいはカテゴリの固定効果（セグメント効果）の総和がゼロになるように制約を加えることで識別可能性を確保します。
本記事では、私が事業会社で日々数千以上のカテゴリデータと格闘してきた実務経験を踏まえ、後者の効果の総和をゼロにする方法を採用しています。ちなみに、先ほどのRでのシミュレーションコードを見返していただくと、セグメントごとの効果がほぼ足してゼロになるように設計している点にもお気づきいただけると思います。
最後の generated quantities ブロックでは、欠損しているセグメントを確率的に復元し、実際の正解データと比較することで F1 スコアを算出し、モデルの推定性能を評価できるようにしています。
data {
int seg_type;
int N;
vector[N] x;
array[N] int seg;
vector[N] y;
int val_N;
vector[val_N] val_x;
array[val_N] int val_seg;
vector[val_N] val_y;
}
parameters {
vector<lower=0>[seg_type] alpha;
simplex[seg_type] theta;
real<lower=0> sigma_beta_x;
vector<lower=0>[seg_type] sigma_beta_seg;
real intercept;
real beta_x;
vector[seg_type] beta_seg_unnormalized;
real<lower=0> sigma_y;
}
transformed parameters {
vector[seg_type] beta_seg;
beta_seg = beta_seg_unnormalized - mean(beta_seg_unnormalized);
}
model {
alpha ~ gamma(1, 1);
theta ~ dirichlet(alpha);
sigma_beta_x ~ inv_gamma(5, 5);
sigma_beta_seg ~ inv_gamma(5, 5);
intercept ~ normal(0, 10);
beta_x ~ normal(0, sigma_beta_x);
beta_seg ~ normal(0, sigma_beta_seg);
for (i in 1:N){
if (seg[i] != 0){
seg[i] ~ categorical(theta);
y[i] ~ normal(intercept + x[i] * beta_x + beta_seg[seg[i]], sigma_y);
} else {
vector[seg_type] case_when;
for (j in 1:seg_type){
case_when[j] = log(theta[j]) + normal_lpdf(y[i] | intercept + x[i] * beta_x + beta_seg[j], sigma_y);
}
target += log_sum_exp(case_when);
}
}
}
generated quantities {
array[val_N] int predicted_seg;
vector[seg_type] F1_score_by_class;
real F1_score_macro;
for (i in 1:val_N){
vector[seg_type] case_when;
for (j in 1:seg_type){
case_when[j] = log(theta[j]) + normal_lpdf(val_y[i] | intercept + val_x[i] * beta_x + beta_seg[j], sigma_y);
}
predicted_seg[i] = categorical_rng(softmax(case_when));
}
for (s in 1:seg_type){
int TP = 0;
int FP = 0;
int FN = 0;
for (i in 1:val_N){
if (val_seg[i] == s && predicted_seg[i] == s){
TP += 1;
}
else if (val_seg[i] != s && predicted_seg[i] == s){
FP += 1;
}
else if (val_seg[i] == s && predicted_seg[i] != s){
FN += 1;
}
}
F1_score_by_class[s] = (2 * TP) * 1.0/(2 * TP + FP + FN);
}
F1_score_macro = mean(F1_score_by_class);
}
モデル推定
まずはモデルをコンパイルして：
m_init <- cmdstanr::cmdstan_model("latent_segment.stan")
変分推論でパラメータを推定します：
> m_estimate <- m_init$variational(
seed = 123,
data = list(
seg_type = 3,
N = nrow(my_df),
x = my_df$x,
seg = my_df$segment_observed,
y = my_df$y,
val_N = nrow(val_my_df),
val_x = val_my_df$x,
val_seg = val_my_df$segment_true,
val_y = val_my_df$y
)
)
------------------------------------------------------------
EXPERIMENTAL ALGORITHM:
This procedure has not been thoroughly tested and may be unstable
or buggy. The interface is subject to change.
------------------------------------------------------------
Gradient evaluation took 0.018211 seconds
1000 transitions using 10 leapfrog steps per transition would take 182.11 seconds.
Adjust your expectations accordingly!
Begin eta adaptation.
Iteration: 1 / 250 [ 0%] (Adaptation)
Iteration: 50 / 250 [ 20%] (Adaptation)
Iteration: 100 / 250 [ 40%] (Adaptation)
Iteration: 150 / 250 [ 60%] (Adaptation)
Iteration: 200 / 250 [ 80%] (Adaptation)
Success! Found best value [eta = 1] earlier than expected.
Begin stochastic gradient ascent.
iter ELBO delta_ELBO_mean delta_ELBO_med notes
100 -262145.232 1.000 1.000
200 -251991.530 0.520 1.000
300 -251600.987 0.347 0.040
400 -251450.319 0.261 0.040
500 -251522.446 0.209 0.002 MEDIAN ELBO CONVERGED
Drawing a sample of size 1000 from the approximate posterior...
COMPLETED.
Finished in 20.5 seconds.
20秒で終わりました！
結果のサマリーを保存します：
m_summary <- m_estimate$summary()
結果の可視化
まずは F1 スコアを確認し、モデルが欠損しているセグメントデータをどの程度正確に復元できているかを見てみましょう：
m_estimate$draws("F1_score_by_class") |>
tibble::as_tibble() |>
`colnames<-`(as.character(c(1:3))) |>
dplyr::mutate(
rid = dplyr::row_number(),
dplyr::across(
dplyr::everything(),
~ as.numeric(.)
)
) |>
tidyr::pivot_longer(!rid, names_to = "type") |>
dplyr::bind_rows(
m_estimate$draws("F1_score_macro") |>
tibble::as_tibble() |>
`colnames<-`("macro") |>
dplyr::mutate(
rid = dplyr::row_number(),
dplyr::across(
dplyr::everything(),
~ as.numeric(.)
)
) |>
tidyr::pivot_longer(!rid, names_to = "type")
) |>
ggplot2::ggplot() +
ggplot2::geom_density(ggplot2::aes(x = value, fill = type), alpha = 0.5) +
ggplot2::labs(
title = "posterior distributions of F1 scores"
)
結果として、すべての F1 スコアが 0.95 を統計的に有意に上回る という非常に良好で完璧に近い性能が得られました！
最後に、セグメントの回帰係数の事後分布を確認します；
m_estimate$draws("beta_seg") |>
tibble::as_tibble() |>
`colnames<-`(as.character(c(1:3))) |>
dplyr::mutate(
rid = dplyr::row_number(),
dplyr::across(
dplyr::everything(),
~ as.numeric(.)
)
) |>
tidyr::pivot_longer(!rid, names_to = "type") |>
ggplot2::ggplot() +
ggplot2::geom_density(ggplot2::aes(x = value, fill = type), alpha = 0.5) +
ggplot2::geom_vline(xintercept = -4.66, linetype = "dashed", color = "red",linewidth = 0.5) +
ggplot2::geom_vline(xintercept = -0.66, linetype = "dashed", color = "red",linewidth = 0.5) +
ggplot2::geom_vline(xintercept = 5.33, linetype = "dashed", color = "red",linewidth = 0.5) +
ggplot2::labs(
title = "posterior distributions of coefficients"
)
破線が正解値を示しています。事後分布の分散はほとんどなく、正解値とほぼ完全に一致しているため、少し見づらいかもしれませんが、こちらも非常に高い精度を達成しています。
結論
本記事では、欠損したセグメント情報を含むデータに対して、潜在変数としてセグメントをモデル化するベイズアプローチを紹介しました。
シミュレーションの結果、モデルは欠損したセグメントラベルを非常に高い精度で復元できることが示され、F1スコアはすべて 0.95 を上回る 結果となりました。また、セグメントごとの回帰係数もほぼ正確に推定できており、事後分布の分散が非常に小さいことから、モデルの安定性と識別性能の高さも確認できました。
このアプローチのポイントは、欠損したラベルを確率的に扱い、その不確実性をベイズ的に取り込むことです。単純に欠損データを削除したり、平均値で置き換えたりするよりも、データの情報を最大限活用できます。
もちろん、本モデルは戦略的・意図的な欠損には対応していませんが、単純な入力忘れや計測の抜け漏れのような「偶発的な欠損」には非常に有効です。事業会社の実務データに限らず、政治学や経済学などさまざまな実証分析の現場で応用可能な手法です。
欠損データに悩まされている方は、ぜひ今回紹介したベイズモデルを試してみてください！
最後に、私たちと一緒に働きたい方はぜひ下記のリンクもご確認ください：