傾向スコアによるマッチングを試す。内容は岩波DS3を、コードは以下を参考。

統計的因果推論(2): 傾向スコア(Propensity Score)の初歩をRで実践してみる - 渋谷駅前で働くデータサイエンティストのブログ

発売日: 2016/06/10
メディア: 単行本（ソフトカバー）

傾向スコアとは

ランダム割当（RCT）となっていないような場合、「処置群になるかならないか」というのは各パネルの特徴（共変量）によって決まる傾向にある。この「施策群になる確率」を傾向スコアと呼ぶ。例えば、CM接触効果を測る場合、TVをよく見ている人ほどCM接触をしやすい（処置群になりやすい）ことから「TVをどれだけ見ているか」ということが共変量のひとつとなる。

そのため、処置群ダミーを目的変数とし、共変量を説明変数としてロジスティク回帰やプロビット回帰などをおこなうことで「処置群になるかどうかの確率（の推定値）」である傾向スコアを求めることができる。

この傾向スコアはあくまで傾向、さらにいえば推定値なので実際に施策がされたかどうかを示しているわけではないが、この傾向スコアによって補正をおこなうことで擬似的にRCT状態を組み込んだモデルを作成することができる。

なお、傾向スコアは共変量によって大きく変わるので、傾向スコアが妥当かどうかはc統計量（ROC曲線の下部面積）で判断でき、0.8以上であれば概ね問題ないとされているらしい。

また、確率を求めるためにロジスティック回帰以外も使うことができるが、機械学習によって求めた場合の確率はやや問題がありキャリブレーション（補正）が必要になるらしい。ここはそのうちちゃんと書きます。

統計的因果推論(4): 機械学習分類器による傾向スコアを調整してみる - 渋谷駅前で働くデータサイエンティストのブログ

ATE, ATT, ATU

母集団に対する施策の定量的な平均効果はATE(Average Treatment Effect:平均処置効果/因果効果)と呼ばれる。

また、処置群における平均処置効果はATT(Average Treatment effenct of the Treated)、対照群における平均処置効果はATU(Average Treatment effect of the Untreated)と呼ぶ。

ATTは「処置群に、施策を行った場合と行わなかった場合の差」となるため、その施策の効果測定をおこなうことができ、ATUは「対象群に、施策を行った場合と行わなかった場合の差」となるので、対照群に対して施策をおこなうかどうかの判断に使うことができる。

ATE, ATT, ATUの求め方

ATTは、「処置群に、施策を行った場合と行わなかった場合の差」から求めることができるが、施策を行った場合は行わなかった場合(あるいはその逆)の結果は得られない。そのため、「強い意味での無視可能性(strong ignorability)」という仮定をおくことで推定をおこなう。

この仮定は、
$p(Y_k | X) = p(Y_k | Z =k, X) k=1,0$
つまり「結果変数Yは共変量Xを条件付けると、介入群・対照群の割り付けZとは無関係（独立）で共変量Xのみに依存する」という仮定。この仮定を置くと、
$E(Y_k) = E_X(E(Y_k | X)) = E_X(E(Y_k | Z=k, X))$
となる。

$E_X(E(Y_k | X))$ は観測値から推定できないが、 $E_X(E(Y_k | Z=k, X))$ は推定可能となるので、観測値だけから $E(Y_k)$ を求めることができるため、 $E(Y_1) - E(Y_0) = ATE$ を求めることができる。

この仮定は例えば「Xのみに基づいて処置群に入れるかを判別する(k)」場合は成り立つし、「Yが高い人ほど処置群になりやすい」ような場合、つまりX以外がzを決めている場合は成り立たない。

ATEを求めるための手法は、処置群と対照群でマッチングさせたペアのY差平均を求める方法の他、岩波DS3では4つほど紹介されている。そのうち「IPW(Inverse Probability Weighting:逆確率重み付け法)」と「DR(Doubly Robust:二重にロバストな推定法)」の2つを考える。

IPWは「各々のサンプルに傾向スコアの逆数（対照群の場合は1から傾向スコアを引いた値）で重み付けをかけてATEを補正したもの」で、DRは「各々のサンプルに対してまず共変量による回帰モデルを作って再当てはめによって目的変数の予測値を出し、これと元の目的変数とで傾向スコアの逆数で加重平均させてATEを補正したもの」とのこと。計算式は岩波DS3のp.75参照。

ATE

サンプルデータlalondeを用いる。

まずは傾向スコアを求める

