はじめに

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

本記事では、異質処置効果モデルを提案してみます。

共変量によって処置効果がどう違うかを推定できるだけでなく、推定のプロセスの可視化までできる手法です。

ぜひ最後まで読んでください！

因果推論の問題意識と既存手法の限界

皆さん、ビジネスやアカデミアで以下のようなことを聞かれた経験はありますか？

営業部のマネージャ：新しい営業施策の処置効果はどこの部署のKPIで計算されたの？

学会の参加者：内戦の経済効果の推定にアメリカのデータも入れてるの？

全然違う観測値を比較してはいけない、というのがこの二つの指摘に共通している前提です。因果推論でいうと共変量をきちんと揃えないと信頼性が問われるということです。

共変量のバランスを揃えるために様々マッチング系の手法が開発されました。一番有名な傾向スコアマッチング（Rosenbaum and Rubin 1983）の他に、計量政治学の世界で開発されたCEM（Coarsened Exact Matching、粗正確マッチング）（Iacus, King and Porro 2012）などがあります。

本記事で紹介する手法は、共変量のマッチング問題をクラスタリング問題に置き換えて、ディリクレ過程というベイズ機械学習で「似た観測値同士」をより精緻に洗い出します。

また、営業施策の場合、年次や等級・役割によって、内戦の効果の場合、経済発展や民主主義の発展などによって、そもそも効果に違いが出る可能性が十分にあります。因果推論でいうと共変量によって異質処置効果が存在しているかもしれないということになります。

異質処置効果の推定が最近流行っている因果機械学習の出番なんですが、解釈が難しいアルゴリズムが多い印象です。

例えば有名な因果フォレスト（Athey, Tibshirani and Wager 2019）は共変量の重要度でモデルの意思決定を可視化できるとされるが、重要度は抽象的にいうと、あくまでも決定木（的なもの）を構築する際に変数がどれほど使われたかを示す指標にすぎないため、「営業年次が低い方が効果が出る」、「民族的多様性が高いの国の方が内戦のダメージが大きい」、などといった意思決定に直結する重要な知見の抽出で利用できません。

本記事で提案する手法は、因果機械学習のように効果の異質性を学習できて、かつ「営業年次が低い方が効果が出る」、「民族的多様性が高いの国の方が内戦のダメージが大きい」、などの知見が比較的簡単に抽出できます。

考え方の概要

本記事で提案する手法は、厳密には処置効果を推定するためのモデルではなく、処置効果推定のための前処理をやってくれるモデルです。

クラスタリングという前処理を施せば、同じクラスター内の観測値同士の共変量は類似しているはずなので、推定結果のモデル依存性（「共変量の二乗項を入れたら処置効果がプラスからマイナスになった！」、など）が軽減されます。これはHo, Imai, King and Stuart（2007）の考えに近いですが、本記事はまず初歩的な概念にフォーカスしたいため、処置効果のモデリングに別の統計モデルを設定するのではなく、とりあえず処置グループと統制グループの平均の差で処置効果を計算します。

前処理としてクラスタリングを行った後、各クラスターごとに平均処置効果を計算しますが、例えばクラスター1内の効果を計算するときは他のクラスターの観測値を全部捨てる、という極端でワイルドなやり方ではなく、ベイズ機械学習モデルが推定したクラスター1に属する事後確率で観測値に重みをつけて平均を計算します。

なぜクラスター1以外の観測値を捨てないのかというと、下記のTVCMの効果測定の事例を一緒に考えてみましょう。

東京都でTVCMを放映する際の効果を測定する際に、大阪府の情報を完全に捨てるのは果たして適切なのでしょうか？東京も大阪も世界的な大都市て、文化と産業構成にある程度差があるものの、適切なウエイトをつければ、東京都の効果を推定する時の参考になるはずです。一方で、大分県になると、CMのドメインによる部分もあるものの、東京都のTVCMの効果の計算の参考にはおよそならないといっても過言ではないかと思われます。

大阪を参考にすることで、実質的に東京が2つ（もしくは1.XXX個）あることになり、サンプル数が増えますので、効果推定の精度も上がるはずです。

