çµ±è¨ˆçš„å› æžœæŽ¨è«–(2): 傾向スコア(Propensity Score)の初歩をRで実践してみる

さて、統計的因果推論についてだらだらと独習していくこのシリーズですが、今回はDonald Rubinが考案したRubinの因果モデルで用いられる、傾向スコア(Propensity Score)を取り上げてみようと思います。「お前岩波DS3で事前に原稿読んで中身は知っているはずじゃないのか」とかいうツッコミはご勘弁ください（笑）。

元々は観察データ（つまりRCTを含む実験データではない）から因果関係を推定するための手法ということで、いかにして観察データに隠れた影響を与える共変量を突き止め、その共変量から及んでくる影響をバランスさせ、真の因果効果を推定するか、というのが主眼でした。つまり、RCTを実践できないような疫学データや社会科学的データに対する適用がメインだったようです。

しかしながら、最近は例えば広告やマーケティングといった「ある程度の介入（処置）*1はできても事実上RCTを徹底することは不可能」な調査データにおいても同様のニーズが増えているとも聞きます。岩波DS3にもそのような事例の紹介が実際にされており、今後も同様の取り組みは広告・マーケティング領域では増えていくものと予想されます。

今回は僕自身の勉強が決定的に不足しているため、基本的にはほぼ全面的に他の資料を参考にしながら備忘録的に内容をまとめたに過ぎない内容に終始しています。そのため大半の数式・定義類は割愛しております。きちんとした解説を読みたいという方は岩波DS3をお求めいただくか、バント効果推定記事をお書きいただいた中村知繁さんのブログ記事をお読みいただくことをお薦めいたします。

岩波データサイエンス Vol.3

そしてRでの実践部分については以下のid:isseing33さん＆id:SAMさんのブログ記事を参考にしています。

http://d.hatena.ne.jp/isseing333/20110511/1305124310

Disclaimer

いつも通り、調べながらの記事なので盛大に間違っている可能性があります。間違いを見つけた方は盛大に突っ込んで下さると有難いです。。。

傾向スコアとは

実は僕は岩波DS3の原稿を校閲している時からずーーーっと分かってなかったんですが、要は「事前に別の説明変数（共変量）で介入群or対照群の割り付けを目的変数としてロジスティック回帰し、そこから得られた「共変量から勘案して理論的に推測される各群に属するであろう確率値」（再当てはめ予測値）を新たな説明変数である傾向スコア』として用意し、その上で本来の目的変数をターゲットとし群割り付け＆傾向スコアを説明変数として改めて回帰する」という因果推論の枠組みの中で用いられるものなんですね。

言い換えると「介入群or対照群の割り付けは外部の人間がある程度恣意的に決めたものに過ぎず、実際には様々な共変量によって『実際にはこう割り付けられるべき』という傾向がその裏側に存在するので、その傾向を『傾向スコア』としてモデルに組み込むことで補正する」ということなのかなと。その上で、本来の目的変数をモデリングして再当てはめした結果が潜在的結果変数という「真の」目的変数になる、ということのようです。

（※ただし本来ならこのような傾向スコアによって補正されるべき偏りはあっては困るわけで、真に理想的にきちんとバランスされたRCTであれば傾向スコアは全て0.5になるということも言えます）

傾向スコアがどのようにして効果を発揮するかという様子を僕なりにイメージして図にしたものが、以下です。

こうして共変量に基づいてロジスティック回帰をかけ、実際に介入群・対照群に対して「共変量から推測される」割り付け確率を「傾向スコア」として算出し、これを説明変数としてモデルに組み込む、ということで（"strong ignorability"の仮定、すなわち「介入群・対照群の割り付けと潜在的結果変数とが共変量を条件付けると独立である」という仮定が成り立つならば）多次元の共変量全てに対して調整を行ったのと同じ結果が得られることになる、そうです。

これにより、多くの共変量からの影響を1つの変数にまとめ上げることになり、目的変数が複数存在するケース（例としてはCMに接触したことで影響を受ける「商品認知度」「利用意向」「実際の購買行動」など）で毎回モデリングのたびに共変量を調整したりせずとも、傾向スコアだけで補正すれば良くなる、というメリットがあるとのことです。下の方で紹介しているCM接触効果データセットも、異なる3つの目的変数に対して全て1つの傾向スコアによる調整で事足りています。

