次元数とクラスター数を自動選定するベイズ投票解析モデル

# はじめに

イデオロギーは古くから政治学の中心になっているのは、もはやいうまでもない。

国内政治において、右派と左派の枠組みは政治学者だけでなく、有権者や政治家同士、投資家などが政治情勢を判断する重要な視座である。

国際政治においても、冷戦期の米ソ対立は、イデオロギー抜きでは説明できない部分が多い。現代でも、アメリカをはじめとする西側諸国と中国の対立は、やはりイデオロギーで解釈した方が妥当である。

このように、投票データから投票者のイデオロギーを推定するのは、政治学において重要なタスクである。

では、我々はこのイデオロギーという概念をいかにデータを用いて定量的に評価するのか。政治学では、心理学などでも活用される項目反応理論を使う研究が多い（Clinton, Jackman and Rivers 2004）。詳細は後述するが、項目反応理論では、法案などの解答項目と、国会議員などの回答者にそれぞれベクトルが付与される。議員Aのベクトルと法案Aのベクトルが近ければ、議員Aは法案Aに賛成票を投じることでより高い効用を得るため、法案Aに賛成する傾向が強くなる。

本研究は、既存の項目反応理論を下記二点で改善し、イデオロギーをより緻密に可視化することに貢献する。

まず、既存手法は一次元のベクトルを利用することが多い。ただ、一という次元数は、あくまでも分析者の仮定に過ぎず、モデルの誤った定式化を招く恐れがある。本研究では、棒折り過程で次元数をノンパラメトリックに特定する。

次に、個々の国会議員のイデオロギーを推定するだけでなく、より広い範囲のイデオロギーレジームも合わせて推定したほうが、より緻密な政治情勢の可視化が期待される。イデオロギーのレジームは、必ずしも所属政党と同一とは限らない。よって、推定されたレジームと実際の所属政党を比較することで、政党を超えた協力関係や、異常な投票行動を取る議員を特定することができ、新たなリサーチクエスチョンの発見や比較政治理論につながる。本研究では、ディリクレ過程混合でレジーム数とそれぞれの議員が属するイデオロギーレジームをノンパラメトリックに推定する。

余談にはなるが、項目反応理論はレコメンドシステムでよく利用される協調フィルタリングと数学的に似ている。本研究の提案手法は協調フィルタリングモデルに簡単に転用できるため、転用する際のビジネス上のメリットも必要に応じて説明する。

ビジネス側と情報工学界隈で機械学習に携わる読者は、議員をユーザー（顧客）、法案をアイテム（商品）にそれぞれ置き換えて考えると問題意識がわかりやすくなる。

# モデルの定式化

ここでは、まずモデルが仮定する効用関数の形を説明する。本研究が主にフォーカスしたい次元数とイデオロギーレジームの特定は、パラメータの事前分布に作用するため、効用関数自体とは関係ない。

## 効用関数

本研究が利用する効用関数は、Clinton, Jackman and Rivers（2004）が提案した効用関数を土台にする。

まず、国会議員Sが法案Bに賛成と反対した際の効用関数をそれぞれ下記の通りに仮定する：

$$
U(senator=S,bill=B,result=Y)= -||senator_{S}-bill_{B,Y}||^2+\epsilon_{S,B,Y}
$$

$$
U(senator=S,bill=B,result=N)= -||senator_{S}-bill_{B,N}||^2+\epsilon_{S,B,N}
$$

- $senator$：議員のイデオロギーを示す無限次元ベクトル
- $bill$：法案のイデオロギーの位置を示す無限次元ベクトル
- $\epsilon$：タイプI極値分布に従う誤差項

よって、国会議員Sが法案Bに賛成票を投じる確率は

$$
P(senator=S,bill=B,result=Y) = P(U(senator=S,bill=B,result=Y)>U(senator=S,bill=B,result=N))
$$

要するに、$U(senator=S,bill=B,result=Y)-U(senator=S,bill=B,result=N)>0$の確率である。

少し左辺を整理すると

```math
\begin{align}
&=-||senator_{S}-bill_{B,Y}||^2+\epsilon_{S,B,Y}-(-||senator_{S}-bill_{B,N}||^2+\epsilon_{S,B,N})\\
&=-senator_{S}^2 - bill_{B,Y}^2 + 2 * senator_{S} * bill_{B,Y} + senator_{S}^2 + bill_{B,N}^2 - 2 * senator_{S} * bill_{B,N} + \epsilon_{S,B,Y} - \epsilon_{S,B,N}
\end{align}
```

ベクトルXの二乗は$X^2=X'X$と約束する。

ここで、$senator_{S}^2$が相殺されるので、残りは

```math
\begin{align}
&=bill_{B,N}^2 - bill_{B,Y}^2 + senator_{S} * 2 * (bill_{B,Y}-bill_{B,N}) + \epsilon_{S,B,Y} - \epsilon_{S,B,N}
\end{align}
```

