傾向スコアを用いた効果検証

　施策の有効性を測定し、その結果を組織や政策の意思決定に活かすことは、ビジネスから社会政策に至るまで幅広い領域で重要なテーマといえます。新商品のキャンペーンや既存顧客へのプロモーション効果を正確に把握するといったマーケティング視点においても、あるいは公共政策の策定や福祉分野など社会的インパクトを評価する場面でも、「この施策は本当に効果があったのか」「競合他社やトレンドなど、他の要因の影響をどこまで排除できるか」を検証するために、近年注目されているのが因果推論の考え方になります。
　因果推論では、「どのようにデータを使って“原因”と“結果”の関係をより厳密に導けるか」を統計的に考えます。理想的にはランダム化比較試験（RCT）を行って効果検証をするのがベストと言われていますが、マーケティング施策や公共政策の実施の現場では予算や時間、運用ルールなどの制約が多く、しばしば観察データのみをもとに効果を推定せざるを得ない状況があります。そこで活用される手法のひとつが、傾向スコア分析です。この手法では、施策を受けたグループと受けていないグループが持つ背景要因（年齢、購買履歴、ウェブサイトの訪問頻度など）を統計的にコントロールし、まるでランダム化実験を行ったかのように効果を推定することができます。
　こちらの資料では、まず因果推論や効果検証の前提や必要性をあらためて整理し、そのうえで傾向スコア分析の基本的な考え方について解説します。さらにLalondeのデータセットを事例として、実際に傾向スコアを活用した効果検証の流れをまとめました。因果推論の基本を押さえたい方や、観察データに基づいた施策評価をより厳密に行いたいときのご参考になるかと思います。

　

　
また分析に用いたコードと資料はこちらのリンク先になります。 [コード(notebook)][資料] 　　
　
スライドのタイトル

  • 効果検証と因果推論について
    　なぜ効果を測るのか / セレクションバイアスと誤った効果測定 / 数式で見るセレクションバイアス / 因果推論の根本問題 / ランダム化比較実験と平均的効果 / ランダム化比較実験の難しさ / 因果推論とは / 因果推論の手法
  • 傾向スコアについて
    　傾向スコアを用いた効果検証の概要 / 効果の種類 / 傾向スコアを用いる場合の仮定 / 傾向スコアとは / 手法1：傾向スコアマッチングについて / 手法2：逆確率重み付け法(IPW)について / 標準化平均差について
  • Lalondeデータの分析
    　分析の問題設定とデータの概要 / データの可視化 / 傾向スコア分析：傾向スコアの計算 /傾向スコアマッチング(CPS1のATTの推定) / IPW(CPS1のATE推定) / IPWでのCPS1のATT推定
  • まとめとその他

以上になります。おかしい点などあれば、ご指摘いただければと思います。

LSA(Latent Semantic Analysis)と特異値分解

　自然言語処理タスクの多くは、文字列→ベクトルまたはベクトル→文字列になります。また前者の代表的なタスクには、文章分類、感情分析、文書検索、情報抽出などがあります。
　またそのための技術として近年は深層学習モデルが中心となってきていますが、古典的な統計的言語モデルも重要な技術です。今回こちらの資料では、その中でもLSA(Latent Semantic Analysis)というトピックモデルおよびそれを用いた文書分類について整理してみました。

　
また分析に用いたコードと資料はこちらのリンク先になります。 [コード(notebook)][資料] 　
　
スライドのタイトル

以上になります。おかしい点などあれば、ご指摘いただければと思います。

回帰不連続デザインに入門2（識別条件と仮定の概要）

因果推論は施策効果を推定する手法になりますが、各手法に通じる大きな考え方として「ランダム化比較実験が行えない状況で、なんとか反事実を作り出し、事実と反事実の差を比べる」というものがあります。
そのため各手法には分析対象となるデータが満たすべき仮定が存在し、その仮定の下で統計学を駆使することで効果推定がされています。 以前の記事でRDDについて整理しましたが、今回の資料ではそのRDDが満たすべき仮定（識別条件）について整理しました。

　
また資料はこちらのリンク先になります。 [資料] 　
　
スライドのタイトル

  • 識別条件と成立しないケース
  • 可視化での判断
  • 検定での判断①
  • 検定での判断②
  • 共変量バランステストでの検定
  • RDDで使えるPythonライブラリ
  • 補足：単回帰分析の偏回帰係数の検定

以上になります。おかしい点などあれば、ご指摘いただければと思います。

モンテカルロDCF法による事業価値の算出

事業価値の評価は、投資判断やM&A、資金調達の際に重要な指標になります。適切な評価により、資本配分の効率化や経営戦略の策定、リスク管理が可能になります。
伝統的な評価方法にDCF法（Discount Cash Flow）があり、これは将来の価値を現在価値に割り引いて算出する方法です。
また拡張した評価方法としてモンテカルロDCF法があり、こちらは乱数生成した変数に基づいて事業価値を算出する方法です。つまり異なるシナリオ下での企業価値の変動範囲を理解できるため、不確実性を考慮した事業価値の評価が可能になります。 経営の意思決定では最悪のシナリオの考慮が重要であるため、モンテカルロDCF法はより良い意思決定のサポートができると考えられます。

こちらの資料ではモンテカルロDCF法の概要に加え、実際に企業の株価データと有価証券報告書に基づいて分析した結果をまとめました。

　

また資料と分析に用いたコードはこちらのリンク先になります。 [コード(notebook)][資料]

　
スライドのタイトル

　
以上になります。おかしい点などあれば、ご指摘いただければと思います。

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

RDD(Regression Discontinuity Design, 回帰不連続デザイン)は因果推論手法の一つであり、施策対象になるかどうかがある一つの基準で決定される場合に用いられる手法です。
例えば購買量に応じたメールマーケティングのように、ある指標に基づいて施策が振り分けられた状況で用いられます。その境界直前・直後の対象者は同質であると考えられることから、内的妥当性の高い手法とされています。
こちらの資料ではRDDの概要とシミュレーションデータを用いた効果推定を整理しました。

　
また資料と分析に用いたコードはこちらのリンク先になります。 [コード(notebook)][資料] 　
　
スライドのタイトル

  • RDD(回帰不連続デザイン)とは？
  • パラメトリックとノンパラメトリック
  • 共分散分析によるATEの推定の問題
  • 局所的な効果推定（ノンパラメトリック分析）
  • バンド幅の最適化
  • シミュレーションデータで分析（分析するデータについて）
  • シミュレーションデータで分析（RDプロットによる可視化）
  • シミュレーションデータで分析（連続性の仮定の診断：分布による確認）
  • シミュレーションデータで分析（潜在結果変数での大域的な平均処置効果）
  • シミュレーションデータで分析（各種手法による大域的な効果推定）
  • シミュレーションデータで分析（局所的な効果推定：ノンパラメトリック分析）

以上になります。おかしい点などあれば、ご指摘いただければと思います。

傾向スコアに入門 / マッチングと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