BUSINESS DATA SCIENCEの続き。

データなどは作者のgitにある。

Business Data Science: Combining Machine Learning and Economics to Optimize, Automate, and Accelerate Business Decisions

作者:Taddy, Matt
発売日: 2019/08/21
メディア: ハードカバー

最近日本語版も出たみたいです

ビジネスデータサイエンスの教科書

作者:マット・タディ
発売日: 2020/07/22
メディア: 単行本

この章はなにか

5章*1では実験環境下、つまりRCTの場合の因果効果を測る話だったが6章はそうでない場合、conditional ignorabilityを用いることで、効果を測りたいdのyへの共変量XでControlをおこなうことで、yへのdだけの効果を取り出す。つまり、共変量Xを用いることで擬似的にRCTと同様の環境にする。

5章・6章(の一部)の話は「効果検証入門」のときにもやっている

knknkn.hatenablog.com

本書ではそこから一歩先に進んで、共変量Xを選ぶのは難しいので自動的にLassoでXを選択する方法を紹介している。また、全体的な平均処置効果Average Treatment Effect(ATE)だけではなくある要素毎での処置効果(例えば、性別毎での処置の違いや年収ごとなど）となるHeterogeneous Treatment Effect(HTE)の紹介もおこなっている。

Conditional Ignorability

Conditional Ignorability(以下CI)は以下の式で表される。

$｛y_i(d) ∀ d｝⊥ d_i | x_i$

これは、共変量xが与えられたときに、効果を検証したい項目d(バイナリの場合は割当)はどのdの場合のyの同時分布とも独立している ことをいう。つまり、dの値は共変量xによってのみ決まり、yによっては決まらない。言い換えると、xが同じならdは同じ値になる。そのため、dとyのストレートな関係を得られるといえるし、RCTと同様になるためセレクションバイアスがかからないともいえる。

書き方を変えると、 $p(d_i | y_i(d), x_i) = p(d_i | x_i)$ および $p( y_i(d) | d, x_i) = p( y_i(d) | x_i)$ となる。

なお、これは「効果検証入門」でConditional Independence Assumption(CIA)が書かれていたが（IgnorabilityとIndependenceの違いはあれど）CIAをd(Z)=1/0(バイナリ)ではなく一般化したもの？

CIA = ${Y_0, Y_1} ⊥ Z_i | X_i$

pira-nino.hatenablog.com

また、これらより強い条件としてStrongly Ignorable Treatment Assignment(強く無視できる割当条件) というものがあり、「どのx(Z)も観測値」+ 「どのiもどのdになる可能性もある」という条件が加わる模様。このあたりは各著者がどういう意図で書き分けてるのかは謎。おおもとのRubinはStrongly Ignorable Treatment Assignmentを置いて効果検証をしているっぽい。

qiita.com

実践

オレンジジュースデータを用いて、値段が売上に与える効果を考える。まずは単回帰をおこなう。

basefit <- lm(log(sales) ~ log(price), data=oj)
coef(basefit)

# (Intercept)  log(price) 
#   10.423422   -1.601307

このときの-1.601307は、値段に対して価格の割当が一定でない、つまり共変量Xが含まれていない。仮に共変量がオレンジジュースのブランドのみだとすると、以下では共変量X=ブランドを加えることで値段の売上に対する割当がランダムとなる。ブランドをA,B,Cとして式で書くと、

$｛Sales_i(Price)｝⊥ Price_i | Brand_i$

$p(Price_i | Sales_i(Price), Brand_i) = p(Price_i | Brand_i)$

となる。つまり、「ブランドさえ同じなら、例えば値段がいくらであろうと、それぞれでの売上の確率分布は同じ」だし、「ブランドさえ同じなら、売上がいくらであろうと値段の確率分布は同じ」ということ。逆に、「ブランド毎で、同じ値段のときの売上の確率分布が変わる」ともいえる。

そのため、ブランドを式に入れてControllすることで、「ブランドによって売上の確率分布が変わる」のを起きなくさせて、値段を通したブランドによる売上への効果を除去した、値段と売上のストレートな関係を観ることができる。

brandfit <- lm(log(sales) ~ brand + log(price), data=oj)
coef(brandfit)

#     (Intercept) brandminute.maid   brandtropicana       log(price) 
#      10.8288216        0.8701747        1.5299428       -3.1386914

You can better contorol for confounder by modeling the treatment assignment process

LTE

CIAを明示的に回帰式に組み込むと以下のLinear Treatment Effect(LTE)で表すことができる。

$y = dγ + x'β + ε, ε | d, x = 0,$ ①
$d = x'τ + ν, ν | x = 0$ ②

dは処置、xはyに対するdの全共変量を指す。x'は yに対するdの共変量の一部(不完全な共変量)を指している。ここで $ε | d, x = 0, ν | x = 0$ がCIAを指している。

つまり、2段目(②)の式 $d = x'τ + ν, ν | x = 0$ によって、x`に足りない全共変量xのdへの影響がνとして表される。
そして、そのx', dを使った1段目(②)の式 $y = dγ + x'β + ε, ε | d, x = 0$ では不完全な共変量x'を使っているが、dはνを用いることで擬似的に全共変量xを含めた値となっていることからγは全共変量xでContorollしたyへの効果を示すことになる。

要するに、全共変量xを求めるのは難しいので、不完全な共変量x'だけで上手くdのyへの効果を抽出するのがLTEとなっている。

ただし、これは共変量xの次元が少ない(次元<<n)の場合を前提としているため高次元xのときは成り立たない。そのため、できる限り重要な共変量をxとして入れて低次元化する必要があるがどれを選んだらいいか試すのは非常に時間がかかるのでlassoとCVを用いることで解決させる。

LTE LassoRegression

LTEをLasso経由でうまく求める。

d = E(d|x) = x`τ としてML(Lasso)で学習をおこない、実際のデータにfittingさせて $\hat{d_i} = x'_i \hat{τ}$ を求める。
$E(y|x,d) = \hat{d}δ + dγ + x'β$ をLassoで求める( $\hat{d}$ はモデルに必要なので除外さえないように罰則項を入れない)

このとき2.で推定された $\hat{τ}$ が効果(treatment effect)になる。

より詳細に書くと、1.で求めた $\hat{d_i} = x'_i \hat{τ}$ と真のdの差分は $d - \hat{d} = \hat{v}$ となる。この $\hat{v}$ はLTEの式②よりx = x'であれば $\hat{v} = 0$ となることから、 $\hat{v}$ は共変量の一部であるx'で表せきれていない共変量xのdへの効果を捉えていると考えられる。
そのため、LTEの式① $E(y | d) = dγ + x'β$ に式② $d = x'τ + ν$ を代入すると、 $(x'τ +ν)γ + x'β$ となり、変形すると $νγ + x'(γτ + β)$ となる。 $(γτ + β) = β'$ とすると $νγ + xβ'$ になり、これとLTEの式①を見比べると、「効果をみたい変数d」の係数と「 $\hat{v}$ = 共変量の一部であるx'で表せきれていない共変量xのdへの効果」の係数はどちらもγとなる。

2.では「 $\hat{d}, x'$ でControllした $d$ の係数γ」は、「 $\hat{v}$ の効果」となる。「 $\hat{v}$ の効果」は「 dの全共変量でコントロールしたときの効果」となることから、2.式より、dの全共変量でコントロールしたときの効果を係数γから求めることができる。

実践

以下の殺人率データを用いて、ケータイ電話の所持率の効果をみる。

github.com

sna <- factor(s, levels=c(NA,levels(s)), exclude=NULL)
x = sparse.model.matrix( ~ t + phone*sna + .^2, data=controls)[,-1] # 共変量x'

# 1.部分
treat <- cv.gamlr(x ,d, lmr=1e-3) # E[d|x] = x'τを作成
dhat <- drop( predict(treat, x, select="min") ) # treatを使ってdhatを求める

# 2.部分
causal <- cv.gamlr(cBind(d,dhat,x),y,free=2,lmr=1e-3)  # dhatのみpenaltyを与えない
coef(causal, select="min")["d",]  # 共変量で統制したdのyへの効果
# => 0.

Sample-Splitting and Orthogonal Machine Learning

LTE LassoではLassoを用いて効果を求めている。Lassoでは通常の方法では信頼区間を求めることができない。しかし、3章で紹介したような方法(※当ブログでは書いてない)でLassoに対してモデル選択後の係数信頼区間をパラメトリックブートストラップとサブサンプリングで求める方法を紹介している。

LTE Lassoでもその方法は可能だが、LTEモデルを用いたLassoに関してより簡易な方法として、データを2分割し、

片方のデータでLassoでモデル選択をおこなう
そのモデルでOLSをおこなう

ことで、OLSが噛むためOLS部分で信頼区間を求めることができる。また、これをCross Validationでおこなうことで汎化的な値を求める。

具体的には、

1. データをK分割する
1-1. k群のデータで、LassoなどのMLを用いて E[d | x]とE[y | x] *2のモデル選択をおこなう
1-2. k群以外のデータで1-1.のモデルを用いたd, y の予実差 $\tilde{d}, \tilde{y}$ を求める
1-3. 上記をk+1, k+2,...Kでもおこなう
1. $\tilde{y} | \tilde{d} = α + \tilde{d} γ$ を求め、 $\hat{γ}$ の信頼区間を求める

LTE LassoRegressionの以下

① d = E[d|x] = x`τ としてML(Lasso)で学習をおこない、実際のデータにfittingさせて $\hat{d_i} = x'_i \hat{τ}$ を求める。
② $E[y|x,d$ = \hat{d}δ + dγ + x'β] をLassoで求める( $\hat{d}$ はモデルに必要なので除外さえないように罰則項を入れない)

と対応させて考えると、1-1. 1-2. は①に当てはまる。そのため、予実差 $\tilde{d}$ は、x(x')だけでは表せきれていない共変量の効果 $\hat{v}$ となる。
また、E[y | x] から作成した $\tilde{y}$ はx(x')だけでは表せきれていないy部分となる。

そのため、 $\tilde{y}$ と $\tilde{d}$ の関係式の $\tilde{d}$ の係数は2. $\tilde{y} | \tilde{d} = α + \tilde{d} γ$ のように効果γとなる。
言い換えると、「yの共変量xで表せきれていない部分」と「dの共変量xで表せきれていない部分」の関係性となるのでそれは効果と等しくなる。

そしてこれはOLSなのでγの信頼区間も合わせて求めることができる、

実践

上記アルゴリズムは以下から読み込める

github.com

source("examples//orthoML.R")
dreg <- function(x,d){ 
    gamlr(x, d, standardize=FALSE, lmr=1e-5) }

yreg <- function(x,d){ 
    gamlr(x, d, family="binomial", standardize=FALSE, lmr=1e-5) }

resids <- orthoPLTE( x=x[,-who], d=x[,who], y=y, 
                dreg=dreg, yreg=yreg, nfold=5)

# fold: 1  2  3  4  5  
# gamma (se) = 0.238736 (0.0212952)


0.238736 + c(-2,2)*0.0212952 # 信頼区間2se
# 0.1961456 0.2813264

2*pnorm(-0.238736/0.0212952) # p値
# => 3.609601e-29

Heterogeneous Treatment Effect

先程まで観ていた効果は、処置有無毎でのAverage Treatment Effect(ATE)となっている。このATEは、例えば男女混合での効果だが男性への効果/女性への効果、のようにあるグループに分けた効果をHeterogeneous Treatment Effect(HTE)と呼ぶ。

ビジネスにおいて、効果が高いHeterogeneousがどれかを把握することは重要なためHTEは非常に役に立つ。

これは、

$E(y_i | x_i, d_i) = α + x'_i β + d_i γ_0 + (d_i × x_i)' γ$

となる。このx'は共変量+ Heterogeneousとなる。 $γ_0$ はbaseとなる効果で、交互作用 $(d_i × x_i)'$ によって、共変量, Heterogeneousそれぞれの効果がベクトルγで現れる。そのため、Heterogeneousの効果は $γ_0 + γ_{Heterogeneous}$ で現れる。

なお、式からわかるようにHeterogeneous毎にサブサンプリングをおこなって通常のLTEをおこなうことでもHeterogeneousの効果を表すことはできる。

また、これより例えばA,Bテストの場合

$HTE_i = E(y|i | x_i, d_B) - E(y|i | x_i, d_A) = (γ_0 + x'_i γ) × (d_B - d_A)$

となる。

実践

データは以下

github.com

まずは、Heterogeneousとしてみたい項目XのNAを無くす。

### linear HTE estimation (Controls chapter)
library(gamlr)
source("examples/naref.R")

## dealing with missingness
# make NA the reference level for all categories
X <- naref(X)

# pull out the numeric variables 
xnum <- X[,sapply(X,class)%in%c("numeric","integer")]
xnum[66:70,]
colSums(is.na(xnum))
# flag missing
xnumna <- apply(is.na(xnum), 2, as.numeric)
xnumna[66:70,]
# impute the missing values
mzimpute <- function(v){ 
    if(mean(v==0,na.rm=TRUE) > 0.5) impt <- 0
    else impt <- mean(v, na.rm=TRUE)
    v[is.na(v)] <- impt
    return(v) }
xnum <- apply(xnum, 2,  mzimpute)
xnum[66:70,]

# replace/add the variables in new data frame 
for(v in colnames(xnum)){
    X[,v] <- xnum[,v]
    X[,paste(v,"NA", sep=".")] <- xnumna[,v] }
X[144:147,]

次に先程のみたい項目X(Heterogeneous)と共変量 P$numhhを結合してx'を作成する。

xhte <- sparse.model.matrix(~., data=cbind(numhh=P$numhh, X))[,-1]
xhte[1:2,1:4]

dim(xhte)
# 6855   64

P %>% head()
# person_id household_id doc_any_12m selected medicaid numhh
# 1         1       100001           0        1        0     1
# 2         2       100002           0        1        1     1
# 5         5       100005           0        1        0     1
# 6         6       100006           1        1        0     1
# 8         8       102094           0        0        0     2
# 9         9       100009           1        0        0     1

次に、x'とdの交互作用項を作成。

dxhte <- P$selected*xhte
colnames(dxhte) <- paste("d",colnames(xhte), sep=".") # 列名をd.x名とする
htedesign <- cBind(xhte,d=P$selected,dxhte)

そして、x'と交互作用項をあわせたhtedesignを用いてLassoを行う。

# include the numhh controls and baseline treatment without penalty 
htefit <- gamlr(x=htedesign, y=P$doc_any_12m, free=c("numhh2","numhh3+","d")) # numhh2, numhh3, d列は罰則項なし
gam <- coef(htefit)[-(1:(ncol(xhte)+1)), ]

このときgamにd.Xの交互作用項の係数をgamに格納している。これをいくつか見てみる。

round(sort(gam)[1:6],4)
round(sort(gam, decreasing=TRUE)[1:6],4)

f:id:chito_ng:20200704190025p:plain

dとなっているものが今回効果をみたいもの(selected)のbase効果となる。ここから、例えばrace_pacific_12mYesであることでd+race_pacific_12mYes = 0.0927+0.0404 = 0.133つまり、13%増加する効果がある。

なお、ATEは [\hat{γ_0} + x'平均 \hat{γ_0} ]で計算できる。

gam["d"] + colMeans(xhte)%*%gam[grep("d.", names(gam))]
# => 0.05616546

その他

今回でBUSINESS DATA SCIENCEの勉強は終わり。他にも自然言語処理とかの章もあるが使わないのでいったん保留。
全体的にビジネスという観点で必要なものをどうやって出すか、ということで実用的だった。

2点だめな点をあげると、Rの書き方が今風(tidyverse)ではないこと、sampledataが散在していること。
何のデータを使うか単位でRファイルがあり、そのファイルは色々な章で登場する。そのため、例えば1,5章でデータXを使う場合1,5章のコードが全てX.Rに記載されている。そして、5章は1章部分のコードの続きという立ち位置なので1章部分を読み込まないと5章部分だけでは動かない。
また、書籍内コードにread_csv部分の記載やどのRファイルを使うか書いてないにも関わらず、ファイル名がわかりづらい。つまり、毎回どのファイルか調べないといけないので糞めんどくさい。

*1:当ブログではまとめてない

*2:E[y | x, d]ではないので注意

まずは蝋の翼から。

学んだことを書きながら確認・整理するためのメモブログ。こういうことなのかな？といったことをふわっと書いたりしていますが、理解が浅いゆえに的はずれなことも多々あると思うのでツッコミ歓迎

BUSINESS DATA SCIENCE 6章 Controls

この章はなにか

Conditional Ignorability

実践

You can better contorol for confounder by modeling the treatment assignment process

LTE

LTE LassoRegression

実践

Sample-Splitting and Orthogonal Machine Learning

実践

Heterogeneous Treatment Effect

実践

その他