Rで始める傾向スコア分析による効果検証

最近因果推論の勉強を始めた管理人です. 因果効果の推定の基礎となるポテンシャルアウトカムフレームワーク, 選択バイアス, 傾向スコア(特に逆確率重み付き推定)についてRでの実装つきでまとめたので記事にしました. 大まかな内容については、以下の本を参考にしています.

参考 - 効果検証入門(安井)

1. 因果効果

対象への働きかけによって, 結果にどの程度影響が及ぶかについて推定したい状況を考えます.
例えば, ある新薬を患者に投与することで症状の改善が見られるかを測りたい状況をイメージするとよいでしょう.
因果推論の分野において, 対象への働きかけを介入や処置と呼び, 介入によって生じた効果を因果効果と呼びます.

1.1 ポテンシャルアウトカムフレームワーク

ポテンシャルアウトカムフレームワーク(potential outcome framework)は, 因果効果を統計的に取り扱う考え方です. 潜在的な観測量(potential outcome)とは, 介入が行われていない場合に観測される結果を指します.

潜在的な観測量を表す確率変数を $Y^{(0)}$ , 介入時の観測量を表す確率変数を $Y^{(1)}$ とします.
この時, 介入の効果は次の量で表されることになります(下図参照).

$\tau = Y^{(1)} - Y^{(0)}$

f:id:kimigayoseishou:20200715225206p:plain — 因果効果

この時, 平均的な因果効果は次のように表されます,

$\hat{\tau} = E[Y^{(1)} - Y^{(0)}$ ]

これは, 平均処置効果(ATE/Average Treatment Effect)と呼ばれる量です.

平均処置効果を推定するには, 十分な数のサンプルを用意し各サンプルに対して $Y^{(1)} - Y^{(0)}$ を観測しそれらの平均値を取るという方法が考えられるでしょう.

しかし, 実際のケースとして各サンプルに対して介入時の結果と介入なしの場合のどちらかしか観測できないことが多く, 上記のような平均処置効果の推定は現実的ではありません(因果推論の根本問題).

そこで, 同一の対象に対して介入時の結果と介入なしの場合の測定をするのではなく,
介入が行われたグループと介入が行われなかったグループに分た各グループ間の比較による因果効果の推定(対照実験)が応用上重要となります.

1.2 選択バイアスと無作為化実験

対照実験によって得られたデータを用いて因果効果を分析する場合, 意図せず介入の仕方が偏ってしまう場合があります.
例えば, 血圧を下げる効果が期待される新薬の治験において, 新薬を投与した患者のグループの投薬前の血圧が他方のグループより高い数値であった場合などがこれにあたります.

このような状況において, 結果へ影響を与えるような性質の偏りが比較したいグループ間に存在することによって因果効果が正しく推定できないことを選択バイアスといいます.
以下, どのようにして選択バイアスが発生するのか説明します.

記号

$N$ : サンプルサイズ
$Z_i (i = 1, \ldots, N)$ : 介入の有無を表す変数(介入ありの場合1, そうでない場合0).
$Y_i (i = 1, \ldots, N)$ : 観測結果
$Y^{(1)}_i (i = 1, \ldots, N)$ : 介入があった場合の観測結果
$Y^{(0)}_i (i = 1, \ldots, N)$ : 介入がない場合の観測結果

介入の効果を推定するのに, 介入を行ったグループの平均値とそうでないグループの平均値の差もしばしば用いられます,

$\hat{\tau}_{naive} = \frac{1}{\sum_{i=1}^N Z_i}\sum_{i = 1}^N Z_iY_i - \frac{1}{\sum_{i=1}^N (1- Z_i)}\sum_{i = 1}^N (1-Z_i)Y_i$

サンプルサイズ $N$ を増やす時, グループ間の平均値の差は以下の量に等しくなります.

$\tau_{naive} = E[Y^{(1)} | Z = 1 ] - E[Y^{(0)} | Z = 0]$

$\tau_{naive}$ と $\tau$ を比較するために式変形を行います,