とは言え、これだけ傾向スコアが共変量に強く依存しているところからも分かるように、共変量の選択が傾向スコアの良し悪しを決めるということも言えます。その選択基準としてc統計量（＝ROCカーブの下側面積）なるものが提唱されていてこれが0.8を超えていればOKというように医学・疫学系では言われているとのことですが、機械学習の文脈だとこれはk-folds CVとかで見た方が良いのではないかなぁという気も。。。この辺の話はまた記事を改めて取りあげようと思います。

なお、岩波DS3のp.74でもコメントされているように、傾向スコアの算出はロジスティック回帰に限る必要はなく、何がしかの機械学習2クラス分類器でも良いようです。陽にクラス分類事後確率を出せるものであればOKと見て良さそうですね*2。

IPW推定量とDR推定量

前回のDID（差分の差分法）でも間接的に触れた話題ですが、基本的にこの手の因果推論で求めたいものというのは端的に言えば「『介入した場合』と『介入しなかった場合』とでの（目的変数への）効果の差」であり、その定量的な大きさです。岩波DS3のp.65で触れられているように、これを因果効果とか平均処置効果(Average Treatment Effect: ATE)と呼びます。因果推論は、このATEがその他の共変量によって歪められるのを防ぐための営みであるとも言えます。

で、傾向スコアで補正したATEを求めるには、岩波DS3で紹介されている範囲では2つ方法があります。一つは逆確率重み付け法(Inverse Probability Weighting: IPW)によるもの、もう一つは二重にロバストな推定法(Doubly Robust: DR)によるものです。

どちらも多分定義式を見るよりはその字面を見た方が分かりやすくて、要はIPWは「各々のサンプルに傾向スコアの逆数（対照群の場合は1から傾向スコアを引いた値）で重み付けをかけてATEを補正したもの」で、DRは「各々のサンプルに対してまず共変量による回帰モデルを作って再当てはめによって目的変数の予測値を出し、これと元の目的変数とで傾向スコアの逆数で加重平均させてATEを補正したもの」です。なお具体的な式は岩波DS3のp.75に全て記載されているのでここでは割愛します。以下のR実行例にも出てきますが、回帰モデルの計算で代用することも可能です。

これによって、例えば岩波DS3のpp.91-100で紹介されているような「共変量（ここでは各リサーチ参加者＝パネルのデモグラフィック情報など）が影響しているせいで常識に反してCMの効果＝ATEがマイナスになって見える」ような現象を補正できるようになる、というわけです。

ちなみに、ATE以外にも「介入群における平均介入効果(Average Treatment Effect on the Treated: ATT)」「対照群における平均介入効果(Average Treatment Effect on the Untreated: ATU)」という概念もあり、広告・マーケティングの文脈で言えば「ターゲティングされた（orされなかった）人々に対してプロモーションを行った場合と行わなかった場合との差の大きさ」を表すものなので、ROI（投資収益率）など効果の大きさを知りたい場合はこちらの方がより実用的かもしれません。

Rでサンプルデータセットに対して傾向スコアを用いた因果推論をやってみる：lalondeデータセット

それでは、実際にRで傾向スコアを用いた因果推論をやってみましょう。まずはド定番でRに最初から入っている傾向スコア実践用のデータセットである"lalonde"を用い、その後岩波DS3で用いられているCM接触データセットを用いることとします。

lalondeデータセットですが、{Matching}パッケージから以下の通り読み込めます。

> library(Matching)
## 
##  Matching (Version 4.9-2, Build Date: 2015-12-25)
##  See http://sekhon.berkeley.edu/matching for additional documentation.
##  Please cite software as:
##   Jasjeet S. Sekhon. 2011. ``Multivariate and Propensity Score Matching
##   Software with Automated Balance Optimization: The Matching package for R.''
##   Journal of Statistical Software, 42(7): 1-52. 
##
> data(lalonde)

パッケージ読み込み時の説明に既に書いてあるように、このパッケージ自体が実は傾向スコアマッチングのためのもので、"lalonde"はその中の代表的なデータセットです。ヘルプを見ると以下のような説明がなされています。

<b>Format</b>

A data frame with 445 observations on the following 12 variables.

age
    age in years.

educ
    years of schooling.

black
    indicator variable for blacks.

hisp
    indicator variable for Hispanics.