分析を行う際は、このように類似している観測値の情報を積極的に活用すべきであり、「他のところの情報をうまく転用する」ことが統計学と機械学習の強さの本当の秘密です。

私は現在会社でデータサイエンティストチームのリーダーを務めていますが、データを切って独立したモデルを大量に作るのは特別な理由がない限りしないでくださいと、いつも部署内外の分析者に言っています。

アルゴリズムの全体フロー

次に、アルゴリズムの全体フローを説明します。

ディリクレ過程モデルでクラスタリング

まずは、ディリクレ過程という柔軟なクラスターモデルで学習データの観測値の共変量をクラスタリングします。ただし、処置フラグ、処置後変数、および結果変数は入れません。

処置フラグを入れないのはただ単に入れていないだけなので、本質的ではないです。今後は別記事で処置フラグを入れてモデルを拡張する方法を紹介します。

処置後変数を入れないのは、入れることで逆に効果の推定にバイアスが発生するからです。Montgomery, Nyhan and Torres（2018）が説明したように、処置後変数のバランスを取ろうとすると、逆に類似していない観測値同士を比較するという逆説的な状況が発生する。

1年後の効果を測定することが目的の治験を考えましょう。薬を投与した一週間後の症状が改善された処置群と統制群を考えてほしいですが、そもそも薬飲んでないのに勝手に病気が治った統制群って、普通の統制群とは違いますよね？元々免疫が強い可能性があり、さらにいうとそもそも病気にかかっておらず、偽陽性だったという可能性もあります。この「変な統制群」と「普通に薬飲んで症状が改善した処置群」と比較するのは極めて不適切で危険です。

最後に、結果変数を入れないのは、モデルのバイアス削減のためです。

Wager and Athey（2018）は決定木を構築するデータと実際に処置効果を推定する目的変数を含むデータを分けることで、バイアスが減ると説明した。

テキストデータを使った因果推論モデルを提案したFong and Grimmer（2016）も、テキストからトピックを推定するためのデータと、実際に処置効果推定に利用するデータを分けることで、因果推論の結論の信憑性がより高くなると説明した。

もちろん、いわゆる処置後の経過状況を分析したいニーズもありますが、その際はZhou and Yamamoto（2023）などの専用のモデル定式化を利用しましょう。

学習データの観測値に所属クラスターの事後確率を計算

クラスタリングモデルを構築できたら、各観測値がそれぞれのクラスターに属する事後確率を求めることができます。やり方の概念はシンプルで、具体的な実装方法はこの記事を参照してください：

理論に関しては、Imai and Tingley（2012）の内容を確認してください。

クラスターごとの処置効果を計算する

最後のステップもシンプルで、各クラスターごとに普通の処置群と統制群の結果変数の重み付けの平均の差を計算すればいいです。重みは各観測値が対象のクラスターに属する事後確率です。

シミュレーションデータ

さて、ここではいよいよ実装スタートです！

モデルの性能を確認するため、シミュレーションデータを作成します。

概念としては、データには5つのクラスターがあり、クラスターごとに共変量の分布と処置効果が異なります。処置はImbens and Rubin（2015）の4章で紹介されたベルヌーイ試行で割り当てられます。

set.seed(12345)
demo_df <- tibble::tibble(
  # クラスターを指定する。モデルには見せない
  latent_group = sample(1:5, size = 5000, replace = TRUE)
) |>
  dplyr::mutate(
    x1 = purrr::map_dbl(
      latent_group,
      \(x){
        # クラスターでx1共変量がサンプルされる分布を決める
        dplyr::case_when(
          x %in% c(1,2) ~ rnorm(1, -2, 0.5),
          x %in% c(3) ~ rnorm(1, 0, 0.3),
          TRUE ~ rnorm(1, 2, 0.2)
        )
      }
    ),
    x2 = purrr::map_dbl(
      latent_group,
      \(x){
        # クラスターでx2共変量がサンプルされる分布を決める
        dplyr::case_when(
          x %in% c(2,5) ~ rnorm(1, -2, 0.3),
          x %in% c(3) ~ rnorm(1, 0, 0.6),
          TRUE ~ rnorm(1, 2, 0.2)
        )
        }
      ),
    # 処置の割り当てはベルヌーイ試行
    treatment = rbinom(dplyr::n(), 1, 0.5),
    error_term = rnorm(dplyr::n(), 0, 1),
    outcome = dplyr::case_when(
      # クラスターで処置効果の生成モデルを決める
      # クラスター1の処置効果は1、クラスター2は-1、クラスター3は10、クラスター4は-5、クラスター5は6
      latent_group == 1 ~ 1 * treatment + error_term,
      latent_group == 2 ~ -1 * treatment + error_term,
      latent_group == 3 ~ 10 * treatment + error_term,
      latent_group == 4 ~ -5 * treatment + error_term,
      TRUE ~ 6 * treatment + error_term
    )
    )

