回帰不連続デザインに入門

会社の勉強会で回帰不連続デザインについて発表する機会があったので、整理した内容をこちらにも挙げておこうと思います。主にこちらの書籍を勉強していきました。
個人的な解釈も多分にあるので、誤りなどあればご指摘いただければと思います。

Amazon.co.jp: 統計的因果推論の理論と実装

傾向スコアに入門 / マッチングとIPTW

会社の勉強会で効果検証入門に取り組みましたが、実装はほとんど手をつけられていなかったり、個人的に気になったこともあったので、傾向スコアの部分を整理していきます。また個人的な解釈も多分にあるので、誤りなどあればご指摘いただければと思います。

https://www.amazon.co.jp/%E5%8A%B9%E6%9E%9C%E6%A4%9C%E8%A8%BC%E5%85%A5%E9%96%80%E3%80%9C%E6%AD%A3%E3%81%97%E3%81%84%E6%AF%94%E8%BC%83%E3%81%AE%E3%81%9F%E3%82%81%E3%81%AE%E5%9B%A0%E6%9E%9C%E6%8E%A8%E8%AB%96-%E8%A8%88%E9%87%8F%E7%B5%8C%E6%B8%88%E5%AD%A6%E3%81%AE%E5%9F%BA%E7%A4%8E-%E5%AE%89%E4%BA%95-%E7%BF%94%E5%A4%AA/dp/4297111179

こちらの書籍はRで実装(以下リンク参照)されているので、今回はRで実装しました。

github.com

また今回整理していたnotebookはこちらに置いてあります。
https://gist.github.com/wakamatsuikuma/ef8cf85b15378b731888fc19c6d9e16f

目次    

因果推論と傾向スコアについて

因果推論について

 因果推論はある処置の効果測定をする手法のことですが、実際はRCT(A/Bテスト)でデザインされたデータでの分析が理想的です。しかしRCTでは多大なコストがかかるため、観察データに対してRCTの結果を近似する様な方法論を提供してくれる因果推論を用いるというのが、ビジネスにおける因果推論の必要性となってきます。

傾向スコアとは

 傾向スコアとは、観察データによる効果をRCTによる効果に近似するためのツールであり、各個体において介入が行われる確率になります。この確率が同一の個体で行った比較はセレクションバイアスが軽減されるという状況を利用して分析を行います。 これは傾向スコアが同一になるようなサンプルの中では、介入がY(0)とは独立に割り振られているという仮定に基づいて、効果を推定することになります。

傾向スコアによる効果推定(その他の因果推論でも)をする上での仮定

 傾向スコアによる効果推定をする上では、強い意味での無視可能な割付けが成立しておく必要があり(これはそのまま、以下の傾向スコアの定理が成立するための条件です)、以下の3点が成立していなければなりません。

① SUTVA(サトヴァ、因果推論する上で重要な仮定)
1. 相互干渉がない: 個体Aに対しての処置が、別の個体Bに対して影響を及ばさない
2. 個体に対する隠れた処置がない: ある処置を受ける個体が、その別の形で受けてはいけない
SUTVAは、ある処置と比較して、別の処置の因果効果を除外するために、外的な情報や背景情報に依拠すると仮定する制約(除外制約)となります。

またSUTVAが成立と仮定した上で、以下の2式が成立します。
② 条件付き独立(無交絡性): 共変量Xを条件としたときに、介入変数Zが潜在的結果変数の組{Y(1), Y(0)} に依存しない(処置の割り付けに影響を与えるのは観測された共変量のみ)

{Y(1), Y(0)} ⊥ Z | X

③ 条件付き正値性: 共変量Xを条件としたときに、どの個体も介入群または対照群に割付けられる確率が0または1でないこと(どの個体も、どちらにも割付けられる可能性がある)

0 < Pr(Z = 1 | X) < 1

この3点が傾向スコアを活用する上での、必要な仮定となります。
また以下が傾向スコアの定理です。これらの仮定から傾向スコアが介入の選択基準とは無関係であるといえるので、擬似的にRCTを再現していることになります。
【傾向スコアの定理】
① バランシング: 介入変数Zと観測された共変量Xは、傾向スコアe(X)が与えられたとき、条件付き独立である(つまり、傾向スコアe(X)が同じ値であれば、介入群と対照群における共変量Xの分布は同じ)。

X ⊥ Z | e(X)

② 条件付き独立性: 傾向スコアe(X)が与えられれば、潜在的結果変数{Y(1), Y(0)}と介入変数Zは条件付き独立である。

{Y(1), Y(0)} ⊥ Z | e(X)

傾向スコアを用いた効果推定

 まず傾向スコアは、介入変数Zを目的変数、共変量を説明変数としたロジスティック回帰で推定されます。また確率として算出できれば良いので、Gradient Boosting Decision Tree(GBDT)も選択として考えられます。

【手法①: 傾向スコアマッチング】
 傾向スコアマッチングは、"(1): 介入された個体"に対して"(2): (1)と同じ傾向スコアを持つ非介入のサンプル"の結果変数を、(1)の非介入時の結果変数と考えます。つまりこの(2)を任意のアルゴリズム(最近傍法、最適ペアマッチング、復元の有無、etc•••)で探し出し、各グループ間の結果変数の差を効果の推定値とします。ただし傾向スコアが高いor低いとこを捨てることになり、解析結果が当てはまる対象集団が限定的になると以下の記事でも述べられていました。