married
    indicator variable for martial status.

nodegr
    indicator variable for high school diploma.

re74
    real earnings in 1974.

re75
    real earnings in 1975.

re78
    real earnings in 1978.

u74
    indicator variable for earnings in 1974 being zero.

u75
    indicator variable for earnings in 1975 being zero.

treat
    an indicator variable for treatment status.

<b>Details</b>

Two demos are provided which use this dataset. The first, DehejiaWahba, replicates one of the models from Dehejia and Wahba (1999). The second demo, AbadieImbens, replicates the models produced by Abadie and Imbens in their Matlab code. Many of these models are found to produce good balance for the Lalonde data.

ここではre78（1978年における実際の収入額）を最終的な目的変数にしましょう。で、まずは傾向スコアを算出します。id:SAMさんの例では全ての共変量を突っ込んでましたが、今回は収入額系は除外してデモグラフィック情報だけを使うことにしてみます。

> model <- glm(treat~age+educ+black+hisp+married+nodegr,lalonde,family=binomial)
> ps <- model$fitted.values
> summary(model)

Call:
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 

ここで得られた傾向スコアpsと群割り付けであるlalonde$treatとを説明変数として、目的変数であるlalonde$re78を線形回帰してみます。

> re.model <- lm(lalonde$re78~lalonde$treat+ps)
> summary(re.model)

Call:
lm(formula = lalonde$re78 ~ lalonde$treat + ps)

Residuals:
   Min     1Q Median     3Q    Max 
 -7173  -4502  -1829   2796  54249 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)     3009.3     1480.5   2.033  0.04269 * 
lalonde$treat   1667.7      643.4   2.592  0.00986 **
ps              3844.1     3540.0   1.086  0.27812   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6578 on 442 degrees of freedom
Multiple R-squared:  0.02044,	Adjusted R-squared:  0.016 
F-statistic: 4.611 on 2 and 442 DF,  p-value: 0.01043

正直言ってあまりモデルの説明力は高くない気がしますが、群割り付けによる影響はp < 0.05で有意かつ、その偏回帰係数は1667.7ということで正の値を示しています。すなわち介入によるre78の向上が見られた、ということになります。

ただし、既に見たようにこれだけでは実際に共変量を考慮して補正したATEを算出することはできません。以下のようにIPW推定量とDR推定量を求めてみましょう。なお、信頼区間の算出はここでは割愛します。

# IPW
> y <- lalonde$re78
> z1 <- lalonde$treat
> (ipwe1 <- sum((z1*y)/ps)/sum(z1/ps))
[1] 6195.964
> (ipwe0 <- sum(((1-z1)*y)/(1-ps))/sum((1-z1)/(1-ps)))
[1] 4563.501
> ipwe1 - ipwe0
[1] 1632.463

# DR
# SAMさんの例にならって関数化してあります
> 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
[1] 6185.915 4566.401
> ret[1] - ret[2]
[1] 1619.514

IPW / DW推定量とも、若干id:SAMさんの結果よりも大きくなりました。また、先に線形回帰で求めた偏回帰係数1667.7よりも、補正済みATE自体の大きさはやや小さくなっていることが分かります。

Rでサンプルデータセットに対して傾向スコアを用いた因果推論をやってみる：岩波DS3のCM接触データセット

データセットの概要が岩波DS3のpp.92-3に載っていますが、要はいわゆるシングルソースパネル(SSP)*3市場調査（マーケティングリサーチ）データです。スマホゲームアプリのプレイ状況が、ある当該期間中に放映したCMに接触（視聴）したか否かによって、どれくらい変化したかを分析することでCMの広告効果を推定しようというのがここでの課題です。内訳としてはCMへの接触の有無、ゲームプレイの有無・回数・秒数、TVの視聴秒数といった行動履歴から、年齢・性別・家族構成・居住地・所得・学歴といったデモグラフィック＆ソシオグラフィック属性を含んでいます。

本文中でも触れられていますが、特にwebによるデジタル化が進行している昨今ではCMに接触するorしないという「結果的に生じる群割り付け」自体が既に様々な環境要因に振り回されることが分かっています。例えば、スマホゲームは高齢者ほどプレイしない傾向が強い一方で高齢者ほどTVはよく視聴しているため、この時点で既に大きなバイアスが存在することが予想されます。また、そもそもこの手のマーケティングリサーチデータによくあるようにゲーム利用者は少数派なので、この点による不均衡データとしてのバイアスも発生します。

