階層棒折り過程ベイズ機械学習モデルで冷戦史を推定してみた：第一報

# はじめに

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

https://dippeople.dip-net.jp/18130/

今回の記事は、2024年夏の[計量・数理政治学会（JSQPS）](https://sites.google.com/view/jsqps/2024-summer-meeting-program)のプレゼン発表でいただいたフィードバックをもとに改善・進化させた研究の第一報です。第一報なので、詳しい先行研究のレビューはまだありません。整理中です。

https://developer.dip-net.co.jp/entry/2024/07/22/JSQPS_2024_Summer_Meeting%E3%81%AB%E3%81%A6%E5%BC%8A%E7%A4%BE%E3%81%AE%E3%83%87%E3%83%BC%E3%82%BF%E3%82%B5%E3%82%A4%E3%82%A8%E3%83%B3%E3%83%86%E3%82%A3%E3%82%B9%E3%83%88%E3%81%8C%E7%99%BA%E8%A1%A8

世界を代表する政治学および関連分野の研究者の皆様と交流する貴重な機会を提供してくださったJSQPSに、心より感謝申し上げます。また、その参加を様々な面でご支援いただいたディップ株式会社、特に私の上長と同僚の皆様にも、深く感謝いたします。

https://dippeople.dip-net.jp/18434/

https://dippeople.dip-net.jp/18731/


2024年夏のJSQPSで発表した内容はこちらの記事から確認できます：

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

前回の研究からの改善・進化なので当たり前ですが、特にモデル定式化のところでかなりの内容が被ります。

具体的な内容としては、**イデオロギーベクトルの次元数が時間と共に変化することを許容する**モデルを開発して、国連総会の投票データに応用しました。

暫定の分析結果としてここで共有したいのは

- 冷戦中より、冷戦後の方が国際社会の争いが単純
- ソ連とそのコアの同盟国は、冷戦中孤立されていた
- アメリカと他の西側諸国は一枚岩ではなかった

では、早速内容に入ります！

# モデル定式化

ここでは、モデルの背景にある国家の選好のミクロ経済学的な基礎づけ（microfoundation）、イデオロギーレジーム、次元数とその変化、そして確率モデル化した結果とその実装方法を説明します。

詳細に入る前に、まず全体像から説明します。

このモデルは、要するに国家と国連の決議案をベクトル化して、その内積から法案に対する賛成・欠席・反対を予測するモデルです。機械学習界隈の読者は、国家に決議案をレコメンドする協調フィルタリングをやっていると理解していただいても構いません。

ベクトル化は簡単です。ただ、何次元がいいでしょうか？さらに、次元数は時間と共にどう変動するのでしょうか？これを解決するために、**年固有の次元の重要度ベクトル**を導入しました。

年固有の次元の重要度ベクトルとは、その年のイシューの重要度のようなものです。

下の図で説明します：

![Screen Shot 2024-09-22 at 17.06.09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/8bb22290-fac4-4fc1-3de0-91b42ba97472.png)


この図では、A国が1951年と1952年で表現するイデオロギーベクトルは異なります。ただ、これはA国のイデオロギーが変動するのではなく、A国のイデオロギーベクトル自体は安定していて、変わったのはその年にイシューであると仮定します。

より具体的な事例を挙げると、A国はB国が大嫌いで、これは第五次元の低さで表現されています。ところが、1951年でも1952年でも、会議でB国関連のことが話題になることが一切なく、第五次元の-3.1がある意味**起動される場がない**、**表現できない**です。ただ、1951年と1952年でB国関連の話が話題になっていないからといって、A国のB国に対する態度がニュートラルになったといっていいのかというと、おそらく違うと思います。

急に話が変わりますが、レコメンドシステムでも同じことが起こります。ウイスキー通販サイトのレコメンドシステムを想像していただきたいです。何かしらの事情で、超有名なバーボンウイスキーの入荷が急に止まりました。これにより、ユーザーが超有名なバーボンウイスキーを選択することがなくなります。ただ、ユーザーが選ばなくなっても、そのバーボンウイスキーの人気度が下がったのではなく、普通に在庫が無くなってユーザーが自身の好みを**起動できない**、**表現できない**だけです。なので、モデリングの時も、ユーザーの好みが変わったのではなく、ユーザーの持つ選好ベクトルに重みがかけられたと定式化した方が、より真の好みを抽出できると思います。

なので、本記事で提案するモデルは、Bailey（2007）や、Bailey, Strezhnev and Voeten（2017）、Martin and Quinn（2002）などの、イデオロギーベクトル・理想点（ideal point）自体が時間と共に変動すると仮定する政治学の先行研究とは大きく異なります。本記事のモデルにおいて、イデオロギーベクトルは時間と共に変動しません。

## 国家の効用関数と選択行動

本研究で利用する国家の効用関数は、Bailey, Strezhnev and Voeten（2017）の定式化を参考にしています。

国家は、**無限次元**のイデオロギーベクトルを持っています。そして、各決議案もまた、**無限次元**のベクトルで表現できます。さらに、年固有の**無限次元**の重要度ベクトルもあります。

まず、Gopalan, Ruiz, Ranganath and Blei（2014）を参考に、イデオロギーベクトルの各要素に、重要度ベクトルで重みをつけて、国家SがY年て表現されるイデオロギーベクトルとします：

$$
表現イデオロギー_{S,Y,d} = イデオロギーベクトル_{S,d} * 年重要度ベクトル_{Y,d};\space d=1から∞まで
$$

次に、表現イデオロギーベクトルと決議案ベクトルの内積を取って、スコアを計算します。

$$
Score_{S,R} = 表現イデオロギーベクトル_{S} '* 決議案ベクトル_{R}
$$

このスコアと各決議案に設定された閾値を下記のように掛け合わせて、国家の投票結果を予測します：

$$
結果_{S,R}=反対\space IF \space スコア_{S.R} < 閾値_{R,1}
$$

$$
結果_{S,R}=棄権\space IF \space 閾値_{R,1} < スコア_{S.R} < 閾値_{R,2}
$$

$$
結果_{S,R}=賛成\space IF \space 閾値_{R,3} < スコア_{S.R}
$$

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

イデオロギーレジームというのは、要するに各国が持つイデオロギーベクトルの事前分布のことです。この分布のことをGと言います。

もちろん、正規分布などを仮定するのも問題ないですが、正規分布のような既知の分布だと、国家のイデオロギーの複雑な様子をうまく捉えられない可能性があります。なので、本研究では、ディリクレ過程を利用して、G自体を推定します。

もちろん、ディリクレ過程を用いて選好をモデリングする先行研究はすでに存在しています。

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

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

このような定式化では、分析対象同士の類似度パターンが限定され、全員が少数のパラメータを共用することになります。特に、政治学的に重要な「A議員に類似した議員」、「A国に類似した国」を比較する際、同一パラメータを共有するグループ内での比較は困難です。もちろん、研究の目的によっては、議員がどのグループに属しているかを知るだけで十分な場合もありますが、議員間の差異をも可視化できた方が、より精密な分析と解釈が可能になるでしょう。

さらに、Subhashis and van der Vaart（2017, p. 102）が指摘するように、ディリクレ過程はその離散的な性質から、確率密度の推定には適していません。イデオロギーや選好の確率密度の推定は、議員や被験者のベクトルが事後分布においてどの程度広がるかを可視化する際に有用です。また、新たに当選した議員がどのように投票するかを予測する際も、有限個の点から予測するよりも、確率密度を用いた方が柔軟な予測が可能です。

### 分布設計

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

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

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

次に、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(0.001, 0.001)
$$

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

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

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

$$
イデオロギーベクトル_{S} \sim Normal(P_{latent, \eta_{S}}, P_{\sigma, \eta_{S}})
$$

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

この一連のサンプリングのプロセス分布を$G$で表します。

## 次元数の変化の定式化

### 棒折り過程の使い方

次元数の推定に関する試みは、先行研究でも多く見られます。初期の研究として、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のみのインスタンスに移行できる可能性があり、モデル性能を維持しつつ、コスト削減と処理速度の向上が期待されます。

また、後述するように、学習時に高次元ベクトルのメモリを削減する必要が生じた場合でも、モデルの保存時に利用されなかった次元を削除することで、少なくとも推論時のインスタンスコストの削減とレスポンス速度の改善が可能です。

### 棒折り過程の経時的変化

ここでは、棒折り過程に経時的変化を許容するため、階層ディリクレ過程（Teh, Jordan, Beal, and Blei, 2006）の考えを導入します。

Teh, Jordan, Beal, and Blei（2006）は、純粋にグループごと（年ごと）にディリクレ（棒折り）過程を設定するだけだと、そもそもグループ間の次元は異なるものになってしまうと指摘しました。その理由は、ディリクレ過程が正規分布のような連続確率分布から離散の値を抽出するため、ディリクレ過程同士の値に重なりが発生する確率は0になると説明しました。

この問題を解決するために、Teh, Jordan, Beal, and Blei（2006）は階層ディリクレ過程を提案しました。

### 分布設計

まず、最初の年の重要度ベクトルを棒折り過程でサンプリングします。具体的には、

$$
\beta \sim Gamma(0.001, 0.001)
$$

をサンプリングし、次に棒折り過程で$年重要度ベクトル_{1,d}$を求める

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

$$
年重要度ベクトル_{1,d} = \delta_{1,d} \prod\limits_{l=1}^{d - 1} (1 - \delta_{1,l})
$$

これをd = 1からd = ∞まで繰り返します。

次に、経時的変化の度合いを示すパラメータをサンプリングします：

$$
\theta\sim Gamma(0.001, 0.001)
$$

そして、Y年の重要度ベクトルを、Y-1年の重要度ベクトルを「事前分布」とする階層棒おり過程からサンプリングします。

$$
\delta_{Y,d} \sim Beta\left( \theta \delta_{Y-1,d}, \theta \left(1 - \sum_{l=1}^{d} \delta_{Y-1,l} \right) \right)
$$

$$
年重要度ベクトル_{Y,d} = \delta_{Y,d} \prod\limits_{l=1}^{d - 1} (1 - \delta_{Y,l})
$$

これもd = 1からd = ∞まで繰り返します。

# モデル実装

実装用のStanコードはこちらです。重要なQOIの事後分布もサンプリングされます。

```c++:dynamic_dimension.stan
functions {
  vector stick_breaking(vector breaks){
    int length = size(breaks) + 1;
    vector[length] result;
    
    result[1] = breaks[1];
    real summed = result[1];
    for (d in 2:(length - 1)) {
      result[d] = (1 - summed) * breaks[d];
      summed += result[d];
    }
    result[length] = 1 - summed;
    
    return result;
  }
  
  real weighted_inner_product(
    vector vec_a, vector vec_b, vector weight
  ){
    vector[size(vec_a)] result_record;
    for (i in 1:size(vec_a)){
      result_record[i] = vec_a[i] * vec_b[i] * weight[i];
    }
    return sum(result_record);
  }
  
  real cosine_similarity(
    vector vec_a, vector vec_b
  ){
    return (vec_a '* vec_b)/sqrt(sum(vec_a.^2) * sum(vec_b.^2));
  }
  
  real entropy(
    vector distribution, real threshold
  ){
    vector[size(distribution)] result;
    
    for (i in 1:size(distribution)){
      if (distribution[i] < threshold){
        result[i] = 0.0;
      } else {
        result[i] = distribution[i] * log(distribution[i]);
      }
    }
    
    return -sum(result);
  }
  
  real partial_sum_lpmf(
    array[] int result,
    
    int start, int end,
    
    array[] int country, array[] int resolution, array[] int year,
    
    array[] vector cutpoints,
    array[] vector dimension_year,
    vector country_propensity, array[] vector country_latent,
    array[] vector resolution_latent
  ){
    vector[end - start + 1] lambda;
    int count = 1;
    for (i in start:end){
      lambda[count] = ordered_logistic_lpmf(
                        result[count] | 
                        country_propensity[country[i]] + 
                        weighted_inner_product(country_latent[country[i]], resolution_latent[resolution[i]], dimension_year[year[i]]), 
                        cutpoints[resolution[i]]
                      );
      count += 1;
    }
    return sum(lambda);
  }
}
data {
  int dimension_type;
  int group_type;
  
  int country_type;
  int resolution_type;
  int year_type;
  int result_type;
  
  array[resolution_type] int resolution_year;
  
  int N;
  array[N] int country;
  array[N] int resolution;
  array[N] int year;
  array[N] int result;
  
  int val_N;
  array[val_N] int val_country;
  array[val_N] int val_resolution;
  array[val_N] int val_year;
  array[val_N] int val_result;
}
parameters {
  real<lower=0> dimension_global_alpha;
  
  real<lower=0> dimension_across_year_alpha;
  array[year_type] vector<lower=0, upper=1>[dimension_type - 1] dimension_year_breaks;
  
  real<lower=0> group_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[group_type - 1] group_breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  vector<lower=0>[group_type] group_sigma;
  array[group_type] vector[dimension_type] group_latent;
  
  vector[country_type] country_propensity;
  array[country_type] vector[dimension_type] country_latent;
  array[resolution_type] vector[dimension_type] resolution_latent;
  
  array[resolution_type] ordered[result_type - 1] cutpoints;
}
transformed parameters {
  array[year_type] simplex[dimension_type] dimension_year;
  simplex[group_type] group;
  
  for (t in 1:year_type){
    dimension_year[t] = stick_breaking(dimension_year_breaks[t]);
  }
  
  group = stick_breaking(group_breaks);
}
model {
  dimension_global_alpha ~ gamma(0.001, 0.001);
  
  dimension_year_breaks[1] ~ beta(1, dimension_global_alpha);
  
  dimension_across_year_alpha ~ gamma(0.001, 0.001);
  
  for (t in 2:year_type){
    for (d in 1:(dimension_type - 1)){
      dimension_year_breaks[t, d] ~ beta(dimension_across_year_alpha * dimension_year[t - 1, d], dimension_across_year_alpha * (1 - sum(dimension_year[t - 1, 1:d])));
    }
  }
  
  group_alpha ~ gamma(0.001, 0.001);
  
  group_breaks ~ beta(1, group_alpha);
  
  country_propensity ~ normal(0, 1);
  
  group_sigma ~ gamma(0.001, 0.001);
  
  for (g in 1:group_type){
    group_latent[g] ~ normal(0, 1);
  }
  
  for (c in 1:country_type){
    vector[group_type] case_vector;
    for (g in 1:group_type){
      case_vector[g] = log(group[g]) + normal_lpdf(country_latent[c] | group_latent[g], group_sigma[g]);
    }
    target += log_sum_exp(case_vector);
  }
  
  for (r in 1:resolution_type){
    resolution_latent[r] ~ normal(0, 1);
  }
  
  int grainsize = 1;
  
  target += reduce_sum(
    partial_sum_lupmf, result,
    
    grainsize,
    
    country, resolution, year,
    
    cutpoints,
    dimension_year,
    country_propensity, country_latent,
    resolution_latent
  );
}
generated quantities {
  array[dimension_type] real drawn_G;
  array[year_type, dimension_type] real drawn_G_weighted_yearly;
  array[year_type, country_type] vector[dimension_type] country_latent_weighted_yearly;
  array[resolution_type] vector[dimension_type] resolution_latent_weighted_yearly;
  array[country_type, country_type] real cos_sim;
  array[country_type] vector[group_type] estimated_eta;
  vector[year_type] entropy_per_year;
  array[val_N] int predicted;
  vector[result_type] F1_score_by_class;
  real F1_score_macro;
  
  {
    int sampled_group = categorical_rng(group);
    drawn_G = normal_rng(group_latent[sampled_group], group_sigma[sampled_group]);
    
    for (t in 1:year_type){
      for (d in 1:dimension_type){
        drawn_G_weighted_yearly[t, d] = drawn_G[d] * dimension_year[t, d];
        for (c in 1:country_type){
          country_latent_weighted_yearly[t, c, d] = country_latent[c, d] * dimension_year[t, d];
        }
      }
    }
  }
  
  for (i in 1:resolution_type){
    for (j in 1:dimension_type){
      resolution_latent_weighted_yearly[i, j] = resolution_latent[i, j] * dimension_year[resolution_year[i], j];
    }
  }
  
  for (c1 in 1:country_type){
    for (c2 in 1:country_type){
      cos_sim[c1, c2] = cosine_similarity(country_latent[c1], country_latent[c2]);
    }
  }
  
  for (c in 1:country_type){
    vector[group_type] case_vector;
    for (g in 1:group_type){
      case_vector[g] = log(group[g]) + normal_lpdf(country_latent[c] | group_latent[g], group_sigma[g]);
    }
    estimated_eta[c] = softmax(case_vector);
  }
  
  for (i in 1:year_type){
    entropy_per_year[i] = entropy(dimension_year[i], 0.00001);
  }
  
  
  for (i in 1:val_N){
    predicted[i] = ordered_logistic_rng(
      country_propensity[val_country[i]] + 
      weighted_inner_product(country_latent[val_country[i]], resolution_latent[val_resolution[i]], dimension_year[val_year[i]]), 
      cutpoints[val_resolution[i]]
      );
  }
  
  for (r in 1:result_type){
    int TP = 0;
    int FP = 0;
    int FN = 0;
    for (i in 1:val_N){
      if (val_result[i] == r && predicted[i] == r){
        TP += 1;
      }
      else if (val_result[i] != r && predicted[i] == r){
        FP += 1;
      }
      else if (val_result[i] == r && predicted[i] != r){
        FN += 1;
      }
    }
    F1_score_by_class[r] = (2 * TP) * 1.0/(2 * TP + FP + FN);
  }
  
  F1_score_macro = mean(F1_score_by_class);
}

```

# モデル推定

今回利用するデータは、こちらからダウンロードできます：

https://dataverse.harvard.edu/dataset.xhtml?persistentId=hdl:1902.1/12379&version=18.0

では、まずデータを読み込んで、マスターテーブル群を作成し、カテゴリデータにIDを付与します：

```R
vote_df <- readr::read_csv("local_data/UNVotes1-72full.csv") |>
  tidyr::drop_na(Countryname, rcid, year.x) |>
  dplyr::filter(!(vote %in% c(8, 9)))

country_master <- vote_df |>
  dplyr::select(Countryname, year.x) |>
  dplyr::distinct() |>
  dplyr::group_by(Countryname) |>
  dplyr::summarise(
    cnt = length(unique(year.x))
  ) |>
  dplyr::ungroup() |>
  dplyr::mutate(
    country_id = dplyr::row_number()
  )

year_master <- vote_df |>
  dplyr::select(year.x) |>
  dplyr::distinct() |>
  dplyr::arrange(year.x) |>
  dplyr::mutate(
    year_id = dplyr::row_number()
  )

resolution_master <- vote_df |>
  dplyr::select(rcid, year.x, descr) |>
  dplyr::distinct() |>
  dplyr::arrange(rcid) |>
  dplyr::mutate(
    resolution_id = dplyr::row_number()
  ) |>
  dplyr::left_join(year_master, by = "year.x")

result_master <- vote_df |>
  dplyr::select(vote) |>
  dplyr::distinct() |>
  dplyr::mutate(
    result_id = dplyr::case_when(
      # yes
      vote == 1 ~ 3,
      # abstain
      vote == 2 ~ 2,
      # no
      vote == 3 ~ 1
    )
  )


vote_df_with_id <- vote_df |>
  dplyr::left_join(country_master, by = "Countryname") |>
  dplyr::left_join(resolution_master |>
                     dplyr::select(rcid, resolution_id), 
                   by = "rcid") |>
  dplyr::left_join(year_master, by = "year.x") |>
  dplyr::left_join(result_master, by = "vote") |>
  tidyr::drop_na(country_id)
```

次に、反対、棄権、賛成から1000ずつの検証データを残して、他のデータを学習データにします：

```R
set.seed(12345)
val_id <- c(
  sample(
    which(vote_df_with_id$result_id == 1),
    1000,
    replace = FALSE
  ),
  sample(
    which(vote_df_with_id$result_id == 2),
    1000,
    replace = FALSE
  ),
  sample(
    which(vote_df_with_id$result_id == 3),
    1000,
    replace = FALSE
  )
)

vote_df_with_id_train <- vote_df_with_id |>
  dplyr::slice(-val_id)

data_list <- list(
  dimension_type = 20,
  group_type = 20,
  
  country_type = nrow(country_master),
  resolution_type = nrow(resolution_master),
  year_type = nrow(year_master),
  result_type = nrow(result_master),
  
  N = nrow(vote_df_with_id_train),
  country = vote_df_with_id_train$country_id,
  resolution = vote_df_with_id_train$resolution_id,
  year = vote_df_with_id_train$year_id,
  result = vote_df_with_id_train$result_id,
  
  val_N = nrow(vote_df_with_id[val_id,]),
  val_country = vote_df_with_id$country_id[val_id],
  val_resolution = vote_df_with_id$resolution_id[val_id],
  val_year = vote_df_with_id$year_id[val_id],
  val_result = vote_df_with_id$result_id[val_id]
)
```

学習データの件数を確認すると：

```R
> nrow(vote_df_with_id_train)
[1] 765501
```

76万になっていますが、それでもStanの変分推論でしたら短時間で推定が終わります。

続いて、モデルをコンパイルして

```R
m_dd_init <- cmdstanr::cmdstan_model("estimation/dynamic_dimension.stan",
                                     cpp_options = list(
                                       stan_threads = TRUE
                                     )
                                     )

```

実際のモデル推定に入ります：

```R
> m_dd_estimate <- m_dd_init$variational(
     seed = 12345,
     threads = 20,
     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.455395 seconds 
1000 transitions using 10 leapfrog steps per transition would take 4553.95 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) 
Iteration: 250 / 250 [100%]  (Adaptation) 
Success! Found best value [eta = 0.1]. 
Begin stochastic gradient ascent. 
  iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
   100      -562325.028             1.000            1.000 
   200      -406138.542             0.692            1.000 
   300      -348548.580             0.517            0.385 
   400      -314513.223             0.415            0.385 
   500      -290846.847             0.348            0.165 
   600      -273179.200             0.301            0.165 
   700      -259876.769             0.265            0.108 
   800      -249813.672             0.237            0.108 
   900      -242047.868             0.214            0.081 
  1000      -235844.362             0.195            0.081 
  1100      -231102.765             0.179            0.065 
  1200      -226923.276             0.166            0.065 
  1300      -223578.950             0.154            0.051 
  1400      -220699.419             0.144            0.051 
  1500      -218327.116             0.135            0.040 
  1600      -216317.484             0.128            0.040 
  1700      -214508.195             0.121            0.032 
  1800      -213306.466             0.114            0.032 
  1900      -211741.406             0.109            0.026 
  2000      -210510.253             0.103            0.026 
  2100      -209557.862             0.099            0.021 
  2200      -208726.108             0.094            0.021 
  2300      -207700.608             0.091            0.018 
  2400      -207206.007             0.087            0.018 
  2500      -206342.186             0.084            0.015 
  2600      -205656.528             0.080            0.015 
  2700      -204923.950             0.078            0.013 
  2800      -204848.466             0.075            0.013 
  2900      -203976.296             0.072            0.011 
  3000      -203630.814             0.070            0.011 
  3100      -202952.362             0.068            0.009   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  2282.9 seconds.
```

38分で終わりました。

結果を保存します：

```R
m_dd_summary <- m_dd_estimate$summary()
```

# 推定結果

## 予測性能

まず、検証データに対する予測性能を確認します：

```R
m_dd_estimate$draws("F1_score_by_class") |>
  tibble::as_tibble() |>
  `colnames<-`(c("nay", "abstain", "yea")) |>
  dplyr::mutate(
    rid = dplyr::row_number(),
    dplyr::across(
      dplyr::everything(),
      ~ as.numeric(.)
    )
  ) |>
  tidyr::pivot_longer(!rid, names_to = "type") |>
  dplyr::bind_rows(
    m_dd_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_scores.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/93c693db-8218-2d29-dab2-b71ccac9027f.png)

棄権クラスに対する予測性能はあまり良くないですが、賛成クラスと反対クラスのF1スコアは0.75付近を超えており、賛成クラスに関しては0.75を統計的に有意に超えていますので、予測性能から見るといいモデルだと言えると思います。

## 次元

### 全体

続いては、年固有の次元の重要度ベクトルを可視化して、国連総会の決議案が発現（reveal）を許す次元の経時的変化を確認します。

```R
m_dd_summary |> 
  dplyr::filter(stringr::str_detect(variable, "dimension_year\\[")) |>
  dplyr::mutate(
    id = variable |>
      purrr::map(
        \(x){
          stringr::str_split(x, "\\[|\\]|,")[[1]][2:3]
        }
      )
  ) |>
  tidyr::unnest_wider(id, names_sep = "_") |>
  dplyr::mutate(
    dplyr::across(dplyr::starts_with("id"),
                  ~ as.integer(.))
  ) |>
  dplyr::left_join(
    year_master, by = c("id_1" = "year_id")
  ) |>
  dplyr::mutate(
    dimension = as.factor(id_2)
  ) |> 
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = year.x, y = mean, color = dimension), alpha = 0.5) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = year.x, ymin = q5, ymax = q95, fill = dimension), alpha = 0.3) + 
  ggplot2::scale_x_continuous(breaks = seq(min(dimension_year$year.x), max(dimension_year$year.x), by = 2)) + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(
    title = "Change of Dimensions' Importance Across Time",
    x = "Year",
    y = "Importance"
  ) + 
  ggplot2::geom_vline(xintercept = 1989, linetype = "dashed", color = "blue", alpha = 0.6, linewidth = 1.5) + 
  ggplot2::annotate("text", x = 1989, y = Inf, label = "End of Cold War", vjust = 2, color = "blue")