(https://riklog.com/research/propensity-score/)

このときの効果推定の計算式は、


\hat{τ} _ {match} = E( E(Y | \hat{P}(X), Z=1) - E(Y | \hat(P)(X), Z=0) | Z=1)

となります。またこの方法で求められる効果は、ATT(Average Treatment effect on Treated)になりますが、マッチングでもATE(Average Treatment effect)を求める方法はあります(この記事では触れません)。

【手法②: 傾向スコアを用いたIPTW】
 IPTWはInverse Probability of Treatment Weighting(逆確率重み付き)推定と言います。
傾向スコアの逆数を重みに使って、全サンプルが介入受けた場合と非介入だった場合の各期待値を算出し、その差分を効果の推定値とします。

 つまりIPTWでは、データを捨てずに効果推定を達成しようとします。そもそも傾向スコアは介入が行われる確率なので、RCTが実現しているデータセットにおいては0.5になるはずであり、そこからずれるということは、何かしらの理由で介入が偏りセレクションバイアスが発生している状態と考えられます。したがって介入群に属しているかつ傾向スコア(介入が行われる確率)が低い個体は、RCTであった場合と比較して、観測した量が少ない個体と考えられます。これは他の個体においても同様と考えられるので、傾向スコアの逆数で重み付けすることで、見かけ上各傾向スコアにおける観測した量を同じにし、セレクションバイアスが無い状態を作り出しています。
 効果推定の計算式は、

 
\hat{τ} _ {IPTW} 
= Y^{(1)} - Y^{(0)}
=
\sum_{i=1}^N\frac{Z _ iY _  i}{\hat{P}(X _ i)} / \sum_{i=1}^N\frac{Z _ i}{\hat{P}(X _ i)}
-
\sum_{i=1}^N\frac{(1-Z _ i)Y _  i}{1-\hat{P}(X _ i)} / \sum_{i=1}^N\frac{(1-Z _ i)}{1-\hat{P}(X _ i)}

となります。データ全体での平均的な効果を推定していることになるので、ATEの推定となります。ただし重み付けのやり方を変えるとATTの推定も可能であり、記事の最後でも触れています。

回帰分析と比較したメリット、デメリット

[回帰分析]

  • OVB(脱落変数バイアス)やSensitivity Analysisといった分析がしやすくなるものがあり、取り組みやすい
  • 共変量の選定が重要であり、入念なモデリングが求められる

[傾向スコア]

  • 介入の決定方法に関する情報についてのヒアリングを行うだけで良く実用的
  • 傾向スコアマッチング: 計算時間がかかる、解析結果が当てはまる対象集団が限定的
  • IPTW: 傾向スコアが0に近いサンプルがあると計算が不安定になりやすい

LaLondeデータセットの分析

 それではこれらの手法を実装していきます。今回の分析での課題となる効果推定は、就労経験を与えることが収入にどのくらい影響するかです。各手法の効果推定の精度を評価するために、まずRCTデータでの効果推定値を実際の効果とします。そしてRCTデータに異なるデータセットを足し、バイアスのあるデータを作成します。このバイアスデータに対して、各手法による効果推定をしていきます。

NSWとCPSのデータ確認

 まずはデータについてです。用いるデータはNSW(National Supported Work)とCPS(Current Population Survey)です。NSWは条件を満たす失業者に対して短期的な就労経験を与えることで就職を助ける試みです。CPSは、NSW外で得られた調査データです。この2つのデータセットを用いて、以下の3つのデータセットを作成します。

  • nswdw_data(NSW): NSWのデータ(RCTにより得られたデータ)
  • cps1_nsw_data(CPS1): NSWの対照群の代わりにCPSを用いたデータ(バイアスのあるデータ)
  • cps3_nsw_data(CPS3): CPSにおいて1976年春時点で雇用されていないと回答しているサンプルに対照群を限定しており、CPS1の対照群から就労経験の必要性が低いと考えられるサンプルを取り除いたデータ

ここから実装です。 まずはパッケージのインストール、ライブラリの読み込み、データの準備です。

# havenパッケージのインストール
install.packages("haven")
install.packages("MatchIt")
install.packages("WeightIt")
install.packages("cobalt")
install.packages("psych")
install.packages("reshape2")

# ライブラリの読み込み
library("tidyverse")
library("haven")
library("broom")
library("MatchIt")
library("WeightIt")
library("cobalt")
library("psych")
library("ggplot2")
library("reshape2")

# データ読み込み
nswdw_data <- read_dta("https://users.nber.org/~rdehejia/data/nsw_dw.dta")
cps1_data <- read_dta("https://users.nber.org/~rdehejia/data/cps_controls.dta")
cps3_data <- read_dta("https://users.nber.org/~rdehejia/data/cps_controls3.dta")

# NSWデータから介入グループだけを取り出して、cps1における介入グループとして扱う
cps1_nsw_data <- nswdw_data %>%
  filter(treat==1) %>%
  dplyr::bind_rows(cps1_data)

# NSWデータから介入グループだけを取り出して、cps3における介入グループとして扱う
cps3_nsw_data <- nswdw_data %>%
  filter(treat==1) %>%
  dplyr::bind_rows(cps3_data)

次にNSWのデータを集計してみます(出力結果はNSWのみに割愛)。だいたいここで色々比較しますが、ここでは変数の確認までとします。
変数としては以下が確認できました。

  • treat: 就労経験を与えたかどうか(介入変数)
  • age: 年齢
  • education: 学歴
  • black, hispanic: 人種
  • married: 結婚しているか
  • nodegree: 学位がないか
  • re78, re75, re74: その年の収入(re78が結果変数)
# データ確認の関数定義
data_check <- function(df){
  print("【形状】:")
  print(dim(df))

  print("【データ型】:")
  print(sapply(df, class))

  print("【欠損量】:")
  print(colSums(is.na(df)))

  print("【ユニーク数】:")
  print(apply(df, 2, function(x) length(unique(x))))

  print("【先頭行】:")
  print(head(df, n=5))

  print("【統計量】")
  summary(df)
}

data_check(nswdw_data)
data_check(cps1_data)
data_check(cps3_data)


[output]
[1] "【形状】:"
[1] 445  11
[1] "【データ型】:"
    data_id       treat         age   education       black    hispanic 
"character"   "numeric"   "numeric"   "numeric"   "numeric"   "numeric" 
    married    nodegree        re74        re75        re78 
  "numeric"   "numeric"   "numeric"   "numeric"   "numeric" 
[1] "【欠損量】:"
  data_id     treat       age education     black  hispanic   married  nodegree 
        0         0         0         0         0         0         0         0 
     re74      re75      re78 
        0         0         0 
[1] "【ユニーク数】:"
  data_id     treat       age education     black  hispanic   married  nodegree 
        1         2        34        14         2         2         2         2 
     re74      re75      re78 
      115       155       308 
[1] "【先頭行】:"
# A tibble: 5 × 11
  data_id      treat   age education black hispanic married nodegree  re74  re75
  <chr>        <dbl> <dbl>     <dbl> <dbl>    <dbl>   <dbl>    <dbl> <dbl> <dbl>
1 Dehejia-Wah…     1    37        11     1        0       1        1     0     0
2 Dehejia-Wah…     1    22         9     0        1       0        1     0     0
3 Dehejia-Wah…     1    30        12     1        0       0        0     0     0
4 Dehejia-Wah…     1    27        11     1        0       0        1     0     0
5 Dehejia-Wah…     1    33         8     1        0       0        1     0     0
# ℹ 1 more variable: re78 <dbl>
[1] "【統計量】"
   data_id              treat             age          education   
 Length:445         Min.   :0.0000   Min.   :17.00   Min.   : 3.0  
 Class :character   1st Qu.:0.0000   1st Qu.:20.00   1st Qu.: 9.0  
 Mode  :character   Median :0.0000   Median :24.00   Median :10.0  
                    Mean   :0.4157   Mean   :25.37   Mean   :10.2  
                    3rd Qu.:1.0000   3rd Qu.:28.00   3rd Qu.:11.0  
                    Max.   :1.0000   Max.   :55.00   Max.   :16.0  
     black           hispanic          married          nodegree    
 Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.000  
 1st Qu.:1.0000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:1.000  
 Median :1.0000   Median :0.00000   Median :0.0000   Median :1.000  
 Mean   :0.8337   Mean   :0.08764   Mean   :0.1685   Mean   :0.782  
 3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:1.000  
 Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.000  
      re74              re75            re78      
 Min.   :    0.0   Min.   :    0   Min.   :    0  
 1st Qu.:    0.0   1st Qu.:    0   1st Qu.:    0  
 Median :    0.0   Median :    0   Median : 3702  
 Mean   : 2102.3   Mean   : 1377   Mean   : 5301  
 3rd Qu.:  824.4   3rd Qu.: 1221   3rd Qu.: 8125  
 Max.   :39570.7   Max.   :25142   Max.   :60308  

またここでは各データの介入変数ごとの結果変数の分布を確認しておきます。

data = nswdw_data
ggplot(data, aes(x=re78, fill=as.factor(treat))) +
  geom_histogram(position = "identity", alpha=0.6, binwidth=2000) +
  ggtitle("distribution of re78") +
  theme(text=element_text(size=20))

data = cps1_nsw_data
ggplot(data, aes(x=re78, fill=as.factor(treat))) +
  geom_histogram(position = "identity", alpha=0.6, binwidth=2000) +
  # scale_y_log10() +
  ggtitle("distribution of re78") +
  theme(text=element_text(size=20))

data = cps3_nsw_data
ggplot(data, aes(x=re78, fill=as.factor(treat))) +
  geom_histogram(position = "identity", alpha=0.6, binwidth=2000) +
  ggtitle("distribution of re78") +
  theme(text=element_text(size=20))

re78の分布

CPS3は就労経験の必要性が低いサンプルを取り除いているためか、NSWと同様に介入群と対照群の分布が似ていることが確認できます。一方CPS1はCPSデータをそのまま追加しているので、対照群のデータがかなり多く、介入群と対照群の分布も異なることが確認できます。

RCTデータの評価

 それではまずRCTデータで、ナイーブな推定と回帰分析による効果推定です。

# RCTデータであるnswdw_dataで回帰分析

# ナイーブな推定: グループ間の平均の比較
tau_nsw_naive <- mean(filter(nswdw_data, treat==1)$re78) - mean(filter(nswdw_data, treat==0)$re78)
paste("nswdw_dataのナイーブな推定効果: ", tau_nsw_naive)

# 共変量なしの回帰分析
nsw_nocov <- lm(data = nswdw_data, 
formula = re78 ~ treat) %>%
  tidy() # t検定などの結果を省略。各変数について、推定効果、標準誤差、t値, p値。
print(nsw_nocov)

# 共変量ありの回帰分析
nsw_cov <- lm(data = nswdw_data, 
formula = re78 ~ treat + age + black + hispanic + married + 
  nodegree + education + re74 + re75) %>%
  tidy() # t検定などの結果を省略
print(nsw_cov)

[output]
'nswdw_dataのナイーブな推定効果:  1794.3423818501'
# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    4555.      408.     11.2  1.15e-25
2 treat          1794.      633.      2.84 4.79e- 3
# A tibble: 10 × 5
   term          estimate std.error statistic p.value
   <chr>            <dbl>     <dbl>     <dbl>   <dbl>
 1 (Intercept)   785.     3375.        0.233  0.816  
 2 treat        1676.      639.        2.62   0.00898
 3 age            55.3      45.3       1.22   0.223  
 4 black       -2160.     1169.       -1.85   0.0654 
 5 hispanic      164.     1549.        0.106  0.916  
 6 married      -139.      880.       -0.158  0.875  
 7 nodegree      -70.7    1004.       -0.0704 0.944  
 8 education     396.      227.        1.74   0.0825 
 9 re74            0.0821    0.0774    1.06   0.289  
10 re75            0.0528    0.135     0.389  0.697  

RCTなデータだけあってか、どれも大体一緒です(この結果からデータセットがRCTと判断できるわけでは無いです、多分できないかと)。共変量ありの回帰分析に注目し、$1,676の収入増加が実際の介入効果とします。統計的な妥当性としても、有意水準0.05 > p値0.00898であることから有意であると判断して良いかと思います。

バイアスデータの評価

ナイーブな推定と回帰分析

 ここからバイアスデータに対しての効果推定です。まずはナイーブな推定を確認すると、

# ナイーブな推定: グループ間の平均の比較
tau_cps1_naive <- mean(filter(cps1_nsw_data, treat==1)$re78) - mean(filter(cps1_nsw_data, treat==0)$re78)
tau_cps3_naive <- mean(filter(cps3_nsw_data, treat==1)$re78) - mean(filter(cps3_nsw_data, treat==0)$re78)

paste("cps1_nsw_dataのナイーブな推定効果: ", tau_cps1_naive)
paste("cps3_nsw_dataのナイーブな推定効果: ", tau_cps3_naive)

[output]
'cps1_nsw_dataのナイーブな推定効果:  -8497.51614813298'
'cps3_nsw_dataのナイーブな推定効果:  -635.026250406911

全然ダメです。就労支援すると収入ダウンという推定結果になっています。データの分布確認では、NSWと似た形であったCPS3でも収入ダウンという結果であり、単純な集計では誤った効果推定をしてしまうことが確認できました。

それでは回帰分析をしてみます。

# 回帰分析
cps1_reg <- lm(data = cps1_nsw_data, 
formula = re78 ~ treat + age + black + hispanic + married + 
  nodegree + education + re74 + re75) %>%
    tidy()

cps3_reg <- lm(data = cps3_nsw_data,
formula = re78 ~ treat + age + black + hispanic + married + 
  nodegree + education + re74 + re75) %>%
    tidy()

print(cps1_reg)
print(cps3_reg)

[output]
# A tibble: 10 × 5
   term        estimate std.error statistic   p.value
   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
 1 (Intercept) 5736.     445.        12.9   8.48e- 38
 2 treat        699.     548.         1.28  2.02e-  1
 3 age         -102.       5.88     -17.3   1.33e- 66
 4 black       -837.     213.        -3.93  8.44e-  5
 5 hispanic    -218.     219.        -0.998 3.18e-  1
 6 married       73.1    142.         0.513 6.08e-  1
 7 nodegree     372.     178.         2.10  3.61e-  2
 8 education    160.      28.6        5.60  2.15e-  8
 9 re74           0.289    0.0121    24.0   1.26e-124
10 re75           0.471    0.0122    38.7   2.28e-313
# A tibble: 10 × 5
   term         estimate std.error statistic     p.value
   <chr>           <dbl>     <dbl>     <dbl>       <dbl>
 1 (Intercept)    66.5   2437.        0.0273 0.978      
 2 treat        1548.     781.        1.98   0.0480     
 3 age            13.0     32.5       0.399  0.690      
 4 black       -1241.     769.       -1.61   0.107      
 5 hispanic      499.     942.        0.530  0.597      
 6 married       407.     695.        0.585  0.559      
 7 nodegree      260.     847.        0.307  0.759      
 8 education     404.     159.        2.54   0.0113     
 9 re74            0.296    0.0583    5.09   0.000000489
10 re75            0.232    0.105     2.21   0.0273     

 CPS3は介入効果が$1548アップと良く推定できているかと思います。ただCPS1は介入効果が$699と、介入効果を低く推定しています。これは介入群が少ないこと および 介入群と対照群でデータの分布(特徴)が異なるためと考えられます。NSWはそもそも介入(職業訓練)の必要性がある人たちを対象にRCTを行っていますが、CPSはそうではありません。そのため、類似サンプル内において介入・非介入の偏りがあったり、サンプル間で効果が非線形であることが考えられます。つまり、効果推定を行うデータセットの分布の乖離が大きいと言ってよいかと思います。この状態のデータセットでは、回帰分析の計算の仕組み上(*1)うまく推定ができません。その介入によって、どういった特徴や前提のあるデータを対象に効果推定したいのかを明確にして、それが可能となるデータセットの作成・準備、分析をする必要があることが窺えます。

傾向スコアの算出と評価(分布とc統計量)

 それではこの記事の本題である傾向スコアに入っていきます。ロジスティック回帰で傾向スコアを算出し、分布を確認してみます。

# glmで傾向スコア算出
ps_model_cps1 <- glm(data = cps1_nsw_data,
formula = treat ~ age + black + hispanic + married + 
  nodegree + education + re74 + re75, 
  family = binomial # ロジスティック回帰を指定
  )
summary(ps_model_cps1)

ps_model_cps3 <- glm(data = cps3_nsw_data,
formula = treat ~ age + black + hispanic + married + 
  nodegree + education + re74 + re75, 
  family = binomial # ロジスティック回帰を指定
  )
summary(ps_model_cps3)

# fitted.valuesで傾向スコア取得
ps_cps1 <- ps_model_cps1$fitted.values
ps_cps3 <- ps_model_cps3$fitted.values

# データフレームに傾向スコア追加
cps1_nsw_data$ps <- c(ps_cps1)
cps3_nsw_data$ps <- c(ps_cps3)

# ヒストグラム描画
ggplot(cps1_nsw_data, aes(x=ps, fill=as.factor(treat))) +
  geom_histogram(position = "identity", alpha=0.6, binwidth=0.05) +
  scale_y_log10() +
  ggtitle("distribution of ps") +
  theme(text=element_text(size=20))

ggplot(cps3_nsw_data, aes(x=ps, fill=as.factor(treat))) +
  geom_histogram(position = "identity", alpha=0.6, binwidth=0.05) +
  # scale_y_log10() +
  ggtitle("distribution of ps") +
  theme(text=element_text(size=20))

傾向スコアの分布

 CPS1は傾向スコアが0.5付近までしか算出されていませんが、CPS3は傾向スコアが高い範囲において介入群に対して対照群が少ないなという印象です。

 次にc統計量(AUC)を確認します。

# c統計量(AUC)
# ROC曲線で評価
ROC1 <- roc(treat ~ ps, data=cps1_nsw_data)
ROC3 <- roc(treat ~ ps, data=cps3_nsw_data)

# AUC出力
print(ROC1)
print(ROC3)

# ROC曲線描画
plot(ROC1)
plot(ROC3)

ROC曲線

 ヒストグラムを見た時の直感とは異なり、CPS1の方が高い結果となりました。CPS1の方が傾向スコアの分布としてはイマイチな感じがありましたが、AUCは高いです。これは対照群のデータ量が介入群と比べて圧倒的に多く、傾向スコアが低い側の精度が良いためかと思います。

傾向スコアマッチングでATTの推定

 それではまずCPS1に対して、RライブラリのMatchItパッケージを使って傾向スコアマッチングを行なっていきます。またglm()で算出した傾向スコアと差異があるかどうかを、各傾向スコアの差の総和が0かどうかで確認しておきます。マッチングのアルゴリズムとしては、わかりやすい非復元の最近傍法で行います。

# glm()で算出した傾向スコアと差があるか確認。

# matchitで傾向スコア算出
m_near_cps1 <- matchit(data = cps1_nsw_data,
formula = treat ~ age + black + hispanic + married + 
  nodegree + education + re74 + re75,
  method="nearest")

cps1_nsw_data$ps_match <- c(m_near_cps1$distance)

# glm()で算出した傾向スコアとの差
sum(abs(cps1_nsw_data$ps - cps1_nsw_data$ps_match))

[output]
0

 glm()での結果と変わらないのでOKです。

 次に得られた傾向スコアが良い傾向スコアかどうかを判定します。これは上記で確認していたc統計量のような指標が用いられることもある様ですが、近年は基本的に共変量のバランスがとれているかどうかが重要との見解が一般的の様です。
 このバランスがとれているかどうかの判断に使われるものが、傾向スコアによる調整前後における共変量ごとの標準化平均差(ASAM: Average Standard Absolute Mean disitance)になります。これはいわゆる効果量ともいえるもので、これが0.1以下(つまり効果が小さい=差異が小さい)であれば、バランスがとれていると判断します。
 このASAMを簡単に可視化できるcobaltパッケージのlove.plotというものがあるので、傾向スコアの分布と併せて可視化します。

# 調整後データ
matched1_data <- match.data(m_near_cps1)

# 傾向スコアのヒストグラム
ggplot(matched1_data, aes(x=ps, fill=as.factor(treat))) +
  geom_histogram(position = "identity", alpha=0.6, binwidth=0.05) +
  scale_y_log10() +
  ggtitle("distribution of ps") +
  theme(text=element_text(size=20))

# バランシングの評価
love.plot(m_near_cps1, threshold = 0.1, binary="std")

CPS1の傾向スコアマッチング後の傾向スコアとバランシングの可視化

 傾向スコアの分布としては、介入群と対照群の分布がよく似ていることがわかります。また重要なバランシングの方ですが、共変量および傾向スコアのASAMが調整後はかなり小さくなっていることがわかります。最近隣法マッチングによる非復元1対1マッチングは批判されている手法ではありますが、バランシングが悪くないのでOKと考えて良いかと思います。 また以下で、一応検算を。

# ASAMの検算用の関数
standardized_mean_error <- function(df, col){
  
  # ベクトルで値を抽出
  values_treat0 <- filter(df, treat==0)[[col]]
  values_treat1 <- filter(df, treat==1)[[col]]

  # サンプルサイズ
  n_0 <- length(values_treat0)
  n_1 <- length(values_treat1)

  # 平均
  mean_0 <- mean(values_treat0)
  mean_1 <- mean(values_treat1)

  # 分散
  s2 <- function(x){ var(x)*(length(x)-1)/length(x) }
  var_0 <- s2(values_treat0)
  var_1 <- s2(values_treat1)

  # 標準平均化誤差
  m_diff <- mean_1 - mean_0
  s_co <- (var_1*(n_1 - 1) + var_0*(n_0 - 1)) / (n_1 + n_0 - 2)

  return ( m_diff / sqrt(s_co) )
}


# 検算
cols_li <- c("age", "black", "hispanic", "married", "nodegree", "education", "re74", "re75")
for (col in cols_li){ 
  print(paste(col, "のASAM"))
  print(paste("調整前: ", standardized_mean_error(cps1_nsw_data, col)))
  print(paste("調整後: ", standardized_mean_error(matched1_data, col)))
}

[output]
[1] "age のASAM"
[1] "調整前:  -0.673045077618078"
[1] "調整後:  0.2435367704803"
[1] "black のASAM"
[1] "調整前:  2.93324709421367"
[1] "調整後:  0.0437130189471938"
[1] "hispanic のASAM"
[1] "調整前:  -0.0486883897274181"
[1] "調整後:  0.0735214622093808"
[1] "married のASAM"
[1] "調整前:  -1.1552824924244"
[1] "調整後:  0.147042924418762"
[1] "nodegree のASAM"
[1] "調整前:  0.903320101340552"
[1] "調整後:  0.0587058294730288"
[1] "education のASAM"
[1] "調整前:  -0.587471521697691"
[1] "調整後:  -0.00913489647323476"
[1] "re74 のASAM"
[1] "調整前:  -1.25103406168653"
[1] "調整後:  0.0462501126863007"
[1] "re75 のASAM"
[1] "調整前:  -1.31388559386733"
[1] "調整後:  0.0834656684538306"

 なぜかlove.plotと合わない結果となりました。いくつかドキュメントなど確認しましたが、数式も見つからなかったので原因はわかりませんでした。ただ調整後の方が小さいので、とりあえず今回はスルーです。わかる方いたら是非教えていただきたいです。。。

 それでは調整後データを用いて効果推定します。傾向スコアマッチング後にも共変量の影響が残っている可能性があるため、解析モデルにも共変量を取り込むことが推奨されているようです(*2)。

# 回帰分析で効果推定(解析モデルにも共変量を使用)
psm1_co_result <- lm(data = matched1_data, 
formula = re78 ~ treat + age + black + hispanic + married + 
  nodegree + education + re74 + re75) %>%
  tidy()

print(psm1_co_result)

[output]
# A tibble: 10 × 5
   term         estimate std.error statistic p.value
   <chr>           <dbl>     <dbl>     <dbl>   <dbl>
 1 (Intercept)  426.      3479.       0.122   0.903 
 2 treat       1625.       721.       2.26    0.0247
 3 age           -9.10      48.3     -0.188   0.851 
 4 black       -746.      1150.      -0.649   0.517 
 5 hispanic    1504.      1933.       0.778   0.437 
 6 married      361.      1062.       0.340   0.734 
 7 nodegree     274.      1104.       0.248   0.804 
 8 education    421.       215.       1.96    0.0510
 9 re74          -0.0110     0.111   -0.0996  0.921 
10 re75           0.312      0.165    1.89    0.0598

 就労経験により$1,625の収入増加効果があることを示しています。RCTの$1,676とかなり近い推定値であり、調整前の$699と比べてかなり改善しました。有意水準0.05 > p値0.0180でもあり、有意な効果とも判断できます。

 CPS3についても傾向スコアマッチングを行ってみます。

# CPS3についても傾向スコアマッチング
m_near_cps3 <- matchit(data = cps3_nsw_data,
formula = treat ~ age + black + hispanic + married + 
  nodegree + education + re74 + re75, 
  method="nearest")

# 調整後データ
matched3_data <- match.data(m_near_cps3)

# 傾向スコア
ggplot(matched3_data, aes(x=ps, fill=as.factor(treat))) +
  geom_histogram(position = "identity", alpha=0.6, binwidth=0.05) +
  # scale_y_log10() +
  ggtitle("distribution of ps") +
  theme(text=element_text(size=20))

# バランシングの評価
love.plot(m_near_cps3, threshold = 0.1, binary="std")

# 回帰分析(解析モデルにも共変量を使用)
psm3_co_result <- lm(data = matched3_data, 
formula = re78 ~ treat + age + black + hispanic + married + 
  nodegree + education + re74 + re75) %>%
  tidy()

print(psm3_co_result)

[output]
# A tibble: 10 × 5
   term          estimate std.error statistic p.value
   <chr>            <dbl>     <dbl>     <dbl>   <dbl>
 1 (Intercept) -2112.      3512.       -0.601 0.548  
 2 treat        1345.       790.        1.70  0.0894 
 3 age             7.80      42.9       0.182 0.856  
 4 black        -469.      1010.       -0.465 0.642  
 5 hispanic     1064.      1300.        0.818 0.414  
 6 married      -158.       986.       -0.160 0.873  
 7 nodegree      923.      1110.        0.832 0.406  
 8 education     602.       224.        2.69  0.00753
 9 re74            0.0264     0.103     0.256 0.798  
10 re75            0.221      0.160     1.38  0.168  

CPS3の傾向スコアマッチング後の傾向スコアとバランシングの可視化

 CPS3は回帰分析と比べて結果が改悪しました。傾向スコアの分布は介入群と対照群で大きく異なっており、共変量のバランシングもうまくできておりません。 これは、

  • 調整前データの対照群がかなり絞られていた
  • 調整前の傾向スコアの高い領域で介入群の方が多かった

ということから傾向スコアの似たサンプルと一致の具合が十分で無かった、つまり良いマッチング相手が無かったためかもしれません。傾向スコアマッチングでのATT推定は、対照群のデータが豊富で多様性のある時の方が良いかもしれないと思いました。

 また補足として、caliperを使った分析結果も確認してみましょう。caliperを設定すれば傾向スコアが任意の差以上となった場合、マッチング無しとできます。なので無理やりマッチングせずに、ある程度比較するのが適切なペアで効果推定できます。(matchitにcaliper=0.25を追加するだけなのでコードは割愛)。

[output]
# A tibble: 10 × 5
   term          estimate std.error statistic p.value
   <chr>            <dbl>     <dbl>     <dbl>   <dbl>
 1 (Intercept) -1716.      4417.      -0.389   0.698 
 2 treat        1887.       962.       1.96    0.0510
 3 age            -2.52      54.6     -0.0462  0.963 
 4 black        -613.      1335.      -0.459   0.647 
 5 hispanic     -987.      1956.      -0.504   0.615 
 6 married       591.      1257.       0.470   0.639 
 7 nodegree     1613.      1413.       1.14    0.255 
 8 education     531.       284.       1.87    0.0628
 9 re74            0.0809     0.120    0.675   0.500 
10 re75            0.143      0.186    0.772   0.441 

CPS3の傾向スコアマッチング後の傾向スコアとバランシングの可視化(caliper=0.25)
 介入により$1,887の収入アップとの効果推定ということで少し改善でしょうか。傾向スコアの分布としては介入群と対照群の分布は近づいており、共変量のバランシングも良好です。データをどんどん削っていくことにはなりますが、caliperは基本的に使用するのが良さそうです。

IPWでATE, ATTの推定

 それでは次にCPS1に対して、WeightItパッケージを用いてIPTWを使っていきます。またWeightItについても、得られる重みがglm()で得られるものと同等か検算してみます。

# CPS1に対して、IPW
weighting1 <- weightit(data=cps1_nsw_data,
formula = treat ~ age + black + hispanic + married + 
  nodegree + education + re74 + re75, 
  method="ps",
  estimand = "ATE")

# データフレームに得られた重みカラムを追加
cps1_nsw_data$weights <- weighting1$weights

# weightitで算出した重みがglm()で算出した傾向スコアと差があるか確認。
# 対照群の確認
print(nrow(filter(cps1_nsw_data, treat==0))
- sum(filter(cps1_nsw_data, treat==0)$ps + 1 / filter(cps1_nsw_data, treat==0)$weights))
# 介入群の確認
print(sum(filter(cps1_nsw_data, treat==1)$ps - 1 / filter(cps1_nsw_data, treat==1)$weights))

[output]
[1] 0
[1] -1.126529e-13

 検算結果はOKかと。

 次に傾向スコアマッチングの時と同様に、バランシングを確認します。またこちらも検算してみます。

# バランシングの評価
love.plot(weighting1, threshold=0.1)

CPS1のIPWによるバランシングの可視化

# IPTWで調整後のASAM
# 傾向スコアから算出される重みのカラムは作成しておく
standardized_mean_error_weighted <- function(df, weights, col, multiple=10){

  # 重みのぶんだけ増やすレコード数のカラム作成(今回は重みの10倍)
  df$weights_multiplied <- cps1_nsw_data[[weights]]*multiple

  # 介入群の調整後のベクトル作成
  df_1 <- filter(df, treat==1)
  weighted_values1 <- vector()
  for (i in 1:nrow(df_1)){
    repeated_value1 <- df_1[[col]][i]
    repeat_num1 <- round(df_1[["weights_multiplied"]][i])
    weighted_values1 <- append(weighted_values1, rep(repeated_value1, repeat_num1))
    }

  # 対照群の調整後のベクトル作成
  df_0 <- filter(df, treat==0)
  weighted_values0 <- vector()
  for (i in 1:nrow(df_0)){
    repeated_value0 <- df_0[[col]][i]
    repeat_num0 <- round(df_0[["weights_multiplied"]][i])
    weighted_values0 <- append(weighted_values0, rep(repeated_value0, repeat_num0))
    }

  # 各群のサンプルサイズ
  n1 <- length(weighted_values1)
  n0 <- length(weighted_values0)

  # 各群の平均値
  m1 <- mean(weighted_values1)
  m0 <- mean(weighted_values0)

  # 各群の不偏分散
  uv <- function(x){ var(x)*(length(x)-1)/length(x) }
  v1 <- uv(weighted_values1)
  v0 <- uv(weighted_values0)

  # 効果量Hedgesのg
  cov <- (v1*(n1-1) + v0*(n0-1)) / (n1+n0-2)
  g <- (m1-m0) / sqrt(cov)

  return (g)

}

# ASAMの検算
cols_li <- c("age", "black", "hispanic", "married", "nodegree", "education", "re74", "re75")
for (col in cols_li){ 
  print(paste(col, "のASAM"))
  print(paste("調整前: ", standardized_mean_error(cps1_nsw_data, col)))
  print(paste("調整後", standardized_mean_error_weighted(cps1_nsw_data, "weights", col)))
}

[output]
[1] "age のASAM"
[1] "調整前:  -0.673045077618078"
[1] "調整後 -0.645065994590819"
[1] "black のASAM"
[1] "調整前:  2.93324709421367"
[1] "調整後 0.544955431150178"
[1] "hispanic のASAM"
[1] "調整前:  -0.0486883897274181"
[1] "調整後 0.0970073665749104"
[1] "married のASAM"
[1] "調整前:  -1.1552824924244"
[1] "調整後 -0.475543303472403"
[1] "nodegree のASAM"
[1] "調整前:  0.903320101340552"
[1] "調整後 0.156113371826795"
[1] "education のASAM"
[1] "調整前:  -0.587471521697691"
[1] "調整後 -0.242858375535186"
[1] "re74 のASAM"
[1] "調整前:  -1.25103406168653"
[1] "調整後 -0.858869692212102"
[1] "re75 のASAM"
[1] "調整前:  -1.31388559386733"
[1] "調整後 -0.786715907260075"

 やはりlove.plotと検算結果が異なります。またlove.plotは調整前のASAMが、matchit()での結果と異なっています。ここは同じになるはずと考えていましたが、何か勘違いをしているんでしょうか。。。 ただ調整前後でも増減の挙動傾向は類似しているのでOKとします。
 結果としてはある程度バランスは取れている傾向ですが、十分ではありません。実際なら推定する意味無いと思いますが、せっかくなので確認してみましょう。

# 回帰分析で効果を推定。(解析モデルにも共変量を使用)
IPTW_result1 <- lm(data=cps1_nsw_data,
formula = re78 ~ treat + age + black + hispanic + married + 
  nodegree + education + re74 + re75,
  weight = weighting1$weights) %>%
    tidy()

print(IPTW_result1)

[output]
# A tibble: 10 × 5
   term         estimate std.error statistic   p.value
   <chr>           <dbl>     <dbl>     <dbl>     <dbl>
 1 (Intercept) 10596.     483.        21.9   3.46e-105
 2 treat       -1506.     133.       -11.3   1.64e- 29
 3 age           -97.1      6.46     -15.0   8.77e- 51
 4 black       -3164.     173.       -18.3   3.68e- 74
 5 hispanic      -78.7    209.        -0.377 7.06e-  1
 6 married     -1313.     136.        -9.68  4.27e- 22
 7 nodegree    -1812.     175.       -10.4   3.92e- 25
 8 education     -35.4     31.9       -1.11  2.67e-  1
 9 re74            0.472    0.0110    43.0   0        
10 re75            0.220    0.0108    20.4   2.54e- 91

 介入の推定効果がマイナスとなっていて、全然効果推定できてませんね。IPTWは介入グループと非介入グループの傾向が違う場合、分析結果が信頼しにくいことが知られているそうです。IPWアルゴリズムがNSWとCPSのデータが混ざったような状況での効果推定となるため、今回推定したかったNSWの実験を再現するような分析となっていないことが原因かと考えられます。データの背景を考えてみると対照群はCPSのデータであり収入が安定している個体が多いため、このような結果となっていると思われます。

 ただIPWでATTを推定する方法もあります。考え方は同じでATT推定の場合は、対照群の分布を介入群に合わせいくイメージで重みが変わるだけです。つまり介入群の重みが1、対照群の重みが(傾向スコア) / (1-傾向スコア)となるだけです(*3)。 実装は以下で、weightit()の引数であるestimandをATTにするだけです。

# cps1_nsw_dataで効果推定
# 対照群の分布を介入群に合わせにいくイメージ。なので重みは、介入群は1、対照群は(傾向スコア)/(1-(傾向スコア))。
weghiting1_att <- weightit(data=cps1_nsw_data,
formula = treat ~ age + black + hispanic + married + 
  nodegree + education + re74 + re75, 
  method="ps",
  estimand="ATT")

# 共変量のバランスを確認
love.plot(weghiting1_att, threshold=0.1)

IPW_result1_att <- lm(data=cps1_nsw_data,
formula=re78 ~ treat + age + black + hispanic + married + 
  nodegree + education + re74 + re75,
  weight=weghiting1_att$weights) %>%
    tidy()

print(IPW_result1_att)

[output]
# A tibble: 10 × 5
   term         estimate std.error statistic  p.value
   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
 1 (Intercept)  2026.     523.        3.88   1.06e- 4
 2 treat        1208.     107.       11.3    2.47e-29
 3 age           -29.2      6.52     -4.48   7.65e- 6
 4 black       -1481.     183.       -8.11   5.32e-16
 5 hispanic      -87.7    281.       -0.312  7.55e- 1
 6 married         5.46   153.        0.0357 9.71e- 1
 7 nodegree       22.5    169.        0.134  8.94e- 1
 8 education     422.      32.3      13.1    9.62e-39
 9 re74            0.101    0.0161    6.27   3.80e-10
10 re75            0.346    0.0232   14.9    1.03e-49

CPS1のIPW(ATT推定)によるバランシングの可視化
 バランシングはかなり良好で、推定結果もかなり改善しました。ただし傾向スコアマッチングと比べると推定結果はいまいちです。やはりデータの背景、算出したい効果に合わせて手法を選ぶことが大事だということでしょうか。データを重みで調整して比べるよりも傾向マッチングで実際に近い個体を比べる方が、実験の再現になるというのは直感的にはわかる気がします。

おわりに

今回は効果検証入門の傾向スコアの部分を整理していきました。

*1 効果検証入門の方で解説されています。

*2 こちらの書籍で紹介されています。 https://www.amazon.co.jp/%E7%B5%B1%E8%A8%88%E7%9A%84%E5%9B%A0%E6%9E%9C%E6%8E%A8%E8%AB%96%E3%81%AE%E7%90%86%E8%AB%96%E3%81%A8%E5%AE%9F%E8%A3%85-Wonderful-R-%E9%AB%98%E6%A9%8B-%E5%B0%86%E5%AE%9C/dp/4320112458www.amazon.co.jp

*3 こちらの記事がわかりやすかったです。 傾向スコアを用いた処置効果推定 #データサイエンス - Qiita

因果探索に入門 / ベイジアンネットワークによる因果探索

 因果分析に興味があるということでこちらの書籍を勉強中。
book.mynavi.jp
 こちらの書籍で紹介されているベイジアンネットワークによる因果探索の手法を整理していきます。また以下の記事も参考にさせていただきました。
統計的因果探索とAI
グラフィカルモデルに基づく因果探索手法の調査

目次

因果探索とベイジアンネットワーク

 因果探索とは、データ項目のつながり(因果グラフ)を推測する手法になります。
 因果探索の手法としてベイジアンネットワークがあり、基本的には
- 有向非巡回グラフ(矢印が巡回せず一方通行)である(DAG: Directed Acyclic Graphと呼ばれる)
- 未観測共通原因(潜在的な交絡因子)が無いこと
というように因果モデルに仮定を置いて
分析を進めるものになります。

ベイジアンネットワークによる因果探索の手法は以下に示すよう大きく3つあります。

[スコアリングによる構造学習]
 さまざまなDAGに対してスコア(※)を求め、最も指標値の良いDAGを選ぶものです。
つまりいろんなパターンのつながりをもったネットワークを考えて(規定して)、それぞれのスコアを算出して比較する手法になります。
(※)観測データに対するネットワークの当てはまりの良さを示すスコア:AIC, BICなど

[条件付き独立性検定による構造学習]
 変数間で独立性検定を各条件付きパターンで行い、スケルトン構造を同定します。その後ルールベースによって部分的な方向づけを繰り返して、スケルトン構造->PDAG(Partially DAG)->DAGと構築していく手法になります。

[ハイブリッド型構造型学習]
 スコアリングと条件付き独立性検定を組み合わせたもので、代表的な手法にMMHCアルゴリズム(Max-Min Hill Climbing algorithm)があります。

 今回は、
- 条件付き独立性検定による構造学習の代表的なアルゴリズムであるPCアルゴリズム
- 条件付き確率表(CPT)とベイズ情報量基準(BIC)を用いたスコアリングによる構造学習
を実装していきます。 また実装していくにあたりpythonライブラリとして主にpgmpyを使っていますが、書籍で紹介されているバージョンが0.1.9でありエラーが出ますので、0.1.21で実装しています(2023.05.13時点最新バージョンは0.1.22)。またv0..で開発段階のため、今後も大きな修正は行われそうなので使用する際は注意となります。

データの説明

 それでは、分析対象の擬似データの構造を示します。 今回も以前の因果推論の記事で用いた構造と同じもので行います。

職場での上司の研修受講が部下の面談満足度にどう影響するのか知りたいという場面を想定します。
双方とも上司の部下育成の熱心さと仕事のパフォーマンスにも影響を受けるとします。加えて、研修受講の効果は上司の部下育成の熱心さが高いと効果も高くなるとします。
( そこそこ妥当、そこそこ違和感というような場面設定ですが)

データの因果関係
今回は上記の様な因果(=矢印)のつながりと方向性を見つけることが目的となります。
(効果の大きさまでは分析しません。この部分に関しては、他の記事でも紹介している因果推論の出番となります。)

生成コードは以下です。

# データ数
num_data = 2000

# x1 : 上司の部下育成の熱心さ。-1から1の一様乱数。
x1 = np.random.uniform(low=-1, high=1, size=num_data)

# x2 : 上司の仕事のパフォーマンス。-1から1の一様乱数。
x2 = np.random.uniform(low=-1, high=1, size=num_data)
# x2 = np.random.choice([-1, -0.5, 0, 0.5, 1], num_data, p=[0.1, 0.2, 0.1, 0.3, 0.3])


# Z : 上司が「上司向け:部下とのキャリア面談のポイント研修」に参加したかどうか
e_z = randn(num_data)  # ノイズの生成
# シグモイド関数により受講するかどうかを熱心さに影響を受ける確率で表現(あくまで確率)
z_prob = expit(7*x1 + 5*x2 + 5*e_z)
Z = np.array([])
for i in range(num_data):
    Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0] # 
    Z = np.append(Z, Z_i)