library(Matching)
data(lalonde)

# model作成
model <- glm(treat ~ age + educ + black + hisp + married + nodegr,lalonde,family=binomial)
# fitさせ傾向スコアpsを求める
ps <- model$fitted.values

summary(model)
# all:
#   glm(formula = treat ~ age + educ + black + hisp + married + nodegr, 
#       family = binomial, data = lalonde)
# 
# Deviance Residuals: 
#   Min       1Q   Median       3Q      Max  
# -1.3941  -0.9933  -0.9237   1.3086   1.6633  
# 
# Coefficients:
#   Estimate Std. Error z value Pr(>|z|)   
# (Intercept)  1.053028   1.047384   1.005  0.31471   
# age          0.005917   0.014267   0.415  0.67833   
# educ        -0.063960   0.071354  -0.896  0.37005   
# black       -0.254369   0.363974  -0.699  0.48464   
# hisp        -0.829159   0.504230  -1.644  0.10009   
# married      0.234241   0.266182   0.880  0.37886   
# nodegr      -0.838552   0.309383  -2.710  0.00672 **
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# (Dispersion parameter for binomial family taken to be 1)
# 
# Null deviance: 604.20  on 444  degrees of freedom
# Residual deviance: 589.43  on 438  degrees of freedom
# AIC: 603.43
# 
# Number of Fisher Scoring iterations: 4

head(ps)
# 1         2         3         4         5         6 
# 0.4279364 0.2572819 0.5519755 0.3580844 0.4118540 0.3809886

c統計量はrms::lrmで求まる。

library(rms)
model_lrm <- lrm(treat ~ age + educ + black + hisp + married + nodegr,lalonde)
model_lrm

# Logistic Regression Model
# 
# lrm(formula = treat ~ age + educ + black + hisp + married + nodegr, 
#     data = lalonde)
# 
# Model Likelihood     Discrimination    Rank Discrim.    
# Ratio Test           Indexes           Indexes       
# Obs           445    LR chi2     14.77    R2       0.044    C       0.603    
# 0            260    d.f.            6    g        0.407    Dxy     0.206    
# 1            185    Pr(> chi2) 0.0221    gr       1.502    gamma   0.207    
# max |deriv| 3e-12                         gp       0.097    tau-a   0.100    
# Brier    0.235

c統計量は0.603で基準値0.8に割と満たしてないがいったんこのままやる。

matchingで傾向スコアマッチングで処置群割当をおこない、グループ平均の差で平均処置効果を測定する。

# 共変量Xにpsを指定しているので傾向スコアが近い人がマッチングされる（傾向スコアマッチング）だが、ageとかを使うとageなどの属性をもとにしたマッチングとなる
#
match_ps <- Match(Y = lalonde$re78, Tr=lalonde$treat, X =ps)
summary(match_ps)

# Estimate...  2302 
# AI SE......  764.97 
# T-stat.....  3.0093 
# p.val......  0.0026187 
# 
# Original number of observations..............  445 
# Original number of treated obs...............  185 
# Matched number of observations...............  185 
# Matched number of observations  (unweighted).  492

傾向スコアマッチング推定だと、ATEは2302。これは対照群/処置群のマッチングペア毎にYの差分を出し、その平均値を出している。

次に、IPW推定を試す。

# IPW
y <- lalonde$re78
z1 <- lalonde$treat

ipwe1 <- sum((z1*y)/ps)/sum(z1/ps)

ipwe0 <- sum(((1-z1)*y)/(1-ps))/sum((1-z1)/(1-ps))

ipwe1 - ipwe0
# 1632.463

ATEのIPW推定量は1632と傾向マッチング推定の2302より小さく出た。また、直接周辺期待値を計算しないで、傾向スコアの逆数で重み付け（補正）をおこなって回帰してtreatの効果をみることでも計算は可能。
この場合は、周辺期待値だけでなく標準誤差なども合わせて求められる。

このときのweightは

$\frac{Z_i}{e_i} + \frac{(1 - Z_i)}{1 - e_i}$

式は以下のp.408

http://www.cardio.med.tohoku.ac.jp/newmember/pdf/ms/hmbr46-399.pdf

y <- lalonde$re78
ivec1 <- lalonde$treat # treatダミー
ivec0 = rep(1, nrow(lalonde))-ivec1 # untreatダミー
ivec <- cbind(ivec1,ivec0)

iestp1 <- z1/ps # treatの場合の逆確率重み付け(untreatの場合0)