まず共変量x1とx2を確認しましょう：

demo_df |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = x1, y = x2, color = as.factor(latent_group))) + 
  ggplot2::labs(color = "隠れクラスター") +
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

こんな感じで5つのクラスターがあり、かつ共変量の分散が不均一なクラスターもあります。モデルは果たしてこのデータから私が設定した処置効果を正しく発見してくれるかを確認しましょう！

モデル定式化

クラスタリング用のモデルはディリクレ過程です。

ディリクレ過程とは、クラスター数と観測値がどのクラスターに属するかを自動で推論する手法である。分析者が明示的にクラスター数を指定しなくてもいいところが強みです。

どのようにクラスター数と観測値がどのクラスターに属するかを推論するのかと言いますと、まず前提として無限のグループが「選ばれうる」が、観測されたデータとモデルの仮定を吟味して、無限のグループの中から有限のグループを自動で選択して実際のクラスタリングを行います（Grimmer 2011）。

詳細はこちらの記事も参考してください：

確率モデル

まずは確率モデルを書きます、詳細の説明は下の方で行います。

全体のディリクレ過程

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

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

全てのクラスターkについて

$$sigma_{x_{1},k}\sim Gamma(0.1, 0.1)$$

$$sigma_{x_{2},k}\sim Gamma(0.1, 0.1)$$

$$latent_{x_{1},k} \sim Normal(0, 5)$$

$$latent_{x_{1},k} \sim Normal(0, 5)$$

観測値nについて

$$\eta_{n} \sim Categorical(z)$$

$$x_{1_{n}} \sim Normal(latent_{x_{1, \eta_{n}}}, sigma_{x_{1},{\eta_{n}}})$$

$$x_{2_{n}} \sim Normal(latent_{x_{2, \eta_{n}}}, sigma_{x_{2},{\eta_{n}}})$$

要するに、まず各クラスターがそれぞれの中心となるx1とx2を選択してから、そのクラスター内のx1とx2のばらつきの範囲を決めるという確率モデルです。

ただ、実際に実装する際は、Stan言語の離散確率変数パラメータに対応しない仕様で、$\eta$を周辺化する必要があります。

Stanでの実装

Stanでの実装コードはこちらです。ファイル名がcem_demo.stanになっているのは、前述のCEM（Iacus, King and Porro 2012）のように乱塊法（blocked randomized experiment）を理想として近似しているからです。CEMと乱塊法の関連の詳細はKing and Nielsen（2019）を参照してください（特にp. 442の4.1辺り）。

cem_demo.stan