# Y : 部下の面談の満足度
## 研修の効果は熱心さで効果が変わってくるとする(介入効果の非線形性)
t = np.zeros(num_data)
for i in range(num_data):
    if x1[i] < 0:
        t[i] = 0.5
    elif x1[i] >= 0 and x1[i] < 0.5:
        t[i] = 0.7
    elif x1[i] >= 0.5:
        t[i] = 1.0
e_y = randn(num_data)
Y1 = 2.0 + t*Z + 0.3*x1 + 0.5*x2 + 0.2*e_y

# テーブルを作成
df = pd.DataFrame({'x1': x1,
                   'x2': x2,
                   'Z': Z,
                   't': t,
                   'Y1': Y1
                   })

df.head()  # 先頭を表示
[output]

        x1        x2      Z      t       Y1
0  -0.616961  0.731925   0.0    0.5    2.112008
1  0.244218   0.424114   1.0    0.7    3.067455
2  -0.124545  0.327542   0.0    0.5    2.022895
3  0.570717   -0.230410  1.0    1.0    3.148009
4  0.559952   -0.641392  0.0    1.0    2.032487

分析パターン①(PCアルゴリズム)

 条件付き独立性検定による構造学習の代表的なアルゴリズムでもあり、制約ベースの因果探索とも位置付けられるものになります。

