R
Exploratory

Exploratoryを使った統計的因果推論


はじめに

今回は、RとExploratoryを使って、傾向スコアを用いた統計的因果推論(Causal Effects)を行ってみようと思います。想定読者層は、統計的因果推論を学び始めた人や私のようなExploratoryユーザーの方です。

この分析記事では統計的因果推論の詳細な説明や数理的な側面については扱いません。ごめんなさい。「Exploratoryで行う統計的因果推論」なので、この後利用するロジスティック回帰分析(Logistic regression)、傾向スコア(Propensity Score)、逆確率重みつき推定法(Inverse probability weighting)、また、これらに関連する諸々の概念も詳しくは扱いません。

数理的な側面を知りたい方は、個人的に知っている範囲の紹介にとどまりますが、星野先生津川先生林先生清水先生KRSK先生などのブログを読んだり、書籍を購入してください。大変わかりやすく難解な統計的因果推論を説明されています。

これ以降で、Exploratoryの使い方含め、間違いがありましたら、暖かくご指摘いただけますと幸いです。


統計的因果推論とは

統計的因果推論とは、統計学者ドナルド・ルービン先生の「ルービンの因果モデル(Rubin causal model)」の枠組みをお借りすると、「もし〜だったら」という潜在的なアウトカム(Potential outcome)という考え方を導入し、統計的因果推論を行います。

つまり、各個人の介入を受けた場合のアウトカム、介入を受けなかった場合のアウトカムが手に入れば、その差が個人レベルでの介入効果(Treatment effect)となり、個人の介入効果を推定でき、そこから、個人の介入効果の平均値を取ることで、集団レベルでの介入効果(=因果効果)を推定できるというものです。言い換えると、「すべての個人が介入を受けた場合」と「すべての個人が介入を受けなかった場合」の差こそが因果効果と考えるわけです。非常にざっくりとした説明ですので、詳細な解説は先程の先生方のブログや書籍を参照ください。

具体的にどんな場面で統計的因果推論が使われているかというと、医学の現場であれば、「右心カテーテル検査が患者の予後に与える影響」を考えるために因果推論が使われてたりします。マーケティングであれば、「TVCMが与える売上への影響」を考えるために利用されていたり、スポーツの分野であれば「野球のバントが得点への効果」を考えるために利用されています。つまり、何らかの介入が与える影響の効果を測定する際に利用されるのが統計的因果推論です。

今回の分析の流れは、ロジスティック回帰で傾向スコアを計算し、逆確率重みつき推定で重みづけた回帰分析で平均因果効果(average treatment effect:ATE)、介入群における処置効果(average treatment effect on the treated:ATT)を推定することで因果効果を求めます。


道具立て

今回使用するデータはLalonde's National Supported Work Demonstration Dataです。これは、Dehejia先生とWahba先生の論文「Causal Effects in Non-Experimental Studies: Reevaluating the Evaluation of Training Programs」で使用されているデータで、1976年に実施されたアメリカの職業訓練プログラム(treat)が、訓練実施後の1978年の年収(re78)をアップさせたかどうか、因果効果があったのかを調べた際に利用されたデータです。{MatchIt}パッケージに内包されているので、そこから利用します。

変数
内容

treat
treatment indicator; 1 if treated in the NSWD, 0 if from the Current Population Survey

age
age, a numeric vector.

educ
years of education, a numeric vector between 0 and 18.

black
a binary vector, 1 if black, 0 otherwise.

hispan
a binary vector, 1 if hispanic, 0 otherwise.

married
a binary vector, 1 if married, 0 otherwise.

nodegree
a binary vector, 1 if no degree, 0 otherwise.

re74
earnings in 1974, a numeric vector.

re75
earnings in 1975, a numeric vector.

re78
earnings in 1978, a numeric vector (outcome variable).


データの探索

ざっくりとデータの状態を確認します。「サマリー」タブから欠損値がないか、外れ値がないかを見たところ、特に問題はなさそうです。445人のうち41%が職業訓練を受けているようです。

causal_infer_-_stat-2.png