$\tau_{naive} \\ = E[Y^{(1)} | Z = 1] - E[Y^{(0)} | Z = 0] \\ = E[Y^{(1)} | Z = 1] - E[Y^{(0)} | Z = 1] + E[Y^{(0)} | Z = 1] - E[Y^{(0)} | Z = 0] \\ = E[Y^{(1)} - Y^{(0)} | Z = 1] + E[Y^{(0)} | Z = 1] - E[Y^{(0)} | Z = 0]$

ここで, 因果効果 $Y^{(1)} - Y^{(0)}$ と介入の有無が独立だと仮定することにより $E[Y^{(1)} - Y^{(0)} | Z = 1] = E[Y^{(1)} - Y^{(0)}] = \hat{\tau}$ となるため, 平均処置効果とグループ間の平均差は $E[Y^{(0)} | Z = 1] - E[Y^{(0)} | Z = 0]$ だけずれることがわかります. この差は偏った介入の仕方に起因するもので選択バイアス(selection bias)と呼ばれています.

例えば, $Z_i$ と $Y^{(0)}_i$ に正の相関があるような場合, $Y^{(0)}_i$ の大きい部分では $Z=1$ の傾向が強くなりますが, そのようなサンプルは観測されません. その結果, 介入なしのグループのサンプルの平均値は真の平均値 $E[Y^{(0)}]$ よりも小さく推定されてしまいます(下図).

f:id:kimigayoseishou:20200715225333p:plain — 選択バイアス

選択バイアスによる影響をなるべく減らすために対象への介入を無作為に割り当てて行う対照実験は 無作為化比較実験(RCT : Randomized Controlled Trial)と呼ばれています.

1.3 Rによる選択バイアスの数値実験

ここでは介入に偏りのあるデータを作成し, 選択バイアスを擬似的に発生させる数値実験をRを用いて行ってみます.

(a)介入が無作為である場合

まずは介入が無作為である場合です. 平均0, 標準偏差0.5の正規分布から200個のサンプルを生成し, これを介入が行われていない場合のデータとします(y0). また介入によって結果は1増えると仮定し, 介入が行われた場合のデータを作成します(y1). これらのデータy0, y1に対して, 無作為に介入を行い, 介入の有無によって2つのグループに分けてその平均値を比較します.

library(dplyr)
library(ggplot2)
set.seed(1)
# 介入なし
y0 <- rnorm(200, 0, 0.5)
# 介入が行われた場合, 観測結果は1増加
y1 <- y0 + 1
# 無作為な介入
z <- sample(c(0, 1), 200, TRUE)

# 介入は無作為になっている
data <- tidyr::tibble(Y = y1 * z + y0 * (1-z), Z = z)

data %>%
    group_by(Z) %>%
    summarise(mean_Y = mean(Y))

# 赤点は各グループの平均, 赤縦線は平均 ± 標準偏差の範囲を示す
ggplot(data, aes(x = as.factor(Z), y = Y)) + 
    geom_point(alpha=.4) +
    stat_summary(
        fun.y = mean,
        fun.ymin = function(x) mean(x) - sd(x), 
        fun.ymax = function(x) mean(x) + sd(x),
        color = "red"
        ) +
    xlab("介入の有無(介入あり Z = 1, 介入なし Z = 0)") + 
    ggtitle('介入ありグループと介入なしグループの平均値の比較(RCT)')

f:id:kimigayoseishou:20200715225417p:plain — RCT　グループごとの平均値

出力結果を確認すると, 介入なしグループの平均値は約-0.02, 介入ありグループの平均値は約1.06となっており, その平均値の差=1.06 - (-0.02) =1.08は真の因果効果=1と近い値になります.

f:id:kimigayoseishou:20200715225446p:plain — RCTプロット

(b)介入に偏りがある場合

次に介入に偏りがある場合を考えます. 以下のコードでは, 「介入なしの場合の値が0より大きい場合は必ず介入」というルールを追加して意図的に介入が偏ったデータを作成します.

