共変量バランス

# はじめに 

前回の投稿ではロジスティック回帰を用いた1変数による傾向スコアマッチングを自作したpython関数にて試みました. 　


https://qiita.com/TSYK11361/items/3e03ca863736cc1875f0

結果としては失敗でしたので、今回はRのmatchit()を用いて得られた傾向スコアマッチングの結果について考察したいと思います.  


著書「効果検証入門」には以下趣旨の記載があります.  
本稿の趣旨はずばりそれです.  

- 共変量のバランスが取れているか可視化により確認することができる.  

# メモ 
下図は、4つのデータでそれぞれ共変量"history"をヒストグラムにし、介入"treatment"の有無で色分けしたものです. 以下の投稿で取り組んだものになります.　　  


https://qiita.com/TSYK11361/items/11e0768f46c24d9002a4


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/facc89e2-3cca-77e0-4a41-69b577de834f.png)


上図の最左のRCTはランダムな介入割当が起こった場合のデータです.  
ご覧のとおり、介入割当と"history"が無関係であることが分かります.  

一方で、上図の左から2番目のbiased_dataでは、介入群(treatment == 1)と非介入群(control: treatment == 0)で"history"の分布の形状が異なることが分かります.  

また、上図の左から3番目のmatched_dataは傾向スコアマッチングにより組成されたデータです.  
ぱっと見、介入群, 非介入群間の"history"の分布の差が、biased_dataのそれよりも小さくなっている、とは見えませんね.  

念の為補足しておくと、傾向スコアはロジスティック関数により得られるとした上で、ロジスティック回帰モデルの係数推定が行われています.  

補足を続けます、下図はpythonによりロジスティック回帰の結果を可視化したものです.  


下図の上段左側は共変量"history", ロジスティック関数の出力をそれぞれx軸, y軸とした座標平面上にデータをプロットしたものです.  
さらに、ロジスティック回帰モデルによる"treatment"の予測結果(1, or 0)に応じて色分けされています.    
一方で下図上段右側は、実際の介入有無(treatment)により色分けしたものです.  


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/42bb5810-26c4-51bc-94ac-c39c8ce1766a.png)

上図をご覧のとおり、共変量"history"の分布は、介入有無によってきれいに分離されていません. 
言い換えると、かなりの重複部分があるため、ロジスティック回帰モデルによる予測精度は一言でいうと低そうです.  

関連して、著書内には次のような文があります  

> 傾向スコアはデータに対する説明力が一定を超えることが重要であるという解釈がされ、c統計量のような指標が一定の値を上回ることが望ましいとされていました。しかし、近年では傾向スコアを利用して重み付けかマッチングを行ったあとのデータにおいて、共変量のバランスがとれているかが重要であるという見解が一般的になってきています（Stuart,2010）。安井 翔太. 効果検証入門〜正しい比較のための因果推論／計量経済学の基礎 (Japanese Edition) (p.166). Kindle 版. 

傾向スコアとは今回の場合ロジスティック関数の出力になります.  
そして、c統計量とは何なのでしょう.  
私にとっては初見のワードですが、きっと望ましくない数値が得られるのだと思います.  



色々調べていると、c統計量とAUCの数値が理論上一致するようです.  
そして、いずれもロジスティック回帰結果から計算するみたいでした.  
ですので、まずはロジスティック回帰を実施してみたいと思います. 
そのための関数としてRにはglm()が用意されていました.  
以下は、回帰モデル(変数名model)を作成するためのコードです. 

```R  
model <- glm(
    treatment ~ history, 
    data = biased_data, 
    family = binomial(link = "logit")
)
```
c統計量を計算するためにはDescToolsというライブラリのCstat()なるものが用意されているようです.  

```R  
install.packages("DescTools")
library("DescTools")
c_statistic =Cstat(model, model$fitted.values > 0.5)
c_statistic
# 0.583676635700113
```
結果は上図のとおり約0.58です.