iestp0 <- (1-z1)/(1-ps) # untreatの場合の逆確率重み付け(treatの場合0)

iestp <- iestp1+iestp0 # 片方は0になるのでtreat,untreatを区別せず逆確率重み付け値

ipwe_gs <- lm(y ~ ivec + 0, weights=iestp) # 逆確率重み付けをした加重最小二乗法。ivec1(treatダミー)=0の場合は0*重み=0になる
summary(ipwe_gs)

# Call:
#   lm(formula = y ~ ivec + 0, weights = iestp)
# 
# Weighted Residuals:
#   Min     1Q Median     3Q    Max 
# -10534  -5886  -2389   4072  90256 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# ivecivec1   6196.0      454.7   13.63   <2e-16 ***
#   ivecivec0   4563.5      454.1   10.05   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 9584 on 443 degrees of freedom
# Multiple R-squared:  0.3929,  Adjusted R-squared:  0.3902 
# F-statistic: 143.4 on 2 and 443 DF,  p-value: < 2.2e-16

6196.0 - 4563.5 = 1632.5となり、直接求めた場合と同じ(丸めているので下2桁は不明だが）。

次にDR推定。

# DR
dre <- function(data, target, treat, ps, formula) {
  n       <- nrow(data)
  y       <- data[target]
  data1   <- data[data[treat]==1,]
  data0   <- data[data[treat]==0,]
  model1  <- lm(formula=formula, data=data1)
  model0  <- lm(formula=formula, data=data0)
  fitted1 <- predict(model1, data)
  fitted0 <- predict(model0, data)
  dre1    <- (1/n)*sum(y+((z1-ps)/ps)*(y-fitted1))
  dre0    <- (1/n)*sum(((1-z1)*y)/(1-ps)+(1-(1-z1)/(1-ps))*fitted0)
  return(c(dre1, dre0))
  }

 ret <- dre(lalonde, "re78", "treat", ps,
               +            re78 ~ age + educ + black + hisp + married + nodegr)
ret
# 6185.915 4566.401

ret[1] - ret[2]
# 1619.514

ATEのDR推定量は1620、IPW推定量の1632、傾向マッチング推定の2302よりも小さい。

処置群の中には処置群なのに傾向スコアが低い人、逆に対照群なのに傾向スコアが高い人がいる。傾向スコアによる逆確率重み付けはこういう状況を補正する。そのため単に処置群/対照群の平均差を出した傾向マッチング推定と、傾向スコアで補正をおこなったDR推定量、IPW推定量で結果が大きく異なっている。

なお、CBPSパッケージに様々な関数があるとのこと。

ATT

引き続き、「処置群に、施策を行った場合と行わなかった場合の差」であるATTを求める。

$ATT = E(Y_1 | Z = 1) - E(Y_0 | Z = 1)$

このときのweightは

$Z_i + \frac{(1 - Z_i) e_i }{1 - e_i}$

式はATE同様に以下のp.409

http://www.cardio.med.tohoku.ac.jp/newmember/pdf/ms/hmbr46-399.pdf

コードが微妙にわかるようなわからんような。。。

model <- glm(treat ~ age + educ + black + hisp + married + nodegr,lalonde,family=binomial)
ps <- model$fitted.values # fitさせ傾向スコアpsを求める

ivec1 <- lalonde$treat  # treatダミー
iestp1_ATT <- ivec1  #treatの場合の重み（そのまま）

ivec0 <- rep(1, nrow(lalonde)) - lalonde$treat # untreatダミー
iestp0_ATT <- ivec0 * ps / (1-ps)  #untreatの場合の重み

ivec <- cbind(ivec1,ivec0)
iestp_ATT <- iestp1_ATT + iestp0_ATT
ipwe_ATT_re74 = lm(lalonde$re74 ~ ivec + 0, weights=iestp_ATT)

summary(ipwe_ATT_re74)
# Call:
#   lm(formula = lalonde$re74 ~ ivec + 0, weights = iestp_ATT)
# 
# Weighted Residuals:
#   Min     1Q Median     3Q    Max 
# -2853  -2096  -1753  -1180  32945 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# ivecivec1   2095.6      363.4   5.766 1.52e-08 ***
#   ivecivec0   2226.0      362.9   6.134 1.90e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 4943 on 443 degrees of freedom
# Multiple R-squared:  0.1379,  Adjusted R-squared:  0.134 
# F-statistic: 35.44 on 2 and 443 DF,  p-value: 5.29e-15