さらに、$bill_{B,N}^2-bill_{B,Y}^2$と$2 * (bill_{B,Y} - bill_{B,N})$は国会議員によって変動しないため、それぞれ法案の人気度$popularity_{B}$と法案ベクトル$bill_{B}$として定義すると、

```math
\begin{align}
&P(senator=S,bill=B,result=Y) = P(popularity_{B} + senator_{S} * bill_{B} + \epsilon_{S,B,Y} - \epsilon_{S,B,N} > 0)
\end{align}
```

に整理できる。

ここで、$\epsilon_{S,B,Y}$と$\epsilon_{S,B,N}$が独立したタイプI極値分布に従うため、$P(senator=S,bill=B,result=Y)$が$popularity_{B} + senator_{S} * bill_{B}$のロジット変換で表せることが知られている（McFadden 1974）。よって、最終的に一つの国会議員と法案のペアの尤度関数は、このように表すことができる

$$
\lambda = popularity_{B} + senator_{S} * bill_{B} \\
$$

$$
P(senator=S,bill=B,result=Y) = logit^{-1}(\lambda)^{result=Y} * (1 - logit^{-1}(\lambda))^{result=N}
$$

$\lambda$の中身を詳しく確認見ると

$$
\lambda = popularity_{B} + \sum\limits_{i = 1}^{\infty} senator_{S,i} * bill_{B,i}
$$

と書けるので、$bill_{B,i}$が共変量（特徴量）、$senator_{S,i}$が傾き（ウェイト）の無限次元線形ロジスティック回帰としても解釈できる（Rudolph, Ruiz, Mandt and Blei, 2016）。

次に、本研究の新規性の中心となる$senator_{S}$と$bill_{B}$の事前分布の説明に移る。

## 事前分布

前述したように、本研究が提案する手法は、次元数とイデオロギーレジームの数とレジーム内の議員をノンパラメトリックに推定する。ノンパラメトリック推定を用いることで、分析者の恣意性を減らすことができ、データに隠された情報を抽出する際に有用である。

ここでは、まず先行研究にある類似したアプローチと本研究のアプローチを比較して、後者の新規性を説明する。

単独の項目を立てるまでもないのでここで説明するが、法案人気度$popularity_{B}$は正規分布からサンプリングされると仮定する

$$
popularity_{B} \sim Normal(0, 1)
$$

### イデオロギーレジーム

#### 先行研究との比較

議員のイデオロギーベクトルをディリクレ過程でモデリングする先行研究は既に存在する。

Spirling and Quinn（2010）は、ディリクレ過程でイギリスの下院議員のイデオロギーベクトルをクラスタリングして、政党内の投票パターンの可視化を行った。心理学でも、人の行動をパターン化するという目的で、Navorro, Griffiths, Steyvers and Lee（2006）がディリクレ過程で個人のパラメータが従う分布をモデリングした。

ただし、Spirling and Quinn（2010）とNavorro, Griffiths, Steyvers and Lee（2006）のモデル定式化では、分析対象者のパラメータが直接ディリクレ過程からサンプリングされるため、全員が少数のパラメータを共用することになる。

先行研究のような定式化だと、全員が少数のパラメータを共用することになるため、分析対象同士の類似度パターンが限定的になる。政治学的に意味のある「A議員と類似している議員」を比較すると、まず同じパラメータを共用する群の中では比較は不可能になる。もちろん、リサーチクエスチョンによっては、議員がどのグループに属するかを知るだけで十分かもしれないが、議員の差も可視化できた方が、結果としてより精緻な分析と解釈が可能になる。

また、Subhashis and van der Vaart（2017, p. 102）の説明のように、ディリクレ過程はその離散性のため、確率密度の推定で利用することは難しい。イデオロギー・選好の確率密度の推定は、国会議員や実験対象者のベクトルが事後分布的にどこまで広がりそうかを可視化する際に有用である。他にも、新しく当選した議員がどのように投票しそうかを予想する際、有限個の点から当てるより、確率密度から当てた方が柔軟性が高い。

ビジネスの世界でも、ユーザー・顧客の選好のクラスタリングは重要であるが、Aセグメントのユーザーは全員同じ選好パラメータを共用するモデルだと、結局Aセグメントの内部でパーソナライズすることはできない。本研究で提案する手法のように、セグメント内にさらに変動を許すモデルにすれば、ユーザーの理解だけでなく、高度にパーソナライズされたコミュニケーションも同時に実現できる。

#### 事前分布設計

イデオロギーレジームの事前分布にディリクレ過程を設定する。

要するに、まずは全体のハイパーパラメーター$\alpha$をガンマ分布からサンプリングする

$$
\alpha \sim Gamma(1, 1)
$$

次に、Shiraito, Lo and Olivella（2023）などの先行研究のように、棒折り過程でディリクレ過程を構築する。無限にあるレジームの中のレジームpについてこのように必要な変数をサンプリングする

$$
\pi_{p}\sim Beta(1, \alpha)
$$