結果として本文pp.93-4及び表2にも記されているように、単純にゲームプレイ秒数を比べると何故かCM接触群の方が非接触群よりも短くなる、という謎の逆転現象が起きてしまっています。これを何とかして正しく評価できないか？というのがこの事例での最大のポイントとなります。

まず元データですが、CSVファイルがGitHubに上がっているのでこちらを使います。

普通にダウンロードして以下のようにインポートしておきましょう。

> d <- read.csv('https://github.com/iwanami-datascience/vol3/raw/master/kato%26hoshino/q_data_x.csv')
> head(d)

試しに表2を再現してみるとこうなります。

> mean(d$gamesecond[d$cm_dummy==1])
[1] 2478.066
> mean(d$gamesecond[d$cm_dummy==0])
[1] 3107.706

確かに、素のままだとCMに接触した方がそうでない場合よりも何故かゲームプレイ秒数が短いという結果になっています。これはいくら何でも直感に反するので、傾向スコアで補正していくということを考えたいと思います。このデータに対するRによる傾向スコアを用いた分析例が既にGitHubに上がっていますが、ここではそれをなぞって写経する形でやってみます。

傾向スコア

まず、傾向スコアを出します。ここではCM接触群か否かが重要なので、例に従ってロジスティック回帰で傾向スコアを求めていきます。

> model <- glm(cm_dummy ~ TVwatch_day + age + sex + marry_dummy + child_dummy + inc + pmoney + area_kanto +area_tokai + area_keihanshin + job_dummy1 + job_dummy2 + job_dummy3 + job_dummy4 + job_dummy5 + job_dummy6 + job_dummy7  + fam_str_dummy1 + fam_str_dummy2 + fam_str_dummy3 + fam_str_dummy4, d, family=binomial)
> ps <- model$fitted.values
> head(ps)
         1          2          3          4          5          6 
0.04649543 0.25479069 0.18673620 0.22727078 0.24095033 0.15785738 
> head(d$cm_dummy)
[1] 0 0 0 0 0 0

これで傾向スコアpsが求まりました。ちなみにc統計量も出すだけなら出せて、{rms}パッケージのlrm関数で以下のように回します。

> library(rms)
> model_lrm <- lrm(cm_dummy ~ TVwatch_day + age + sex + marry_dummy + child_dummy + inc + pmoney + area_kanto +area_tokai + area_keihanshin + job_dummy1 + job_dummy2 + job_dummy3 + job_dummy4 + job_dummy5 + job_dummy6 + job_dummy7  + fam_str_dummy1 + fam_str_dummy2 + fam_str_dummy3 + fam_str_dummy4, d)
> model_lrm

Logistic Regression Model

lrm(formula = cm_dummy ~ TVwatch_day + age + sex + marry_dummy + 
    child_dummy + inc + pmoney + area_kanto + area_tokai + area_keihanshin + 
    job_dummy1 + job_dummy2 + job_dummy3 + job_dummy4 + job_dummy5 + 
    job_dummy6 + job_dummy7 + fam_str_dummy1 + fam_str_dummy2 + 
    fam_str_dummy3 + fam_str_dummy4, data = d)
                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       
Obs         10000    LR chi2    2726.06    R2       0.321    C       0.792    
 0           5856    d.f.            21    g        1.509    Dxy     0.583    
 1           4144    Pr(> chi2) <0.0001    gr       4.521    gamma   0.585    
max |deriv| 9e-06                          gp       0.275    tau-a   0.283    
                                           Brier    0.183                     

                Coef    S.E.   Wald Z Pr(>|Z|)