data {
  int N;
  int K;
  
  array[N] real x1;
  array[N] real x2;
  
  int N_treated;
  int N_control;
  
  array[N] int treatment;
  array[N] real outcome;
}
parameters {
  real<lower=0> d_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[K - 1] breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  
  array[K] real<lower=0> x1_sigma;
  array[K] real x1_latent;
  array[K] real<lower=0> x2_sigma;
  array[K] real x2_latent;
}
transformed parameters {
  simplex[K] 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 (k in 2:(K - 1)) {
      z[k] = (1 - sum) * breaks[k];
      sum += z[k];
    }
    z[K] = 1 - sum;
  }
}
model {
  x1_sigma ~ gamma(0.1, 0.1);
  x1_latent ~ normal(0, 5);
  x2_sigma ~ gamma(0.1, 0.1);
  x2_latent ~ normal(0, 5);
  
  for (i in 1:N){
    vector[K] case_vector;
    for (k in 1:K){
      case_vector[k] = log(z[k]) + 
                       normal_lpdf(x1[i] | x1_latent[k], x1_sigma[k]) + 
                       normal_lpdf(x2[i] | x2_latent[k], x2_sigma[k]);
    }
    target += log_sum_exp(case_vector);
  }
}
generated quantities {
  array[N] vector[K] eta;
  vector[K] group_ate;

  for (i in 1:N){
    vector[K] this_eta;
    for (k in 1:K){
      this_eta[k] = z[k] * exp(
                         normal_lpdf(x1[i] | x1_latent[k], x1_sigma[k]) + 
                         normal_lpdf(x2[i] | x2_latent[k], x2_sigma[k])
      );
    }
    eta[i] = this_eta./sum(this_eta);
  }
  for (k in 1:K){
    vector[N_treated] this_treated;
    vector[N_treated] this_treated_eta;
    vector[N_control] this_control;
    vector[N_control] this_control_eta;
    int count_treated = 1;
    int count_control = 1;
    for (i in 1:N){
      if (treatment[i] == 1){
        this_treated[count_treated] = outcome[i] * eta[i, k];
        this_treated_eta[count_treated] = eta[i, k];
        count_treated += 1;
      }
      else {
        this_control[count_control] = outcome[i] * eta[i, k];
        this_control_eta[count_control] = eta[i, k];
        count_control += 1;
      }
    group_ate[k] = sum(this_treated)/sum(this_treated_eta) - sum(this_control)/sum(this_control_eta);
    }
  }
}

モデル推定

ここではStanに投入するための処理を実施して、実際のモデル推定を行いましょう！

demo_list <- list(
  N = nrow(demo_df),
  # わざと多めに最大クラスター数を設定して、クラスター数を発見する性能を確認する
  K = 10,
  
  x1 = demo_df$x1,
  x2 = demo_df$x2,
  
  N_treated = sum(demo_df$treatment == 1),
  N_control = sum(demo_df$treatment == 0),
  
  treatment = demo_df$treatment,
  outcome = demo_df$outcome
)


m_demo_cem_init <- cmdstanr::cmdstan_model("cem_demo.stan")