$$
p_{p} = \pi_{p} \prod\limits_{l=1}^{p - 1} (1 - \pi_{l})
$$

また、レジームpの中心イデオロギーベクトルとそのばらつきは下記のようにサンプリングされる

$$
P_{latent,p} \sim Normal(0,1)
$$

$$
P_{\sigma,p} \sim Gamma(1, 1)
$$

国会議員Sのイデオロギーベクトルをサンプリングする際は、まず上記の棒折り過程で構築した無限次元の確率ベクトルをパラメーターとするカテゴリ分布から、インデックスをサンプリングする

$$
\eta_{S} \sim Categorical(p)
$$

次に、サンプリングされた$\eta_{s}$を基に、イデオロギーベクトルをサンプリングする

$$
senator_{S} \sim Normal(P_{latent, \eta_{S}}, P_{\sigma, \eta_{S}})
$$

ここでわかるように、$senator_{S}$はそれぞれ正規分布からサンプリングされるため、何個かのクラスター付近で分布することが想定されるが、Spirling and Quinn（2010）とNavorro, Griffiths, Steyvers and Lee（2006）のようにサンプリングの結果のベクトルが他の議員と完全に被ることはない。

わかりやすさのため、上記の

$$
\eta_{S} \sim Categorical(p)
$$

$$
senator_{S} \sim Normal(P_{latent, \eta_{S}}, P_{\sigma, \eta_{S}})
$$

分布を$G$で表す。

### 次元数

#### 先行研究との比較

本研究は、無限次元ベクトル$bill_{B}$の大半の要素をゼロにすることで、次元数を特定する。

先行研究でも次元数を推定する試みが見られる。初期では、Poole and Rosenthal（1991）次元数を推定する重要性を説き、次元数を増やすことでモデルの分類性能がどう変化していくのかを可視化した。より最近の研究でいうと、Kim, Londregan and Ratkovic（2018）はベイズLassoを事前分布として設定することで、一部の次元の値をゼロにする。他にも、遺伝子発現のファクターモデルでの応用例になるが、正規分布の分散を次元数の増加に従ってゼロに近づけさせる定式化もある（Bhattacharya and Dunson 2011）。

本研究では、レコメンドシステム分野の研究で提案された、棒折り過程で最初の数次元以外をゼロにする定式化を利用する（Gopalan, Ruiz, Ranganath and Blei 2014）。

ビジネスの世界では、ユーザーの選好を何次元で表せるかの知識は、ユーザー理解を超えたバリューがある。

例えば、もともとユーザーの好みを500次元でモデリングしていたが、実は38次元で表現できることがわかったら、まずユーザー情報を格納するために必要なメモリとストレージが減る。

さらに、K次元ベクトル同士の内積の計算量は$O(K)$である。そのため、500次元から38次元にベクトルを圧縮できたら、計算量は7.6%になり、場合によってGPU付きのインスタンスからCPUのみのインスタンスに切り替えることも可能になる。これによって、モデル性能を担保する上で、大幅なコスト削減と速度向上が期待される。

後述のように、学習時にやむをえず高次元ベクトルのメモリを抑えなくてはいけなくても、モデルに利用されなかった次元はモデルファイル保存時に削除することができるため、少なくとも推論時のインスタンスコストの削減と返却速度の改善ができる。

#### 事前分布設計

イデオロギーレジームの事前分布と似ているが、Gopalan, Ruiz, Ranganath, Blei（2014, p. 277）のように、まず全体のハイパーパラメータ$\beta$をガンマ分布からサンプリングする

$$
\beta \sim gamma(1, 1)
$$

次に棒折り過程で$bill_{B}$と内積を取る$dimension_{d}$を求める

$$
\delta_{d}\sim Beta(1, \beta)
$$

$$
dimension_{d} = \delta_{d} \prod\limits_{l=1}^{d - 1} (1 - \delta_{l})
$$

最後に、法案ベクトル$bill_{B}$の「標準化前」の値$bill\space unnormalized_{B}$を正規分布からサンプリングして、その結果を$dimension$と内積を取ることで$bill_{B}$を求める

$$
bill\space unnormalized_{B} \sim Normal(0, 1)
$$

$$
bill_{B} = bill\space unnormalized_{B} '* dimension
$$

本研究の表記でいうと、Gopalan, Ruiz, Ranganath, Blei（2014）の定式化では、全ての$bill_{B}$に異なる$dimension$を設定するが、本研究の場合$dimension$は全$bill_{B}$共通と仮定する。

# モデル実装

上記のモデルをStan言語で実装したコードは下記の通りである（Stan Development Team, 2024）。ただ、並列処理を実行するため、```reduce_sum```を利用した。また、Stan言語は$\eta_{S}$のような離散確率変数をパラメータとして使うことを容認しないため、周辺化を実施する。