Intercept       -1.7709 0.2615  -6.77 <0.0001 
TVwatch_day      0.0001 0.0000  31.99 <0.0001 
age             -0.0026 0.0029  -0.89 0.3750  
sex              0.0006 0.0647   0.01 0.9927  
marry_dummy     -0.0781 0.0856  -0.91 0.3610  
child_dummy      0.3142 0.0743   4.23 <0.0001 
inc             -0.0005 0.0002  -2.92 0.0035  
pmoney           0.0119 0.0077   1.54 0.1227  
area_kanto       0.4050 0.0790   5.12 <0.0001 
area_tokai      -0.7233 0.0765  -9.45 <0.0001 
area_keihanshin -2.0420 0.0756 -27.01 <0.0001 
job_dummy1       0.1752 0.1562   1.12 0.2619  
job_dummy2       0.1651 0.1677   0.98 0.3250  
job_dummy3       0.5399 0.1604   3.37 0.0008  
job_dummy4       0.3604 0.2433   1.48 0.1386  
job_dummy5       0.6414 0.1520   4.22 <0.0001 
job_dummy6       0.2848 0.1581   1.80 0.0717  
job_dummy7       0.1540 0.1834   0.84 0.4013  
fam_str_dummy1   0.7640 0.2045   3.74 0.0002  
fam_str_dummy2   1.0033 0.2176   4.61 <0.0001 
fam_str_dummy3   0.6137 0.2021   3.04 0.0024  
fam_str_dummy4   0.1799 0.2216   0.81 0.4167  

本文中にも出ていますが、c統計量は0.792ということで概ね問題なしということになりました。もっともこれを（例えば）ランダムフォレストとかでやってみたらどうなるのかなぁという気はしますが。。。

ゲームをプレイしたか否かへのCM接触の影響

GitHubの例ではこの後色々分析がケース分けされていきますが、一旦ここではCM接触がゲームをプレイしたか否かに影響したかを見てみましょう。

> summary(glm(as.factor(d$gamedummy)~d$cm_dummy+ps, family=binomial))

Call:
glm(formula = as.factor(d$gamedummy) ~ d$cm_dummy + ps, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.4535  -0.4086  -0.3903  -0.3691   2.3977  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -2.35915    0.07406 -31.854  < 2e-16 ***
d$cm_dummy   0.17658    0.08907   1.982  0.04744 *  
ps          -0.59756    0.18533  -3.224  0.00126 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 5277.3  on 9999  degrees of freedom
Residual deviance: 5266.5  on 9997  degrees of freedom
AIC: 5272.5

Number of Fisher Scoring iterations: 5

CMに接触すると、ゲーム（アプリ）をプレイするという因果効果があることが分かります*4。以下、ゲームプレイ秒数にターゲットを絞って分析を進めていきます。その他のゲームプレイ有無・ゲームプレイ回数のATT推定については、岩波DS3本書並びにGitHubのR実行例をご覧ください。

ATEの推定

まず、岩波DS3の例に従って最初にATEを推定します。ゲームプレイ秒数のATEは以下のようにして求まります。ちなみにid:SAMさんの記事でも微妙に触れられていますが、目的変数と介入ダミー変数とを傾向スコアで重み付け線形回帰（ただし切片を抜く）した時の偏回帰係数としても、ATEは求まります。こちらの方が標準偏差も同時に算出できるので便利ですね。

> ivec1 <- d$cm_dummy # Treated group
> ivec0 <- rep(1, nrow(d))-ivec1 # Untreated group
> ivec <- cbind(ivec1,ivec0)
> iestp1 <- (ivec1/ps)
> iestp0 <- (ivec0/(1-ps))
> iestp <- iestp1+iestp0
> ipwe_gs <- lm(d$gamesecond ~ ivec+0, weights=iestp)
> summary(ipwe_gs)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-25199  -5328  -3601  -3011 531077 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
ivecivec1   4143.3      277.3   14.94   <2e-16 ***
ivecivec0   2639.4      262.4   10.06   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 27220 on 9998 degrees of freedom
Multiple R-squared:  0.03143,	Adjusted R-squared:  0.03124 
F-statistic: 162.2 on 2 and 9998 DF,  p-value: < 2.2e-16

岩波DS3の本文表4に示されているのと同じように、CM接触群のゲームプレイ時間が4143.3秒、非接触群では2639.4秒というATEが得られました。

ATTの推定

一方、本文でも指摘されているようにマーケティング上重要なのはむしろ「CMに接触したらどれくらいゲームプレイ時間が伸びるか」というATTの方なので、こちらも求めてみます。

> iestp1_ATT <- ivec1
> iestp0_ATT <- ivec0*ps/(1-ps)
> iestp_ATT <- iestp1_ATT+iestp0_ATT
> ipwe_ATT_gs = lm(d$gamesecond ~ ivec+0, weights=iestp_ATT)
> summary(ipwe_ATT_gs)

Call:
lm(formula = d$gamesecond ~ ivec + 0, weights = iestp_ATT)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-19751  -2478  -1918  -1143 220618 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
ivecivec1   2478.1      227.4  10.897   <2e-16 ***
ivecivec0   2080.3      209.0   9.952   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14640 on 9998 degrees of freedom
Multiple R-squared:  0.02132,	Adjusted R-squared:  0.02112 
F-statistic: 108.9 on 2 and 9998 DF,  p-value: < 2.2e-16

CM接触群におけるゲームプレイ秒数に対するATTは約400秒ぐらいのプラス、という結果に至りました。

傾向スコアでCM接触効果を調整した上で、パネル属性が与える影響を推定する

次に、傾向スコアで調整した線形回帰モデルを推定し、その結果から「どのようなデモグラフィック属性がCMへの接触の有無に関係なくゲームプレイ秒数に影響を受けないか」を特定することを目指します。

> idx1 <- which(d$cm_dummy==1)
> idx0 <- which(d$cm_dummy==0)
> ME_tr_gs <- lm(gamesecond ~ child_dummy + area_kanto +area_tokai + area_keihanshin + T + F1 + F2 + F3 + M1 + M2 ,weights=(1/ps[idx1]) ,data=d[d$cm_dummy==1,])
> ME_un_gs <- lm(gamesecond ~ child_dummy + area_kanto +area_tokai + area_keihanshin + T + F1 + F2 + F3 + M1 + M2,weights=(1/(1-ps[idx0])), data=d[d$cm_dummy==0,])
> summary(ME_tr_gs)

Call:
lm(formula = gamesecond ~ child_dummy + area_kanto + area_tokai + 
    area_keihanshin + T + F1 + F2 + F3 + M1 + M2, data = d[d$cm_dummy == 
    1, ], weights = (1/ps[idx1]))

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-64781 -10218   -718   1582 470814 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)        498.6      819.6   0.608  0.54298    
child_dummy      -1938.9      666.3  -2.910  0.00363 ** 
area_kanto       -2729.2     1083.7  -2.518  0.01182 *  
area_tokai       -3030.9      996.9  -3.040  0.00238 ** 
area_keihanshin   7833.9      807.3   9.704  < 2e-16 ***
T                 6413.3     2904.3   2.208  0.02728 *  
F1               -1042.6     1221.1  -0.854  0.39326    
F2                1142.6     1047.0   1.091  0.27521    
F3                -265.1     1539.0  -0.172  0.86326    
M1                1247.3     1122.2   1.111  0.26646    
M2                9624.8      929.9  10.350  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30060 on 4133 degrees of freedom
Multiple R-squared:  0.07586,	Adjusted R-squared:  0.07362 
F-statistic: 33.93 on 10 and 4133 DF,  p-value: < 2.2e-16