> m_demo_cem_estimate <- m_demo_cem_init$variational(
     seed = 12345,
     iter = 50000,
     data = demo_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.007347 seconds 
1000 transitions using 10 leapfrog steps per transition would take 73.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       -17286.389             1.000            1.000 
   200       -13254.972             0.652            1.000 
   300       -10429.013             0.525            0.304 
   400       -10372.743             0.395            0.304 
   500       -10351.234             0.317            0.271 
   600       -10404.807             0.265            0.271 
   700       -11145.974             0.236            0.066 
   800       -10306.938             0.217            0.081 
   900       -10330.912             0.193            0.066 
  1000       -10292.896             0.174            0.066 
  1100       -10373.761             0.159            0.008   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  124.2 seconds.

2分で終わりましたー

最後は結果のsummaryを保存しましょう：

m_demo_cem_summary <- m_demo_cem_estimate$summary()

結果の可視化

クラスター数の精度

まずは推定されたクラスター数の全体の割合から確認します：

> m_demo_cem_summary |> dplyr::filter(stringr::str_detect(variable, "z"))
# A tibble: 10 × 7
   variable      mean    median       sd       mad         q5      q95
   <chr>        <dbl>     <dbl>    <dbl>     <dbl>      <dbl>    <dbl>
 1 z[1]     0.199     0.199     0.00497  0.00484   0.190      0.207   
 2 z[2]     0.201     0.201     0.00608  0.00588   0.191      0.211   
 3 z[3]     0.000277  0.000162  0.000340 0.000149  0.0000275  0.000898
 4 z[4]     0.196     0.196     0.00540  0.00520   0.187      0.205   
 5 z[5]     0.000170  0.000112  0.000208 0.0000906 0.0000256  0.000491
 6 z[6]     0.000221  0.000133  0.000297 0.000122  0.0000272  0.000670
 7 z[7]     0.200     0.200     0.00548  0.00549   0.191      0.209   
 8 z[8]     0.204     0.203     0.00578  0.00572   0.194      0.213   
 9 z[9]     0.0000956 0.0000488 0.000148 0.0000564 0.00000287 0.000356
10 z[10]    0.000101  0.0000491 0.000155 0.0000568 0.00000299 0.000384

上のコードのコメントに書いているように、実際が5つのクラスターしかないですが、クラスター数自体を推定する性能を確認するため、わざと最大クラスター数10をモデルに渡しました。

結果として、クラスター1、2、4、7、8がそれぞれおよそ20%を占めていて、他のクラスターの割合は0に近いです。これは私が設定した、5つのクラスターがそれぞれ20%を占めるという状況をうまく学習できたかと思います。

処置効果の推定精度

続いては処置効果の推定結果です：

> m_demo_cem_summary |> dplyr::filter(stringr::str_detect(variable, "group_ate"))
# A tibble: 10 × 7
   variable         mean median         sd      mad     q5    q95
   <chr>           <dbl>  <dbl>      <dbl>    <dbl>  <dbl>  <dbl>
 1 group_ate[1]    9.93   9.93   0.00125    0.00122  9.93   9.93 
 2 group_ate[2]    6.15   6.15   0.000121   0        6.15   6.15 
 3 group_ate[3]  NaN     NA     NA         NA       NA     NA    
 4 group_ate[4]   -1.08  -1.08   0.00305    0.00308 -1.08  -1.07 
 5 group_ate[5]  NaN     NA     NA         NA       NA     NA    
 6 group_ate[6]  NaN     NA     NA         NA       NA     NA    
 7 group_ate[7]   -5.09  -5.09   0.0000974  0       -5.09  -5.09 
 8 group_ate[8]    0.976  0.976  0.00334    0.00328  0.971  0.982
 9 group_ate[9]  NaN     NA     NA         NA       NA     NA    
10 group_ate[10] NaN     NA     NA         NA       NA     NA

クラスター1の処置効果は9.93になっていて、シミュレーションのクラスター3の処置効果（10）に近い数字です。

クラスター2の処置効果は6.15になっていて、シミュレーションのクラスター5の処置効果（6）に近い数字です。

クラスター4の処置効果は-1.08になっていて、シミュレーションのクラスター2の処置効果（-1）に近い数字です。

クラスター7の処置効果は-5.09になっていて、シミュレーションのクラスター4の処置効果（-5）に近い数字です。

クラスター8の処置効果は0.976になっていて、シミュレーションのクラスター1の処置効果（1）に近い数字です。

クラスターの番号は本質的ではないので、処置効果もうまく学習できたと思います。

一部の「利用されていない」クラスターの処置効果がおそらく重み付け平均の分母が小さすぎてNaNになっていますが、これは今後Stanのgenerated quantitiesブロック内の例外処理で改善したいと思います。

ここで一点強調したいのは、summaryのtibbleに格納された信用区間（q5とq95）は、頻度主義的なサンプリングの誤差ではなく、ディリクレ過程のクラスタリング時の不確実性から派生された不確実性です。なので結果変数自体のばらつきは処置効果の信用区間に反映されていません。

推定の意思決定の可視化

さて、最後は各クラスターの$\eta$（因果効果計算の重み）を大きさで可視化することで、モデルがどのような比較をして処置効果を計算したのかを可視化しましょう：

demo_df |>
  dplyr::bind_cols(
    eta = eta[,1]
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = x1, y = x2, color = as.factor(treatment), size = eta)) + 
  ggplot2::scale_color_manual(values = c("red", "blue")) + 
  ggplot2::labs(color = "処置フラグ", size = "因果効果計算の重み") + 
  ggplot2::ggtitle("クラスター1の処置効果計算ロジックの可視化") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

このように、真ん中のクラスターの処置効果を計算する際は、真ん中の観測値の重みが最も高く、他のクラスターの重みはどれも低く、情報はあまり計算で利用されません。ただ、一つ着目していただきたいのは、左下の真ん中のクラスターと被っているところの観測値は、少し大きめな点になっています。これは前述したように、「類似していれば、他のクラスターの観測値もある程度利用する」というこの手法の強みです。