1. 変数間の独立性検定を各条件付きパターンで繰り返して、スケルトン構造を同定
2. オリエンテーションルールで部分的に方向づけを繰り返してPDAG(Partially DAG)をつくる
3. 2を繰り返してDAGを推定

 独立性検定とは統計的仮説検定の1つであり、背理法により"独立でない" or "独立でないとはいえない"を判定する手法です。 統計WEBの記事などが参考になるかと思います。 https://bellcurve.jp/statistics/course/9309.html https://bellcurve.jp/statistics/course/9496.html

 今回は、"独立でない""関連があるので、因果関係にある"という様に解釈して分析を進めていく手法となります。 また今回は検定統計量として、以下の数式で表されるカイ二乗統計量を用います。

 
χ ^ 2=\sum_{i=1}^r\sum_{i=1}^c\frac{(n _ {ij}-E _ {ij}) ^ 2}{E _ {ij}}

( r: 度数分布表の行数、 c: 度数分布表の列数、 n _ {ij}:  i行目 j列目の度数、 E _ {ij}:  i行目 j列目の推定期待度数)
また
 E _ {ij}=\frac{n _ i \times n _ j}{N}
( n _ i:  i行目のデータ数、 n _ j:  j行目のデータ数、 N: 全データ数)

データの離散化(ビン分割)

 PCアルゴリズムは独立性の検定を行うため連続値のデータは扱えないので、まず離散値に変換します。 今回は以下のように、pandasのcut関数で閾値での分割を行います。

df_bin = df.copy()
del df_bin["t"]

# x1:部下育成への熱心さ
df_bin["x1"], x_bins = pd.cut(df["x1"], 5, labels=[1, 2, 3, 4, 5], retbins=True)

# x2:仕事のパフォーマンス
df_bin["x2"], x_bins = pd.cut(df["x2"], 5, labels=[1, 2, 3, 4, 5], retbins=True)

# Z:上司が「上司向け:部下とのキャリア面談のポイント研修」に参加したかどうか
df_bin["Z"], Z_bins = pd.cut(df["Z"], 2, labels=["0", "1"], retbins=True)

# Y1:部下の面談の満足度
df_bin["Y1"], Y_bins = pd.cut(df["Y1"], 5, labels=[1, 2, 3, 4, 5], retbins=True)

# データ型変更
df_bin["x1"] = df_bin["x1"].astype("int")
df_bin["x2"] = df_bin["x2"].astype("int")
df_bin["Z"] = df_bin["Z"].astype("int")
df_bin["Y1"] = df_bin["Y1"].astype("int")

# 確認
df_bin.head()
[output]

   x1 x2 Z Y1
0  1  5  0  2
1  4  4  1  4
2  3  4  0  2
3  4  2  1  4
4  4  1  0  2

離散化できました。また各ビンのデータ量も確認しておきます。

for col_name in df_bin.columns:
  print(f"カラム名:{col_name}")
  display(df_bin[col_name].value_counts())
  print("=====")
カラム名:x1
1    452
5    414
3    410
4    381
2    343
Name: x1, dtype: int64
=====
カラム名:x2
4    432
3    409
1    400
5    393
2    366
Name: x2, dtype: int64
=====
カラム名:Z
1    1009
0     991
Name: Z, dtype: int64
=====
カラム名:Y1
2    693
3    545
4    469
1    189
5    104
Name: Y1, dtype: int64

閾値で分割したため、Y1ではデータ数に偏りが出ていますね。x1の場合分けによるtの値とrandn()の正規分布でのノイズの影響が出ている様です。

データの準備もできたので、PCアルゴリズムによる因果探索をしていきます。
ちなみに分析パターン①は以下のコードのみで結果が出ますが、中身を見ていくために今回は条件付き独立性検定とオリエンテーションルールの手順を一つずつ実行していきます。

# 以下のコードでPCアルゴリズムが実行できる
from pgmpy.estimators import PC
from pgmpy.base import DAG
from pgmpy.independencies import Independencies

est = PC(df_bin)
skel, seperating_sets = est.build_skeleton(ci_test="chi_square", max_cond_vars=6, significance_level=0.05)
pdag = est.skeleton_to_pdag(skel, seperating_sets)
model = pdag.to_dag()
print("DAG edges:        ", model.edges())

条件付き独立性検定

1. 0次の独立性検定
 まずは0次の独立性の検定です。0次なので、2つの変数間にある変数0個を条件付きにした際の検定ですのでそのまま行います。カイ二乗検定なので、以下のコードで一つづつ実行します。

# こちらは今回のパラメータ設定です。
from pgmpy.estimators.CITests import chi_square # 独立性の検定手法としてカイ二乗検定を選択

# パラメーターの設定
boolean_flg = True # Trueだと検定結果を返し、Falseだと検定統計量とp値を返す
level_value = 0.05 # 有意水準
# 4変数なので6通り
print(chi_square(X="x1", Y="Z", Z=[], data=df_bin, boolean= boolean_flg, significance_level=level_value))
print(chi_square(X="x1", Y="x2", Z=[], data=df_bin, boolean= boolean_flg, significance_level=level_value))
print(chi_square(X="x1", Y="Y1", Z=[], data=df_bin, boolean= boolean_flg, significance_level=level_value))
print(chi_square(X="Z", Y="x2", Z=[], data=df_bin, boolean= boolean_flg, significance_level=level_value))
print(chi_square(X="Z", Y="Y1", Z=[], data=df_bin, boolean=boolean_flg, significance_level=level_value))
print(chi_square(X="x2", Y="Y1", Z=[], data=df_bin, boolean=boolean_flg, significance_level=level_value))
[output]
False
True
False
False
False
False

x1 - x2間のみが独立でないとはいえないという結果なので、この部分だけエッジがない構造を以下の様に描くことができます。

分析パターン①: 0次の独立性検定の結果

2. 1次の独立性検定
 引き続いて1次の独立性の検定です。1次なので、2つの変数間にある変数1個を条件付きにした際の検定です。また条件付きにする変数は検定対象の変数と関連があるものから選びます。 例えば変数x1についての1次の独立性の検定は、以下のコードで実行します。

# 変数x1について: 2変数(Z, Y1)と関連があったので、1*2=2通り

print(chi_square(X="x1", Y="Z", Z=["Y1"], data=df_bin, 
                 boolean= boolean_flg, significance_level=level_value))
print(chi_square(X="x1", Y="Y1", Z=["Z"], data=df_bin, 
                 boolean= boolean_flg, significance_level=level_value))
[output]
True
False

x1とZが独立でないとはいえないという結果なので、x1 - Zのエッジが消去されて以下の様に構造が描けます。

分析パターン①: x1についての1次の独立性検定の結果

省略しますが他の変数についても検定をしていきます。結果としては、どの検定でも独立とは判定されなかったので1次の独立性検定の結果としては上記になります。

3. 2次の独立性検定
 続いて2次の独立性の検定です。1次の独立性検定と同様の説明ですが、2つの変数間にある変数2個を条件付きにした際の検定です。また条件付きにする変数は検定対象の変数と関連があるものから選びます。 2次ということから3変数以上と関連がある変数についてしか検定できないので、検定対象はY1のみとなります。以下のコードで実行します。

# 変数Y1について: 4変数(x1, Z, x2, Y2)なので、3*4=12通り
print(chi_square(X="Y1", Y="x1", Z=["Z", "x2"], data=df_bin, 
                 boolean= boolean_flg, significance_level=level_value))
print("=====")
print(chi_square(X="Y1", Y="Z", Z=["x1", "x2"], data=df_bin, 
                 boolean= boolean_flg, significance_level=level_value))
print("=====")
print(chi_square(X="Y1", Y="x2", Z=["x1", "Z"], data=df_bin, 
                 boolean= boolean_flg, significance_level=level_value))
[output]
False
=====
False
=====
False

いずれも独立でないという結果なので、2次の独立性検定で消去されるエッジはありませんでした。

オリエンテーションルールによる方向づけ

 独立性の検定で得られた構造がDAGのベースとなるスケルトンが得られました。ここから方向づけのためのオリエンテーションフェーズです。

1. V字構造に着目してスケルトンからPDAGへ変換
 まずはスケルトンをPDAGへ変換していくために、はじめにV字構造の部分に対して方向づけを行います。V字構造では下記の図の様に4パターンが考えられます。因果関係を絞るためにCI(A, C | B)を検定しますが、検定結果が独立でないと判定された場合のみ、パターン1の因果構造とできます。

V次構造に対する方向づけのパターン

ケルトンからx1-Y1-Zとx1-Y1-x2がV字構造と判断できるので、これらについて独立性の検定をしていきます。 以下のコードで実行します。

# v字構造での方向づけ

# x1-Y1-Z
print(chi_square(X="x1", Y="Z", Z=["Y1"], data=df_bin, 
                 boolean= boolean_flg, significance_level=level_value))

# x1-Y1-x2
print(chi_square(X="x1", Y="x2", Z=["Y1"], data=df_bin, 
                 boolean= boolean_flg, significance_level=level_value))
[output]
True
False

x1-Y1-Zは独立でないという結果なので、因果関係は以下のようになります。

分析パターン①: オリエンテーションルール-V字構造についてのスケルトンからPDAGへの変換

