はじめに

こんにちは、事業会社で働いているデータサイエンティストです。

本記事では、ベイズ階層埋め込みの手法を提案して、アメリカの上院の投票データから、議員のイデオロギーベクトルの推定と可視化を行います。

結果としては、モデルが設定した事前分布を突破するほど、投票傾向が共和党に近い民主党議員を見つけました。

さて、誰でしょうか？ぜひ最後まで読んでください！

データ説明

今回の記事で利用するのは、第109回アメリカ合衆国上院（2005 ~ 2006）の投票結果データです。

早速中身を確認しましょう。

data("s109", package = "pscl")

投票結果の行列はs109$votesに格納されています。

> s109$votes[1:10, 1:10]
                 1-1 1-2 1-3 1-4 1-5 1-6 1-7 1-8 1-9 1-10
BUSH (R USA)       9   1   1   9   9   9   9   9   1    1
SESSIONS (R AL)    6   1   1   1   1   6   6   6   1    1
SHELBY (R AL)      9   1   1   1   1   6   6   6   1    1
MURKOWSKI (R AK)   9   1   1   9   1   6   6   6   1    1
STEVENS (R AK)     6   1   1   1   1   6   6   6   1    1
KYL (R AZ)         9   1   1   1   1   6   6   6   1    1
MCCAIN (R AZ)      9   1   1   1   1   6   6   6   1    1
PRYOR (D AR)       6   1   1   1   6   1   1   1   6    1
LINCOLN (D AR)     6   1   6   1   1   1   6   1   1    1
BOXER (D CA)       1   6   6   1   6   1   1   1   6    1

次に、このデータをデータベース形式に整形しましょう。

ちなみに投票結果データのコーディングルールは

> s109$codes
$yea
[1] 1 2 3

$nay
[1] 4 5 6

$notInLegis
[1] 0

$missing
[1] 7 8 9

です。

rc_df <- s109$votes |> 
  as.data.frame() |>
  tibble::rownames_to_column(var = "senator") |>
  tibble::tibble() |>
  tidyr::pivot_longer(!senator, names_to = "bill", values_to = "result_cd") |>
  dplyr::filter(result_cd %in% c(1,2,3,4,5,6)) |>
  dplyr::mutate(
    result = dplyr::case_when(
      result_cd %in% c(1,2,3) ~ 1,
      TRUE ~ 0
    )
  )

中身はこんな感じです：

> rc_df
# A tibble: 62,857 × 4
   senator      bill  result_cd result
   <chr>        <chr>     <dbl>  <dbl>
 1 BUSH (R USA) 1-2           1      1
 2 BUSH (R USA) 1-3           1      1
 3 BUSH (R USA) 1-9           1      1
 4 BUSH (R USA) 1-10          1      1
 5 BUSH (R USA) 1-19          6      0
 6 BUSH (R USA) 1-52          6      0
 7 BUSH (R USA) 1-83          6      0
 8 BUSH (R USA) 1-84          6      0
 9 BUSH (R USA) 1-87          1      1
10 BUSH (R USA) 1-104         1      1
# ℹ 62,847 more rows
# ℹ Use `print(n = ...)` to see more rows

次に、マスターテーブルを作成して投票結果テーブルにジョインします。

senator_master <- rc_df |>
  dplyr::select(senator) |>
  dplyr::distinct() |>
  dplyr::arrange(senator) |>
  dplyr::mutate(
    senator_id = dplyr::row_number()
  ) |>
  dplyr::mutate(
    party_id = dplyr::case_when(
      # 議員名カラムの命名規則は「苗字（所属政党 州名）」なので、「(R」が入っていてば共和党（Republican）になります。
      # Dは民主党（Democrat）です。
      stringr::str_detect(senator, "\\(R") == TRUE ~ 1,
      TRUE ~ 2
    )
  )

bill_master <- rc_df |>
  dplyr::select(bill) |>
  dplyr::distinct() |>
  dplyr::arrange(bill) |>
  dplyr::mutate(
    bill_id = dplyr::row_number()
  )