このシミュレーションデータの場合、処置効果は乱暴に（？）設定されたんですが、現実世界では、確率分布が近いクラスターは処置効果も近いはずです。このように少し異なる観測値から情報を借りることで、実質サンプル数を増やしてより精度の高い推定が期待されます。

全てのクラスターを可視化すると長くなるので、最後は実際にモデルに利用されたクラスター2とモデルに利用されていないクラスター3を可視化しましょう：

demo_df |>
  dplyr::bind_cols(
    eta = eta[,2]
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = x1, y = x2, color = as.factor(treatment), size = eta)) + 
  ggplot2::scale_color_manual(values = c("red", "blue")) + 
  ggplot2::labs(color = "処置フラグ", size = "因果効果計算の重み") + 
  ggplot2::ggtitle("クラスター2の処置効果計算ロジックの可視化") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

クラスター2は定性的にはクラスター1と変わらないですね！

クラスター3がどうなっているかというと：

demo_df |>
  dplyr::bind_cols(
    eta = eta[,3]
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = x1, y = x2, color = as.factor(treatment), size = eta)) + 
  ggplot2::scale_color_manual(values = c("red", "blue")) + 
  ggplot2::labs(color = "処置フラグ", size = "因果効果計算の重み") + 
  ggplot2::ggtitle("クラスター3の処置効果計算ロジックの可視化") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

どの観測値の重みもほぼ0です（点の大きさが示す数字の違いに気をつけていただきたいです）。

結論

いかがでしたか？

この手法を使えばクラスターごとに存在する異質処置効果を推定できるだけでなく、効果の推定で具体的にはどの観測値を一番参考にしているのかも簡単に可視化できます。

また、クラスタリングは確率モデルを使っているため、文字や都道府県、職業、国などのカテゴリデータや電話数などのカウントデータにも柔軟に対応できます。

カテゴリデータの場合はx1とx2の正規分布をカテゴリ分布、カウントデータの場合はポワソン分布などにすれば問題ないです。

皆さんはぜひこのモデルをご自身の研究か業務に活かしてみてください！

参考文献

Athey, Susan, Julie Tibshirani, and Stefan Wager. "Generalized Random Forests." The Annals of Statistics 47. 2 (2019): 1148-1178.

Fong, Christian, and Justin Grimmer. "Discovery of treatments from text corpora." Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers). 2016.

Grimmer, Justin. "An introduction to Bayesian inference via variational approximations." Political Analysis 19.1 (2011): 32-47.

Ho, Daniel E., Kosuke Imai, Gary King and Elizabeth A. Stuart. "Matching as nonparametric preprocessing for reducing model dependence in parametric causal inference." Political analysis 15.3 (2007): 199-236.

Iacus, Stefano M., Gary King, and Giuseppe Porro. "Causal inference without balance checking: Coarsened exact matching." Political analysis 20.1 (2012): 1-24.

Imai, Kosuke, and Dustin Tingley. "A statistical method for empirical testing of competing theories." American Journal of Political Science 56.1 (2012): 218-236.

Imbens, Guido W., and Donald B. Rubin. Causal inference in statistics, social, and biomedical sciences. Cambridge university press, 2015.

King, Gary, and Richard Nielsen. "Why propensity scores should not be used for matching." Political analysis 27.4 (2019): 435-454.

Montgomery, Jacob M., Brendan Nyhan, and Michelle Torres. "How conditioning on posttreatment variables can ruin your experiment and what to do about it." American Journal of Political Science 62.3 (2018): 760-775.

Rosenbaum, Paul R., and Donald B. Rubin. "The central role of the propensity score in observational studies for causal effects." Biometrika 70.1 (1983): 41-55.

Wager, Stefan, and Susan Athey. "Estimation and inference of heterogeneous treatment effects using random forests." Journal of the American Statistical Association 113.523 (2018): 1228-1242.

Zhou, Xiang, and Teppei Yamamoto. "Tracing causal paths from experimental and observational data." The Journal of Politics 85.1 (2023): 250-265.