論理の流刑地

地獄の底を、爆笑しながら闊歩する

X→Yの効果と交互作用のあるZを見つける

統計学 R

Introduction

p-hackingやasterisk-seekingなんて言葉も人口に膾炙している通り、調査データから「有意な関連」を見出すこと自体が目的化しがちなのが経験科学における定量的研究である。
HARKing, p-hacking, asterisk-seekingを助長している学術界 - 講義のページへようこそ
↑参考：京大関口先生（経営学）のブログ

わかりやすいp-hackingというのは、X→Yの関連性を有意にする（"asteriskをつける"）ためにデータや変数を加工したり手法を勘案したりすることである。
これは独立変数と従属変数という2変数間の問いである。

もう少し進んだ問いとして、経験科学では「交互作用の問い」が出てくることがある。
X→Yの効果の調節弁(moderator)となるようなZを見つける、という問いである。
（例：通塾が学力に与える効果は性別によって異なるか、TVCMが販売を促進する効果は年齢層によって異なるか、など）
そこでは、Zの水準によってX→Yの効果の大きさ・強さが変容することを明らかにすることが"発見"となる。

じゃあこちらの種類の問いが俎上にあるときに、p-hackingするにはどうすればいいかというと、
「X→Yの効果のバラつき」を評価基準として、

  • 大きいバラつきをもたらしてくれるようなZがどの変数であるか
  • Zのどの値の水準でデータを分割したときバラつきが大きくなるか

を教えてくれる手法が欲しいわけですが、それはすでにあるよね、という話をする。

How To

やりかたを書く

Data

別にデータはなんでもいいのですが、試しにここではpanelrパッケージに含まれるWageDataを使う。

library(panelr)
data(WageData)
head( WageData ,2)
# exp wks occ ind south smsa ms fem union ed blk
# 1   3  32   0   0     1    0  1   0     0  9   0
# 2   4  43   0   0     1    0  1   0     0  9   0
# lwage t id
# 1 5.56068 1  1
# 2 5.72031 2  1

データの説明(CRANのref. manualから引用)は以下の通り

  • wks：週間労働時間
  • lwage：対数変換所得
  • union：組合所属
  • ms：既婚ダミー
  • occ: ホワイトカラーダミー
  • ind: 製造業ダミー
  • south:南部居住ダミー
  • smsa：大都市部居住ダミー
  • fem: 女性ダミー
  • blk：アフリカ系ダミー
  • ed:　教育年数
  • exp: 職場における経験年数

本データは縦断データだが、本稿ではその性質は無視してすすめる。

説明をわかりやすくするために、四大卒か否かの変数(univ)を作成して、本記事における処置変数とする。また、分析対象は男性に絞る。

wage_df <- WageData %>>%
  mutate( univ = if_else( ed >= 16  , 1, 0 ))%>>%
  filter( fem == 0 )

Pre Analysis

今回のケースでXは学歴（大卒か否か）で、Yは（対数）収入である。
まず基本モデルとして、経験年数、労働時間、学歴のみ回帰したモデルを推定する。

mdl0 <- lwage ~ exp + wks
mdl1 <- lwage ~ exp + wks + univ

lm0 <- lm(mdl0 ,data = wage_df)
lm1 <- lm(mdl1 ,data = wage_df)
summary(lm1)
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept) 6.2036136  0.0641618  96.687  < 2e-16 ***
#   exp         0.0097034  0.0005979  16.229  < 2e-16 ***
#   wks         0.0048399  0.0013119   3.689 0.000228 ***
#   univ        0.3714632  0.0146922  25.283  < 2e-16 ***

職場での経験年数や労働時間を統制しても、exp(0.3714) = 1.447倍大卒者のほうが平均的に稼いでいるとわかる。
これがoverallでみたときのX→Yの効果であるが、さて、問題はこの効果の調整弁となるZをデータから見つけることとなる。

How to Find : Causal Tree

occ, ind, south, smsa, ms , union , blkやその組み合わせから、効果の調整弁となるZを見つけていく。

ここでは、causalTreeパッケージ*1のhonest.causalTree()を使う。

準備として、処置変数「以外」を投入したモデルを推定して、残差を推定しておく。
この残差所得を使うことで、交絡による処置群/統制群間のbaseline bias（cf. Morgan and Winship(2014)）は除去されている

resid_fml <- lwage ~ exp + wks + ind + south + smsa + ms + union + blk
resid_lm <- lm( resid_fml , data = wage_df)
wage_df$lwage_resid <- residuals( resid_lm)

次にCausal Tree（因果木）を用いて、Zをみつけていく。
Causal TreeはX→Yの効果のバラつきを最大化するようにデータを分割していく木を生成する手法であるが、開発者による解説にもある通り、分割に用いられる基準は多数用意されている(TOT/fit-H/fit-A/CT-H/CT-A/tstats)。
ここでは、CT-Hを開発者が想定していないちょっとトリッキーな方法で使うことで、Zを探索する