2. PDAGからの部分的な方向づけ
 できたPDAGに対して、以下の図を参考にさらにオリエンテーションルールに従った方向づけを行います。(成り立つ理由は理解できていませんが、PCアルゴリズムの原著に示されている様です。)

PDAGからの部分的な方向づけのパターン

上記のパターンを適用すると、Y1->Zの方向をつけることができ、さらにそこからx2->Zも方向付けることができます。

分析パターン②: PDAGからの方向づけの結果

 これで因果構造を描くことは出来ました。しかし真の構造とはかなり違ったものとなり、良い分析はできておりません。そもそも今回の場面設定ではZ: 研修の受講がY1: 部下の面談の満足度にどの程度影響しているかを知りたいというのが因果分析を行うモチベーションとなるはずです。ですので、Y1 -> Zという因果構造になったことからも分析結果に対して懸念を持つことになるかと思います。 そこで分析の方法を変更してみます。

分析パターン②(条件付き独立性検定 -> スコアリングによる構造学習)

 それでは分析の改善のため、データの離散化の方法を変更します。
分析パターン①では、離散化を閾値で等間隔になる様に分割を行なったため、ビン間にデータ量の偏りがありました。PCアルゴリズムは独立性の検定により、エッジの有無判定、スケルトンからPDAGへの変換を行うので離散化の方法による影響、つまりデータの偏りや量は注目するべき点になります。そこで分析パターン②では、データ量が偏らない様に分割していきます。 データの生成は同じなので、説明は割愛します。

ビンの分割方法の変更

以下のコードで実行します。pandasのqcut関数を用います。

df_bin = df.copy()
del df_bin["t"]

# x1:部下育成への熱心さ
df_bin["x1"], x_bins = pd.qcut(df["x1"], 5, labels=[1, 2, 3, 4, 5], retbins=True)

# x2:仕事のパフォーマンス
df_bin["x2"], x_bins = pd.qcut(df["x2"], 5, labels=[1, 2, 3, 4, 5], retbins=True)

# Z:上司が「上司向け:部下とのキャリア面談のポイント研修」に参加したかどうか
df_bin["Z"], Z_bins = pd.cut(df["Z"], 2, labels=["0", "1"], retbins=True)

# Y:部下の面談の満足度
df_bin["Y1"], Y_bins = pd.qcut(df["Y1"], 5, labels=[1, 2, 3, 4, 5], retbins=True)

# データ型変更。
df_bin["x1"] = df_bin["x1"].astype("int")
df_bin["x2"] = df_bin["x2"].astype("int")
df_bin["Z"] = df_bin["Z"].astype("int")
df_bin["Y1"] = df_bin["Y1"].astype("int")

# 確認
df_bin.head()
[output]
   x1  x2  Z  Y1
0  2   5   0  3
1  4   4   1  5
2  3   4   0  2
3  4   2   1  5
4  4   1   0  2
for col_name in df_bin.columns:
  print(f"カラム名:{col_name}")
  display(df_bin[col_name].value_counts())
  print("=====")
[output]
カラム名:x1
2    400
4    400
3    400
5    400
1    400
Name: x1, dtype: int64
=====
カラム名:x2
5    400
4    400
2    400
1    400
3    400
Name: x2, dtype: int64
=====
カラム名:Z
1    1009
0     991
Name: Z, dtype: int64
=====
カラム名:Y1
3    400
5    400
2    400
4    400
1    400
Name: Y1, dtype: int64
=====

偏りなくビン分割できています。それではこちらのデータで同様に条件付き独立性の検定を行なっています。

条件付き独立性検定

 こちらも分析パターン①と同様に、0次、1次、2次と条件付き独立性の検定を行なっていくだけなので割愛します。 ただ得られた構造は分析パターン①と異なり、下記のように真の構造と同じエッジが残りました。 結果の流れとしては、0次で得られた構造からエッジが無くなることなく2次までの検定が完了しました。

分析パターン②: 0次、1次、2次の独立性検定で得られたスケルト

ただこの構造ではV字構造が無いためPDAGへの変換ができず、オリエンテーションルールが適用できません。そこでスコアリングによる構造学習で方向づけを考えていきます。

スコアリングによる構造学習

 スコアリングによる構造学習は、

(再喝)
さまざまなDAGに対してスコア(※)を求め、最も指標値の良いDAGを選ぶ。
(※)観測データに対するネットワークの当てはまりの良さを示すスコア:AIC, BICなど

という手法です。今回はBICを指標とします。
手順としては、

1. あるモデル(= 因果構造)を規定する
2. 観測データから条件付き確率表(CPT)を推定する
3. 得られた推定したCPTから各変数の各サンプルの確率値とパラメーター数が得られる
4. 各モデルのBICが求められるので比較して、最も良好な値のモデルを採用する

となります。
 はじめにモデルを規定して指標を比較するので、変数が多くなるとDAGのパターンが爆発的に増加してしまうのが欠点になります(そのため厳密性を犠牲にして計算量を減らす手法もいくつか提案されている様です)。

BICについて
 BICが以下の式で計算されます。
 BIC _ {m} = -2l _ {m} (θ|X) + k _ {m} (logN)
mはモデルを表しており、対象としているベイジアンネットワークのDAGになります。 Xは観測されたデータ、 k _ {m}はモデルのパラメータ数で自由度に対応する様な値、 Nは観測されたデータ数になります。
また l _ {m} (θ|X)はそのモデルと規定した場合の確率値から求められる対数尤度であり、比例する値として以下の様にも計算されます。

 l _ {m} (θ|X) \propto \sum_{i}\sum_{j}\sum_{k}(N _ {ijk})log\hat{θ} _ {i,j,k}

 iは変数、 jは変数iについてのとある条件パターン、 kは変数iのとる値)
負の対数尤度に、過剰な適合を考慮するためのペナルティとしてパラメータ数とデータ数を掛け合わせたものを足しているので、値が小さいほど良いモデルと考える指標になります。

 それでは分析を進めるために、まず因果構造を考えます。今回はZ: 研修受講がY1: 部下の面談の満足度にどう影響するかを分析することが目標になります。ですので、Z -> Y1という因果方向を前提としていくつかのモデルを考えます。DAGであることから循環しないモデルを考えると、以下の9パターンが考えられます。

因果構造の候補のパターン

以下のコードで各モデル(因果構造)におけるBICを計算します。pgmpyでのBICの計算値は上記で示したものに-0.5をかけたものになっているので、値が大きいほど良いモデルと考えます。

# スコアリング構造学習で判断する
from pgmpy.models import BayesianNetwork
from pgmpy.estimators import BicScore

model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('x1', 'Y1'), ('Z', 'x2'), ('x2', 'Y1')])  
print(BicScore(df_bin).score(model))

model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('x1', 'Y1'), ('Z', 'x2'), ('Y1', 'x2')])  
print(BicScore(df_bin).score(model))

model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('Y1', 'x1'), ('x2', 'Z'), ('x2', 'Y1')])  
print(BicScore(df_bin).score(model))

model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('Y1', 'x1'), ('Z', 'x2'), ('x2', 'Y1')])  
print(BicScore(df_bin).score(model))

model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('Y1', 'x1'), ('Z', 'x2'), ('Y1', 'x2')])  
print(BicScore(df_bin).score(model))

model = BayesianNetwork([('Z', 'Y1'), ('Z', 'x1'), ('Y1', 'x1'), ('x2', 'Z'), ('x2', 'Y1')])  
print(BicScore(df_bin).score(model))

model = BayesianNetwork([('Z', 'Y1'), ('x1', 'Z'), ('x1', 'Y1'), ('Z', 'x2'), ('x2', 'Y1')])  
print(BicScore(df_bin).score(model))

model = BayesianNetwork([('Z', 'Y1'), ('x1', 'Z'), ('x1', 'Y1'), ('Z', 'x2'), ('Y1', 'x2')])  
print(BicScore(df_bin).score(model))

model = BayesianNetwork([('Z', 'Y1'), ('x1', 'Z'), ('x1', 'Y1'), ('x2', 'Z'), ('x2', 'Y1')])  
print(BicScore(df_bin).score(model))
[output]
-9652.730917876488
-9468.571910366429
-9468.571910366429
-9468.571910366429
-9468.571910366429
-9468.571910366429
-9652.730917876488
-9468.571910366429
-9675.163065400724

出力結果より、真のモデルは一番値が小さく分析結果としては一番当てはまらないモデルという結果になってしまいました。。。。。

おわりに

 今回は因果探索への入門として、ベイジアンネットワークを用いた手法を実装していきました。 手法の原理もある程度とばしながらですが一通りの実装はできたと思います。データの分割方法の変更、PCアルゴリズムとスコアリングによる構造学習との併用も実験的に行えて勉強になりました。ただ4変数しかないにもかかわらず真の構造を推定できず、因果探索自体が難しいタスクという感触でした。データからの結果だけで無く、ドメイン知識による因果関係の推定も重要になってきそうです。また手法などは色々と提案されている様ですが、データの分割方法などの前処理の方法論は特に無いようですので、試しながら肌感覚を養っていけたらと思います。

(追記 注意点)
 便利に思われるPCアルゴリズムですが、検定する変数の順番が異なればエッジが変わることもあります。これについては私も実験して確認し、著者の方に問い合わせたところ(以下のリンク)事実の様です。またこれに対して、

どうするのが正解なのか、きちんと解説している文献を見たことがありません。

とのことです。これは非常に大事な注意点かと思います。
https://github.com/YutaroOgawa/causal_book/issues/39

因果推論に入門( 処置効果が非線形なデータで ) / Meta-Learnersの実装・比較

因果分析に興味があるということでこちらの書籍を勉強中。

book.mynavi.jp

こちらの書籍で紹介されている因果推論の手法を整理していきます。
処置効果が線形な場合は一度整理したので、今回は処置効果が非線形場合に、適用可能な手法の整理となります。

目次

Meta-Learnersとは

以前の記事で整理した手法である回帰分析、IPTW法、DR法は因果関係にある変数間の処置効果が線形の場合を想定していました。
それに対してMeta-Learnersとは非線形な因果関係および処置効果を扱える因果推論手法であり、いわゆるCATE(条件付き処置効果)を推定する手法となります。
ただし樹木モデル( ランダムフォレスト、勾配ブースティング木など )やニューラルネットワークといった分析モデルを指すのでは無く、処置効果を求めるために"分析モデルを用いた解き方・使い方"のようなものとなります。
ちなみに今回は分析モデルにはランダムフォレストを用います。

データの説明

以下の構造をもつデータで因果推論を行います。
職場での上司の研修受講が部下の面談満足度に影響するという場面を想定します。
また双方とも上司の部下育成の熱心さと仕事のパフォーマンスにも影響を受けるとします。
加えて、研修受講の効果は上司の部下育成の熱心さが高いと効果も高くなるとします。
( そこそこ妥当、そこそこ違和感というような場面設定ですが気にせず )

データの因果関係

とはいっても、これだけではイメージしきれない部分もあると思うので、
データの生成コードと可視化結果を載せます。
非線形な処置効果を表しているのは、部下の面談の満足度のコードになります。
熱心さの場合分けをif文で処理し、研修受講したかどうかにかかる係数を指定しています。
やっていることは、熱心さが-1 から 1の値を取る中、
-1以上0未満では0.5、
0以上0.5未満では0.7、
0.5以上1以下では1
の処置効果となるようにしています。
ですのでおおよそ平均的には0.675の処置効果となります。

import random
import numpy as np

# 乱数のシードを固定
random.seed(9999)
np.random.seed(9999)

# Scipy 平均0、分散1に正規化(標準化)関数
import scipy.stats

# シグモイド関数をimport
from scipy.special import expit

# 標準正規分布の生成用
from numpy.random import *

# グラフの描画用
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns

# その他
import pandas as pd

# データ数
num_data = 500

# enthusiam: 上司の部下育成の熱心さ。-1から1の一様乱数。
x1 = np.random.uniform(low=-1, high=1, size=num_data)

# performance: 上司の仕事のパフォーマンス。-1から1の一様乱数。
x2 = np.random.uniform(low=-1, high=1, size=num_data)
# x2 = np.random.choice([-1, -0.5, 0, 0.5, 1], num_data, p=[0.1, 0.2, 0.1, 0.3, 0.3])

# train_attend: 上司が「上司向け:部下とのキャリア面談のポイント研修」に参加したかどうか
e_z = randn(num_data)  # ノイズの生成
# シグモイド関数により受講するかどうかを熱心さに影響を受ける確率で表現(あくまで確率)
z_prob = expit(5*x1 + 5*x2 + 5*e_z)
Z = np.array([])
for i in range(num_data):
    Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0] # 
    Z = np.append(Z, Z_i)

# stsf_level: 部下の面談の満足度
## 研修の効果は熱心さで効果が変わってくるとする(介入効果の非線形性)
t = np.zeros(num_data)
for i in range(num_data):
    if x1[i] < 0:
        t[i] = 0.5
    elif x1[i] >= 0 and x1[i] < 0.5:
        t[i] = 0.7
    elif x1[i] >= 0.5:
        t[i] = 1.0
e_y = randn(num_data)
Y = 2.0 + t*Z + 0.3*x1 + 0.3*x2 + 0.2*e_y

# テーブルを作成
df_B = pd.DataFrame({'enthusiam': x1,
                   'performance': x2,
                   'train_attend': Z,
                   'effect': t,
                   'stsf_level': Y
                   })

df_B.head()  # 先頭を表示
[output]
# データの先頭行
    enthusiam  performance  train_attend  effect  stsf_level
0   0.646778   -0.201372        1.0         1.0    2.968777
1   -0.562795  0.574169         1.0         0.5    3.018853
2   -0.937078  0.979220         1.0         0.5    2.252320
3   0.091852   -0.512450        0.0         0.7    2.098059
4   -0.691089  -0.040364        0.0         0.5    1.863024

データの可視化

データをseabornのpairplotで簡単に可視化します。
満足度(sfsf_level)に対して、上司の熱心さ(enthusiam)、上司のパフォーマンス(performance)、研修を受講したかどうか(train_attend)は正の相関が見られます。
また処置効果(effect、受講->満足度)は熱心さに影響を受けて階段状になっており、非線形な影響を受けていることがわかります。ただし勿論ですが、この部分の可視化は現実では確認することはできません。

sns.pairplot(df_B)

各変数間の可視化