```

![dimension_change.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/3c4c4183-62a5-27b1-0c3a-f267a7c497d8.png)

この図からまず言えるのは、冷戦の前後の年重要度ベクトルの多様性の違いです。冷戦後、基本的には二つの次元が安定して各国の国連総会決議案投票行動を説明できていますが、冷戦前は言葉を選ばずにいうとカオスです。

国際政治理論の攻撃的リアリズムを代表する学者、ジョン・ミアシャイマー（John Joseph Mearsheimer）は、冷戦終結直後の1990年の論文で、冷戦の終結はヨーロッパの不安定化につながると予言しました（Mearsheimer, 1990）。

ミアシャイマーは論文の中で歴史を振り返り、多極体制は二極体制より不安定かしがちと主張しました。米ソ二極体制が冷戦期の安定をもたらし、冷戦終結による国際社会の多極体制への回帰は、不安定なヨーロッパ情勢につながると主張しました。

ネオリアリズムは元々国連のような国際機関の意義を疑問視しており（Waltz, 2000）、そもそも国連での投票行動は見る価値もないと主張するかもしれないが、あえて冷戦終結は世界情勢の不安定化につながるという予言を国連での投票行動に適用すると、ミアシャイマーの想定はモデル推定結果と真逆で、二極体制下の冷戦の方が投票行動を説明するための次元数が多くし、冷戦後に初めて安定しました。

もちろん、そもそも多極体制ではなくアメリカによる一極体制が冷戦後の実態なので、極の数自体の予想が間違っているという説明もありうると思いますが、いずれにせよ、安定性に関する予想はモデルの推定結果と異なる事実に変わりはありません。

では、具体的には次元の多様性は時間と共にどのように変動してきたのか？

これを厳密に測定するため、棒折り過程は確率分布であることを利用して、確率分布の「雑多さ」を数値化するエントロピーを利用します。年固有の次元の重要度ベクトルのエントロピーの事後分布を求めるだけなので概念は単純です。

Stanのコードの中で既にサンプリングしているため、これを取り出して可視化しましょう：

```R
m_dd_summary |> 
  dplyr::filter(stringr::str_detect(variable, "entropy_per_year\\[")) |>
  dplyr::bind_cols(
    year = year_master$year.x
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = year, y = mean)) +
  ggplot2::geom_ribbon(ggplot2::aes(x = year, ymax = q95, ymin = q5), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::labs(
    title = "Change of Entropy of Stick Breaking Process",
    x = "Year",
    y = "Entropy"
  ) + 
  ggplot2::scale_x_continuous(breaks = seq(min(dimension_year$year.x), max(dimension_year$year.x), by = 2)) + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::geom_vline(xintercept = 1989, linetype = "dashed", color = "red", alpha = 0.6, linewidth = 1.5) + 
  ggplot2::annotate("text", x = 1989, y = Inf, label = "End of Cold War", vjust = 2, color = "red")
```

![entropy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/5e65385a-6735-f00f-6a5f-41311e1bb069.png)

確率分布の雑多さを示す指標なので、低ければ低いほど、内容が単純ということになります。図から分かるように、冷戦後、特に1995年以降、エントロピーは常に下降傾向にあり、2000年代以降は国連創設以来最も低い値をほぼ毎年更新しています。

### 個別

上の重要度ベクトルの図を見ればわかるように、推定時に無限の（実際は無限の近似として20）の次元を設定しましたが、実際は一部の次元しかモデルに利用されていません。ここでは、実際にモデルに利用された最初の6つの次元を確認しましょう：

```R
purrr::map(
  1:6,
  \(x){
    m_dd_summary |> 
      dplyr::filter(stringr::str_detect(variable, "dimension_year\\[")) |>
      dplyr::mutate(
        id = variable |>
          purrr::map(
            \(x){
              stringr::str_split(x, "\\[|\\]|,")[[1]][2:3]
            }
          )
      ) |>
      tidyr::unnest_wider(id, names_sep = "_") |>
      dplyr::mutate(
        dplyr::across(dplyr::starts_with("id"),
                      ~ as.integer(.))
      ) |>
      dplyr::left_join(
        year_master, by = c("id_1" = "year_id")
      ) |>
      dplyr::mutate(
        dimension = as.factor(id_2)
      ) |> 
      dplyr::filter(dimension == x) |>
      ggplot2::ggplot() + 
      ggplot2::geom_line(ggplot2::aes(x = year.x, y = mean), color = "blue", alpha = 0.5) + 
      ggplot2::geom_ribbon(ggplot2::aes(x = year.x, ymin = q5, ymax = q95), fill = "blue", alpha = 0.3) + 
      ggplot2::scale_x_continuous(breaks = seq(min(dimension_year$year.x), max(dimension_year$year.x), by = 10)) + 
      ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
      ggplot2::labs(
        title = stringr::str_c("dimension ", x),
        x = "Year",
        y = "Importance"
      ) + 
      ggplot2::geom_vline(xintercept = 1989, linetype = "dashed", color = "red", alpha = 0.6, linewidth = 1.5) + 
      ggplot2::annotate("text", x = 1989, y = Inf, label = "End of Cold War", vjust = 2, color = "red") + 
      ggplot2::ylim(c(0, 1))
  }
) |>
  do.call(what = gridExtra::grid.arrange, args = _)
```

![per_dimension_lines.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/a2971f20-c99d-b7f5-889d-af1a2d91a3f0.png)

続いては、各次元の意味をより明確化するため、各次元の代表的な決議案を確認します。まずは内容のローデータから判断しましょう。決議案を当該次元の大きさで並べるだけです。まず第一次元から確認しましょう：


```R
> m_dd_summary |> 
     dplyr::filter(stringr::str_detect(variable, "resolution_latent")) |>
     dplyr::pull(mean) |>
     matrix(ncol = 20) |>
     as.data.frame() |>
     tibble::as_tibble() |>
     dplyr::bind_cols(resolution_id = resolution_master$resolution_id) |>
     dplyr::left_join(resolution_content_master, by = "resolution_id") |>
     dplyr::arrange(V1) |>
     dplyr::select(descr) |>
     head(10)
# A tibble: 10 × 1
   descr                                                                                                                                                                                  
   <chr>                                                                                                                                                                                  
 1 Situation of human rights and fundamental freedoms in Chile                                                                                                                            
 2 Situation of human rights in the Sudan                                                                                                                                                 
 3 Situation of human rights in the Islamic Republic of Iran : resolution / adopted by the General Assemly                                                                                
 4 Situation of human rights in the Islamic Republic of Iran : resolution / adopted by the General Assembly                                                                               
 5 TO REQUEST CONSIDERATION OF THE IDEA OF ELABORATING A SECOND OPTION PROTOCOL TO THE INTERNATIONAL COVENANT ON CIVIL AND POLITICAL RIGHTS, AIMING AT THE ABOLITION OF THE DEATH PENALTY.
 6 Situation of human rights in Turkmenistan : resolution / adopted by the General Assemly                                                                                                
 7 Situation of human rights in the Islamic Republic of Iran : resolution / adopted by the General Assembly                                                                               
 8 Situation of human rights in the Islamic Republic of Iran : resolution / adopted by the General Assembly                                                                               
 9 RES DECIDING TO ABOLISH PER CAPITA CEILING PRINCIPLE IN FORMULATION & ESTAB OF RATES OF ASSESSMENT, COMMENCING WITH SCALE FOR 1977-79 TRIENNIUM.                                       
10 Elaboration of a 2nd optional protocol to the International Covenant on Civil and Political Rights, aiming at the abolition of the death penalty  
```

かなりわかりやすいですが、人権系です。第一次元が冷戦終結直前から増加傾向にあり、冷戦後は最も重要度の高い次元になっていることから見ると、冷戦後の国連は人権のことは資本主義と共産主義のようなイデオロギーではなく、人権をより多く語るようになったと言えます。ただ、気に入らない国家は人権を守っていないと批判することも多いので、人権をより多く語るようになったことは人権を本気で重視するようになったとは限らないと強調したいです。

第二次元にどんなものが入っているかを確認します：

```R
m_dd_summary |> 
     dplyr::filter(stringr::str_detect(variable, "resolution_latent")) |>
     dplyr::pull(mean) |>
     matrix(ncol = 20) |>
     as.data.frame() |>
     tibble::as_tibble() |>
     dplyr::bind_cols(resolution_id = resolution_master$resolution_id) |>
     dplyr::left_join(resolution_content_master, by = "resolution_id") |>
     dplyr::arrange(V2) |>
     dplyr::select(descr) |>
     head(10)
# A tibble: 10 × 1
   descr                                                                                                                                                                                   
   <chr>                                                                                                                                                                                   
 1 TO ADOPT 4TH COMM. DRAFT RESOLUTION II (A/4672) REQUESTING BELGIUM TO REVOKE ITS SUSPENSION OF THE POWERS OF THE MWAMI OF RUANDA AND DECIDING TO CONDUCT A U.N.-SUPERVISED REFERENDUM T…
 2 TO APPROVE A RESOLUTION STRONGLY CONDEMNING THE COLLUSION BY THE GOVERNMENTS OF CERTAIN WESTERN COUNTRIES, AND IN PARTICULAR THOSE OF THE UNITED STATES, THE FEDERAL REPUBLIC OF GERMAN…
 3 TO DENOUNCE THE USE BY DEVELOPED COUNTRIES OF ECONOMIC COERCION SUCH AS BLOCKADES AND EMBARGOES AGAINST DEVELOPING NATIONS.                                                             
 4 TO RECOMMEND THAT THE MEMBER STATES CONCERNED AND THE COMPETENT INTERNATIONAL ORGANIZATIONS SHOULD, AS A MATTER OF URGENCY, GIVE DUE CONSIDERATION TO THE FORMULATION OF POLICIES WITH …
 5 TO CALL FOR THE SECRETARY-GENERAL TO CREATE AN INTER-AGENCY GROUP OF REPRESENTATIVES FROM APPROPRIATE ORGANIZATIONS WITHIN THE UNITED NATIONS TO STUDY THE ISSUE OF REVERSE TRANSFER OF…
 6 TO APPROVE A RESOLUTION REQUESTING THE COMMITTEE ON CONTRIBUTIONS TO SUBMIT TO THE GENERAL ASSEMBLY A THOROUGH STUDY ON ALTERNATIVE METHODS TO ASSESS THE REAL CAPACITY TO PAY OF MEMBE…
 7 strongly condemns the continuing and increasing collaboration of Israel with the racist regime of South Africa, especially in the economic, military and nuclear fields;                
 8 TO ADOPT WHOLE 4TH COMMITTEE DRAFT RESOLUTION (A/6089) URGING U.K. TO CEASE MILITARY REPRESSIONS IN ADEN.                                                                               
 9 TO APPROVE A RESOLUTION URGING THE DEVELOPED COUNTRIES TO EXTEND FULL COOPERATION TO ENSURE THE SUCCESSFUL CONCLUSION OF THE CURRENT MULTILATERAL TRADE NEGOTIATIONS, GIVING FULL REGAR…
10 TO REQUEST A CONFERENCE TO CONVENE ON THE REVERSE TRANSFER OF TECHNOLOGY.  
```

これはおそらく第三世界関連次元ですね。

第三次元は：

```R
m_dd_summary |> 
     dplyr::filter(stringr::str_detect(variable, "resolution_latent")) |>
     dplyr::pull(mean) |>
     matrix(ncol = 20) |>
     as.data.frame() |>
     tibble::as_tibble() |>
     dplyr::bind_cols(resolution_id = resolution_master$resolution_id) |>
     dplyr::left_join(resolution_content_master, by = "resolution_id") |>
     dplyr::arrange(V3) |>
     dplyr::select(descr) |>
     head(10)
# A tibble: 10 × 1
   descr                                                                                                                                                                                   
   <chr>                                                                                                                                                                                   
 1 "TO AMEND TO RES (A/L.650/REV.1) CALLING ON THE SYG TO REACTIVATE THE MISSION OF THE SPECIAL REPRESENTATIVE OF THE SYG TO THE MIDDLE EAST, PURSUANT TO S/RES/242(1967)."                
 2 "Nuclear weapons prohibition convention (draft text)"                                                                                                                                   
 3 "TO AMEND TO OPERATIVE PARA 1 OF RES (A/L.650/REV.1) DELETING WORDS \\AND THAT, CONSEQUENTLY, TERRITORIES THUS OCCUPIED MUST BE RESTORED.\\\\\"\t1\t0\t0\t0\t1\t0"                      
 4 "TO AMEND TO RES (A/L.650/REV.1) REPLACING OP PARA 4 & 5 NITH A NEW OP PARA."                                                                                                           
 5 "Reducing nuclear danger"                                                                                                                                                               
 6 "Convention on the Prohibition of the Use of Nuclear Weapons"                                                                                                                           
 7 "Convention on the Prohibition of the Use of Nuclear Weapons"                                                                                                                           
 8 "Convention on the Prohibition of the Use of Nuclear Weapons"                                                                                                                           
 9 "Nuclear disarmament"                                                                                                                                                                   
10 "TO AMEND TO RES (A/L.650/REV.1) CALLING ON THE PARTIES TO RESUME IMMEDIATELY THE CONVERSATIONS UNDER AUSPICES OF THE SPECIAL REPRESENTATIVE WITH A VIEW TO CONCLUDING A PEACE AGREEMEN…
```

核兵器廃絶次元ですね。冷戦前からゆっくりと上昇し、冷戦後は人権次元に続いて二番目に重要な次元になっています。

第四次元は：

```R
m_dd_summary |> 
     dplyr::filter(stringr::str_detect(variable, "resolution_latent")) |>
     dplyr::pull(mean) |>
     matrix(ncol = 20) |>
     as.data.frame() |>
     tibble::as_tibble() |>
     dplyr::bind_cols(resolution_id = resolution_master$resolution_id) |>
     dplyr::left_join(resolution_content_master, by = "resolution_id") |>
     dplyr::arrange(V4) |>
     dplyr::select(descr) |>
     head(10)
# A tibble: 10 × 1
   descr                                                                                                                                                                        
   <chr>                                                                                                                                                                        
 1 TO ADOPT PARAGRAPH 2 OF A GENERAL COMM. RECOMMENDATION (A/4214), DECIDING NOT TO CONSIDER ANY FURTHER PROPOSALS ON CHINA UN REPRESENTATION DURING CURRENT SESSION.           
 2 TO ADOPT PARAGRAPH 1 OF A GENERAL COMM. RECOMMENDATION (A/4214), REJECTING INDIAN REQUEST FOR AN ITEM ON CHINA REPRESENTATION IN THE UN, IN THE AGENDA.                      
 3 TO ADOPT WHOLE GENERAL COMMITTEE RECOMMENDATION ON                                                                                                                           
 4 TO ADOPT A GENERAL COMM. DRAFT RESOL. (A/3670) ON CHINESE REPRESENTATION.                                                                                                    
 5 TO ADOPT PARAGRAPH 1 OF THE GENERAL COMM. RECOMMENDATION (A/3926), SAID PARAGRAPH REJECTING INDIAN REQUEST.                                                                  
 6 TO ADOPT PARAGRAPH 1 OF THE GENERAL COMM. DRAFT RESOL. (A/3670), SAID PARAGRAPH REJECTING INDIA-PROPOSED AGENDA ITEM ON REPRESENTATION OF CHINA IN THE UN.                   
 7 TO ADOPT PARAGRAPH 1 OF THE GENERAL COMM. RECOMMENDATION (A/4520) SAID PARAGRAPH REJECTING U.S.S.R. REQUEST FOR AN ITEM ON REPRESENTATION OF CHINA IN THE UN, IN THE AGENDA. 
 8 TO ADOPT PARAGRAPH 2 OF THE GENERAL COMM. DRAFT RESOL. (A/3670), DECIDING NOT TO CONSIDER ANY PROPOSALS TO UNSEAT NATIONALIST CHINA UN REPRESENTATION DURING CURRENT SESSION.
 9 TO ADOPT CREDENTIAL COMMITTEE DRAFT RESOLUTION (A/4574) ACCEPTING THE CREDENTIALS OF THE REPUBLIC OF CONGO (LEOPOLDVILLE) REPRESENTATIVES.                                   
10 TO ADOPT GENERAL COMMITTEE RECOMMENDATION (A/3926) ON CHINA REPRESENTATION.   
```

あとで紹介する第六次元とかなり内容が被り、両方とも中国代表権問題関連の決議案が上位に出てきていますが、重要度の高さを維持する期間から判断すると、第四次元は冷戦次元と解釈した方が妥当そうです。

第五次元は：

```R
m_dd_summary |> 
     dplyr::filter(stringr::str_detect(variable, "resolution_latent")) |>
     dplyr::pull(mean) |>
     matrix(ncol = 20) |>
     as.data.frame() |>
     tibble::as_tibble() |>
     dplyr::bind_cols(resolution_id = resolution_master$resolution_id) |>
     dplyr::left_join(resolution_content_master, by = "resolution_id") |>
     dplyr::arrange(V5) |>
     dplyr::select(descr) |>
     head(10)
# A tibble: 10 × 1
   descr                                                                                                                                                                                   
   <chr>                                                                                                                                                                                   
 1 PROPOSAL THAT RIGHT TO SPEAK IN DEBATE ON SUBJECT SHLD BE LIMITED TO SINGLE INTERVENTION, WITHOUT PREJUDICE TO PROVISIONS OF RULES OF PROCEDURE RE EXERCISE OF RIGHTS OF REPLY.         
 2 TO COMMEND THE SPECIAL COMMITTEE TO INVESTIGATE ISRAELI PRACTICES AFFECTING THE HUMAN RIGHTS OF THE POPULATION OF THE OCCUPIED TERRITORIES FOR ITS EFFORTS IN PERFORMING THE TASKS ASSI…
 3 TO REQUEST THE SECURITY COUNCIL TO CONSIDER AND ADOPT THE NECESARY RESOLUTIONS AND MEASURES IN ORDER TO ENABLE THE PALESTINIAN PEOPLE TO EXERCISE ITS INALIENABLE NATIONAL RIGHTS IN AC…
 4 TO AFFIRM THAT THE EARLY RESUMPTION OF THE PEACE CONFERENCE ON THE MIDDLE EAST WITH THE PARTICIPATION OF ALL PARTIES CONCERNED, INCLUDING THE PALESTINIAN LIBERATION ORGANIZATION, IN A…
 5 TO CONDEMN THE MASSIVE, DELIBERATE DESTRUCTION OF QUNEITRA PERPETRATED DURING ISRAELI OCCUPATION AND PRIOR TO THE WITHDRAWL OF ISRAELI FORCES FROM THAT CITY IN 1974.                   
 6 RES, I.A., COMMENDING SPECIAL COMTE FOR ITS WORK; CALLING UPON ISRAEL TO DESIST IMMEDIATELY FROM CHANGING PHYSICAL CHARACTER, DEMOGRAPHIC COMPOSITION, ETC., OF SUCH TERRITORIES DECLAR…
 7 TO CONDEMN THE CONTINUING AND INCREASING COLLABORATION BY ISRAEL WITH THE RACIST REGIME OF SOUTH AFRICA AS A FLAGRANT VIOLATION OF THE RESOLUTIONS OF THE UNITED NATIONS AND AS AN ENCO…
 8 TO CONDEMN THE INTENSIFIED ACTIVITIES OF THOSE FOREIGN ECONOMIC, FINANCIAL AND OTHER INTERESTS WHICH CONTINUE TO EXPLOIT THE NATURAL AND HUMAN RESOURCES OF THE COLONIAL TERRITORIES AN…
 9 RES, I.A., AFFIRMING RIGHT OF ARAB STATES & PEOPLES UNDER FOREIGN OCCUPATION TO PERMANENT SOVEREIGNTY; REAFFIRMING ISRAELI MEASURES TO EXPLOIT SUCH RESOURCES TO BE ILLEGAL; & AFFIRMIN…
10 TO PROCLAIM THAT THE RACIST REGIME OF SOUTH AFRICA IS ILLEGITIMATE AND HAS NO RIGHT TO REPRESENT THE PEOPLE OF SOUTH AFRICA AND TO REAFFIRM THAT THE NATIONAL LIBERATION MOVEMENTS RECO…
```

今後の更なる調査・分析の中心になりますが、内容から判断すると、おそらく南アフリカか第四次中東戦争関連の特殊情勢かと思います。なぜこの時期に各国の投票行動が「狂い」、冷戦次元と第三世界次元で説明できなくなったのかを深掘りしていきたいと思います。

最後の第六次元が国連設立当初の会員資格周りを表していると思われます：
```R
m_dd_summary |> 
     dplyr::filter(stringr::str_detect(variable, "resolution_latent")) |>
     dplyr::pull(mean) |>
     matrix(ncol = 20) |>
     as.data.frame() |>
     tibble::as_tibble() |>
     dplyr::bind_cols(resolution_id = resolution_master$resolution_id) |>
     dplyr::left_join(resolution_content_master, by = "resolution_id") |>
     dplyr::arrange(V6) |>
     dplyr::select(descr) |>
     head(10)
# A tibble: 10 × 1
   descr                                                                                                                                                                        
   <chr>                                                                                                                                                                        
 1 TO ADOPT PARAGRAPH 2 OF A GENERAL COMM. RECOMMENDATION (A/4214), DECIDING NOT TO CONSIDER ANY FURTHER PROPOSALS ON CHINA UN REPRESENTATION DURING CURRENT SESSION.           
 2 TO ADOPT PARAGRAPH 1 OF A GENERAL COMM. RECOMMENDATION (A/4214), REJECTING INDIAN REQUEST FOR AN ITEM ON CHINA REPRESENTATION IN THE UN, IN THE AGENDA.                      
 3 TO ADOPT WHOLE GENERAL COMMITTEE RECOMMENDATION ON                                                                                                                           
 4 TO ADOPT A GENERAL COMM. DRAFT RESOL. (A/3670) ON CHINESE REPRESENTATION.                                                                                                    
 5 TO ADOPT PARAGRAPH 1 OF THE GENERAL COMM. RECOMMENDATION (A/3926), SAID PARAGRAPH REJECTING INDIAN REQUEST.                                                                  
 6 TO ADOPT PARAGRAPH 1 OF THE GENERAL COMM. DRAFT RESOL. (A/3670), SAID PARAGRAPH REJECTING INDIA-PROPOSED AGENDA ITEM ON REPRESENTATION OF CHINA IN THE UN.                   
 7 TO ADOPT PARAGRAPH 1 OF THE GENERAL COMM. RECOMMENDATION (A/4520) SAID PARAGRAPH REJECTING U.S.S.R. REQUEST FOR AN ITEM ON REPRESENTATION OF CHINA IN THE UN, IN THE AGENDA. 
 8 TO ADOPT PARAGRAPH 2 OF THE GENERAL COMM. DRAFT RESOL. (A/3670), DECIDING NOT TO CONSIDER ANY PROPOSALS TO UNSEAT NATIONALIST CHINA UN REPRESENTATION DURING CURRENT SESSION.
 9 TO ADOPT CREDENTIAL COMMITTEE DRAFT RESOLUTION (A/4574) ACCEPTING THE CREDENTIALS OF THE REPUBLIC OF CONGO (LEOPOLDVILLE) REPRESENTATIVES.                                   
10 TO ADOPT GENERAL COMMITTEE RECOMMENDATION (A/3926) ON CHINA REPRESENTATION.   
```

以上がローデータに基づく説明でした。ローデータだけだと、ではそれぞれの次元に特殊なワードはなんなのか、といったような研究者が分析結果を解釈する上で必要な情報を可視化できません。よって、ここではTaddy（2015）が提案した多項分散回帰で、各次元の代表的なワードを抽出します。要するに、決議案のデータの単語の出現回数を従属変数にして、年固有の次元の重要度ベクトルで重み付けした決議案ベクトルを独立変数にして、係数の絶対値が強い単語を確認します。

モデルはこのように推定します：

```R
resolution_dfm <- resolution_master |>
  dplyr::pull(descr) |>
  stringr::str_remove_all("[^a-zA-Z ]") |>
  quanteda::tokens() |>
  quanteda::tokens_remove(quanteda::stopwords()) |>
  quanteda::tokens_ngrams(1:2) |>
  quanteda::tokens_tolower() |>
  quanteda::dfm() |>
  quanteda::dfm_trim(min_termfreq = 5)

resolution_latent_weighted_yearly_matrix <- m_dd_summary |> 
  dplyr::filter(stringr::str_detect(variable, "resolution_latent_weighted_yearly\\[")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 20)

cl <- parallel::makeCluster(16)

m_dmr_resolution <- distrom::dmr(
  covars = resolution_latent_weighted_yearly_matrix,
  counts = resolution_dfm,
  cl = cl,
  standardize = FALSE,
  verb = TRUE
)

parallel::stopCluster(cl)
```

係数の絶対値が大きい単語は下記の通りです：

```R
> m_dmr_resolution |>
     coef() |>
     as.matrix() |>
     as.data.frame() |>
     tibble::rownames_to_column(var = "dimension") |>
     tibble::tibble() |>
     tidyr::pivot_longer(!dimension, names_to = "vocab") |>
     dplyr::filter(dimension != "intercept" & dimension %in% c(1:6)) |>
     split(~ dimension) |>
     purrr::map(
       \(df){
             df |>
                 dplyr::mutate(
                     value = abs(value)
                 ) |>
                 dplyr::arrange(-value) |>
                 dplyr::select(vocab) |>
                 dplyr::slice_head(n = 20)
         }
     ) |>
     dplyr::bind_cols() |>
     `colnames<-`(stringr::str_c("dimension_", 1:6)) |>
     print(n = 20)
New names:
• `vocab` -> `vocab...1`
• `vocab` -> `vocab...2`
• `vocab` -> `vocab...3`
• `vocab` -> `vocab...4`
• `vocab` -> `vocab...5`
• `vocab` -> `vocab...6`
# A tibble: 20 × 6
   dimension_1             dimension_2              dimension_3                dimension_4              dimension_5              dimension_6          
   <chr>                   <chr>                    <chr>                      <chr>                    <chr>                    <chr>                
 1 rights_sudan            israel_south             res_alrev                  development_organization res_alrev                nepal_amendment      
 2 rights_cuba             res_alrev                nationals                  appropriations_year      exercise_palestinian     council_admission    
 3 democratic_republic     rights_sudan             disarmament_general        moratorium               condemn_israels          palestine_question   
 4 moratorium              vote_operative           reducing_nuclear           budget_appropriations    moratorium               al_general           
 5 moratorium_use          possible_offer           nuclear_danger             appropriations_biennium  continuing_increasing    adopt_trusteeship    
 6 use_death               offer_necessary          legality                   biennium_amount          increasing_collaboration approved_th          
 7 penalty_resolution      resources_meet           legality_threat            appropriations           korea_continue           reject               
 8 rights_democratic       meet_needs               danger_resolution          nepal_amendment          israel_racist            consider_proposals   
 9 armaments_resolution    near_ea                  arrangements_strengthening expel                    crime_humanity           section_ad           
10 rights_iraq             us_motion                reducing                   un_expel                 rights_population        weapons_testing      
11 transparency_armaments  comprehensive_just       jerusalem_resolution       financing_appropriations collaboration_israel     nepal                
12 penalty                 exercise_palestinian     security_nonnuclearweapon  expel_forthwith          exploit                  add_corr             
13 death_penalty           conflict_middle          rights_islamic             biennium                 al_corr                  corr_corr            
14 promotion_democratic    continuing_increasing    us_motion                  institute_training       regime_south             corr_rev             
15 democratic_equitable    increasing_collaboration taking_forward             training_research        resources_territories    committee_report     
16 international_order     collaboration_israel     forward_multilateral       proposal_change          vii_charter              hungary_agenda       
17 equitable_international israel_racist            multilateral_nuclear       change_representation    offer                    league_nations       
18 extrajudicial           rights_cuba              palestine_department       korea_continue           pream_para               republic_congo       
19 extrajudicial_summary   core                     disarmament_negotiations   deciding_recognize       assigned                 recommending_security
20 executions              reaffirms_conviction     disarmament_prevention     human_settlements        performing               permanent_members    
```

これ決議案のローデータを見て判断した次元の内容の妥当性がよりわかりやすくなったのではないかと思います。

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

最後は、推定されたイデオロギーベクトルで各国の投票行動を可視化します。具体的には、t-SNE（Van der Maaten and Hinton, 2008）で各国のイデオロギーベクトルの事後分布の平均とイデオロギーベクトルの全体の事前分布で散布図を作ります：

```R
set.seed(12345)
country_latent_and_G_tsne <- m_dd_summary |> 
  dplyr::filter(stringr::str_detect(variable, "country_latent\\[")) |>
  dplyr::pull(mean) |>
  matrix(ncol = 20) |>
  rbind(
    m_dd_estimate$draws("drawn_G") |> 
      as.data.frame() |> 
      as.matrix()
  ) |>
  Rtsne::Rtsne()

tsne_df <- tibble::tibble(
  type = c(rep("country", nrow(country_master)), rep("sample", 1000)),
  name = c(country_master$Countryname, rep("sample", 1000)),
  x1 = country_latent_and_G_tsne$Y[,1],
  x2 = country_latent_and_G_tsne$Y[,2]
)

g_tsne <- tsne_df |>
  dplyr::filter(type == "country") |>
  ggplot2::ggplot() +
  ggplot2::geom_text(ggplot2::aes(x = x1, y = x2, label = name), color = ggplot2::alpha("blue", 0.5)) + 
  ggplot2::geom_density_2d(data = tsne_df |>
                             dplyr::filter(type == "sample"),
                           ggplot2::aes(x = x1, y = x2),
                           color = ggplot2::alpha("black", 0.2),
                           show.legend = FALSE
  ) +
  ggplot2::xlim(c(-22, 22)) + 
  ggplot2::ylim(c(-22, 22)) + 
  ggplot2::labs(title = "tsne of the country latent vectors and G")

plotly::ggplotly(g_tsne)
```

![Screen Shot 2024-09-22 at 21.12.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/61e98cb2-6f2c-ce31-42f3-e49d3d712717.png)

等高線がG分布で、青い文字が各国のイデオロギーベクトルの位置です。少しわかりにくいですが、まず最初に注目していただきたいのは、上のG分布の範囲からはみ出した国です。拡大しますね：

![Screen Shot 2024-09-22 at 21.15.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/78d0c9a6-86e6-1e77-9f13-04b509ad0e4c.png)

共産圏のコアの国ですね。これは何を意味しているのかというと、共産圏のコアの国のイデオロギーベクトルは、他の国と違いすぎて、G分布の観点からいうと「統計的に有意に」孤立していると言えます。

次に、下の方を拡大しましょう：

![Screen Shot 2024-09-22 at 21.17.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/5daf7b03-6f8e-26e5-1152-f789ae70bd4d.png)

いわゆる西側諸国です。

ただし、ここでも注目していただきたいのは、アメリカ、イスラエル、オーストラリア、カナダと他の西側諸国の距離です。これは今後厳密に測定する方法を考案したいですが、第二次中東戦争（スエズ戦争）に代表されるように、西側諸国にもそれぞれの思惑があり、外交政策においても必ずしも一枚岩ではないことを示していると思います。

最後に、地図で各国の次元の分布を可視化します。ただし、比較軸をわかりやすくするため、アメリカの値が正になるように調整します。具体的には、アメリカの値が負になっている次元に対しては、全体に-1をかけます。

```R
1:6 |>
  purrr::map(
    \(i){
      target_df <- m_dd_summary |> 
        dplyr::filter(stringr::str_detect(variable, "country_latent\\[")) |> 
        dplyr::mutate(
          id = variable |>
            purrr::map(
              \(x){
                stringr::str_split(x, "\\[|\\]|,")[[1]][2:3]
              }
            )
        ) |>
        tidyr::unnest_wider(id, names_sep = "_") |>
        dplyr::mutate(
          dplyr::across(dplyr::starts_with("id"),
                        ~ as.integer(.))
        ) |>
        dplyr::filter(id_2 == i) |>
        dplyr::left_join(
          country_master, by = c("id_1" = "country_id")
        ) |>
        dplyr::mutate(
          Countryname = dplyr::case_when(
            Countryname == "United States of America" ~ "USA",
            Countryname == "Russian Federation" ~ "Russia",
            Countryname == "Taiwan, Province of China" ~ "Taiwan",
            Countryname == "United Kingdom of Great Britain and Northern Ireland" ~ "UK",
            Countryname == "Democratic People's Republic of Korea" ~ "North Korea",
            Countryname == "Democratic Republic of the Congo" ~ "Democratic Republic of Congo",
            Countryname == "Côte D'Ivoire" ~ "Ivory Coast",
            Countryname == "Lao People's Democratic Republic" ~ "Laos",
            Countryname == "Syrian Arab Republic" ~ "Syria",
            Countryname == "Viet Nam" ~ "Vietnam",
            Countryname == "Federal Republic of Germany" ~ "Germany",
            Countryname == "German Democratic Republic" ~ "Germany",
            Countryname == "Iran (Islamic Republic of)" ~ "Iran",
            Countryname == "Gambia (Islamic Republic of the)" ~ "Gambia",
            Countryname == "Bolivia (Plurinational State of)" ~ "Bolivia",
            Countryname == "Venezuela (Bolivarian Republic of)" ~ "Venezuela",
            Countryname == "Micronesia (Federated States of)" ~ "Micronesia",
            Countryname == "Republic of Moldova" ~ "Moldova",
            Countryname == "Cabo Verde" ~ "Cape Verde",
            Countryname == "Myanmar" ~ "Burma",
            Countryname == "Republic of Korea" ~ "South Korea",
            Countryname == "Palestine, State of" ~ "Palestine",
            Countryname == "United Republic of Tanzania" ~ "Tanzania",
            TRUE ~ Countryname
          )
        )
      if (target_df$mean[which(target_df$Countryname == "USA")] < 0){
        target_df$mean <- -1 * target_df$mean
      }
      
     g <- ggplot2::map_data("world") |>
       dplyr::left_join(
          target_df, by = c("region" = "Countryname")
        ) |>
       ggplot2::ggplot() + 
       ggplot2::geom_polygon(ggplot2::aes(x = long, y = lat, group = group, fill = mean)) + 
       ggplot2::scale_fill_gradient2(
         name = "",
         low = "red",
         mid = "white",
         high = "blue",
         midpoint = 0,
         labels = NULL
       ) + 
       ggplot2::xlab("") +
       ggplot2::ylab("") + 
       ggplot2::ggtitle(stringr::str_c("dimension ", i)) + 
       ggplot2::guides(fill = "none") + 
       ggplot2::theme(
         axis.text.x = ggplot2::element_blank(),
         axis.text.y = ggplot2::element_blank()
       )
     
     return(g)
    }
  ) |>
  do.call(what = gridExtra::grid.arrange, args = _)
```

![country_latent_map.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2257087/1a67b5d6-ce52-8c7b-3156-06cb26bc8cd2.png)

青がアメリカと同じ方向、赤がアメリカと逆方向と理解していただければ読みやすいのかと思います。

この地図でまず着目していただきたいのは、第三次元の核兵器廃絶次元です。西欧諸国がアメリカと逆の色になっているのはこの次元だけです。Mearsheimer（1990, 1994-1995）が予言したような戦争や軍事的不安定にならなくても、リアリズムの説明を借りると、ソ連による軍事的な脅威が消えた冷戦後では、西欧諸国はもはやアメリカの（核）軍事力を必要としないため、アメリカとの政策選好が距離遠くなると解釈できます。

次は第四次元の冷戦次元です。南米の色はアメリカより濃くなっているのは深掘りする価値がありそうです。先行研究を参照していない状態で仮説を述べさせていただくと、おそらく自国の力でソ連の影響に対抗できない南米諸国が国連という国際機関に頼らざるを得ないのかと思われます。また、旧ソ連以外で負（赤）になっている国がないのも興味深いです。

# 結論

いかがでしたか？まだ分析の途中なので、全ての内容を説明できたというわけではありません。StanでサンプリングしたQOIの事後分布も一部しか説明できていません。

今後は、表現イデオロギーベクトルのコサイン類似度の経時的変化や各国が属するクラスターなどの未解説のQOIだけでなく、

- 各次元の意味の深掘り
- 国際政治学の各理論がどのように冷戦の終結を説明したかのまとめと本モデルとの比較
- 国際政治史、特に冷戦史を踏まえた年重要度ベクトルの変化の説明

などの内容にも取り組んでいきたいと思います！

最後に、私たちと働いて、データサイエンスの力で社会を改善したい方はこちらをご確認ください：

https://www.wantedly.com/projects/1021866

# 参考文献

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

Bailey, Michael A. "Comparable preference estimates across time and institutions for the court, congress, and presidency." American Journal of Political Science 51.3 (2007): 433-448.

Bailey, Michael A., Anton Strezhnev, and Erik Voeten. "Estimating dynamic state preferences from United Nations voting data." Journal of Conflict Resolution 61.2 (2017): 430-456.

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.

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.

Mearsheimer, John J. “Back to the Future: Instability in Europe after the Cold War.” International Security, vol. 15, no. 1, 1990, pp. 5–56.

Mearsheimer, John J. “The False Promise of International Institutions.” International Security, vol. 19, no. 3, 1994-1995, pp. 5–49.

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.

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.

Taddy, Matt. "Distributed multinomial regression." The Annals of Applied Statistics 9.3 (2015): 1394-1414.

Teh, Yee Whye, Michael I. Jordan, Matthew J. Beal, and David M. Blei. "Hierarchical Dirichlet processes." Journal of the American Statistical Association 101.476 (2006): 1566-1581.

Van der Maaten, Laurens, and Geoffrey Hinton. "Visualizing data using t-SNE." Journal of machine learning research 9.11 (2008).

Waltz, Kenneth N. "Structural realism after the Cold War." International security 25.1 (2000): 5-41.