ct_fml <- lwage ~ exp + wks + ind + south + smsa + ms + union + blk #formula
ct_alpha <- 0.7
ct_res <- honest.causalTree(ct_fml , data = wage_df , treatment = wage_df$univ , est_data = wage_df,est_treatment = wage_df$univ ,split.Rule="CT",split.Honest = T, HonestSampleSize =  nrow(wage_df), split.alpha = ct_alpha , split.Bucket = T , bucketNum = 20 , cv.option ="CT", cv.Honest = T, cv.alpha= ct_alpha)

※細かい引数の意味などについては、上記の開発者(S. Athey氏)の出している解説資料を読んでほしい。

トリッキーな方法というのは、dataとest_dataに全く同じデータを想定している、という点である。
元々の目的として、Honest推定においては、葉の分割に使う（共変量の）訓練用データと葉の内部でXのYへの効果を推定するための推定用データを分割している。
これは、機械学習の手法ならではの予測性能の汎化を企図したものである。
※ちなみにこの手法の元論文でウリにされているオリジナルな点でもある。

しかし、今回の目的はただ目の前にあるデータにおけるX→Yの調整弁たる変数Zを見つけることにあるので、
推定データにも訓練用データにも同一のデータを指定するのである。

さて、本題に戻ると、さきほどのデータを枝切りして*2その詳細を見てみると、
以下のようになっている

pruned_ct  <-  prune_rpart( ct_res, nsplits=6) #分割数で枝切り
pruned_ct
f:id:ronri_rukeichi:20210618141818p:plain
Causal Treeの分割の様子

上のsplitの情報からは、

  • 製造業であるか否かが一つ目の大きな分水嶺になっている（製造業のほうが大卒のメリットが大きい）
  • 製造業の場合、大都市か否かが大卒の収益に関わっている（製造業だと、大都市でないほうが大卒のメリットが大きい）
  • 製造業でない場合、労働時間が大きなカギとなる（非製造業だと労働時間が長い場合に大卒のメリットが大きくなる）

ということが分かる。

じゃあ、という感じで学歴の主効果と交互作用項を追加したモデルを推定してみよう。

add_fml1  <- lwage ~ exp + wks + univ + ind + south + smsa + ms + union + blk + ind:univ #formula
add_fml2 <- update( add_fml1 , .~.+ind:univ + ind:univ:wks + univ:smsa + univ:smsa:ind)
add_lm1 <- lm( add_fml1 , data = wage_df)
add_lm2 <- lm( add_fml2 , data = wage_df)
summary(add_lm1)
summary(add_lm2)

出てきたモデルの交互作用項目をみると、しっかりと（0.1%水準で）「有意」になっている。

f:id:ronri_rukeichi:20210618150814p:plain

別にこのcausalTreeによる探索のプロセスを、論文なり何なりに書く必要はないわけで、
はじめからこの交互作用を想定していたような問題設定や仮説を呈示して、論文を書くことは可能だろう。*3
これは学術内部の規範からして褒められた分析（あるいは執筆）フローではないのだけど、
剽窃や改竄のような、研究倫理上の問題として罰せられるようなものでもない*4。
一丁上がり、である。

Conclusion

一部の経験科学*5においては、汎化性能を持った予測モデルを構築することではなく、
今目の前にある調査データを分析して、「そのデータから何を引きだせるか」が全てになりがちである。
やや誇張した言い方をすれば、「いま、ここにあるデータ*6にいかに語らせるか」が全てになってしまう。

探索的な分析と仮説ドリブンの分析を厳密に分けるってのは難しい、みたいな話は『社会科学の方法論争』でもなされていたテーマでもあった。

社会科学の方法論争 [原著第2版]: 多様な分析道具と共通の基準

社会科学の方法論争 [原著第2版]: 多様な分析道具と共通の基準

Amazon

deductionでもinductionでもないabductionがあるよね、というのは古くから言われ続けてきたことではあるので、
何がデータに対して「誠実」であるのかも、一筋縄ではいかないものだな、というのが個人的な感想ではある。
まぁ自分は研究者ではないから、そこまで考える必要はないのかもですが。


www.youtube.com

Enjoy!!

*1:CRANには上がっていないので、install_githubしましょう

*2:自作関数を使った

*3:たとえばこの場合だったら、「非大都市圏における製造業と学歴主義」みたいなことをあたかもはじめから考えていたようにテーマにしちゃえばいいのだ。

*4:まぁそもそも「それ楽しいの？」っていうのはありますが

*5:政治学とか労働経済学とか想定してます。もちろんこれらの領域でもそうでない場合があることも理解はしているんだけど。

*6:リクルートのキャッチコピーっぽい