また上司が研修を受講したかどうかでの部下の面談の満足度のヒストグラムも確認します。 研修の受講自体は良い影響が出ていそうとわかります。
しかし各分布の平均値をみると、満足度の差が1.0近くもあります(全体平均では0.675程度のはず)。
部下育成に熱心な上司は、そもそも研修を受講しやすいというセレクションバイアス(比較している集団の潜在的な傾向の違いによるバイアス)がかかっています(となるようなデータの生成をしています)。
また当然ですがこれらの可視化は、単変量や2変量を確認するものなので、他の変量による非線形な影響にも気付けません。

# 研修を受けたかどうかで購入量のヒストグラムを比較
df_B[df_B["train_attend"]==1]["stsf_level"].hist(color="r", alpha=0.3, label="train_attend=1")
df_B[df_B["train_attend"]==0]["stsf_level"].hist(color="b", alpha=0.3, label="train_attend=0")
plt.show()

上司が研修を受けたかどうか別の部下の面談の満足度のヒストグラム

[output]
enthusiam       0.130303
performance     0.205564
train_attend    1.000000
effect          0.715789
stsf_level      2.838375
dtype: float64
======
enthusiam      -0.223002
performance    -0.167072
train_attend    0.000000
effect          0.597233
stsf_level      1.883290
dtype: float64

傾向スコアによるIPTWでわかること、わからないこと

それではまず傾向スコアによるIPTW(逆確率重み付け法)です。 以前の記事でも整理したので割愛しますが、介入有り・無しそれぞれの期待値を傾向スコアの調整・推定し、その差分を処置効果とするものです。

以下で、実行します。

from sklearn.linear_model import LogisticRegression

# 傾向スコアを求めるために、原因変数を目的変数としてモデリング

# 説明変数
X = df_B[["enthusiam", "performance"]]

# 目的変数(原因変数)
Z = df_B["train_attend"]

# 学習(回帰の実施)
model = LogisticRegression()
model.fit(X,Z)

# 平均処置効果 ATEを計算
Y = df_B["stsf_level"]
ATE_i = Y / Z_pred[:, 1] * Z - Y / Z_pred[:, 0] * (1-Z)
ATE = 1 / len(Y) * ATE_i.sum()
print("推定したATE", ATE)
[output]
推定したATE 0.687

推定したATE(平均処置効果)は0.687と結構良い推定結果になっており、セレクションバイアスも考慮できていそうです。
しかしIPTWでは全体の平均的な処置効果を推定できるだけであり、 非線形な影響を確認できません。

Meta-Learnersの実装

そこでここからがMeta-Learnersの実装になります。
今回は4つの手法を比較していきます。
再喝ですが、分析モデルにはランダムフォレストを用いていきます。

T-Learner

T-Learnerは、原因変数をモデル構築には使用せず、集団ごとの分割の基準のみに使います。
Tは2=Twoの略らしく、2つのモデルを構築します。
手順としては以下です。
1. まず介入有り無しでデータを分ける
2. それぞれのデータで、結果変数と共変量(原因変数(=介入)以外)でモデル構築
3. それぞれのモデルで、介入有り無し分割前の共変量(原因変数(=介入)以外)で結果変数の推定
 介入有りの推定値:mu_1
 介入無しの推定値:mu_2
4. 各サンプルにおける各モデルの推定値の差の平均がATE
サンプルごとの推定処置効果(各サンプルにおける各モデルの推定値の差) = mu_1(i) - mu_2(i)
ATE(平均) = (サンプルごとの推定処置効果).mean()

また以下のようにATT(処置群における平均処置効果)とATU(対象群における平均処置効果)も推定できます。
ATT = ( 介入有り結果変数 - 介入無しモデル(介入有り共変量) ).mean()
ATU = ( 介入有りモデル(介入無し共変量) - 介入無し結果変数).mean()

それでは以下のコードで、T-Learnerを実装し、推定処置効果の可視化まで行います。
model0とmodel1という様に、2つのモデルを構築しています。

# ライブラリをインポート
from sklearn.ensemble import RandomForestRegressor

# 介入の有りと無しでデータを分割
df_0 = df_B[df_B["train_attend"]==0]
df_1 = df_B[df_B["train_attend"]==1]

# それぞれでモデル構築(目的変数: 結果変数、説明変数: 共変量)
model0 = RandomForestRegressor(max_depth=4)
M_0 = model0.fit( df_0[["enthusiam", "performance"]], df_0["stsf_level"]) # 介入を受けてないモデル

model1 = RandomForestRegressor(max_depth=4)
M_1 = model1.fit( df_1[["enthusiam", "performance"]], df_1["stsf_level"]) # 介入を受けたモデル

# 各モデルで分割前の全データで結果変数を推定
mu_0 = M_0.predict(df_B[["enthusiam", "performance"]])
mu_1 = M_1.predict(df_B[["enthusiam", "performance"]])

# サンプルごとの処置効果の推定とATEの算出
ATE_i = mu_1 - mu_0
ATE = ATE_i.mean()
print(f"ATE: {ATE}")

# 推定した処置効果を描画
t_estimated = ATE_i
plt.scatter(df_B["enthusiam"], df_B["effect"], c="k", alpha=1, label="true_effect") # 正解を描画
plt.scatter(df_B["enthusiam"], t_estimated, c="r", alpha=0.3, label="estimated_effect") # 推定値を描画
plt.legend()
plt.show()
[output]
ATE: 0.6874744694414802

推定したATEは0.687とIPTWと同様の性能を示しました。
ただMeta-Learnersはサンプルごとの処置効果も推定でき、以下が可視化結果になります。
綺麗に正解に沿っている訳ではありませんが、非線形性を表現できており、熱心さが非線形に効いていることは充分に判断できます。

T-Learnerによる熱心さに対する推定処置効果(研修受講->満足度)

またもう一つの共変量である上司のパフォーマンスに対する処置効果の可視化です。
こちらは比較的一様な分布になっていることから、同じ共変量であっても熱心さの大小が処置効果に大きく影響しており、パフォーマンスはさほど(もしくは全く)影響していないことが判断できます。

T-Learnerによるパフォーマンスに対する推定処置効果(研修受講->満足度)

S-Learner

T-LearnerだけでもMeta-Learnersのメリットを感じられたと思います。
ここからは他の手法についても眺めるように見ていきます。

まずS-Learnerです。SはSingleの頭文字の様で、T-Learnerに対する感じですね。
モデルが1つということで、T-Learnerと異なり、原因変数をモデル構築に使用します。
手順としては以下です。
1. 結果変数と共変量+原因変数の全データでモデル構築
2. 介入が全て有り or 無しとして置き換えたデータセットを作成
  X_0 : 介入を全て無しに置き換えた原因変数+共変量
  X_1 : 介入を全て有りに置き換えた原因変数+共変量
3. サンプルごとの推定処置効果 = モデル(X_1) - モデル(X_0)
4. ATE = (サンプルごとの推定処置効果).mean()

以下で可視化まで実装します。

# ライブラリをインポート
from sklearn.ensemble import RandomForestRegressor

# モデルを構築
model = RandomForestRegressor(max_depth=4)
X = df_B[["enthusiam", "performance", "train_attend"]]
model.fit(X, df_B["stsf_level"])

# 目的変数を各処置の場合を1つのモデルから求めて、ATEを計算
X_0 = X.copy()
X_0["train_attend"] = 0 

X_1 = X.copy()
X_1["train_attend"] = 1

ATE_i = model.predict(X_1) - model.predict(X_0)
ATE = ATE_i.mean()
print(f"ATE: {ATE}")

# 推定した処置効果を描画
t_estimated = ATE_i
plt.scatter(df_B["enthusiam"], df_B["effect"], c="k", alpha=0.3, label="true_effect") # 正解を描画
plt.scatter(df_B["enthusiam"], t_estimated, c="r", alpha=0.3, label="estimated_effect") # 推定値を描画
[output]
ATE: 0.706241672477917

推定したATEは0.706と性能は若干低い結果でした。ただ良し悪しに再現性は無さそうな程度感かなとも思います。
以下に可視化結果です。T-Learnerとの差はあまりわかりませんが、S-Learnerも非線形性を表現できており、熱心さが非線形に効いていることは充分に判断できるかと思います。

S-Learnerによる熱心さに対する推定処置効果(研修受講->満足度)

X-Learner

X-Learnerは、T-Learnerの結果を傾向スコアを用いて補正する手法となります。
気持ちとしては、自分のデータを使わずに構築されたモデルでの推定結果に重きを置くことで、トータルとしての推定精度を上げたいというものの様です(そこにさらに傾向スコアでの補正をしているイメージかなと認識してます)。
手順としては以下です。
1. まず介入有り無しでデータを分ける
2. それぞれのデータセットで、結果変数と共変量(原因変数(=介入)以外)でモデル構築
3. それぞれのデータセットとモデル推定値の各サンプルにおいて、ITTとITUを求める
 ITU = 介入有りモデル(介入無し共変量) - 介入無し結果変数
 ITT = 介入有り結果変数 - 介入無しモデル(介入有り共変量)
  (ATU = ITU.mean(),、ATT = ITT.mean())
4. 各データセットの共変量を説明変数、上記ITTおよひITUを目的変数として、モデル構築
   ITUモデル = model.fit(介入無し共変量, ITU)
   ITTモデル = model.fit(介入有り共変量, ITT)
5. 上記の各モデルを全データの共変量を使って、ITT(hat)とITE(hat)を推定
   ITU(hat) = ITUモデル(全データの共変量)
   ITT(hat) = ITTモデル(全データの共変量)
6. 原因変数と共変量で傾向スコア( g(x) : 介入有りの確率)算出
7. 各効果の推定値を傾向スコアで調整
 (サンプルごとの推定した処置効果) = g(x) * ITU(hat) + ( 1 - g(x) ) * ITT(hat)
 ATE = (サンプルごとの推定した処置効果).mean()

以下で可視化まで実装します。
大まかには前半はT-Learnerを行い、そこから傾向スコアの算出、処置効果の補正という流れです。

# ライブラリをインポート
from sklearn.ensemble import RandomForestRegressor

# 介入の有りと無しでデータを分割
df_0 = df_B[df_B["train_attend"]==0]
df_1 = df_B[df_B["train_attend"]==1]

# それぞれでモデル構築(目的変数: 結果変数、説明変数: 共変量)
model0 = RandomForestRegressor(max_depth=4)
M_0 = model0.fit( df_0[["enthusiam", "performance"]], df_0["stsf_level"]) # 介入を受けてないモデル

model1 = RandomForestRegressor(max_depth=4)
M_1 = model1.fit( df_1[["enthusiam", "performance"]], df_1["stsf_level"]) # 介入を受けたモデル

# ITUとITTを求める
tau_0 = M_1.predict(df_0[["enthusiam", "performance"]]) - df_0["stsf_level"]
tau_1 = df_1["stsf_level"] - M_0.predict(df_1[["enthusiam", "performance"]])

# ITUとITTの推定モデルを構築
model_ITU = RandomForestRegressor(max_depth=4)
M_ITU = model_ITU.fit(df_0[["enthusiam", "performance"]], tau_0)

model_ITT = RandomForestRegressor(max_depth=4)
M_ITT = model_ITT.fit(df_1[["enthusiam", "performance"]], tau_1)

# 傾向スコアを求める
from sklearn.linear_model import LogisticRegression

# 説明変数: 共変量
X = df_B[["enthusiam", "performance"]]

# 目的変数: 結果変数
Z = df_B["train_attend"]

# モデルを構築
model = LogisticRegression()
model.fit(X, Z)

# 傾向スコア算出
Z_proba = model.predict_proba(X)

# それぞれのモデルで全データの処置効果を推定し、傾向スコアで調整
ITU_pred = M_ITU.predict(df_B[["enthusiam", "performance"]])
ITT_pred = M_ITT.predict(df_B[["enthusiam", "performance"]])
ATE_i = Z_proba[:, 1]*ITU_pred + Z_proba[:, 0]*ITT_pred
ATE = ATE_i.mean()
print(f"ATE: {ATE}")

# 推定した処置効果を描画
t_estimated = ATE_i
plt.scatter(df_B["enthusiam"], df_B["effect"], c="k", alpha=0.3, label="estimated_TE") # 正解を描画
plt.scatter(df_B["enthusiam"], t_estimated, c="r", alpha=0.3, label="estimated_TE") # 推定値を描画
[output]
ATE: 0.6757366813565363

推定したATEは0.676と良好な結果でした。
以下に可視化結果です。他の手法と比べて振れ幅も狭くなり、より良好な推定をしている様に見えます。

X-Learnerによる熱心さに対する推定処置効果(研修受講->満足度)
傾向スコアでの補正前の推定値とも比べました。傾向スコアによる補正前について、推定ITUは熱心さが大きい範囲で性能が良く・分散が小さい、推定ITTは熱心さが小さい範囲で性能が良く・分散が小さくなっています。そして傾向スコアでの補正で良いとこどりができていることが確認できます。
(ただなぜこうなるかは考察しきれませんでした。処置効果が大きい範囲は反実仮想の値を推定する際の精度によっては分散が大きくなるイメージはつくんですが、、、)
傾向スコアによる補正前後_X-Learnerによる熱心さに対する推定処置効果(研修受講->満足度)

DR-Learner

DR-Learnerは、傾向スコアを反実仮想の推定に利用し、より性能の良い機械学習モデルを構築する手法の様です(正直腹落ちしておらず、結果変数の調整に傾向スコアと反実仮想を用いているイメージです)。

X-Learnerとの違いは、傾向スコアの使い方が異なります。

X-Learnerは、モデルはよくある流れで構築して、推定したITTとITUに傾向スコアを用いて効果を推定しています。
1. T-Learner的に算出したITTおよびITUを目的変数に、共変量を説明変数としてモデル構築
2. そのモデルに再度共変量を入力して、(サンプルごとの推定した処置効果)を推定
3. 推定した処置効果を介入有り無しで逆側の傾向スコアで加重平均して、介入効果を求める

DR-Learnerは、モデル構築に用いる結果変数の調整に用い、その調整された結果変数をもとにしたITTとITUを各サンプルのITEとして、モデル構築および推定を行っています。
1. T-Learner的にITTとITUを算出する際に、結果変数側を傾向スコアと推定結果変数で調整して算出する(精度良く推定しながら、推定誤差も調整することで、個別の誤差を少しでも緩和)
2. 算出したITTとITUを目的変数に、共変量を説明変数にしてモデル構築
3. そのモデルに再度共変量を入力して、推定処置効果を求める