次に介入群とコントロール群に分けて平均値の差を確認します。ブランチを切って進めます。画像を貼るまでもないと思いますので、言葉で説明します。まずgroup_by(treat)でグループ化し、summarise_all(funs(mean))で平均値を計算します。そこから、treat以外の列をgather(key, value, age:u75)でwide型からLong型に変換し、valueround(value,2)をかけて、spread(key = treat, cal = value)で、Long型からwide型に変換し、平均値の差をmutate()で計算します。


R

# Steps to produce the output

exploratory::read_delim_file("/Users/name/Desktop/lalonede.csv" , ",", quote = "\"", skip = 0 , col_names = TRUE , na = c('','NA') , locale=readr::locale(encoding = "UTF-8", decimal_mark = "."), trim_ws = TRUE , progress = FALSE) %>%
readr::type_convert() %>%
exploratory::clean_data_frame() %>%
group_by(treat) %>%
summarize_all(funs(mean)) %>%
gather(key, value, age:u75, na.rm = TRUE, convert = TRUE) %>%
mutate(value = round(value, digits = 2)) %>%
spread(treat, value) %>%
rename(not_treat = `0`, treat = `1`) %>%
mutate(diff = treat-not_treat)

計算結果を見てみると、訓練実施後の1978年の年収(re78)はすでに、介入群のほうが高いようです。しかし、この表を見てわかるように、群間で各共変量で差があるのがわかります。こうなると、treatの介入によって、年収がアップしたのか、hispだから年収が低くなっているのか、効果の切り分けができません。(※もう少し分かりやすく差がでてくれればよかったのに)

causal_infer_-_stat.png


傾向スコアの算出

まずはロジスティック回帰分析を使って傾向スコアを求めます。ロジスティック回帰のアウトカムはtreatで、共変量にはre78を除いた他の項目を使います。

傾向スコアは、ランダム化比較試験(Randomizes Controlled Trial)が様々な理由によって実施できない場合に力を発揮します。つまり、傾向スコアは、今回のような観察研究を行った際に発生する魑魅魍魎とした交絡に対して、共変量を調整して因果効果を推定するために用いられるバランス調整の統計手法です。

スクリーンショット 2019-05-18 3.43.28.png

そして、ここで一つ「ブランチ」を切ってROCを計算します。理由は、計算した傾向スコアが介入(treat)されるかどうかを正しく識別できているかを判断するためです。

今回の因果推論の方法の他にも、傾向スコアをマッチングに使用し、因果効果を推定する方法もあります。その場合にAUCが0.6程度しかない場合、マッチングを行うと、似たような分布からサンプルサイズを減らすだけで、反対に、AUCが0.9以上のような場合、マッチングを行うと、似たような分布からのオーバーラップが少なく、マッチングできる数が減ってしまい、サンプルサイズを極端に減ってしまうなど、傾向スコアマッチングを行う際の判断材料として利用されます。

スクリーンショット_2019-05-18_3_46_24(2).png

予測された確率の列には、ロジスティック回帰で求めたPredicted_probabilityを入れ、実際の値の列には、treatを入れます。

causal_infer_-_cStat.png

ROCからAUCを可視化します。見ての通り、treatの分類があまりうまく行われていません。

causal_infer_-_cStat.png

実際の数値で確認する場合は、ブランチの前に戻り、「アナリティクス」タブから同様の設定でロジスティック回帰分析を実行すれば、「サマリ」タブからAUCが確認できます。0.62なので、あまり良くないですね。慣習的には0.8くらいを目指せば良いそうです。なので、モデルを改善する必要があります。今回はチュートリアルなので、このままで行きます。

causal_infer_-_lalonde-2.png


逆確率による重み付け

逆確率重みつき推定で重みづけた回帰分析で因果効果を推定するために、重みを計算します。今回は、平均因果効果(ATE)と介入群における処置効果(ATT)を推定するので、2つの重みを傾向スコアから作ります。ざっくりとしたイメージとしては、重みをつけることで、介入群とコントロール群のベースライン特性のバランスを揃えるという感じです。

