ディリクレ過程埋め込みモデルが勝手に議員の政党を発見してくれた件

# はじめに

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

https://ainow.ai/dx-magazine/?article=article-988

これは投票データで国会議員のイデオロギーを分析するシリーズの三本目の記事です。

二本目はこちらです：

https://qiita.com/Gotoubun_taiwan/items/c8198f81b94747df78dd

今までの記事で利用したモデルは、政党の数（共和党の民主党なので2）と議員の所属政党を明示的に知ることができたが、今回の記事では

- 政党の数
- 議員の所属政党

を投票データだけで推定するモデルを提案します。

本記事は、[今までの記事](https://qiita.com/Gotoubun_taiwan/items/c8198f81b94747df78dd)のコード・実行結果と解説を前提としています。

Dockerの記事はもう少し待ってください、、、

# モデル

## レジームの説明

言葉で説明しますと、今回の記事で提案するモデルと[一本目の記事](https://qiita.com/Gotoubun_taiwan/items/7d761159777bca174f74)のモデルの違いは、「所属政党」の数と議員の「所属政党」を自動で割り出せるところです。

ただし、「所属政党」という表現は少し比喩的で、正しくは議員のイデオロギーベクトルを教師なし的にクラスタリングして、議員がどのクラスに所属しているのかを推定するモデルです。これは政党というよりはイデオロギーのグループ・レジームとして解釈した方が良さそうなので、この記事ではレジームという表現を使います。

また、レジームが全体で何個あるかはモデルに明示的に教えておらず、理論上は無限にあるレジームからモデルが自動で利用するレジームをピックアップする感じです。これをディリクレ過程というノンパラメトリックなベイズモデリングの技法で実現します。

以上で説明したモデルの性質を踏まえて、このモデルを**ディリクレ過程レジーム埋め込みモデル**（**D**irichlet **P**rocess **R**egime **E**mbedding、**DPRE**）と命名したいと思います。

## 強みと政治学的な意義

同じ政党に所属しても、一枚岩とは限りません。これは基本的にはどの組織でもそうです。なので、分析者が議員の所属政党を明示的にモデルに渡すと、政党の中に存在しているミニ政党や、所属政党が違うものの投票行動で協力している超党派グループの割り出しなど、比較政治学的に意味のあるquantities of interestの抽出が難しくなるかもしれません。

もちろん、明示的に所属政党を渡しても、ベイズモデルがエビデンスを踏まえてよしなに仮説を棄却してくれるが、どうしても学習の柔軟性が下がります。

そこで、イデオロギーレジームの数も所属レジームも仮定しないディリクレ過程を利用すれば、政党という括りでの分析ではわからない傾向が明らかになることが期待されます。

最終的にイデオロギーレジームが政党と完全一致しても、これはこれでモデルが再発見できるほど、政党内の投票行動に一貫性があることを示すことで、政治学的な意義を発揮します。

## 定式化

概ね一本目の記事のままになりますが、今回の記事で提案するモデルの数理的定式化はこんな感じです：

全ての法案billについて

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

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

埋め込みレジーム全体について

$$
d\\_alpha \sim Gamma(1,1)
$$

$$
z \sim Stick\space Breaking(d\\_alpha)
$$

無限にある埋め込みレジームpについて

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

全ての議員sについて

$$
\eta_{s} \sim Multinomial(z)
$$

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

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

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

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

ディリクレ過程はstick breaking過程で構築します。

zは無限次元の単体（simplex）で、議員全体の中でレジームp所属の議員が占める割合が格納されています。$\eta$はそれぞれの議員の所属レジームを表すフラグです。

続いてはStanでの実装なんですが、数理モデルとの相違点をまとめると：

- 現代のコンピュータは無限次元のデータを受け付けません。なので、数理モデルで無限次元と表現されているパラメーターは、全て高次元のパラメーターで近似されます。
- Stanは$\eta$のような離散変数を受け付けません。なので、$\eta$は周辺除去されて、事後分布を計算する際にその期待値が再現されます。

## Stanでの実装

今回のモデルは離散変数を周辺化する処理がありますので、詳しく実装時の考え方に興味ある方はこちらの記事をご覧ください：

https://qiita.com/Gotoubun_taiwan/items/b40f3fdd25291b8aaf6f

コードはこんな感じです：


```cpp: senate_dper.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 P;
  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 {
  real<lower=0> d_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[P - 1] breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  vector[bill_type] bill_popularity;
  array[bill_type] vector[K] bill_latent;
  
  array[P] vector[K] P_latent;
  
  vector<lower=0>[senator_type] senator_sigma;
  array[senator_type] vector[K] senator_latent;
}
transformed parameters {
  simplex[P] z; // 
  // stick-breaking representationの変換開始
  // https://discourse.mc-stan.org/t/better-way-of-modelling-stick-breaking-process/2530/2 を参考
  z[1] = breaks[1];
  real sum = z[1];
  for (p in 2:(P - 1)) {
    z[p] = (1 - sum) * breaks[p];
    sum += z[p];
  }
  z[P] = 1 - sum;
}
model {
  d_alpha ~ gamma(1, 1);
  
  breaks ~ beta(1, d_alpha);
  
  for (p in 1:P){
    P_latent[p] ~ std_normal();
  }
  
  bill_popularity ~ std_normal();
  
  for (b in 1:bill_type){
    bill_latent[b] ~ std_normal();
  }
  
  senator_sigma ~ gamma(0.01, 0.01);
  
  for (s in 1:senator_type){
    vector[P] case_vector;
    for (p in 1:P){
      case_vector[p] = log(z[p]) + normal_lupdf(senator_latent[s] | P_latent[p], senator_sigma[s]);
    }
    target += log_sum_exp(case_vector);
  }
  
  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] vector[P] eta;
  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 (s in 1:senator_type){
      vector[P] this_eta;
      for (p in 1:P){
        this_eta[p] = z[p] * exp(normal_lpdf(senator_latent[s] | P_latent[p], senator_sigma[s]));
      }
      eta[s] = this_eta./sum(this_eta);
    }
    
    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));
      }
    }
  }
}

```

全体的には一本目の記事のモデルの特に尤度を定義する部分以外であまり変わらず、離散変数を周辺化する処理と周辺化されてしまった$\eta$の事後分布を出すところが増えただけです。

# モデル推定

## 前処理

念のため前処理のコードも全部含めたんですが、詳細は一本目の記事の説明を参照してください。

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

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
    )
  )

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(
      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") 

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

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)

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]
)