以下が手順となります。
1. T-learner的に機械学習モデルを構築する
2. 結果変数を傾向スコアと1の機械学習モデルによる推定した結果変数で調整する
3. 潜在的結果変数も1の機械学習モデルから求め、介入有り無しでそれぞれのITEを求める
4. 3で求められたITEと共変量でモデル構築
5. 共変量で介入効果を推定する

# ライブラリをインポート
from sklearn.ensemble import RandomForestRegressor

# T-Learnerで処置のケースごとにモデル構築
# 介入の有りと無しでデータを分割
df_0 = df_B[df_B["train_attend"]==0]
df_1 = df_B[df_B["train_attend"]==1]

# それぞれでモデル構築(目的変数: 結果変数、説明変数: 共変量)
model0 = RandomForestRegressor(max_depth=4)
M_0 = model0.fit( df_0[["enthusiam", "performance"]], df_0["stsf_level"]) # 介入を受けてないモデル

model1 = RandomForestRegressor(max_depth=4)
M_1 = model1.fit( df_1[["enthusiam", "performance"]], df_1["stsf_level"]) # 介入を受けたモデル

# 傾向スコアを求めるモデル構築
from sklearn.linear_model import LogisticRegression

X = df_B[["enthusiam", "performance"]]
Z = df_B["train_attend"]
model_ps = LogisticRegression()
model_ps.fit(X, Z)

g_x = LogisticRegression().fit(X, Z)
g_x_val = g_x.predict_proba(X)

# DR法に基づく推定
# 処置群について
Y_1 = (
    M_1.predict(df_1[["enthusiam", "performance"]])
    + (df_1["stsf_level"] - M_1.predict(df_1[["enthusiam", "performance"]]))
    /g_x.predict_proba(df_1[["enthusiam", "performance"]])[:, 1]
    )
df_1["ITE"] = Y_1 - M_0.predict(df_1[["enthusiam", "performance"]])

# 非処置群について
Y_0 = (
    M_0.predict(df_0[["enthusiam", "performance"]]) 
    + (df_0["stsf_level"] - M_0.predict(df_0[["enthusiam", "performance"]]))
    /g_x.predict_proba(df_0[["enthusiam", "performance"]])[:, 0]
)
df_0["ITE"] = M_1.predict(df_0[["enthusiam", "performance"]]) - Y_0

 # データを結合
df_DR = pd.concat([df_0, df_1])
display(df_DR.head())

# ITEに対してモデル構築し、サンプルごとの処置効果を推定する
model_DR = RandomForestRegressor(max_depth=4)
M_DR = model_DR.fit(df_DR[["enthusiam", "performance"]], df_DR["ITE"])
ATE_i = M_DR.predict(df_DR[["enthusiam", "performance"]])
ATE = ATE_i.mean()
print(f"ATE: {ATE}")

# 推定した処置効果を描画
t_estimated = ATE_i
plt.scatter(df_DR["enthusiam"], df_DR["effect"], c="k", alpha=1, label="true_effect") # 正解を描画
plt.scatter(df_DR["enthusiam"], t_estimated, c="r", alpha=0.3, label="estimated_effect") # 推定値を描画
plt.legend()
plt.show()
[output]
ATE: 0.6705706151150326

推定したATEは0.671と良好な結果でした。 以下に可視化結果です。大きく外れてしまうサンプルもある様な結果でした。
y軸を拡大した可視化も見てみると、他手法と比べて各熱心さの範囲でx軸に並行な推定結果となっており、処置効果の傾向にはより近づこうとしている印象です。

DR-Learnerによる熱心さに対する推定処置効果(研修受講-&gt;満足度)
y軸拡大_DR-Learnerによる熱心さに対する推定処置効果(研修受講->満足度)

モデルのハイパラ調整(X-Learner、木の最大深さと葉ノードに属する最小サンプル数)

Meta-Learnerは使用する分析モデルは何でも良いので、実際はモデル比較やハイパーパラメータ調整を行うのだろうかと思います。
ということでハイパーパラーメータ調整を行い、より精度の高い機械学習モデルを用いると結果がどう変わるのか確認しました。
- 用いる手法: X-Learner
- 調整するハイパーパラメータ: 木の最大深さ、葉ノードに属する最小サンプル数
とします。
ハイパーパラメータ調整を行う箇所は、T-Learnerでの2つのモデル構築とITUとITTの推定モデルの構築の部分です。
またハイパーパラメータの最適化は、評価指標をMAE(平均絶対誤差)とした交差検証による平均値比較でのグリッドサーチを行います。

以下に実装コードです。
交差検証とグリッドサーチは多用するので関数を用意しておきます。

from sklearn.model_selection import KFold
from sklearn.metrics import mean_absolute_error

# 交差検証の評価指標の結果を返す。評価指標はMAE(絶対平均誤差)で行う。
def cv_mse(X, y, model, split_num = 5, shuffle = True):
  cv = KFold(n_splits=split_num, random_state=0, shuffle=shuffle)
  mae_list = []
  for train_index, test_index in cv.split(X):
      X_train, X_test = X.iloc[train_index, :], X.iloc[test_index, :]
      y_train, y_test = y.iloc[train_index], y.iloc[test_index]

      model.fit(X_train, y_train)
      y_pred = model.predict(X_test)

      mae = np.round(mean_absolute_error(y_test, y_pred), 3)
      mae_list.append(mae)

  return mae_list

# 交差検証とグリッドサーチをして、最も良いパラメータを返す
def rf_cv_grid_search(X, y, depth=range(1, 9), leaf=range(1, 11)): 
  best_val_score = np.inf  # best_scoreを初期化
  best_params = {}   # best_paramを初期化
  for i in depth:
    for j in leaf:
      model = RandomForestRegressor(
          max_depth=i, 
          min_samples_leaf=j, 
          random_state=0)
      score_array = np.array(cv_mse(X, y, model))
      mean_score = score_array.mean()
      
      if mean_score < best_val_score:
        best_val_score = mean_score
        best_params = {'max_depth': i, 'min_samples_leaf': j}

  print('best_val_score: {:0.3f}'.format(best_val_score))
  print('best_param: ', best_params)
  
  return best_params

まず今回使ってきたmax_depth=4のランダムフォレストでの交差検証による精度を確認します。

# ライブラリをインポート
from sklearn.ensemble import RandomForestRegressor

# 介入の有りと無しでデータを分割
df_0 = df_B[df_B["train_attend"]==0]
df_1 = df_B[df_B["train_attend"]==1]

print("【T-Learnerの精度評価_MSE】")
print("==介入無し==")
score_TL_0 = cv_mse(df_0[["enthusiam", "performance"]], df_0["stsf_level"],
       RandomForestRegressor(max_depth=4))
print(f"mse: {score_TL_0}")
print(f"mean_mse: {round(np.array(score_TL_0).mean(), 3)}")

print("==介入有り==")
score_TL_1 = cv_mse(df_1[["enthusiam", "performance"]], df_1["stsf_level"],
       RandomForestRegressor(max_depth=4))
print(f"mse: {score_TL_1}")
print(f"mean_mse: {round(np.array(score_TL_1).mean(), 3)}")

# それぞれでモデル構築(目的変数: 結果変数、説明変数: 共変量)
model0 = RandomForestRegressor(max_depth=4)
M_0 = model0.fit( df_0[["enthusiam", "performance"]], df_0["stsf_level"]) # 介入を受けてないモデル

model1 = RandomForestRegressor(max_depth=4)
M_1 = model1.fit( df_1[["enthusiam", "performance"]], df_1["stsf_level"]) # 介入を受けたモデル

# ITUとITTを求める
tau_0 = M_1.predict(df_0[["enthusiam", "performance"]]) - df_0["stsf_level"]
tau_1 = df_1["stsf_level"] - M_0.predict(df_1[["enthusiam", "performance"]])

# 推論時には異なる集団の推論も実施するためのモデルであるが評価ができないので、同じ集団内で評価する
print("【ITUとITTのモデル構築の精度評価_MSE】")

print("==ITU==")
score_ITU = cv_mse(df_0[["enthusiam", "performance"]], tau_0,
       RandomForestRegressor(max_depth=4))
print(f"mse: {score_ITU}")
print(f"mean_mse: {round(np.array(score_ITU).mean(), 3)}")

print("==ITT==")
score_ITT = cv_mse(df_1[["enthusiam", "performance"]], tau_1,
       RandomForestRegressor(max_depth=4))
print(f"mse: {score_ITT}")
print(f"mean_mse: {round(np.array(score_ITT).mean(), 3)}")
[output]
【T-Learnerの精度評価_MSE】
==介入無し==
mse: [0.151, 0.174, 0.144, 0.167, 0.171]
mean_mse: 0.161
==介入有り==
mse: [0.164, 0.192, 0.164, 0.17, 0.143]
mean_mse: 0.167
【ITUとITTのモデル構築の精度評価_MSE】
==ITU==
mse: [0.163, 0.168, 0.13, 0.168, 0.173]
mean_mse: 0.16
==ITT==
mse: [0.173, 0.189, 0.181, 0.163, 0.151]
mean_mse: 0.171

以下が、グリッドサーチを用いたモデルで順次処理していく実装です。

# ライブラリをインポート
from sklearn.ensemble import RandomForestRegressor

# 介入の有りと無しでデータを分割
df_0 = df_B[df_B["train_attend"]==0]
df_1 = df_B[df_B["train_attend"]==1]

# 介入無しのモデルのハイパラのグリッドサーチ
print("==介入無しモデルハイパラ調整の結果==")
best_params_0 = rf_cv_grid_search(df_0[["enthusiam", "performance"]], df_0["stsf_level"])

# 最も良かった精度のパラメータで介入無しのモデル構築
model0 = RandomForestRegressor(**best_params_0)
M_0 = model0.fit(df_0[["enthusiam", "performance"]], df_0["stsf_level"]) 

# 介入有りのモデルのハイパラのグリッドサーチ
print("==介入有りモデルハイパラ調整の結果==")
best_params_1 = rf_cv_grid_search(df_1[["enthusiam", "performance"]], df_1["stsf_level"])

# 最も良かった精度のパラメータで介入有りのモデル構築
model1 = RandomForestRegressor(**best_params_1)
M_1 = model1.fit(df_1[["enthusiam", "performance"]], df_1["stsf_level"]) 

# ITUとITTを求める
tau_0 = M_1.predict(df_0[["enthusiam", "performance"]]) - df_0["stsf_level"]
tau_1 = df_1["stsf_level"] - M_0.predict(df_1[["enthusiam", "performance"]])

# ITUのモデルのハイパラのグリッドサーチ
print("==ITUモデルのハイパラ調整の結果==")
best_params_ITU = rf_cv_grid_search(df_0[["enthusiam", "performance"]], tau_0)

# ITUの推定モデルを構築
model_ITU = RandomForestRegressor(**best_params_ITU)
M_ITU = model_ITU.fit(df_0[["enthusiam", "performance"]], tau_0)

# ITTのモデルのハイパラのグリッドサーチ
print("==ITTモデルのハイパラ調整の結果==")
best_params_ITT = rf_cv_grid_search(df_1[["enthusiam", "performance"]], tau_1)

# ITTの推定モデルを構築
model_ITT = RandomForestRegressor(**best_params_ITT)
M_ITT = model_ITT.fit(df_1[["enthusiam", "performance"]], tau_1)

# 傾向スコアを求める
from sklearn.linear_model import LogisticRegression

# 説明変数: 共変量
X = df_B[["enthusiam", "performance"]]

# 目的変数: 結果変数
Z = df_B["train_attend"]

# モデルを構築
model = LogisticRegression()
model.fit(X, Z)

# 傾向スコア算出
Z_proba = model.predict_proba(X)

# それぞれのモデルで全データの処置効果を推定し、傾向スコアで調整
ITU_pred = M_ITU.predict(df_B[["enthusiam", "performance"]])
ITT_pred = M_ITT.predict(df_B[["enthusiam", "performance"]])
ATE_i = Z_proba[:, 1]*ITU_pred + Z_proba[:, 0]*ITT_pred
ATE = ATE_i.mean()
print("==ATE==")
print(f"ATE: {ATE}")

# 推定した処置効果を描画
t_estimated = ATE_i
plt.scatter(df_B["enthusiam"], df_B["effect"], c="k", alpha=0.3, label="true_ITE") # 正解を描画
plt.scatter(df_B["enthusiam"], ITU_pred, c="b", alpha=0.3, label="estimated_ITU")
plt.scatter(df_B["enthusiam"], ITT_pred, c="g", alpha=0.3, label="estimated_ITT")
plt.scatter(df_B["enthusiam"], t_estimated, c="r", alpha=0.3, label="estimated_effect") # 推定値を描画
plt.legend()
plt.show()
[output]
==介入無しモデルハイパラ調整の結果==
best_val_score: 0.159
best_param:  {'max_depth': 4, 'min_samples_leaf': 4}
==介入有りモデルハイパラ調整の結果==
best_val_score: 0.164
best_param:  {'max_depth': 4, 'min_samples_leaf': 7}
==ITUモデルのハイパラ調整の結果==
best_val_score: 0.159
best_param:  {'max_depth': 3, 'min_samples_leaf': 3}
==ITTモデルのハイパラ調整の結果==
best_val_score: 0.165
best_param:  {'max_depth': 2, 'min_samples_leaf': 1}
==ATE==
ATE: 0.6706889541529172

X-Learner(ハイパラ調整)による熱心さに対する推定処置効果(研修受講->満足度)

精度に関しては、若干良くなりました。可視化結果についても、正解に近いところでハイパーパラメータ調整前よりも推定範囲が狭くなっていて良好になっているように見受けられます。ただ階段状に推定できていません(=熱心さが0.0~0.5の範囲が特に)。効果の大きさが切り替わって遷移している部分は難しいとかあるのでしょうか?
境界付近(熱心さが0または0.5あたり)は推定が難しくなるのはそれはそうだと思うので、両端に境界があると推定がうまくいかないイメージももつことはできますが。

おわりに

主にMeta-Learnersと呼ばれる、非線形な因果関係および処置効果を扱える因果推論手法を整理・実装しました。Meta-Learnersを用いることで、共変量の処置効果への非線形な影響を捉えられることが確認できました。また各手法で少しづつ挙動の違いがあり、今回のケースではX-Learnerが良い手法だなという印象でした。 しかしX-Learnerでの推定ITTおよびITUの挙動の違い、DR-Learnerでの大きく推定が外れるサンプルがあること、Meta-Learnersでも階段状にはうまく推定できていないところはうまく説明できそうになかったです。このあたりはさらに集計で深掘りをしたり、そもそも論文を読んでみたりすることで色々わかるのかもしれません。
今回は入門ということで、とりあえず実装できたことに満足です。