> summary(ME_un_gs)

Call:
lm(formula = gamesecond ~ child_dummy + area_kanto + area_tokai + 
    area_keihanshin + T + F1 + F2 + F3 + M1 + M2, data = d[d$cm_dummy == 
    0, ], weights = (1/(1 - ps[idx0])))

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-17693  -4250  -2398   -932 411626 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       631.86     596.66   1.059 0.289648    
child_dummy      1827.26     511.62   3.572 0.000358 ***
area_kanto        -28.27     846.67  -0.033 0.973368    
area_tokai         41.52     646.05   0.064 0.948758    
area_keihanshin -1091.39     607.28  -1.797 0.072358 .  
T                7141.70    2363.25   3.022 0.002522 ** 
F1                382.83     858.52   0.446 0.655671    
F2                944.73     874.04   1.081 0.279797    
F3               -282.52    1035.28  -0.273 0.784942    
M1               7164.78     847.54   8.454  < 2e-16 ***
M2                981.14     699.32   1.403 0.160673    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23780 on 5845 degrees of freedom
Multiple R-squared:  0.01878,	Adjusted R-squared:  0.0171 
F-statistic: 11.19 on 10 and 5845 DF,  p-value: < 2.2e-16

岩波DS3の表6と同じ結果になっているはずです。この結果からは「子供の有無」(child_dummy)の効果がCM接触群では負に、CM非接触群では正に作用していることが分かります。つまり「子供がいる層」にはCMは影響を与えない（何故ならCMに接触しているにもかかわらず本来想定されるプラスの効果に達せず、逆に接触しなかったにもかかわらず本来想定されるビハインドがない）ことが推測されます。