最後の```generated quantities```ブロックでは、モデルの検証データに対する予測精度の評価指標や$\eta_{S}$、議員同士のイデオロギーベクトルのコサイン類似度、全体の平均イデオロギーベクトルとの距離、$G$の事後分布なども併せてサンプリングする。

```cpp: senate_double_dpre.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[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> dimension_alpha;
  vector<lower=0, upper=1>[K - 1] dimension_breaks;
  
  real<lower=0> regime_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[P - 1] regime_breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  vector[bill_type] bill_popularity;
  array[bill_type] vector[K] bill_latent_unnormalized;
  
  array[P] vector[K] P_latent;
  vector<lower=0>[P] P_sigma;
  
  array[senator_type] vector[K] senator_latent;
}
transformed parameters {
  simplex[K] dimension;
  simplex[P] regime; 
  array[bill_type] vector[K] bill_latent;
  {
    // stick-breaking representationの変換開始
    // https://discourse.mc-stan.org/t/better-way-of-modelling-stick-breaking-process/2530/2 を参考
    dimension[1] = dimension_breaks[1];
    real sum_1 = dimension[1];
    for (k in 2:(K - 1)) {
      dimension[k] = (1 - sum_1) * dimension_breaks[k];
      sum_1 += dimension[k];
    }
    dimension[K] = 1 - sum_1;
    // stick-breaking representationの変換開始
    // https://discourse.mc-stan.org/t/better-way-of-modelling-stick-breaking-process/2530/2 を参考
    regime[1] = regime_breaks[1];
    real sum_2 = regime[1];
    for (p in 2:(P - 1)) {
      regime[p] = (1 - sum_2) * regime_breaks[p];
      sum_2 += regime[p];
    }
    regime[P] = 1 - sum_2;
    
    
    for (b in 1:bill_type){
      for (k in 1:K){
        bill_latent[b,k] = bill_latent_unnormalized[b,k] * dimension[k];
      }
    }
  }
}
model {
  dimension_alpha ~ gamma(1, 1);
  dimension_breaks ~ beta(1, dimension_alpha);
  
  regime_alpha ~ gamma(1, 1);
  
  regime_breaks ~ beta(1, regime_alpha);
  
  for (p in 1:P){
    P_latent[p] ~ std_normal();
    P_sigma[p] ~ gamma(1, 1);
  }
  
  for (b in 1:bill_type){
    bill_popularity[b] ~ std_normal();
    bill_latent_unnormalized[b] ~ std_normal();
  }
  
  for (s in 1:senator_type){
    vector[P] case_vector;
    for (p in 1:P){
      case_vector[p] = log(regime[p]) + normal_lupdf(senator_latent[s] | P_latent[p], (P_sigma[p]));
    }
    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;
  vector[K] senator_mean_latent;
  vector[senator_type] distance_from_mean_latent;
  array[K] real posterior_distribution_of_dirichlet_process_mixture;
  {
    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] = log(regime[p]) + normal_lpdf(senator_latent[s] | P_latent[p], P_sigma[p]);
      }
      eta[s] = softmax(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));
      }
    }
    
    for (k in 1:K){
      vector[senator_type] dimension_k_collector;
      for (s in 1:senator_type){
        dimension_k_collector[s] = senator_latent[s,k];
      }
      senator_mean_latent[k] = mean(dimension_k_collector);
    }
    
    for (s in 1:senator_type){
      distance_from_mean_latent[s] = sqrt(sum((senator_latent[s] - senator_mean_latent)^2));
    }
    
    int sampled_P = categorical_rng(regime);
    posterior_distribution_of_dirichlet_process_mixture = normal_rng(P_latent[sampled_P], (P_sigma[sampled_P]));
  }
}

```

# 前処理

本研究は実証分析のデータとして、アメリカ上院の2005年から2006年までの投票データを利用する。

本研究は理論の方で無限次元のベクトルと無限次元のクラスターを仮定すると説明したが、現代の計算機は有限のものしか扱えない。そのため、無限大の代わりに、最大で考慮すべき次元数とクラスター数を高めに設定する。ここでは、最大で考慮すべき次元数を10、最大で考慮すべきクラスター数を30とする。

```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(
      # s109$codesに従ってyeaを示す1,2,3に場合に1を付与する
      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)

data_double_dpre_list <- list(
  N = nrow(rc_df_with_id[train_id,]),
  P = 30,
  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]
)

```

# モデル推定

モデル推定では、変分推論を利用する（Kucukelbir, Tran, Ranganath, Gelman and Blei, 2017）。

```R
> m_senate_double_dpre_init <- cmdstanr::cmdstan_model("estimation/senate_double_dpre.stan",
                                                      cpp_options = list(
                                                          stan_threads = TRUE
                                                      )
 )
> m_senate_double_dpre_estimate <- m_senate_double_dpre_init$variational(
     seed = 12345,
     threads = 6,
     iter = 50000,
     data = data_double_dpre_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.00716 seconds 
1000 transitions using 10 leapfrog steps per transition would take 71.6 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       -14528.229             1.000            1.000 
   200       -13995.097             0.519            1.000 
   300       -13840.500             0.350            0.038 
   400       -13766.735             0.264            0.038 
   500       -13675.984             0.212            0.011 
   600       -13609.310             0.178            0.011 
   700       -13586.733             0.153            0.007   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  15.3 seconds.
> m_senate_double_dpre_summary <- m_senate_double_dpre_estimate$summary()
```