# 意図的に選択バイアスを発生させる
data_biased <- tidyr::tibble(Z = ifelse(y0 > 0, 1, z)) %>%
    mutate(Y = y1 * Z + y0 * (1-Z))

data_biased %>%
    group_by(Z) %>%
    summarise(mean_Y = mean(Y))

# 赤点は各グループの平均, 赤縦線は平均 ± 標準偏差の範囲を示す
ggplot(data_biased, aes(x = as.factor(Z), y = Y)) + 
    geom_point(alpha=.4) +
    stat_summary(
        fun.y = mean,
        fun.ymin = function(x) mean(x) - sd(x), 
        fun.ymax = function(x) mean(x) + sd(x),
        color = "red"
        ) +
    xlab("介入の有無(介入あり Z = 1, 介入なし Z = 0)") + 
    ggtitle('介入ありグループと介入なしグループの平均値の比較(非RCT)')

f:id:kimigayoseishou:20200715225522p:plain — 非RCT 平均値

出力結果を確認すると, 介入なしグループの平均値は約-0.34, 介入ありグループの平均値は約1.17となっており, その平均値の差=1.17 - (-0.34) = 1.51は真の因果効果=1と比較して大きく推定されてしまいます. このデータではy0 > 0のデータについて必ず介入を行ったことでy0 > 0のサンプルが欠測となり, その結果介入なしのグループ平均値は真の平均値 $E[Y^{(0)}$ =0]よりも小さく推定されていることがわかります.

f:id:kimigayoseishou:20200715225608p:plain — 非RCT プロット

上記の例をまとめると以下のようになります.

RCTの場合は, グループ毎の平均値の差によって平均処置効果を推定できる.
非RCTの場合は, 選択バイアスによりグループ毎の平均値の差は平均処置効果と解離することがある.

特にグループ毎の平均値を用いて因果効果の推定を行う場合は, 介入の仕方に偏りがないか注意を払う必要があります.

2 傾向スコア

傾向スコア(propensity score)は, 無作為でない介入が行われたサンプルからでも因果効果を推定できる手法として広く用いられています. 以下, 傾向スコアについて最小限の事項をまとめています.

詳細については,

などなどを参照されたし....

2.1 バランシングスコアと傾向スコア

以下傾向スコアを導入するにあたり, いくつか用語の準備をします.

記号

$X$ : 説明変数
$Y$ : 応答変数
$Y^{(1)}$ : 介入が行われた場合の結果
$Y^{(0)}$ : 介入が行われな買った場合の結果
$Z$ : 介入の有無を表す変数(介入ありの場合1, そうでない場合0)

定義-バランシングスコア(balancing score)

ベクトル値関数 $b(X)$ がバランシングスコア(balancing score)であるとは, 次の条件が成り立つことをいう

$X \perp Z | b(X).$

定義-傾向スコア(propensity score)

$x$ の関数 $P:x \mapsto p(Z = 1|X =x)$ を傾向スコア(propensity score)と定義する.

定義-強く無視できる割り当て(strongly ignorable treatment assignment)

介入の割り当てが確率変数 $W$ のもとで強く無視できる割り当て(strongly ignorable treatment assignment)であるとは, 次の条件が満たされることをいう

$(Y^{(1)}, Y^{(0)}) \perp Z | W.$

さらに $X$ に関して強く無視できる割り当てである場合, 単に強く無視できる割り当てと呼ぶ.

以上の定義のもとで次が成り立ちます.

定理
1. 傾向スコアはバランシングスコアである.
2. 介入の割り当てが強く無視できる割り当て => 任意のバランシングスコアのもとで強く無視できる割り当てになっている.

証明

こちらのスライドを一部参考にさせていただきました.

また, $p(Z=1 | P=e)$ を $X$ に関して周辺化し計算を行う. $p(Z = 1|P = e) \\ = \int p(Z = 1, X=x|P = e) dx \\ = \int p(Z = 1| X=x, P = e) p(X = x|P=e)dx \\ = \int_{P(x) = e} p(Z = 1, | X=x, P = e) p(X = x|P=e)dx \\ = \int_{P(x) = e} p(Z = 1, | X=x) p(X = x|P=e)dx \\ = \int_{P(x) = e} P(x) p(X = x|P=P(x_2))dx \\ = e \int_{P(x) = e} p(X = x|P=e)dx \\ = e$