data_dper_list <- list(
  N = nrow(rc_df_with_id[train_id,]),
  P = 5,
  senator_type = nrow(senator_master),
  bill_type = nrow(bill_master),
  K = 10,
  
  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]
)

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



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

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

```

ここでは一本目の記事で紹介したモデルでも、本記事で紹介するモデルでも、埋め込み次元数Kは10に設定し、本記事で紹介するモデルにしかない最大レジーム数（実は無限大にしたいがコンピュータにそんなことができないのでやむをえず大きい値を設定するやつです）Pは5に設定しました。

## モデル推定

では、早速モデルの推定に入りましょう！

```r
> 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.005447 seconds 
1000 transitions using 10 leapfrog steps per transition would take 54.47 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.2 seconds. 
```

```r
> m_senate_dper_estimate <- m_senate_dper_init$variational(
     seed = 12345,
     threads = 6,
     iter = 50000,
     init = 20,
     data = data_dper_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.005675 seconds 
1000 transitions using 10 leapfrog steps per transition would take 56.75 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  -4607148015032416768.000             1.000            1.000 
   200  -47410052693459552.000            48.588           96.177 
   300  -587857039245164.875            58.942           79.649 
   400  -6239207484071.819            67.511           93.220 
   500  -36074836884.399            88.399           93.220 
   600   -323683824.166            92.075           96.177 
   700    -11239140.149            82.893           93.220 
   800      -145448.795            82.065           93.220 
   900       -23457.360            73.524           79.649 
  1000       -19661.522            66.191           79.649 
  1100       -18257.642            60.181           76.272   MAY BE DIVERGING... INSPECT ELBO 
  1200       -17498.721            55.169           76.272   MAY BE DIVERGING... INSPECT ELBO 
  1300       -16825.604            50.929           27.800   MAY BE DIVERGING... INSPECT ELBO 
  1400       -16321.371            47.293           27.800   MAY BE DIVERGING... INSPECT ELBO 
  1500       -15862.618            44.142            5.201   MAY BE DIVERGING... INSPECT ELBO 
  1600       -15400.431            41.385            5.201   MAY BE DIVERGING... INSPECT ELBO 
  1700       -15089.516            38.952            1.000   MAY BE DIVERGING... INSPECT ELBO 
  1800       -14836.410            36.789            1.000   MAY BE DIVERGING... INSPECT ELBO 
  1900       -14610.979            34.854            0.193   MAY BE DIVERGING... INSPECT ELBO 
  2000       -14420.575            33.111            0.193   MAY BE DIVERGING... INSPECT ELBO 
  2100       -14249.703            31.535            0.077   MAY BE DIVERGING... INSPECT ELBO 
  2200       -14099.469            30.102            0.077   MAY BE DIVERGING... INSPECT ELBO 
  2300       -13981.346            28.794            0.043   MAY BE DIVERGING... INSPECT ELBO 
  2400       -13902.582            27.594            0.043   MAY BE DIVERGING... INSPECT ELBO 
  2500       -13844.914            26.491            0.040   MAY BE DIVERGING... INSPECT ELBO 
  2600       -13797.300            25.472            0.040   MAY BE DIVERGING... INSPECT ELBO 
  2700       -13784.665            24.529            0.031   MAY BE DIVERGING... INSPECT ELBO 
  2800       -13767.241            23.653            0.031   MAY BE DIVERGING... INSPECT ELBO 
  2900       -13727.187            22.837            0.030   MAY BE DIVERGING... INSPECT ELBO 
  3000       -13737.419            22.076            0.030   MAY BE DIVERGING... INSPECT ELBO 
  3100       -13729.230            21.364            0.029   MAY BE DIVERGING... INSPECT ELBO 
  3200       -13715.785            20.696            0.029   MAY BE DIVERGING... INSPECT ELBO 
  3300       -13705.381            20.069            0.021   MAY BE DIVERGING... INSPECT ELBO 
  3400       -13714.721            19.479            0.021   MAY BE DIVERGING... INSPECT ELBO 
  3500       -13698.703            18.922            0.017   MAY BE DIVERGING... INSPECT ELBO 
  3600       -13692.291            18.397            0.017   MAY BE DIVERGING... INSPECT ELBO 
  3700       -13674.169            17.900            0.015   MAY BE DIVERGING... INSPECT ELBO 
  3800       -13691.249            17.429            0.015   MAY BE DIVERGING... INSPECT ELBO 
  3900       -13679.782            16.982            0.013   MAY BE DIVERGING... INSPECT ELBO 
  4000       -13668.253            16.557            0.013   MAY BE DIVERGING... INSPECT ELBO 
  4100       -13665.016            16.153            0.012   MAY BE DIVERGING... INSPECT ELBO 
  4200       -13659.933            15.769            0.012   MAY BE DIVERGING... INSPECT ELBO 
  4300       -13681.971            15.402            0.011   MAY BE DIVERGING... INSPECT ELBO 
  4400       -13642.981            15.052            0.011   MAY BE DIVERGING... INSPECT ELBO 
  4500       -13650.686            14.718            0.008   MEDIAN ELBO CONVERGED   MAY BE DIVERGING... INSPECT ELBO 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  37.8 seconds.
```

今回の記事で提案するDPREモデルの所要時間の方が長いですが、それでもかなり早いですね。まあ、複雑になったから仕方ないと思います。読者がより大規模世論調査やECサイトのトランザクションデータなどの複雑なデータにDPREを転用する際は、AWSかGCPでのGPU利用を是非検討してみてください。

最後に結果のサマリーを保存します。

```r
m_senate_emb_summary <- m_senate_emb_estimate$summary()

m_senate_dper_summary <- m_senate_dper_estimate$summary()
```

## F1スコアの事後分布の比較

さて、まずこの二つのモデルは本質的には分類器なので、F1スコアの事後分布で精度比較します。

```r
m_senate_emb_estimate$draws("F1_score") |> 
  tibble::as_tibble() |>
  dplyr::mutate(
    F1_score = as.numeric(F1_score),
    model = "party"
    ) |>
  dplyr::bind_rows(
    m_senate_dper_estimate$draws("F1_score") |> 
      tibble::as_tibble() |>
      dplyr::mutate(
        F1_score = as.numeric(F1_score),
        model = "DPRE"
      )
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_density(ggplot2::aes(x = F1_score, fill = model)) + 
  ggplot2::ggtitle("F1 scoreの事後分布") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3") + 
  ggplot2::scale_fill_manual(
    values = c("party" = ggplot2::alpha("blue", 0.3), 
               "DPRE" = ggplot2::alpha("red", 0.3))
  )
```

ただし、party（政党）モデルが一本目の記事のモデル、DPREが今回のモデルです。

![compare_f1_scores.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/30a90ab7-e0f8-f970-edff-df0d013b76c9.png)

F1スコアの事後分布を比較すると、DPREモデルの方が精度的に優れていますが、事後分布がかなり重なっており、両者の違いは統計学的に有意ではないといえます。よって、どれがより優れているかは統計学的に不明です。

ただ、イデオロギーのレジームを自動で発見したのに、分析者が明示的に所属政党を提供するモデルと変わらない精度を達成できるのはすごいですね。

## t-SNE

次に、t-SNEを可視化します：

```r
P_dper_latent <- m_senate_dper_summary |>
  dplyr::filter(stringr::str_detect(variable, "P_latent")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 10)

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

set.seed(12345)
tsne_dper_result <- Rtsne::Rtsne(
  rbind(P_dper_latent, senator_dper_latent)
)

g_dper_tsne <- tsne_dper_result$Y |>
  as.data.frame() |>
  dplyr::tibble() |>
  dplyr::bind_cols(
    name = c("1", "2", "3", "4", "5", senator_master$senator),
    party = c(rep("latent", 5), c("R", "D")[senator_master$party_id]),
    label_size = c(rep(1.2, 5), 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),
      "latent\ngroup" = ggplot2::alpha("black", 0.5)
    )
  ) + 
  ggplot2::guides(size = "none") + 
  ggplot2::xlab("dimension 1") + 
  ggplot2::ylab("dimension 2") + 
  ggplot2::ggtitle("Visualization of Estimated Ideology Vector")

g_dper_tsne
```

![dpre_tsne.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/11b62e3f-0c76-954a-eb34-24f431791727.png)

この図の見方を説明しますと、色がついている文字が議員のt-SNEの解析結果で、青い文字は民主党員で、赤い文字は共和党員です。そして図の中に入っている数字はモデルが推定したイデオロギーレジームの代表的なイデオロギーベクトルです。イデオロギーレジームは全部で5つ推定されたんですが（数理モデル的には無限個推定しないといけないですが、、、）、1と2以外では実質的にモデルに使われず、同じところに溜まっています。

おお！見事に民主党と共和党で分かれました！しかもx = 5, y = 5のところを見てください。共和党の議員のイデオロギーベクトルが集まっているところに、青い文字があります。お馴染みの民主党の上院議員のネルソンさんです。

一本目の記事の繰り返しになりますが、

https://en.wikipedia.org/wiki/Ben_Nelson

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

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

そもそもドメイン知識的に共和党にかなり近い民主党員です。

なのでt-SNEでの可視化に関して言えば、一本目の記事の結果をほぼ再現できているといえます。

## 離散パラメータの再現

続いて、$\eta$を再現する前に、まずzを確認しましょう：

```r
> m_senate_dper_summary |>
     dplyr::filter(stringr::str_detect(variable, "z"))
# A tibble: 5 × 7
  variable     mean    median      sd       mad           q5     q95
  <chr>       <dbl>     <dbl>   <dbl>     <dbl>        <dbl>   <dbl>
1 z[1]     0.547    0.545     0.0515  0.0533    0.465        0.631  
2 z[2]     0.447    0.449     0.0534  0.0548    0.364        0.529  
3 z[3]     0.00429  0.000735  0.0141  0.00101   0.0000181    0.0175 
4 z[4]     0.00107  0.0000967 0.00381 0.000141  0.000000566  0.00478
5 z[5]     0.000515 0.0000233 0.00336 0.0000342 0.0000000774 0.00179
```

これは何を意味しているのかというと、全体の議員の中で、1番目のイデオロギーレジームになっている人が54.7%（z[1]）いて、2番目のイデオロギーレジームになっている人が44.7%（z[2]）いるということを示しています。そして3番目から5番目はほぼゼロです。

t-SNEの解析結果を踏まえると、1番目のレジームは共和党のイデオロギーレジームで、2番目のレジームは民主党のイデオロギーレジームになります。要するに、全体の中の共和党的な人が54%、民主党的な人が44%をそれぞれ占めていると解釈できます。

では、可視化のために、zをパラメータとする多項分布から抽出された$\eta$の事後分布をまず少し整理しましょう：

```r
eta_dper_df <- m_senate_dper_summary |>
  dplyr::filter(stringr::str_detect(variable, "eta")) |>
  dplyr::mutate(
    senator_id_list = purrr::map(
      variable, \(x){
        # eta[X,Y]のXは議員ID、YはイデオロギーレジームのID
        stringr::str_split(x, "\\[|\\]|,")[[1]][2:3] |> 
          as.integer()
      }
    )
  ) |>
  tidyr::unnest_wider(col = senator_id_list, names_sep = "_") |>
  dplyr::rename(
    senator_id = senator_id_list_1,
    latent_group_id = senator_id_list_2
  ) |>
  dplyr::left_join(senator_master, by = "senator_id")
```

まずはバイデンさんの$\eta$を確認しましょう。

```r
eta_dper_df |>
  dplyr::filter(senator == "BIDEN (D DE)") |> 
  ggplot2::ggplot() + 
  ggplot2::geom_bar(ggplot2::aes(x = as.factor(latent_group_id), y = mean), 
                    stat = "identity", fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::geom_errorbar(
    ggplot2::aes(x = as.factor(latent_group_id), ymin = q5, ymax = q95),
    width = 0.2
  ) + 
  ggplot2::labs(
    title = "Posterior Distribution of eta\nJoe Biden (D DE)",
    x = "ideology regime",
    y = "distribution")
```

![biden eta.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/2957523d-53c9-95fb-d52a-23bd9c3da4c8.png)

バイデンさんは完全に民主党的な（2番目の）イデオロギーレジームに従っていることがわかりました。黒い水平線が信用区間なんですが、もう見えないくらい不確実性がないです。

では不確実性がある議員として、チェイフィーさんの$\eta$を確認しましょう：

```r
eta_dper_df |>
  dplyr::filter(senator == "CHAFEE (R RI)") |> 
  ggplot2::ggplot() + 
  ggplot2::geom_bar(ggplot2::aes(x = as.factor(latent_group_id), y = mean), 
                    stat = "identity", fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::geom_errorbar(
    ggplot2::aes(x = as.factor(latent_group_id), ymin = q5, ymax = q95),
    width = 0.2
  ) + 
  ggplot2::labs(
    title = "Posterior Distribution of eta\nLincoln Chafee (R RI)",
    x = "ideology regime",
    y = "distribution")
```

![chafee_eta.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/1ccbe240-2e02-a076-af61-7cdb6154a2e0.png)

平均でいうと共和党の（1番目の）イデオロギーレジームに従っていますが、信用区間がはっきり出ており、上限でいうと、逆の民主党の（2番目の）イデオロギーレジームに従う確率も60%を超えます。

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

https://en.wikipedia.org/wiki/Lincoln_Chafee


> Descended from a long line of moderate, center-right New England Republicans, Chafee's stances became increasingly liberal, more so than his father's positions had been.
> （ニューイングランド州共和党の穏健派中道右派の家系に生まれたチェフィー議員は、次第にその父親よりもリベラルな立場をとるようになった）

ドメイン知識的に、共和党の中でもかなり民主党よりのイデオロギーに近い人です。

最後に、お馴染みのネルソンさんの$\eta$を確認しましょう：

```r
eta_dper_df |>
  dplyr::filter(senator == "NELSON (D NE)") |> 
  ggplot2::ggplot() + 
  ggplot2::geom_bar(ggplot2::aes(x = as.factor(latent_group_id), y = mean), 
                    stat = "identity", fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::geom_errorbar(
    ggplot2::aes(x = as.factor(latent_group_id), ymin = q5, ymax = q95),
    width = 0.2
  ) + 
  ggplot2::labs(
    title = "Posterior Distribution of eta",
    x = "ideology regime",
    y = "distribution")
```

![nelson_eta.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/7cbd7da6-344d-c1b5-b9d9-cdc167735144.png)

不確実性もありますが、完全に共和党ですね。民主党の議員なのに民主党の（2番目の）イデオロギーレジームに従う確率の下限は衝撃のゼロです。

## コサイン類似度

まず、一本目の記事の方法でコサイン類似度可視化用のデータを整形します：

```r
cos_sim_dper_df <- m_senate_dper_summary |>
  dplyr::filter(stringr::str_detect(variable, "cos_sim")) |>
  dplyr::mutate(
    senator_id_list = purrr::map(
      variable, \(x){
        # cos_sim[X,Y]のXとYはそれぞれ上院議員マスターのIDに該当する
        stringr::str_split(x, "\\[|\\]|,")[[1]][2:3] |> 
          as.integer()
      }
    )
  ) |>
  tidyr::unnest_wider(col = senator_id_list, names_sep = "_") |>
  dplyr::left_join(senator_master, by = c("senator_id_list_1" = "senator_id")) |>
  dplyr::left_join(senator_master, by = c("senator_id_list_2" = "senator_id")) |>
  dplyr::rename(
    party_id_1 = party_id.x,
    party_id_2 = party_id.y,
    senator_1 = senator.x,
    senator_2 = senator.y
  ) |>
  dplyr::mutate(
    party_1 = dplyr::case_when(
      party_id_1 == 1 ~ "Republican",
      TRUE ~ "Democrat"
    ),
    party_2 = dplyr::case_when(
      party_id_2 == 1 ~ "Republican",
      TRUE ~ "Democrat"
    )
  )
```

まずはオバマさんから確認しましょう。

```r
cos_sim_dper_df |> 
  dplyr::filter(senator_1 == "OBAMA (D IL)") |> 
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = mean, y = reorder(senator_2, mean), color = party_2)) + 
  ggplot2::geom_errorbarh(ggplot2::aes(xmin = q5, xmax = q95, y = senator_2, color = party_2), position = "dodge", height = 0.2) + 
  ggplot2::scale_color_manual(values = c("Democrat" = "blue", 
                                         "Republican" = "red")) + 
  ggplot2::theme(axis.text.y = ggplot2::element_text(size = 4)) + 
  ggplot2::labs(
    title = "Cosine Similarity of Ideology Vector\nBarack Obama (D IL)",
    x = "Cosine Similarity",
    y = "Senator Name",
    color = "Party")
```

![dpre_obama_cos_sim.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/bbb91c1e-32ad-46b5-ac76-0f9d87b0213b.png)

同じように民主党員との類似度が高く、共和党員との類似度が低い構造になっています。

バイデンさんも同じです：

```r
cos_sim_dper_df |> 
  dplyr::filter(senator_1 == "BIDEN (D DE)") |> 
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = mean, y = reorder(senator_2, mean), color = party_2)) + 
  ggplot2::geom_errorbarh(ggplot2::aes(xmin = q5, xmax = q95, y = senator_2, color = party_2), position = "dodge", height = 0.2) + 
  ggplot2::scale_color_manual(values = c("Democrat" = "blue", 
                                         "Republican" = "red")) + 
  ggplot2::theme(axis.text.y = ggplot2::element_text(size = 4)) + 
  ggplot2::labs(
    title = "Cosine Similarity of Ideology Vector\nJoe Biden (D DE)",
    x = "Cosine Similarity",
    y = "Senator Name",
    color = "Party")