rc_df_with_id <- rc_df |>
  dplyr::left_join(senator_master, by = "senator") |>
  dplyr::left_join(bill_master, by = "bill")

これでモデルに入れるためのデータの準備ができましたので、続いてはモデルの説明に入りたいと思います。

モデル説明

では、まずベイズモデルの定式化からはじめます。

全ての法案billについて

$$
法案人気度_{bill} \sim Normal(0, 1)
$$

$$
法案ベクトル_{bill} \sim Normal(0, 1)
$$

全ての政党pについて

$$
政党イデオロギーベクトル_{p} \sim Normal(0, 1)
$$

全ての議員sについて

$$
sigma_{s} \sim Gamma(0.01, 0.01)
$$

$$
議員イデオロギーベクトル_{s} \sim Normal(政党イデオロギーベクトル_{政党_{s}}, sigma_{s})
$$

全ての投票結果(議員と法案の組み合わせ)iについて

$$
賛成フラグ \sim Bernoulli(logit(法案人気度_{法案_{i}} + 法案ベクトル_{法案_{i}}' * 議員イデオロギーベクトル_{議員_{i}}))
$$

言葉で複雑なところを説明しますと、まず議員が法案に賛成するかどうかは、法案の人気度（誰も反対しないような法案など）に法案ベクトルと議員イデオロギーベクトルの内積を足した値で決定されると仮定します。

法案の人気度が高いか法案ベクトルと議員イデオロギーベクトルの向きが近い場合に賛成票が投じられやすいです。

次に議員のイデオロギーベクトルは、その所属政党のイデオロギーベクトルを平均とする正規分布に従うと仮定されるので、その所属政党のイデオロギーベクトル付近に推移します。

これは議員のイデオロギーベクトルは、その所属政党のイデオロギーに近いのであろうというドメイン知識をベイズモデルに伝える方法です。

ただ、データから十分なエビデンスがあれば、議員のイデオロギーベクトルの分散パラメータの値を高く設定することで、その所属政党のイデオロギーを突破することも許容します。

このような柔軟性こそ、ベイズ機械学習がアカデミアとビジネスの世界で価値を発揮できる理由だと思います。

Stanのコードはこちらです：

senate_emb.stan

functions {
  real partial_sum_lpmf(
    array[] int result,
    int start, int end,
    
    array[] int senator,
    array[] int bill,
    
    array[] vector senator_latent,
    vector bill_popularity,
    array[] vector bill_latent
  ){
    vector[end - start + 1] lambda;
    int count = 1;
    for (i in start:end){
      lambda[count] = bill_popularity[bill[i]] + bill_latent[bill[i]] '* senator_latent[senator[i]];
      count += 1;
    }
    return (
      bernoulli_logit_lupmf(result | lambda)
    );
  }
}
data {
  int N;
  int party_type;
  int senator_type;
  int bill_type;
  int K;
  
  // 上院議員マスター
  array[senator_type] int party_id;
  
  //　投票データ
  array[N] int senator;
  array[N] int bill;
  array[N] int result;
  
  //　検証データ
  int val_N;
  array[val_N] int val_senator;
  array[val_N] int val_bill;
  array[val_N] int val_result;
}
parameters {
  vector[bill_type] bill_popularity;
  array[bill_type] vector[K] bill_latent;
  
  array[party_type] vector[K] party_latent;
  
  vector<lower=0>[senator_type] senator_sigma;
  array[senator_type] vector[K] senator_latent;
}
model {
  bill_popularity ~ std_normal();
  
  for (b in 1:bill_type){
    bill_latent[b] ~ std_normal();
  }
  
  for (p in 1:party_type){
    party_latent[p] ~ std_normal();
  }
  
  senator_sigma ~ gamma(0.01, 0.01);
  
  for (s in 1:senator_type){
    senator_latent[s] ~ normal(party_latent[party_id[s]], senator_sigma[s]);
  }
  
  int grainsize = 1;
  
  target += reduce_sum(
    partial_sum_lupmf, result, 
    
    grainsize,
    
    senator,
    bill,
    
    senator_latent,
    bill_popularity,
    bill_latent
  );
}
generated quantities {
  real F1_score;
  array[senator_type, senator_type] real cos_sim;
  {
    array[val_N] int predicted;
    int TP = 0;
    int FP = 0;
    int FN = 0;
    for (i in 1:val_N){
      predicted[i] = bernoulli_logit_rng(bill_popularity[val_bill[i]] + bill_latent[val_bill[i]] '* senator_latent[val_senator[i]]);
      if (val_result[i] == 1 && predicted[i] == 1){
        TP += 1;
      }
      else if (val_result[i] == 0 && predicted[i] == 1){
        FP += 1;
      }
      else if (val_result[i] == 1 && predicted[i] == 0){
        FN += 1;
      }
    }
    F1_score = (2 * TP) * 1.0/(2 * TP + FP + FN);
    
    for (i in 1:senator_type){
      for (j in 1:senator_type){
        cos_sim[i, j] = (senator_latent[i] '* senator_latent[j])/sqrt(sum(senator_latent[i].^2) * sum(senator_latent[j].^2));
      }
    }
  }
}

モデルファイルの中には議員のイデオロギーベクトルのコサイン類似度を計算するコードもありますが、これは今後のdocker関連の記事のための計算なので、本記事では利用しません。

モデル推定

モデルを推定する前に、学習データと検証データに分ける必要があります。

set.seed(12345)
train_id <- sample(seq_len(nrow(rc_df_with_id)), nrow(rc_df_with_id) * 0.8)

# データは全体的に1（賛成）が多いので、検証データでは1と0（反対）が同じ数になるように整理します。

rc_df_with_id_test_all <- rc_df_with_id[-train_id,]

rc_df_with_id_test_all_zero_id <- which(rc_df_with_id_test_all$result == 0)

# 反対の数（少ない方）と同じ数の1の観測値のidを抽出

rc_df_with_id_test_all_one_id <- sample(
  which(rc_df_with_id_test_all$result == 1),
  length(which(rc_df_with_id_test_all$result == 0))
)

rc_df_with_id_test_balanced <- rc_df_with_id_test_all[c(rc_df_with_id_test_all_zero_id, rc_df_with_id_test_all_one_id),]

data_list <- list(
  N = nrow(rc_df_with_id[train_id,]),
  party_type = 2,
  senator_type = nrow(senator_master),
  bill_type = nrow(bill_master),
  K = 10,
  
  party_id = senator_master$party_id,
  
  senator = rc_df_with_id$senator_id[train_id],
  bill =  rc_df_with_id$bill_id[train_id],
  result = rc_df_with_id$result[train_id],
  
  val_N = nrow(rc_df_with_id[-train_id,]),
  val_senator = rc_df_with_id$senator_id[-train_id],
  val_bill =  rc_df_with_id$bill_id[-train_id],
  val_result = rc_df_with_id$result[-train_id]
)

次に、Stanをコンパイルして、データをモデルに入れます：

m_senate_emb_init <- cmdstanr::cmdstan_model("estimation/senate_emb.stan",
                                             cpp_options = list(
                                               stan_threads = TRUE
                                             )
                                             )

並列処理と変分推論のおかげでモデルの推定はすぐ終わりました。

> m_senate_emb_estimate <- m_senate_emb_init$variational(
     seed = 12345,
     threads = 6,
     iter = 50000,
     data = data_list
 )
------------------------------------------------------------ 
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.005608 seconds 
1000 transitions using 10 leapfrog steps per transition would take 56.08 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       -17458.033             1.000            1.000 
   200       -15150.353             0.576            1.000 
   300       -14414.824             0.401            0.152 
   400       -14101.542             0.306            0.152 
   500       -13943.512             0.247            0.051 
   600       -13860.808             0.207            0.051 
   700       -13815.225             0.178            0.022 
   800       -13768.956             0.156            0.022 
   900       -13755.614             0.139            0.011 
  1000       -13714.420             0.125            0.011 
  1100       -13716.733             0.114            0.006   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  14.0 seconds.

最後に、結果ファイルを抽出します。

m_senate_emb_summary <- m_senate_emb_estimate$summary()

結果の可視化

精度

まず、F1スコアの事後分布を確認しましょう。

m_senate_emb_estimate$draws("F1_score") |> 
  tibble::as_tibble() |>
  dplyr::mutate(F1_score = as.numeric(F1_score)) |>
  ggplot2::ggplot() + 
  ggplot2::geom_density(ggplot2::aes(x = F1_score),
                        fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::ggtitle("F1 scoreの事後分布") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

0.9くらいでなかなかいい精度ですね。

イデオロギーベクトル

さて、一番面白いところに入ります。

まず、政党と議員のイデオロギーベクトルを抽出して、t-SNEで処理します：

party_latent <- m_senate_emb_summary |>
  dplyr::filter(stringr::str_detect(variable, "party_latent")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 10)

senator_latent <- m_senate_emb_summary |>
  dplyr::filter(stringr::str_detect(variable, "senator_latent")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 10)

tsne_result <- Rtsne::Rtsne(
  rbind(party_latent, senator_latent)
  )

最後に、描画用のラベルなどを付けて、ggplot2で可視化します：

g_tsne <- tsne_result$Y |>
  as.data.frame() |>
  dplyr::tibble() |>
  dplyr::bind_cols(
    name = c("R", "D", senator_master$senator),
    party = c("party", "party", c("R", "D")[senator_master$party_id]),
    label_size = c(rep(1.2, 2), rep(1, nrow(senator_master)))
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_text(ggplot2::aes(x = V1, y = V2, label = name, color = party, size = label_size)) + 
  ggplot2:: scale_color_manual(
    values = c(
      "D" = ggplot2::alpha("blue", 0.5),
      "R" = ggplot2::alpha("red", 0.5),
      "party" = "yellow"
    )
  ) + 
  ggplot2::guides(size = "none") + 
  ggplot2::xlab("dimension 1") + 
  ggplot2::ylab("dimension 2") + 
  ggplot2::ggtitle("Visualization of Estimated Ideology Vector")

g_tsne

この図では、青い文字は民主党の上院議員のイデオロギーベクトル、赤い文字は共和党の上院議員のイデオロギーベクトルをそれぞれ表しています。

まず、ベイズモデルの設定通り、民主党の上院議員のイデオロギーベクトルは民主党のイデオロギーベクトル（左上の黄色いD）付近で推移し、共和党の上院議員のイデオロギーベクトルも共和党のイデオロギーベクトル（右下の黄色いR）付近で推移してい、、、

うん？

！？

あれ？なんか右下に青い文字が、、、

そこは共和党で赤いはずなのに、、、

plotlyで少し拡大してみましょう。

plotly::ggplotly(g_tsne)

なるほど、この方のようですね：

英語のウィキペディアのページにはこんな記述があります。

Nelson was one of the most conservative Democrats in the Senate, frequently voting against his party.
（Nelson議員は上院の中で最も保守的な民主党員で、よくその所属政党に反する形で投票していた。）

なるほど、アメリカの国内政治にそこまで詳しくないですが、そもそもドメイン知識的に共和党にかなり近い民主党員のようです。

これをベイズモデルはしっかり検出して、事前分布の形で教えた「議員はその所属政党とイデオロギーが近いぞ！」という仮説を見事に「棄却」し、Nelson議員のイデオロギーベクトルを共和党付近で推定できました。

結論

このように、ドメイン知識を教えながら、エビデンスが十分だったらよしなにドメイン知識を「棄却」して、柔軟にデータから情報を学習するのが、ベイズモデルの強みです。僕はいつも勝手にこれをソフトなルールベースとビジネス側に説明しています。

皆さんもぜひアカデミアとビジネスの世界で、ベイズモデルを活用してください！