次に、可視化に必要なパラメータを保存する

```R
P_double_dpre_latent <- m_senate_double_dpre_summary |>
  dplyr::filter(stringr::str_detect(variable, "P_latent\\[")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 10)

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

bill_double_dpre_latent <- m_senate_double_dpre_summary |>
  dplyr::filter(stringr::str_detect(variable, "bill_latent\\[")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 10)

eta_double_dpre_df <- m_senate_double_dpre_summary |>
  dplyr::filter(stringr::str_detect(variable, "eta")) |>
  dplyr::mutate(
    senator_id_list = purrr::map(
      variable, \(x){
        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")

```

# 結果の可視化

## 精度確認

検証データに対する分類制度に関しては、筆者が以前記事で紹介した、次元数を固定するモデルより「統計的に有意に」低い結果になっているが、定性的には両方とも客観的に高い精度を達している。

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

```R
m_senate_double_dpre_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")
```

![f1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/cef56fa4-9b3f-a08a-aef1-216b31a44786.png)

詳細は後述するが、今回のモデルは2次元だけで以前の記事の10次元のモデルと実質的に同じクラスの精度を達成できたのは、以前のモデルでは次元数が高すぎたせいで、モデルがデータの形を丸暗記してしまった可能性を示唆する。

## 次元数

まず、何件かの$bill_{B}$の棒グラフを確認する

```R
bill_double_dpre_latent[1,] |> 
  barplot(main = "bill 1", xlab = "dimension", ylab = "value", col = ggplot2::alpha("blue", 0.3))
```

![b1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/1fea7c26-b020-a25f-969e-4c5b8665ad39.png)

```R
bill_double_dpre_latent[12,] |> 
  barplot(main = "bill 12", xlab = "dimension", ylab = "value", col = ggplot2::alpha("blue", 0.3))
```

![b12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/70a67096-62bf-058a-5e72-b2ea955be49d.png)

```R
bill_double_dpre_latent[51,] |> 
  barplot(main = "bill 51", xlab = "dimension", ylab = "value", col = ggplot2::alpha("blue", 0.3))
```

![b51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/d401133f-f92e-3ee3-35aa-e98d4173de97.png)

ここでわかるように、一番目と五番目の次元以外はほぼゼロになっている。

$dimension$の事後分布も併せて確認すると

```R
m_senate_double_dpre_summary |>
  dplyr::filter(stringr::str_detect(variable, "dimension\\[")) |>
  dplyr::bind_cols(
    dimension = 1:10
  ) |>
  ggplot2::ggplot() +
  ggplot2::geom_bar(ggplot2::aes(x = as.factor(dimension), y = mean), 
                    stat = "identity", fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::geom_errorbar(
    ggplot2::aes(x = as.factor(dimension), ymin = q5, ymax = q95),
    width = 0.2
  ) + 
  ggplot2::labs(
    title = "Posterior Distribution of dimension",
    x = "dimension id",
    y = "distribution")
    
```

![g_dimension.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/c404502f-8ee8-9039-eec7-6c303b5339c7.png)

ここからもわかるように、モデルは一番目の次元と五番目の次元以外の利用を棄却している。しかも信用区間も短いため、二つの次元でアメリカ上院の2005年から2006年までの投票データを説明できるのは、いわゆる「統計的に有意な」結果といえよう。

## クラスター数

次に、クラスター数の可視化に入る。

まず、$p$の事後分布から確認する

```R
m_senate_double_dpre_summary |>
    dplyr::filter(stringr::str_detect(variable, "regime\\[")) |>
    dplyr::bind_cols(
        dimension = 1:30
    ) |>
    ggplot2::ggplot() +
    ggplot2::geom_bar(ggplot2::aes(x = as.factor(dimension), y = mean), 
                      stat = "identity", fill = ggplot2::alpha("blue", 0.3)) + 
    ggplot2::geom_errorbar(
        ggplot2::aes(x = as.factor(dimension), ymin = q5, ymax = q95),
        width = 0.2
    ) + 
    ggplot2::labs(
        title = "Posterior Distribution of p",
        x = "ideology regime id",
        y = "distribution")
```

![g_p.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/bcd10fe8-71c8-66f2-f558-9861c33ef110.png)

ほとんどの上院議員のイデオロギーベクトルが、イデオロギーレジーム2と6から抽出されたことがわかる。

では、何人かの議員の$\eta$を確認しよう

```R
eta_double_dpre_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 \nBIDEN (D DE)",
        x = "ideology regime id",
        y = "distribution")
```