```

![dpre_biden_cos_sim.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/2b26d6a8-44d5-4034-7120-a60a02c56a98.png)

次はチェフィーさんのコサイン類似度を確認しましょう：

```r
cos_sim_dper_df |> 
  dplyr::filter(senator_1 == "CHAFEE (R RI)") |> 
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = mean, y = reorder(senator_2, mean), color = party_2)) + 
  ggplot2::geom_errorbarh(ggplot2::aes(xmin = q5, xmax = q95, y = senator_2, color = party_2), position = "dodge", height = 0.2) + 
  ggplot2::scale_color_manual(values = c("Democrat" = "blue", 
                                         "Republican" = "red")) + 
  ggplot2::theme(axis.text.y = ggplot2::element_text(size = 4)) + 
  ggplot2::labs(
    title = "Cosine Similarity of Ideology Vector\nLincoln Chafee (R RI)",
    x = "Cosine Similarity",
    y = "Senator Name",
    color = "Party")
```

![dpre_chafee_cos_sim.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/1e05dccd-9900-fc9e-5699-6a8e12d21abf.png)

共和党員と民主党員が混ざっている感じですね。

最後はネルソンさんを確認しましょう：

```r

cos_sim_dper_df |> 
  dplyr::filter(senator_1 == "NELSON (D NE)") |> 
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = mean, y = reorder(senator_2, mean), color = party_2)) + 
  ggplot2::geom_errorbarh(ggplot2::aes(xmin = q5, xmax = q95, y = senator_2, color = party_2), position = "dodge", height = 0.2) + 
  ggplot2::scale_color_manual(values = c("Democrat" = "blue", 
                                         "Republican" = "red")) + 
  ggplot2::theme(axis.text.y = ggplot2::element_text(size = 4)) + 
  ggplot2::labs(
    title = "Cosine Similarity of Ideology Vector\nBen Nelson (D NE)",
    x = "Cosine Similarity",
    y = "Senator Name",
    color = "Party")
```

![dpre_nelson_cos_sim.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/49c860d2-5f57-68e1-fea4-227b4a343018.png)

チェフィーさんほど混ざっていないですが、それにしても民主党なのに、コサイン類似度の高い議員の共和党率が高いですね。

# 結論

いかがでしたか？

このように、ベイズ機械学習モデルは、高い精度を達成できるだけでなく、状況に応じてドメイン知識を付与したり、ドメイン知識を発見させたりして、アカデミアとビジネスで使えるquantities of interestを複数抽出してくれる、柔軟性の高いモデルです。

皆さんも是非、ご自身のアカデミアか（と？）ビジネスで直面する課題にベイズ機械学習を活用してみてください！