Cstatの引数には回帰モデル(変数名model)と、分類結果(予測)を投入するみたいでした.  
model = glm()そのものには1, or 0の分類結果が保持されていないように見えたので、Bingに聞いたところ、glm()$fitted.values > 0.5 のようにせよとのアドバイスがあり、上記のようにしました.  
(0.5が閾値なのだと理解


さらに、c統計量やAUCはROC曲線なるものと密接な関係を持つようです.  
ですのでROC曲線も作ってみます.  
以下は、ほぼググったもののコピペなので、中身がよく理解できていないRのコードです.  


```R 
install.packages("ROCR")
library("ROCR")
roc <- model %>% 
  predict(type = "response") %>% 
  ROCR::prediction(labels = biased_data$treatment == 1) %>% 
  performance("tpr", "fpr") 

df_roc <- data.frame(tpr = roc@y.values[[1]], fpr = roc@x.values[[1]])

ggplot(df_roc, aes(x = fpr, y = tpr)) + 
  geom_line(color = "red") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  coord_fixed() +
  labs(x = "False Positive Rate", y = "True Positive Rate")

```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/a49cfdf3-31ba-3509-b3aa-612d0eb98577.png)

上図のグラフは赤い線が45度線(破線)に近づくほど、対象のモデルの精度が良くない、ということになるそうです. 
逆に精度が極めて良い場合は、原点からy軸方向に直進し、突き当たりを右折するようなプロットになるようです.  
 
更に、AUCの値の範囲は0.5〜1.0となり、45度線と一致するときは0.5、つまり0.5に近いほど悪い. 
ということになります.  


上述のコードの流れでAUCを計算する方法がありました. 
以下はそのコードです.  
どなたかがインターネットを通じて公開してくれていることに感謝します.  

```R  
auc <- model %>% 
  predict(type = "response") %>% 
  ROCR::prediction(labels = biased_data$treatment == 1) %>% 
  performance("auc")

auc@y.values[[1]]  
# 0.583676635700112
```  
コードの仕組みは良くわからないものの、先ほどのc統計量≒0.58と一致しました.  
どちらかというと1.0よりも0.5に近いのでこのロジスティック回帰モデルの推定精度はあまり良くないということになりますね.  
0.58が非常に悪いc統計量の実現値(?)またはAUCということなのであれば、直感的にもそのとおりだと思います.  

さて、次に共変量のバランスについてです.  
こちらは、決まった可視化の方法があり、cobaltというパッケージで簡単に記述できるみたいでした.  
以下は著書で紹介されていたコードです.  

```R  
install.packages("cobalt")
# (9) 共変量のバランスを確認
## ライブラリの読み込み 
library("cobalt")

## マッチングしたデータでの共変量バランス
love.plot(m_near, 
            threshold = .1)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/96ca500f-9ad9-1422-8222-77599c0e97d3.png)

上図では、傾向スコアマッチングの結果が青い点(Adjusted)で表されています.  
図の見方としては、共変量の標準化平均差が、xy軸でいうところのx = 0 上に整列するのが望ましいとのことでした.  
(distanceが何なのかは不明  

上図から、共変量"history"の標準化平均差は、傾向スコアマッチングの結果0に近づくどころか、Unadjustedの場合よりも高い値になり、平均差が拡大していることが分かります.  

下図は再掲となりますが、確かにグラフのタイトルでいうところのBiasedよりもMatchItの方が介入, 非介入での分布のばらつきが大きいようにも見えます.  
そう見えるのは私の主観ですが、客観的な数値であるc統計量の数値(或いはAUC)では、分類したい2群がほぼ同じであることを警告していますし、標準化平均差を定型化された方法でグラフに可視化すると、傾向スコアマッチングによるバイアスの補正が、少なくともある観点では成功していないことが客観的に分かります.  

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1276668/facc89e2-3cca-77e0-4a41-69b577de834f.png)

最後に、著書内の件のくだりを再び見てみます.  



> 傾向スコアはデータに対する説明力が一定を超えることが重要であるという解釈がされ、c統計量のような指標が一定の値を上回ることが望ましいとされていました。しかし、近年では傾向スコアを利用して重み付けかマッチングを行ったあとのデータにおいて、共変量のバランスがとれているかが重要であるという見解が一般的になってきています（Stuart,2010）。安井 翔太. 効果検証入門〜正しい比較のための因果推論／計量経済学の基礎 (Japanese Edition) (p.166). Kindle 版. 

普通に読むと、c統計量のような指標が一定の値を上回らなくても、共変量のバランスが取れていれば良いという主張だと受け止めることができます.  
今回の1変量による傾向スコアマッチングでは、残念ながらc統計量も共変量のバランスも残念な結果でした.  

# まとめ
- 分類タスクの予測精度が悪いのでは？という疑いに関しては、回帰に使用した1特徴量について2群のヒストグラムを描画するのも良いがc統計量やAUC、そしてROCという標準化された指標計算や可視化が可能.  
- c統計量が0.5に近い1変量による傾向スコアマッチングでは、共変量バランスが元データよりも悪化する事例があった.