![g_biden.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/89e7d054-74bc-505b-730a-6c265a40b65c.png)

バイデンのイデオロギーベクトルがほぼ100%の確率でイデオロギーレジーム2からサンプリングされた。

```R
eta_double_dpre_df |>
    dplyr::filter(senator == "OBAMA (D IL)") |> 
    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 \nOBAMA (D IL)",
        x = "ideology regime id",
        y = "distribution")
```

![g_obama.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/9e634e8c-fbbd-1e40-3252-8fcd280a932e.png)

オバマも同様である。後でイデオロギーベクトルを見る際にまた説明するが、イデオロギーレジーム2は民主党のイデオロギーを示しているといえる。

続いて、共和党側を確認する

```R
eta_double_dpre_df |>
    dplyr::filter(senator == "MCCAIN (R AZ)") |> 
    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 \nMCCAIN (R AZ)",
        x = "ideology regime id",
        y = "distribution")
```

![g_mccain.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/4b84f518-bedd-31b2-0bc4-a59f596138b4.png)

McCainのイデオロギーベクトルがほぼ90%くらいの確率でイデオロギーレジーム6からサンプリングされた。

他の共和党の議員も概ね同じ結果で、イデオロギーレジーム6からイデオロギーベクトルがサンプリングされた確率が一番高い。

Craigもそうである。

```R
eta_double_dpre_df |>
    dplyr::filter(senator == "CRAIG (R ID)") |> 
    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 \nCRAIG (R ID)",
        x = "ideology regime id",
        y = "distribution")
```

![g_craig.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/c90a0eb0-e7eb-cf4b-6ee4-4403d1f2768f.png)

このように、イデオロギーレジーム2が民主党、6が共和党を示していて、全ての議員が自身の所属政党のイデオロギーレジームからサンプリングされたイデオロギーベクトルを持つのかというと、そうとは限らない。

例えば、Nelsonの場合、民主党藉にもかかわらず、共和党を示すイデオロギーレジーム6からイデオロギーベクトルがサンプリングされた確率が一番高い

```R
eta_double_dpre_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 \nNELSON (D NE)",
        x = "ideology regime id",
        y = "distribution")
```

![g_nelson.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/62fab4f9-e73c-c557-07d5-65b806e5aeb4.png)

Nelsonの$\eta$の事後分布の2番と6番の差はいわゆる「統計的有意」ではないが、少なくとも民主党のイデオロギーレジームから来たのか、それとも共和党のイデオロギーレジームから来たのかの判断がつかないくらい、民主党との共通点が低く、共和党的な投票行動をしているといえる。

## イデオロギーベクトル

最後に、上院議員のイデオロギーベクトルの1番目の次元と5番目の次元（実際にモデルに利用された次元）、並びに更新された事前分布$G$を可視化する。$G$はイデオロギーベクトルの分布として解釈できる。

```R
posterior_distribution_of_dirichlet_process_mixture <- m_senate_double_dpre_estimate$draws("posterior_distribution_of_dirichlet_process_mixture") |>
  tibble::as_tibble() |> 
  dplyr::select(dplyr::matches("\\[1\\]|\\[5\\]")) |>
  `colnames<-`(c("x1", "x2")) |>
  dplyr::mutate(
    x1 = as.numeric(x1),
    x2 = as.numeric(x2)
  )

tibble::tibble(
  x1 = c(P_double_dpre_latent[,1], senator_double_dpre_latent[,1]),
  x2 = c(P_double_dpre_latent[,5], senator_double_dpre_latent[,5]),
  name = c(as.character(1:30), senator_master$senator),
  party = c(rep("latent", 30), c("R", "D")[senator_master$party_id]),
  label_size = c(rep(1, 30), rep(0,8, nrow(senator_master)))
) |>
  ggplot2::ggplot() + 
  ggplot2::geom_text(ggplot2::aes(x = x1, y = x2, label = name, color = party, size = label_size)) + 
  ggplot2::geom_density_2d(data = posterior_distribution_of_dirichlet_process_mixture, 
                           ggplot2::aes(x = x1, y = x2),
                           color = ggplot2::alpha("black", 0.2),
                           show.legend = FALSE
                           ) +
  ggplot2::scale_color_manual(
    values = c(
      "D" = ggplot2::alpha("blue", 0.6),
      "R" = ggplot2::alpha("red", 0.6),
      "latent\ngroup" = ggplot2::alpha("black", 0.2)
    )
  ) + 
  ggplot2::guides(size = "none") + 
  ggplot2::xlab("dimension 1") + 
  ggplot2::ylab("dimension 5") + 
  ggplot2::ggtitle("Visualization of Estimated Ideology Vector")
```

![ideology_and_G.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/645b595f-7a5a-981e-232f-eab5fa8cf1f4.png)

この図からわかるように、5番目の次元は、おそらく民主党・共和党の軸に該当する。また、y軸の5と-5のところにそれぞれ$G$の山がある。よって、クラスタリングするディリクレ過程は、政党のイデオロギーを軸に議員を分類したと判断できる。