「子供がいる層」を除外した場合のATTの推定

そこで駄目押しとして、「子供がいる層」を除外してもう一度ATTを計算し直します。

> d2 <- subset(d, d$cm_dummy==0 | d$child_dummy==0)
> ivec1_ltd <- d2$cm_dummy
> ivec0_ltd <- rep(1, nrow(d2))-ivec1_ltd
> model2 <- glm(cm_dummy ~ TVwatch_day + age + sex + marry_dummy + child_dummy + inc + pmoney  + area_kanto +area_tokai + area_keihanshin + job_dummy1 + job_dummy2 + job_dummy3 + job_dummy4 + job_dummy5 + job_dummy6 + job_dummy7  + fam_str_dummy1 + fam_str_dummy2 + fam_str_dummy3 + fam_str_dummy4, family=binomial, data = d2)
> ps2 <- model2$fitted
> 
> ivec_ltd <- cbind(ivec1_ltd,ivec0_ltd)
> iestp1_ltd <- ivec1_ltd
> iestp0_ltd <- ivec0_ltd*ps2/(1-ps2)
> iestp_ltd <- iestp1_ltd+iestp0_ltd
> ipwe_ltd_gs <- lm(d2$gamesecond ~ ivec_ltd+0, weights=iestp_ltd)
> summary(ipwe_ltd_gs)

Call:
lm(formula = d2$gamesecond ~ ivec_ltd + 0, weights = iestp_ltd)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-11639  -2430  -1136      0 226619 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
ivec_ltdivec1_ltd   2430.0      253.6   9.583  < 2e-16 ***
ivec_ltdivec0_ltd   1880.3      241.0   7.801 6.87e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12380 on 8239 degrees of freedom
Multiple R-squared:  0.0182,	Adjusted R-squared:  0.01796 
F-statistic: 76.35 on 2 and 8239 DF,  p-value: < 2.2e-16

「子供がいる層」を除外したところ、CM接触群のゲーププレイ秒数は2430.0、非接触群は1880.3という結果になりました。これは全ての層で算出した値よりそれぞれ僅かに-48, -200ぐらいの差しかなく、接触群でのゲームプレイ秒数はほとんど変わらない上にむしろ2群の差は広がっています。

これにより、この例では「子供がいる層へのCMのターゲティングを取り止めてもCM接触群に与えられるゲームプレイ秒数増大効果はむしろ大きくなる」ことが言えるわけで、CM出稿予算をより「子供がいない層」へと多く振り分けることでさらにCM効果が高まることすら期待されます。

このように、傾向スコアを用いて適切にマーケティングリサーチデータを分析することで、冒頭に見たような不適切な推論に惑わされることなく、より正確かつ効果的な分析結果を得ることができるというわけです。

感想

ということで、遅れ馳せながら傾向スコアについて学んでみました。実はCM接触効果データセットと全く同様の逆転現象を昔の職場で見かけたことがあるので、あの時もこれを使った方が良かったのかなぁと今更ながら反省しております（汗）。

そして傾向スコアを求めるモデリングのところなんですが、ここは機械学習でも何でもOKという話なので、それこそ例えばランダムフォレストとか非線形分類器を使うとか、L1正則化で汎化性能を上げるとか、クラス重み付けで不均衡データを補正するとか、やってみたらいいのかなぁと思ったのでした。

後は超絶細かいことなんですが、岩波DS3のGitHubのRスクリプトが僕個人の好みとは違うので適当に書き直しました（笑）。この辺は単なる好みの問題なので、何卒ご容赦を。。。

*1:僕個人はこれまで「処置」という語は使ってこなかったので、以後「介入」に基本的には統一しています

*2:ただし原理的にはそんな非線形分類器的なものまで投入しなければいけない時点で何かがまずいとも言えそう

*3:マーケティングリサーチ会社が有償で募集したパネル＝ボランティアに様々なデバイスを用いて行動を記録＆デモグラフィック・ソシオグラフィック属性情報を提供してもらうことで、広告などのマーケティングの効果を測定するための枠組みです

*4:ただしこの計算は不均衡データによるものなので、クラス重み付けをかけると実はp = 0.1864となって有意でなくなるという。。。