よって, 示すべき等式の左辺 $=ep(X=x'|P = e) = p(Z = 1|P = e)p(X=x'|P = e)=$ 示すべき等式の右辺となり結果が従う.

2 $p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0)|b = b') = p(Z=z|b = b')p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|b = b')$ を示せば良い. 以下のように $X$ について周辺化を行うことで, 示すことができる. $p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0)|b = b') \\ = \int p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0), X = x|b = b') dx \\ = \int p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x, b = b') p(X = x|b=b')dx \\ = \int_{b(x) = b'} p(Z=z, (Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(X = x|b=b')dx \\ = \int_{b(x) = b'} p(Z=z|X=x)p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(X = x|b=b')dx \quad (\because (Y^{(1)}, Y^{(0)}) \perp Z | X)\\ = \int_{b(x) = b'} p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(Z = z, X = x|b=b')dx \\ = \int_{b(x) = b'} p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(Z = z|b=b')p(X=x|b=b')dx \quad (\because X \perp Z |b)\\ = p(Z = z|b=b') \int_{b(x) = b'} p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|X = x) p(X=x|b=b')dx \\ = p(Z = z|b=b') \int p((Y^{(1)}, Y^{(0)}) = (y_1, y_0), X = x|b = b')dx \\ = p(Z = z|b=b')p((Y^{(1)}, Y^{(0)}) = (y_1, y_0)|b = b').$

以上の結果より, 介入の割り当てが強く無視できる割り当てであるとき, $E[Y^{(i)}|Z=i, P]\ (i=0, 1)$ は真の平均値 $E[Y^{(i)}]$ の不偏推定量になっていることが次のようにして導かれます,

$E[Y^{(i)}] = E[E[Y^{(i)}|P]] = E[E[Y^{(i)}|Z=i, P]] \quad (\because Z \perp Y^{(i)} | P).$

注目すべきは $Y^{(i)}|Z=i, P$ は欠測とはならないため, (傾向スコアの推定を行った上で)サンプルから推定が行える点です.

以上をまとめると, 傾向スコアを用いて平均処置効果を推定する大まかな流れは以下のようになります.

介入の割り当てが強く無視できる割り当てとみなせるような説明変数 $X$ の選定を行う.
(GLM等のモデルを用いた)傾向スコア $P(x) = p(Z = 1 | X =x)$ の推定.
各傾向スコア $P$ に対して, $E[Y^{(i)}|Z=i, P]$ をサンプルから推定し, 真の平均値 $E[Y^{(i)}]$ を推定する.
3で推定した値の差をとった量を平均処置効果の推定量とする.

2.2 逆確率重み付き推定

逆確率重み付き推定(IPW : Inverse Probability Weighting)は傾向スコアを利用した真の平均値 $E[Y^{(i)}]$ の推定法の一つです.

$N$ をサンプルサイズ, $X_i$ をサンプル $i$ における説明変数, $Y_i(i = 1, \ldots, N)$ をサンプル $i$ における観測値, $Z_i(i = 1. \ldots, N)$ をサンプル $i$ における介入の有無とする. サンプル $i$ の傾向スコアを $P_i = p(Z = 1|X = X_i)$ とするとき, $\hat{Y}^{(i)}_N(i=0, 1)$ を

$\hat{Y}^{(1)}_N = \sum_{j=1} ^N \frac{Y_iZ_i}{P_i} \bigr/ \sum_{j=1}^N \frac{Z_i}{P_i} \\ \hat{Y}^{(0)}_N = \sum_{j=1} ^N \frac{Y_i(1-Z_i)}{1-P_i} \bigr/ \sum_{j=1}^N \frac{1-Z_i}{1-P_i}$

と定義します.

命題-(逆確率重み付き推定)
$\hat{Y}^{(i)}_N(i=0, 1)$ は $E[Y^{(i)}](i=0, 1)$ の一致推定量である.

証明

$\hat{Y}^{(i)}_N \rightarrow E[Y^{(i)}](N \rightarrow \infty)$ を示せば良い. $i=1$ の場合のみ示す.

$E[YZ|P] \\ = E[(ZY^{(1)} + (1-Z)Y^{(0)})Z | P] \\ = E[ZY^{(1)}| P] \\ = E[Z| P]E[Y^{(1)}|P] \quad ( \because Y^{(1)} \perp Z | P) \\ = P E[Y^{(1)}|P]$

であるので, 大数の法則より

$\frac{1}{N}\sum_{j=1}^N\frac{Y_iZ_i}{P_i} \\ \rightarrow E \biggl[ \frac{YZ}{P} \biggr] \quad (N \rightarrow \infty)\\ = E \biggl[ \frac{E[YZ | P]}{P} \biggr] \\ = E \biggl[ \frac{PE[Y^{(1)}|P] }{P} \biggr] \\ = E [E[Y^{(1)}|P] ] \\ = E[Y^{(1)}] \\$

となる. 再び大数の法則より,

$\frac{1}{N}\sum_{j=1}^N \frac{Z_i}{P_i} \\ \rightarrow E \biggl[ \frac{Z}{P} \biggr] \quad (N \rightarrow \infty)\\ = E \biggl[ \frac{E[Z|P]}{P} \biggr] \\ = E \biggl[ \frac{P}{P} \biggr] = 1$

以上より, $\hat{Y}^{(i)}_N = \frac{ \frac{1}{N}\sum_{j=1}^N\frac{Y_iZ_i}{P_i}}{\frac{1}{N}\sum_{j=1}^N \frac{Z_i}{P_i}} \rightarrow \frac{E[Y^{(1)}]}{1} = E[Y^{(1)}] \quad (N \rightarrow \infty).$ となるから, 確かに $\hat{Y}^{(1)}$ は $E[Y^{(1)}]$ の一致推定量である.

逆確率重み付き推定の数値実験

ここでは, 擬似データを用いて逆確率重み付き推定による因果効果の推定の数値実験を行います.

あるアパレル小売店Aでは, 販売促進のためにクーポン券の配布を行いました. そこでクーポンの利用によってクーポン利用可能期間中のユーザーの購買額が増加したかを分析したいとします.

$X_i$ : ユーザー $i$ の衣料品に対する平均支出額
$Y^{(0)}_i$ : ユーザー $i$ のクーポンがない場合の購入金額
$Y^{(1)}_i$ : ユーザー $i$ のクーポンを利用した場合の購入金額
$Y_i$ : ユーザー $i$ の実際の購入金額
$Z_i$ : クーポン利用の有無(クーポンを利用した場合1, そうでない場合0)

クーポンがない場合の購入金額は衣料品に対する平均支出額を説明変数とした線形モデルで表せるとし、さらにクーポン利用時は $\lambda$ だけ購入金額が増加すると仮定します.

$Y^{(0)}_i \sim {\sf Normal}(\alpha_0 + \alpha_1 X_i, \tau^{2}) \quad (\tau {\text ＞}0) \\ Y^{(1)}_i = Y^{(0)}_i + \lambda.$

次にクーポン利用についてです. ここでは, クーポンの利用者は非利用者よりも衣料品支出額が多いという仮定をおきます.

クーポン非利用者の衣料品に対する平均支出額は平均 $\mu_0$ , 分散 $\sigma^{2}$ の正規分布に, クーポン利用者の衣料品に対する平均支出額は平均 $\mu_1$ , 分散 $\sigma^{2}$ の正規分布に, それぞれ従うとします

$X|Z=j \sim {\sf Normal}(\mu_j, \sigma^{2})\quad (\mu_j {\text ＞} 0, \sigma^{2} {\text ＞} 0, j= 0, 1),$

以上のモデルにしたがって200人分(うちクーポン利用者100人)の擬似データを作成します. 各パラメータは次のようにしています,

$\alpha_0=0.1, \alpha_1 = 0.05, \tau = 0.05, \lambda = 0.2, \mu_1 = 10, \mu_0 = 6, \sigma = 2$

# 擬似データの作成
library(dplyr)
set.seed(1)
alpha0 <- .1
alpha1 <- .05
tau <- .1
lambda <- .2
N <- 200L
mu1 <- 10
mu0 <- 6
sigma <- 2
# クーポン利用者の衣料品に対する平均支出額
x1 <- rnorm(N/2, mu1, sigma)
# クーポン非利用者の衣料品に対する平均支出額
x0 <- rnorm(N/2, mu0, sigma)
# 全ユーザーの衣料品に対する平均支出額
x <- c(x0, x1)

# クーポン非利用時の衣料品に対する平均支出額
y0 <- alpha0 + alpha1 * x + rnorm(N, 0, tau)
y1 <- y0 + lambda
z <- c(rep(0L, N/2), rep(1L, N/2))

# 200件の購買データ
purchase <- tidyr::tibble(Z = z, Y = z * y1 + (1-z) * y0, X = x)

purchase %>%
    ggplot(aes(x = Y, fill = Z)) +
    geom_histogram(alpha = .6, position = "identity") +
    ggtitle('実際の購入金額の分布')

f:id:kimigayoseishou:20200715225733p:plain — 購買金額のヒストグラム

まず, 介入ありグループと介入なしグループで平均値の差を確認します.

lm(data = purchase, formula = Y ~ Z) %>%
    broom::tidy()

f:id:kimigayoseishou:20200715225656p:plain — グループ毎平均値

出力結果のterm=Zの行のestimateの値(=0.42)は, グループ毎の平均値による平均処置効果の推定値になります. この値は真の平均処置効果 $\lambda = 0.2$ の約2.1倍の数値で, 選択バイアスを含んだ値となっています.

そこで選択バイアスの影響を除去するため, 逆確率重み付き推定による平均処置効果の推定値を行います.

まずクーポンの利用率 $p(Z = 1) = 0.5$ として, 傾向スコアを計算すると,

$P(x) \\ = p(Z = 1|X=x) \\ = \frac{p(X = x|Z =1)p(Z=1)}{p(X = x|Z =0)p(Z=0) + p(X = x|Z =1)p(Z=1)} \\ = \frac{p(X = x|Z =1)/p(X = x|Z =0)}{1 + p(X = x|Z =1)/p(X = x|Z =0)} \\ = \frac{e^{\beta_0 + \beta_1x}}{1 + e^{\beta_0 + \beta_1x}}$

となり傾向スコアはロジスティック関数で表されます. ただし, ここで

$\beta_0 = - \frac{\mu_1 ^{2}-\mu_0 ^{2}}{2\sigma ^{2}}, \beta_1 = \frac{\mu_1-\mu_0}{\sigma^{2} }$

としています.

注意今の状況下では, $\beta_0, \beta_1$ は $\mu_0, \mu_1, \sigma$ からカンニングできてしまいますが, 実際のデータ分析時にはデータから推定する必要があります.

この式から傾向スコアの算出を行います.

# 傾向スコアの計算
beta0 <- -.5 * (mu1^2-mu0^2)/sigma^2
beta1 <- (mu1 - mu0) / sigma^2
prop_score <- 1 / (1 + exp(-(beta0 + beta1 * purchase$X)))

# 逆確率重み付き推定
# クーポン利用時の購買金額推定値
y1_estim  <- sum(purchase$Y * purchase$Z / prop_score) / sum(purchase$Z / prop_score)

# クーポン非利用時の購買金額推定値
y0_estim <- sum(purchase$Y * (1 - purchase$Z) / (1 - prop_score)) / sum((1 - purchase$Z) / (1 - prop_score))

# 逆確率重み付き推定によるクーポン利用時/非利用時の購買金額推定値の差が, 
# 平均的な因果効果の推定値となる
y1_estim - y0_estim