因果推論に入門( 処置効果が線形なデータで )

因果分析に興味があるということでこちらの書籍を勉強中。

book.mynavi.jp

こちらの書籍で紹介されている因果推論の手法を整理していきます。 今回は処置効果が線形な場合に、適用可能な手法の整理となります。

目次

因果推論とは

処置の効果測定をする手法のことです。まず因果関係というのは、ある変数を変化させたときに着目している変数が変化する関係にあります。(ある変数が大きいときに、着目している変数が大きかったり小さかったりしているのは相関関係です。) そしてこのときにどれくらい変化するのか=どれくらい効果があるのかを推定するのが、因果推論になります。

実際に因果推論する場合は因果探索やドメイン知識を活用し、因果ダイアグラムなどを作成して変数間の因果関係を示す必要があります。ただ今回は因果推論の手法の整理ということで以下のような擬似データを作成し、分析していきます。

データの説明

以下の構造を持つデータで因果推論を行います。
テレビCMが商品の購入量に影響するという場面を想定しています。
"商品の購入量: Y"が"テレビCMを見た: Z"という処置の影響を受けるとし、このときどれくらいの効果があるかを推定します。 また商品をどのくらい購入するかとテレビCMを見るかどうかは年齢と性別にも影響を受けるとします。

データの因果関係

このとき、
購入量を結果変数(処置の結果となる変数)、
CMを見たかどうかを原因変数(処置の原因となる変数)、
年齢、性別を共変量(因果関係に影響を及ぼす変数)
と呼びます。

データは以下で生成します。 処置効果が線形という点ですが、"商品の購入量: Y"を"テレビCMを見た: Z"の一次式(今回は係数を10として)で生成している点になります。

import random
import numpy as np

# 乱数のシードを固定
random.seed(9999)
np.random.seed(9999)

# シグモイド関数をimport
from scipy.special import expit

# 標準正規分布の生成用
from numpy.random import *

# グラフの描画用
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns

# その他
import pandas as pd

# 乱数のシードを固定
random.seed(9999)
np.random.seed(9999)

# データ数
num_data = 500

# 原因変数以外の説明変数( 共変量 )
# 年齢
x_1 = randint(15, 76, num_data) 

# 性別
x_2 = randint(0, 2, num_data)

# 原因変数
# シグモイド関数に入れる部分。テレビCMを見る可能性を確率変数として生成。(大きいほど見る確率は高くなる)
e_z = randn(num_data) # ノイズ(バイアス項)の生成
z_base = x_1 + (1-x_2)*10 - 40 + 5*e_z

# テレビCMを見る確率を計算 = シグモイド関数を計算
z_prob = expit(0.1*z_base)

# テレビCMを見たかどうかの変数(0は見ていない、1は見た)
Z = np.array([])
for i in range(num_data):
  Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0] # 0または1になる確率に対して、ランダムに選択している
  Z = np.append(Z, Z_i)

# 目的変数
# 各効果も指定する。年齢が高いほど減少、男性の方が増加、テレビCMを見ていると増加。
e_y = randn(num_data) # ノイズ(バイアス項)の生成
Y = -x_1 + 30*x_2 + 10*Z + 80 + 10*e_y

# テーブルを作成
df_A = pd.DataFrame({'age': x_1,
                   'gender': x_2,
                   'CM': Z,
                   'pur_vol': Y,
                   })

df_A.head(5)

[output]
# データの先頭行
age  gender  CM   pur_vol
37     0     1.0  56.691968
64     1     1.0  66.464076
62     0     1.0  14.130339
60     1     1.0  73.538865
69     1     1.0  40.870089

データの可視化

データを可視化します。

まずseabornのpairplotで簡単に見てみます。
購入量に対して、年齢、性別は傾向があります。しかしCMの効果については、CMを見た(1)がCMを見ていない(0)に対して、上にも下にもプロットされており、効果があるのかどうかわかりません。

sns.pairplot(df_B)

各変数間の可視化

またCMを見たかどうかで購入量のヒストグラムを比較した図が以下です。
CMを見ていないほうが購入量が右側に寄っており、可視化だけではCMが逆効果になっていそうと判断してしまいそうです。
また各カラムの平均を見ると年齢と性別の平均に乖離があります。

# CMを見たかどうかで購入量のヒストグラムを比較
df_A[df_A["CM"]==1]["pur_vol"].hist(color="r", alpha=0.3, label="CM=1")
df_A[df_A["CM"]==0]["pur_vol"].hist(color="b", alpha=0.3, label="CM=0")
plt.legend()
plt.show()

print(df_A[df_A["CM"]==1].mean())
print("======")
print(df_A[df_A["CM"]==0].mean())

CMを見たかどうか別の購入量のヒストグラム

[output]
==CMを見た集団の各カラムの平均==
age        51.702454
gender      0.463190
CM          1.000000
pur_vol    52.210321
dtype: float64
==CMを見ていない集団の各カラムの平均==
age        31.781609
gender      0.557471
CM          0.000000
pur_vol    64.855318

回帰分析

まずは回帰分析です。今回は結果変数が連続値なので、線形回帰を用います。
回帰モデルによる因果推論は比較的分かりやすいです。
目的変数として結果変数を、説明変数として考量すべき変数(原因変数、共変量)を用いたモデルを構築し、そのときの係数が効果の量になります。 以下で実行します。

# 説明変数
X = df_A[["age", "gender", "CM"]]

# 被説明変数(目的変数)
y = df_A["pur_vol"]

# 学習(回帰の実施)
model = LinearRegression()
model.fit(X, y)

# 回帰した結果の係数を出力
print("係数:", model.coef_)

[output]
係数: [-1.02228972 29.62534334 10.51298688]

各係数は、年齢: -1.0、性別: 29.6、CM: 10.5となりました。
データ生成時が、年齢: -1、性別: 30、CM: 10なのでそれなりに良い結果かと。
また各変数の影響が一括で出せるのもお手軽で良いかなと思いました。

傾向スコアによるIPTW(Inverse Probability of Treatment Weighting : 逆確率重み付け法 )

次は傾向スコアを用いたIPTWになります。 こちらは傾向スコア(ある属性情報に応じた処置を受ける確率)を用いて、ATE(Average Treatment effect, 平均処置効果)と呼ばれる推定した処置効果を算出する手法です。
手順としては、
1. 原因変数を共変量でモデリング
2. 共変量の元での傾向スコアを算出
3. ATEを調整化公式で変形した数式に傾向スコアを代入し推定したえ処置効果を算出
となります 傾向スコアは、今回はロジスティック回帰で求めます。
また3について、調整化公式とはdoオペレータで表した数式を、doオペレータ無しで表現する方法であり、
こちらと条件付き確率から同時確率へ変換で、ATEは以下のように変形できます。
【調整化公式】

 
P(Y=y | do(Z=z))

 
= \sum_{x=1}^n P(Y=y|Z=z, X=x) P(X=x)

 
= \sum_{x=1}^n \frac{P(Y=y, Z=z, X-x)}{P(Z=z, X=x)}

【ATEの式変形】
 ATE

 = E(Y1|do(Z=1)) - E(Y0|do(Z=0))

 
= \frac{1}{N} \sum_{x=1}^n\frac{yi}{P(Z=1|X=xi)}Zi - \frac{1}{N} \sum_{x=1}^n\frac{yi}{P(Z=0|X=xi)}(1-Zi)

IPTWを以下で実行します。

# 傾向スコアを求めるために、原因変数を目的変数としてモデリング
# 説明変数
X = df_A[["age", "gender"]]

# 目的変数(原因変数)
Z = df_A["CM"]

# 学習(回帰の実施)
model = LogisticRegression()
model.fit(X,Z)

# 係数を出力
print("--係数--")
print("係数beta:", model.coef_)
print("係数alpha:", model.intercept_)

# 各サンプルの傾向スコア
Z_pred = model.predict_proba(X)
print("--傾向スコアの推論値--")
print(Z_pred[0:5])  # 推定結果
print("--正解--")
print(Z[0:5])  # 正解

# 処置の効果を推定
Y = df_A["pur_vol"]
ATE_i = Y / Z_pred[:, 1] * Z - Y / Z_pred[:, 0] * (1-Z) # 式変形後
ATE = 1 / len(Y) * ATE_i.sum() # 平均処置効果 ATEを計算
print("推定したATE: ", ATE)

[output]
--係数--
係数beta: [[ 0.09254961 -0.67589802]]
係数alpha: [-2.84351622]
--傾向スコアの推論値--
[[0.35874405 0.64125595]
 [0.08288635 0.91711365]
 [0.05242315 0.94757685]
 [0.11572386 0.88427614]
 [0.05383419 0.94616581]]
--正解--
0    1.0
1    1.0
2    1.0
3    1.0
4    1.0

推定したATE 13.269226060428515

まずシグモイド関数の係数について正解: 0.1, -1, -3に対して0.09、-0.68、-2.84という結果でした。
傾向スコアは出力した分に関しては、1に対してのスコアが高くそこそこ妥当な結果かと思います。 処置効果の推定値は回帰分析よりは、正解から離れた結果となりました。  

以下に傾向スコアとのヒストグラムを描画します。

# 傾向スコアのヒストグラムを確認
temp_df = pd.merge(pd.DataFrame(Z_pred).reset_index(), df_A["CM"].reset_index()) # サンプル内で処置に対応する傾向スコアのテーブル作成
pd.concat([temp_df[temp_df["CM"]==0][0], temp_df[temp_df["CM"]==1][1]], ignore_index=True).hist() # 傾向スコアのヒストグラムを描画

# 各サンプルのTE(ITE)のヒストグラム
pd.DataFrame(ATE_i).hist()

IPTW法における傾向スコアのヒストグラム
IPTW法における各サンプルに対して推定された処置効果のヒストグラム

IPTW法は傾向スコアで割っているので、傾向スコアが0に近いと計算が不安定になります。
ヒストグラムを見ると小さい値が確認されるので、このあたりでそこそこの推定結果になったことが伺えます。 実際に各サンプルに対する処置効果の推定値のヒストグラムは-600~1000と広範な値をとっていました。

DR(Doubly Robust法)

最後はIPTW法と回帰分析を組み合わせたDoubly Robust(DR法)になります。
こちらはIPTW法では考慮できていなかった各サンプルの反実仮想を、回帰分析などで推定し計算に加える手法です。

手順としては、
1. 反実仮想の推定値を計算: 目的変数として結果変数を、説明変数として原因変数と共変量を用いてモデル構築し、原因変数が各場合における(今回は0または1)推定値を出力する
2. 傾向スコアを求める
3. 下記の式に従ってサンプルごとの推定処置効果を計算し、ATEを求める

1については、購入量に対して回帰モデルを構築して、推論時にCMを見たかどうかについて0、1それぞれにおいて予測して反実仮想の値としてやります。
3については、1で求めた値を使ってIPTWで用いた前半の項と後半の項をそれぞれ加えてやります。
処置の推定効果をそれぞれ
前半の項を、

 
\frac{yi}{P(Z=1|X=xi)}Zi → \frac{yi}{P(Z=1|X=xi)}Zi+(1-\frac{Zi}{P(Z=1|X=xi)})\hat{Y1}^i

後半の項を、

 
\frac{yi}{P(Z=0|X=xi)}(1-Zi) → \frac{yi}{P(Z=0|X=xi)}(1-Zi)+(1-\frac{(1-Zi)}{P(Z=0|X=xi)})\hat{Y0}^i

として差を求めます。
DR法を以下で実行します。

# 反実仮想の推定値(購入量)を求める

# 説明変数
X = df_A[["age", "gender", "CM"]]

# 目的変数
y = df_A["pur_vol"]

# モデリング
model = LinearRegression()
model.fit(X, y)

# 介入Z=0(テレビCMを見ていない)の場合の購入量を予測
X_0 = X.copy()
X_0["CM"] = 0
Y_0 = model.predict(X_0)

# 介入Z=1(テレビCMを見た)の場合の購入量を予測
X_1 = X.copy()
X_1["CM"] = 1
Y_1 = model.predict(X_1)

# 傾向スコアを求める
# 説明変数
X = df_A[["age", "gender"]]

# 被説明変数(目的変数)
Z = df_A["CM"]

# 回帰の実施
model = LogisticRegression()
model.fit(X, Z)

# 傾向スコアを求める
Z_pre = model.predict_proba(X)

# ATEを求める
Y = df_A["pur_vol"]

ATE_1_i = Y/Z_pre[:, 1]*Z + (1-Z/Z_pre[:, 1])*Y_1
ATE_0_i = Y/Z_pre[:, 0]*(1-Z) + (1-(1-Z)/Z_pre[:, 0])*Y_0
pred_eff = ATE_1_i-ATE_0_i
ATE = 1/len(Y)*pred_eff.sum()
print("推定したATE", ATE)

[output]
推定したATE 11.578048448712533

推定した処置効果は、11.6とIPTWの結果(13.2)より10に近く良好となりましたが、
線形回帰の方が良好な結果となりました こちらも各サンプルの処置効果の推定値を確認してみましょう。

DR法における各サンプルに対して推定した処置効果のヒストグラム

IPTW法と比べてとる値の範囲は狭くなりましたが、-200や500くらいの絶対値が大きな値が見られます。

各手法の結果まとめ

以下に今回の結果とデータ量を変えてみた(500 -> 250)結果をまとめました。

データ量  正解 回帰分析 IPTW 法 DR 法
250 10 7.3 11.6 9.1
500 10 10.5 13.3 11.6

データ量が少なると線形回帰は結果が悪化し、傾向スコアを用いた手法は結果が良くなりました。 とはいってもこのあたりはよく分かりません。(線形回帰の結果が悪化するのはイメージつきますが)

おわりに

今回は「つくりながら学ぶ! Pythonによる因果分析」を参考に因果推論の各手法を簡単に試しましたが、これらの手法でもまだ色々バリエーションはありそうです。例えば傾向スコアはロジスティック回帰ではなく、ニューラルネットワークや樹木系モデル(ランダムフォレスト、勾配ブースティング系)でも求めることができるよう(できそう)です (ここまで書いて思いましたが、シグモイド関数のモデル係数はあまり良くなかった)。 また今回は処置効果が線形なデータを使用しましたが、現実の課題は交互作用があったり非線形であったりすることがほとんどです。
どんなものが実務で使えるものかわかりませんがとりあえず勉強続けます。