【R】傾向スコアマッチングによる因果推論の基礎

傾向スコアマッチングによる因果推論について調べたので備忘録を残しておきます。

研究デザイン

参考書籍の『調査観察データの統計科学―因果推論・選択バイアス・データ融合』では研究デザインを大きく2つに分類・定義している。

  • 実験研究 (experimental study) : 原因となる変数を研究者が操作し, 結果となる変数がどのように変化するかを調査
  • 観察研究 (observational study) : 調査観察研究/相関研究。無作為割り当てを伴わない研究

相関分析, 連関分析が観察研究であるのに対して, 因果分析では実験研究であることが望ましい。[1]

原因となる変数を直接操作できる工学分野と異なり, 社会学や疫学では倫理的な理由などにより直接操作できなかったり無作為割り当てが行えない場合がある。無作為割り当てが行えないということは, 興味がある変数について2群に分けた時に, 様々な要因がコントロールされていないことになる。

  • 実験群 (experimental group） : 無作為割り当てが行なわれている実験研究において特別な条件を与えた (介入された) 群
  • 処置群 (treatment group) : 自然実験, 観察研究において何らかの介入が行われた群
  • 対照群 (control group） : 特別な条件を与えない (介入されなかった) 群

しかし, 無作為割り当てが行えないような場合でも因果効果を推論する手法がある。

ルービンの因果モデル

潜在的結果変数 (potential outcomes) とは, 独立変数が取り得る値の数と同じ数だけ存在する 「仮想的な従属変数」 のことである。
介入有無のインディケータ変数 Ｚ に依らず, 観測できるのが Y[0] と Y[1] のどちらかであっても共に潜在的結果変数と言う。 つまり, 観測できる潜在的結果変数, 観測できない潜在的結果変数という言い方をする。

例えば, ある教育を受けた/受けないの2つの状態において, 本来は両方とも値が存在するが現実では一方しか観測されない。起きた世界と起きなかった世界を同時に観測することは不可能であり, これが因果推論の根本的問題となる。(反事実的依存性) [2]

しかし, ここで本来あり得た結果を考えるのがルービンの因果モデルの特徴である。あくまで実際には起きていない結果を仮想的に考えるので, 反実仮想モデル (counterfactual model) とも言われる。

無作為割り当て可能な場合の因果効果

無作為割り当てを行える場合は因果効果を処置群 Y[1] の期待値と対照群 Y[0] の期待値の差で表せる。[3] まず, これを確認してみる。

条件付き確率は

    \begin{eqnarray*} P(Y=y|X=x) = \frac{P(Y=y,X=x)}{P(X=x)} \end{eqnarray*}

条件付き期待値は条件付き確率を用いて

    \begin{eqnarray*} E(Y|X=x)\\ =y_1\cdot P(Y=y_1|X=x)+...+y_n\cdot P(Y=y_n|X=x)\\ =\sum_{i=1}^{n}y_i\frac{P(Y=y_i,X=x)}{P(X=x)} \end{eqnarray*}

また, X, Y が独立な場合の条件付き期待値は以下となる。

    \begin{eqnarray*} E[Y|X]=E[Y]\\ E[X|Y]=E[X] \end{eqnarray*}

介入の割り当てを確率変数 Z とすると潜在的結果変数は以下で表せる。 Z を与えると Y は観測できる潜在的結果変数 Y[1], または観測できない潜在的結果変数 Y[0] となる。

    \begin{eqnarray*} Y=\begin{cases} Y[1], Z=1 \\ Y[0], Z=0 \end{cases}\\ =ZY[1]+(1-Z)Y[0] \end{eqnarray*}

処置群の平均と対照群の平均の差は条件付き期待値の差である。ここで, 割り当て Z が無作為であると仮定すると Z と Y は独立となり, ルービンの因果効果の平均処置効果 (average treatment effect; ATE) [5] と等価となる。

    \begin{eqnarray*} E[Y[1]|Z=1]-E[Y[0]|Z=0]\\ =E[Y[1] - Y[0]]\\ =E[Y[1]] - E[Y[0]] \end{eqnarray*}

次に, 無作為割り当てが行えない場合に因果効果を推論する手法を紹介する。

傾向スコア

傾向スコア (propensity score) は共変量 [4] の値を x, 割り当て変数の値を z とするときに群1に割り当てられる確率 e である。(式は参考書籍を参照)
共変量を用いて因果効果以外の効果を除去することを共変量調整という。様々な共変量を傾向スコアという1つの変数に縮約することができ, これを基準として因果効果を推定する。

傾向スコアを用いた共変量調整は, 前提条件である 「強く無視できる割り当て」 条件が満たされていれば全ての共変量を用いて調整を行ったのと同じだけ偏りを減少させることができる。そして, 傾向スコアの値が等しい処置群と対照群のペアをマッチングすると, 差の平均はルービンの因果効果の不偏推定量となる。

rosenbaum によると, 傾向スコアを用いて因果効果を推定する方法は下記などがある。

  • マッチング
  • 層別解析
  • 共分散分析

マッチング・層別解析は似たものを比較する手法, SEM, 共分散分析は回帰モデルによって共変量の影響を調整する手法である。傾向スコア分析には下記の利点がある。

  • 処置群対象者に対応する対照群対象者が見つからない問題が解決される
  • 共変量と従属変数のモデル設定を行わなくてもよい
  • モデルの誤設定に強い(傾向スコア分析は回帰モデルよりロバストな結果)

一方で, マッチング・層別解析による傾向スコア分析には標準誤差や周辺期待値が求められないことや群毎のサンプルサイズに偏りがある場合に無駄なデータがでてしまうなどのいくつか問題点がある。詳しくは参考書籍を参照。

傾向スコアマッチング

傾向スコアマッチング (Propensity Score Matching) は, 処置群と対照群の2つの群で傾向スコアが等しい対象者をペアにして, その期待値の差をもって因果効果の推定値とする。

R の {Matching} はマッチング機能を提供する。ロジスティック回帰で傾向スコアを計算し Matching::Match() でマッチングを行う流れとなる。

lalonde datasetは, 1976年の米国職業訓練プログラムを受けた群/受けなかった群において, 1978年時の収入にどの程度影響したかに関するデータで今回はこれを用いる。

> help(lalonde)

lalonde {Matching}	R Documentation
Lalonde Dataset

Description

Dataset used by Dehejia and Wahba (1999) to evaluate propensity score matching.

Usage

data(lalonde)
Format

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.

ロジスティック回帰で傾向スコアを計算。

model.glm <- glm(treat ~., data = lalonde[,-9], family = binomial)
summary(model.glm)

高卒以上であるか否かのP値が最も小さくなった。

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.4884  -0.9934  -0.8708   1.2242   1.7403

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.622e+00  1.102e+00   1.472   0.1410
age          8.276e-03  1.452e-02   0.570   0.5687
educ        -8.282e-02  7.230e-02  -1.145   0.2520
black       -2.216e-01  3.684e-01  -0.601   0.5476
hisp        -8.557e-01  5.128e-01  -1.669   0.0952 .
married      1.960e-01  2.784e-01   0.704   0.4813
nodegr      -8.981e-01  3.146e-01  -2.855   0.0043 **
re74        -4.466e-05  3.010e-05  -1.483   0.1380
re75         2.924e-05  4.788e-05   0.611   0.5414
u74         -1.927e-01  3.765e-01  -0.512   0.6088
u75         -3.369e-01  3.213e-01  -1.048   0.2945
---
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: 584.26  on 434  degrees of freedom
AIC: 606.26

Number of Fisher Scoring iterations: 4

傾向スコアを用いてマッチングを行う。従属変数に78年の収入, 独立変数に処置変数, Xに傾向スコアの推定値 (ロジスティック回帰モデルによる推定値)を指定して Match() を実行する。estimand では, 平均処置効果 (average treatment effect; ATE) または 処置群の平均処置効果 (average treatment on treated; ATT), 対照群の平均処置効果 (average treatment effect for the controls; ATC) を指定できる。

mout <- Match(
    Y = lalonde$re78,
    Tr = lalonde$treat,
    X = model.glm$fitted,
  )

処置群 185人 に対して, 対照群は 260人 の中から 173人 のサンプルを使っている。つまり, 重複して使われているサンプルが存在している。
マッチングの結果の Pair は index.treated, index.controlを使って確認できる。

# check pair
lalonde$propensity_score <- model.glm$fitted
treated_group <- lalonde[mout$index.treated,]
control_group <- lalonde[mout$index.control,]
pair <- cbind(treated_group, control_group)

傾向スコアを用いたマッチングによる因果効果の推定量は 2138となった。

Estimate...  2138.6
AI SE......  797.76
T-stat.....  2.6807
p.val......  0.0073468

Original number of observations..............  445
Original number of treated obs...............  185
Matched number of observations...............  185
Matched number of observations  (unweighted).  322

caliper matching を TRUE に指定すると, 最近傍マッチングを行った場合にある特定以上の距離以上になる時はマッチングしないようにできる。

また, MatchBalance() は, マッチングによる共変量調整の結果によって, 共変量の分布が処置群と対照群の2つの群間で近づいているかどうかのチェックできる。

MatchBalance(
    treat ~ .,
    match.out = mout,
    nboots = 1000,
    data = lalonde
  )

下記のように, 各群の共変量の前後の差がわかり共変量調整が確認できる。

***** (V1) age *****
                       Before Matching 	 	 After Matching
mean treatment........     25.816 	 	     25.816 
mean control..........     25.054 	 	      25.87 
std mean diff.........     10.655 	 	   -0.75043 

...

各共変量に対して調整前と調整後での処置群と対照群の平均値の差の計算や検定を行うことができる。

Before Matching Minimum p.value: 0.0020368 
Variable Name(s): nodegr  Number(s): 6 

After Matching Minimum p.value: 0.061 
Variable Name(s): re75  Number(s): 8

IPWE

先に少し触れた傾向スコアを用いたマッチング・層別解析の問題点から, 傾向スコアの逆数を用いて目的変数の値に重み付けを行う IPW推定量 （Inverse Probability Weighting Estimator; IPWE） が提案された。 IPWEでは処置群と対照群の周辺期待値と標準誤差が求めることができる。

傾向スコアを計算するところまでは同様。

# Inverse Probability Weighting Estimator
ivec1 <- lalonde$treat
ivec2 <- rep(1, nrow(lalonde)) - ivec1

ivec <- cbind(ivec1, ivec2)

iestp1 <- (ivec1/model.glm$fitted) * (length(ivec1)/sum(ivec1))
iestp2 <- (ivec2/(1-model.glm$fitted)) * (length(ivec2)/sum(ivec2))

iestp <- iestp1 + iestp2
head(iestp)
 #        1         2         3         4         5         6
 # 6.124465 10.588840  4.526984  7.320260  6.053291  6.745951

ipwe <- lm(lalonde$re78 ~ ivec -1, weight=iestp, data=lalonde)
summary(ipwe)

因果効果の推定量は, 共変量調整した場合の回帰係数の差で 6213 – 4571 = 1642となった。

Weighted Residuals:
   Min     1Q Median     3Q    Max
-17241  -7932  -3538   5383 129799

Coefficients:
          Estimate Std. Error t value Pr(>|t|)
ivecivec1     6213        435  14.284   <2e-16 ***
ivecivec2     4571        514   8.894   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14200 on 443 degrees of freedom
Multiple R-squared:  0.3899,	Adjusted R-squared:  0.3872
F-statistic: 141.6 on 2 and 443 DF,  p-value: < 2.2e-16

おわりに

参考書籍は『 調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)』(1-3章) です。

Code は GitHub にあります。本記事における間違いを見つけられましたら, コメントを頂ければと思います。


[1] 数学セミナー増刊統計学ガイダンスの part3 相関・連関分析から因果分析へ…黒田正博
[2] 相関と因果について考える：統計的因果推論、その(不)可能性の中心
[3] The Central Role of the Propensity Score in Observational Studies for Causal Effects [rosenbaum and rubin, 1983]
[4] 共変量とは従属変数と独立変数の関係を明確にするために, 統制する必要がある変数。従属変数と独立変数の両方に関係がある。統制変数, または医学分野では交絡変数とも呼ばれる。
[5] Chapter 5 Matching Estimators of Causal Effects
[6] Rで学ぶ 傾向スコア解析入門 – 無作為割り当てが出来ない時の因果効果推定 –