ここで強調しておきたいのは、ほとんどの議員が自身の所属政党の側にある中で、反対側にいるNelsonをピンポイントで可視化できたのも、このモデルの性能の良さと政治学における有用性を示す。

次に、1番目の次元を確認する。$G$の二つの山の頂点はx軸においてほぼ垂直線になっているので、1番目の次元は政党とあまり関係のない投票傾向だと思われる。

1番目の次元の値が最も高い議員と最も低い議員をそれぞれ抽出する

```R
> senator_master |> dplyr::bind_cols(e = senator_double_dpre_latent[,1]) |> dplyr::arrange(-e)
# A tibble: 102 × 4
   senator         senator_id party_id      e
   <chr>                <int>    <dbl>  <dbl>
 1 COBURN (R OK)           23        1  0.479
 2 SUNUNU (R NH)           95        1  0.415
 3 MCCAIN (R AZ)           70        1  0.144
 4 GREGG (R NH)            47        1 -0.190
 5 KYL (R AZ)              60        1 -0.301
 6 ENSIGN (R NV)           40        1 -0.310
 7 INHOFE (R OK)           52        1 -0.318
 8 FEINGOLD (D WI)         42        2 -0.626
 9 DEMINT (R SC)           33        1 -0.784
10 SESSIONS (R AL)         88        1 -1.01 
# ℹ 92 more rows
# ℹ Use `print(n = ...)` to see more rows
> senator_master |> dplyr::bind_cols(e = senator_double_dpre_latent[,1]) |> dplyr::arrange(e)
# A tibble: 102 × 4
   senator            senator_id party_id     e
   <chr>                   <int>    <dbl> <dbl>
 1 PRYOR (D AR)               79        2 -9.92
 2 LINCOLN (D AR)             66        2 -9.37
 3 COCHRAN (R MS)             24        1 -8.17
 4 MIKULSKI (D MD)            73        2 -7.80
 5 BENNETT (R UT)              7        1 -7.74
 6 ROCKEFELLER (D WV)         83        2 -7.56
 7 STEVENS (R AK)             94        1 -7.03
 8 MURKOWSKI (R AK)           74        1 -6.82
 9 DOMENICI (R NM)            37        1 -6.72
10 MURRAY (D WA)              75        2 -6.65
# ℹ 92 more rows
# ℹ Use `print(n = ...)` to see more rows
```

