会社の勉強会で効果検証入門に取り組みましたが、実装はほとんど手をつけられていなかったり、個人的に気になったこともあったので、傾向スコアの部分を整理していきます。また個人的な解釈も多分にあるので、誤りなどあればご指摘いただければと思います。
こちらの書籍はRで実装(以下リンク参照)されているので、今回はRで実装しました。
また今回整理していた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)} に依存しない(処置の割り付けに影響を与えるのは観測された共変量のみ)
③ 条件付き正値性: 共変量Xを条件としたときに、どの個体も介入群または対照群に割付けられる確率が0または1でないこと(どの個体も、どちらにも割付けられる可能性がある)
この3点が傾向スコアを活用する上での、必要な仮定となります。
また以下が傾向スコアの定理です。これらの仮定から傾向スコアが介入の選択基準とは無関係であるといえるので、擬似的にRCTを再現していることになります。
【傾向スコアの定理】
① バランシング: 介入変数Zと観測された共変量Xは、傾向スコアe(X)が与えられたとき、条件付き独立である(つまり、傾向スコアe(X)が同じ値であれば、介入群と対照群における共変量Xの分布は同じ)。
② 条件付き独立性: 傾向スコアe(X)が与えられれば、潜在的結果変数{Y(1), Y(0)}と介入変数Zは条件付き独立である。
傾向スコアを用いた効果推定
まず傾向スコアは、介入変数Zを目的変数、共変量を説明変数としたロジスティック回帰で推定されます。また確率として算出できれば良いので、Gradient Boosting Decision Tree(GBDT)も選択として考えられます。
【手法①: 傾向スコアマッチング】
傾向スコアマッチングは、"(1): 介入された個体"に対して"(2): (1)と同じ傾向スコアを持つ非介入のサンプル"の結果変数を、(1)の非介入時の結果変数と考えます。つまりこの(2)を任意のアルゴリズム(最近傍法、最適ペアマッチング、復元の有無、etc•••)で探し出し、各グループ間の結果変数の差を効果の推定値とします。ただし傾向スコアが高いor低いとこを捨てることになり、解析結果が当てはまる対象集団が限定的になると以下の記事でも述べられていました。
(https://riklog.com/research/propensity-score/)
このときの効果推定の計算式は、
となります。またこの方法で求められる効果は、ATT(Average Treatment effect on Treated)になりますが、マッチングでもATE(Average Treatment effect)を求める方法はあります(この記事では触れません)。
【手法②: 傾向スコアを用いたIPTW】
IPTWはInverse Probability of Treatment Weighting(逆確率重み付き)推定と言います。
傾向スコアの逆数を重みに使って、全サンプルが介入受けた場合と非介入だった場合の各期待値を算出し、その差分を効果の推定値とします。
つまりIPTWでは、データを捨てずに効果推定を達成しようとします。そもそも傾向スコアは介入が行われる確率なので、RCTが実現しているデータセットにおいては0.5になるはずであり、そこからずれるということは、何かしらの理由で介入が偏りセレクションバイアスが発生している状態と考えられます。したがって介入群に属しているかつ傾向スコア(介入が行われる確率)が低い個体は、RCTであった場合と比較して、観測した量が少ない個体と考えられます。これは他の個体においても同様と考えられるので、傾向スコアの逆数で重み付けすることで、見かけ上各傾向スコアにおける観測した量を同じにし、セレクションバイアスが無い状態を作り出しています。
効果推定の計算式は、
となります。データ全体での平均的な効果を推定していることになるので、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))
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)
ヒストグラムを見た時の直感とは異なり、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")
傾向スコアの分布としては、介入群と対照群の分布がよく似ていることがわかります。また重要なバランシングの方ですが、共変量および傾向スコアの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は回帰分析と比べて結果が改悪しました。傾向スコアの分布は介入群と対照群で大きく異なっており、共変量のバランシングもうまくできておりません。 これは、
- 調整前データの対照群がかなり絞られていた
- 調整前の傾向スコアの高い領域で介入群の方が多かった
ということから傾向スコアの似たサンプルと一致の具合が十分で無かった、つまり良いマッチング相手が無かったためかもしれません。傾向スコアマッチングでの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
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)
# 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
おわりに
今回は効果検証入門の傾向スコアの部分を整理していきました。
*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