傾向スコアが高い個人は、介入を受ける確率が高く、介入群に多く存在し、コントロール群には少ない一方で、傾向スコアが低い個人は、介入を受ける確率が低く、介入群に少なく、コントロール群には多く存在しているはずです。そのため、重みをかけることによって、両群のバランスを整えるのが、逆確率により重み付け法です。重みの計算式は下記のとおりです。

重み
介入群
コントロール群

ATE Weight
1/傾向スコア
1/(1-傾向スコア)

ATT Weight
1
傾向スコア/(1-傾向スコア)

causal_infer_-_lalonde.png


因果効果の推定

回帰分析を行ってアメリカの職業訓練プログラム(treat)が、訓練実施後の1978年の年収(re78)をアップさせることにつながったのかを分析していきます。まず、計算結果を分かりやすく表示させるために、意図的に新しい変数not_treatを作成しておきます。後の回帰分析でre78 ~ treat + I(1 - treat) + 0と記述しても同じです。

causal_infer_-_lalonde.png

準備が整ったので、ブランチをATE用とATT用に分けます。実際の分析では、使用用途に合わせて、平均因果効果(ATE)と介入群における処置効果(ATT)の計算を行います。


平均因果効果(ATE)

まずは平均因果効果(ATE)の効果を見てみます。GLMから下記のように設定します。Weight Vectorに先程計算した重みipw_ateを入力して、実行します。

causal_infer_-_IPW_ATE.png

結果を確認すると、treatのパラメタが6212.99となっており、not_treatが4571.41となっています。つまり、6212.99-4571.41=$1641.58が職業訓練プログラム(treat)の因果効果と言えます。つまり、職業訓練を受けたことで、$1674.5分、年収がアップしたことがわかりました。

causal_infer_-_IPW_ATE-2.png


介入群における処置効果(ATT)

次に、介入群における処置効果(ATT)の効果を見てみます。GLMから下記のように設定します。Weight Vectorに先程計算した重みipw_attを入力して、実行します。

causal_infer_-_IPW_ATT.png

結果を確認すると、treatのパラメタが6349.15となっており、not_treatが4594.59となっています。つまり、6349.15-4594.59=$1754.56が職業訓練プログラム(treat)の因果効果と言えます。つまり、職業訓練を受けたことで、$1754.56分、年収がアップしたことがわかりました。

causal_infer_-_IPW_ATT.png


おわりに

統計的因果推論は難しいですが、使えるようになると実務でも非常に役立つと思います(私自身、使いこなせてないですが…泣)。星野先生津川先生林先生清水先生KRSK先生などのブログを読んだり、書籍は本当にわかりやすいので、ぜひ一読されることをおすすめします。


参考文献

調査観察データの統計科学―因果推論・選択バイアス・データ融合


Rスクリプト


R

# Set libPaths.

.libPaths("/Users/name/.exploratory/R/3.5")

# Load required packages.
library(janitor)
library(lubridate)
library(hms)
library(tidyr)
library(stringr)
library(readr)
library(forcats)
library(RcppRoll)
library(dplyr)
library(tibble)
library(exploratory)
library(bit64)

# Steps to produce the output
exploratory::read_delim_file("/Users/name/Desktop/lalonede.csv" , ",", quote = "\"", skip = 0 , col_names = TRUE , na = c('','NA') , locale=readr::locale(encoding = "UTF-8", decimal_mark = "."), trim_ws = TRUE , progress = FALSE) %>%
readr::type_convert() %>%
exploratory::clean_data_frame() %>%
build_lr(treat ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74 + u75) %>%
prediction_binary(data = "training") %>%
select(-(predicted_value:standardised_residuals), -(conf_low:predicted_label)) %>%
rename(ps = predicted_probability) %>%
mutate(ipw_ate = if_else(treat == 1, 1/ps, 1/(1-ps)),
ipw_att = if_else(treat == 1, 1 , ps/(1-ps)))%>%
select(-(age:re75), -ps, -u74, -u75) %>%
mutate(not_treat = 1 - treat) %>%
build_glm(re78 ~ treat + not_treat + 0, family = gaussian(link = "identity"), weights = ipw_att)
# build_glm(re78 ~ treat + not_treat + 0, family = gaussian(link = "identity"), weights = ipw_ate)