正しい比較のための因果推論/計量経済学の基礎」のRソースコードを、Pythonで(ほぼ)再現しましたhttps://github.com/nekoumei/cibook-python本記事では、主にRではライブラリどーん!で済むけどPythonではそうはいかない部分の解説をします書籍の紹介https://www.amazon.co.jp/dp/B0834JN23Y上記Amazonに目次が載っているのでそれを見るのが早い気がしますが。正確な意思決定を行うためにどうやってバイアスを取り除くか?に焦点を当てて種々の因果推論の手法(傾向スコア/DiD/RDDなど)をRソースコードによる実装とともに紹介されています。全体を通して、現実問題の効果検証に因果推論を活用するにはどうすればいいか?という観点で書かれており、非常に実用的だなーと感じました。Notebookを作成しました。本Python実装にあたって、グラフの可視化ライブラリは一部を除いてplotly.expressを使用しています。Pagesをご確認ください。(images配下のhtmlを表示しています)plotly.expressに関しては、日本語だと下記記事がとても参考になります。令和時代のPython作図ライブラリのデファクトスタンダードPlotlyExpressの基本的な描き方まとめPythonで書くうえでの主要なポイント回帰分析statsmodelsのOLSおよびWLSを利用しています。results.summary().tables[1]nonrct_mreg_coefここで、Rのlm関数との違いが2つあります。1lmではカテゴリカル変数のダミー変数化を自動で行っているようです。pd.get_dummies()を利用しています。2lmではバイアス項が自動で追加されます。こちらも手動で追加します。実装内で利用されているデータセットの中には、experimentdatarのようにRパッケージとして公開されていたり、RData形式で公開されているデータセットがあります。これらをPythonのみで読み込む方法はないです。(あったら教えてください)コメント欄にて、upuraさんにご教示いただきました!https://qiita.com/nekoumei/items/648726e89d05cba6f432#comment-0ea9751e3f01b27b0adb.rdaファイルは、rdataというパッケージで、Rなしで読み込みが可能です。DataFrameへの変換を行う場合は下記のとおりです。r['vouchers']ch2_voucher.ipynbに記載しているとおり、Rパッケージのインストール等は事前にRの対話環境で行っています。ch3_lalonde.ipynbにて使用するデータセットは.dtaファイルです。こちらは、pandasのread_stata()によって読み込むことができます。pd.read_stata('https://users.nber.org/~rdehejia/data/nsw_dw.dta')傾向スコアマッチング(最近傍マッチング)傾向スコアマッチングのマッチング手法はいくつかありますが、書籍ではMatchItの最近傍マッチングを使っているようです。Pythonにはちょうどいいライブラリはなさそうだったので愚直に実装しています。print('2群のマッチングしかできません。1]で表現してください。0をgroup2とする。y.loc[matched_indices]treatment群の1点と最も傾向スコアが近いcontrol群1点をマッチングさせる、を反復しています。その際、距離がstdの0.2倍に収まるペアのみ抽出するようしきい値を設けています。このあたりの詳細は傾向スコアの概念とその実践が詳しかったです。ch3_pscore.ipynbにてMatchItによるマッチング結果との比較をしていますが、完全に一致はしていません。結論が変わらない且つ共変量のバランスも良い感じ(下記参照)のため概ね良しとしています。love.plot()のような便利ライブラリはないので自力で可視化しています。(images/ch3_plot1.htmlより)逆確率重み付き推定(IPW)IPWは実装がシンプルなのがいいところですね。こちらもWeightItと比較していますが、概ね正しそうです。print('2群のマッチングしかできません。1]で表現してください。0をgroup2とする。weightsCausalImpactch4_did.ipynbに記載しているとおり、下記2つのライブラリを使い、比較しています。https://github.com/tcassou/causal_impactoutputがキレイなのでいつもはdafitiのcausalimpactを使っているんですが、推定結果が書籍(Rの本家causalimpact)に近いのはtcassouのcausal_impactでした。両者とも、statsmodelsの状態空間モデルで推定してるみたいですが、、ちょっと実装の違いはよくわからなかったです。どちらも、R実装と比較して推定誤差がめちゃめちゃ小さくなってるんですよね。共変量が多いのがあまりよくないのかもしれない。dafiti/causalimpactのplottcassou/causal_impactのplot回帰不連続デザイン(RDD)Rだとrddtoolsでサクっと実行できますが、Pythonだとそうはいかないです。rddtoolsのrdd_reg_lmで使われている回帰式を確認します。余談ですが、私はてっきりRDDって、cut-off値の左と右でそれぞれ回帰モデルつくって、cut-off値での推定値の差分を取るんだと思ってたんですが、1つの回帰式で表すんですね。意味的には同じなのかな。以前なら0の値を持つバイナリ変数です。Dのcoefが効果量として確認したい値です。cはcut-off値です。以上を踏まえて、実装します。実装にあたってはrddtoolsのソースコードも参照しています。X)多項式回帰を行うための前処理として、scikit-learnのPolynomialFeaturesを利用しています。(images/ch5_plot2_3.htmlより)良いかんじに推定できました。Notebookを見ていただければわかりますが、効果量の推定も書籍とほぼ一致しています。(RDestimate)できませんでした、、、冒頭の「Pythonで(ほぼ)再現しました」のほぼの部分です。手法を用いて、最適なバンド幅を選択するんですが、そのバンド幅の推定がどうやってるのかよくわかんなかったです。ch5_rdd.ipynbでは、雑にバンド幅を変えたときのMSEが最適なバンド幅を見つけるように実装してみましたが、あまりうまく行っているように見えないです。(ch5_plot4.htmlより)うーんってかんじですね。純粋にバンド幅を狭めて推定するとバンド幅外は予測できなくなるのはそれはそうな気がするんですが、実際どうやるんですかね。cut-off近傍の推定にだけ興味があるからあまり気にしなくていいのかな。おわりに理解や実装が間違っていたり、なんかおかしいところとかがあったらぜひ教えてほしいです。upLoginユーザーは見つかりませんでした。