確かに上位でも下位でも民主党と共和党の議員が混在している。まず、上位の議員の共通点を確認すると、Coburn、Sununu、McCainはこの[記事](https://www.freedomworks.org/press/freedomworks-salutes-senate-fiscal-watch-team/)の説明のように、「財政監視チーム」（Fiscal Watch Team）に所属する議員で、財政支出に対しては保守的な思想を持つことで知られている。

Grimmer（2010）の研究は、McCainとDemintが提案し、利益供給系の法案を一時延期するEarmark Reformに、有権者への利益供給に依存する議員が反対したことを示した。そのMcCainとDemintが1番目の次元の値が最も高い議員のトップ10に入る。

1番目の次元が最も低い議員の方を確認する。本記事で利用するデータに含まれていないが、Pryorが[オバマケア](https://www.healthreformvotes.org/congress/300080)に賛成した。これは政府の財政支出の増加を容認する態度で、「財政監視チーム」らとは真逆の思想である。

また、[ポリティコ](https://www.politico.com/story/2014/05/earmark-debate-congress-106470)の記事によると、同じく1番目の次元の値が低い共和党のCochranは、ミシシッピ州に利益をもたらすことで有名になっていると説明する。同記事は、財政政策の思想が共和党を分断させているとも指摘している。

よって、より詳細な検証が必要だが、財政政策の思想が、所属政党のイデオロギーだけでは説明しきれない5番目の次元の正体の可能性が高い。

Poole and Rosenthal（1991）は、政党の次に強い次元は市民権と人種問題と説明したが、本研究の結果だと財政政策の思想の次元が二番目に強いという結果になっている。これは市民権と人権問題などが落ち着いたからか、それとも2005年から2006年の特殊な状況なのかは、今後の研究では分析対象とする期間を伸ばすことで検証したい。

## 極端な議員の洗い出し

最後に、```generated quantities```ブロックで計算した、議員と平均イデオロギーベクトルのユークリッド距離を可視化する

```R
m_senate_double_dpre_summary |>
    dplyr::filter(stringr::str_detect(variable, "distance_from_mean_latent\\[")) |> 
    dplyr::bind_cols(senator = senator_master$senator,
                     party = ifelse(senator_master$party_id == 1, "Republican", "Democrat")) |>
    ggplot2::ggplot() +
    ggplot2::geom_point(ggplot2::aes(x = mean, y = reorder(senator, mean), color = party)) + 
    ggplot2::geom_errorbarh(ggplot2::aes(xmin = q5, xmax = q95, y = senator, color = party), 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 = "Distance from the Mean Ideology Vector",
        x = "Distance",
        y = "Senator Name",
        color = "Party")
```

![dist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/e2cc175f-f66a-e90c-1cf2-d23218aa6751.png)

このように、イデオロギーベクトルの平均からの距離と所属政党の関係性はあまりなく、かつ離散的に遠い議員もいない。よって、少なくともデータの期間においては、極端な立場を持つ上院議員はいなかったといえよう。

# 結論

本研究は、棒折り過程で次元数を、ディリクレ過程混合でクラスター数をノンパラメトリックベイズで推定する項目反応理論モデルを提案した。

ノンパラメトリックベイズは分析者の細かい仮説を必要としないため、仮説や理論の発見で政治学とビジネスに貢献できるが、その出てきた仮説の意義を検証する必要がある。今後の研究では、アメリカの上院議員などの財政政策の態度をより深く追う予定である。

また、本研究で利用するデータの期間が短く、しかもモデルに時系列的な要素を取り入れていないため、モデルが利用する次元は時間と共に変動するのかを検証できない。

モデルが利用する次元とは議会での議論の軸にもなるため、その経時的変化を可視化するのは、政治学的に有意義のはいうまでもない。よって、今後はMartin and Quinn（2002）などの先行研究を参考に、モデルに時系列的な要素を取り入れる方法を検討する。

解釈可能性の観点から、本研究が提案したモデルは法案内容をある種ブラックボックス化した影響で、結果の説明が難しい問題がある。ここに関しては、Kim, Londregan and Ratkovic（2018）やVafa, Suresh and Blei（2020）、Lauderdale and Clark（2014）のように、テキストデータを取り入れるモデル拡張で対応する。

最後に、細かいところになるが、そもそもこの複雑なノンパラメトリックベイズは数学的に存在しうる・定義しうるかより厳密に担保する根拠を示したい。

# 参考文献

Bhattacharya, Anirban, and David B. Dunson. "Sparse Bayesian infinite factor models." Biometrika 98.2 (2011): 291-306.

Clinton, Joshua, Simon Jackman, and Douglas Rivers. "The statistical analysis of roll call data." American Political Science Review 98.2 (2004): 355-370.

Ghosal, Subhashis, and Aad W. van der Vaart. Fundamentals of nonparametric Bayesian inference. Vol. 44. Cambridge University Press, 2017.

Gopalan, Prem, Francisco J. R. Ruiz, Rajesh Ranganath and David M. Blei. "Bayesian nonparametric poisson factorization for recommendation systems." Artificial Intelligence and Statistics. PMLR, 2014.

Grimmer, Justin. "A Bayesian hierarchical topic model for political texts: Measuring expressed agendas in Senate press releases." Political Analysis 18.1 (2010): 1-35.

Kim, In Song, John Londregan, and Marc Ratkovic. "Estimating spatial preferences from votes and text." Political Analysis 26.2 (2018): 210-229.

Kucukelbir, Alp, Dustin Tran, Rajesh Ranganath, Andrew Gelman and David Blei. "Automatic differentiation variational inference." Journal of machine learning research 18.14 (2017): 1-45.

Lauderdale, Benjamin E., and Tom S. Clark. "Scaling politically meaningful dimensions using texts and votes." American Journal of Political Science 58.3 (2014): 754-771.

Martin, Andrew D., and Kevin M. Quinn. "Dynamic ideal point estimation via Markov chain Monte Carlo for the US Supreme Court, 1953–1999." Political analysis 10.2 (2002): 134-153.

McFadden, Daniel. "Conditional Logit Analysis of Qualitative Choice Behavior." Frontiers in Econometrics, Academic Press, 1974.

Navarro, Daniel J., Thomas L. Griffiths, Mark Steyvers and Michael D. Lee. "Modeling individual differences using Dirichlet processes." Journal of Mathematical Psychology 50.2 (2006): 101-122.

Poole, Keith T., and Howard Rosenthal. "Patterns of congressional voting." American journal of political science (1991): 228-278.

Rudolph, Maja, Francisco J. R. Ruiz, Stephan Mandt and David M. Blei. "Exponential family embeddings." Advances in Neural Information Processing Systems 29 (2016).

Shiraito, Yuki, James Lo, and Santiago Olivella. "A Nonparametric Bayesian Model for Detecting Differential Item Functioning: An Application to Political Representation in the US." Political Analysis 31.3 (2023): 430-447.

Spirling, Arthur, and Kevin Quinn. "Identifying intraparty voting blocs in the UK House of Commons." Journal of the American Statistical Association 105.490 (2010): 447-457.

Stan Development Team. 2024. Stan Modeling Language Users Guide and Reference Manual,2.34.1. https://mc-stan.org

Vafa, Keyon, Suresh Naidu, and David M. Blei. "Text-based ideal points." arXiv preprint arXiv:2005.04232 (2020).