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

| 0件のコメント

今回, 参考にさせて頂いたのは調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)の, 1章-3章の内容です。

本投稿における間違いは私の理解不足に依るものなので, 本書の内容の正確性とは一切関係ありません。

研究デザイン

研究デザインを大きく2つに分類・定義している。

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

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

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

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

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

ルービンの因果モデル

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

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

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

傾向スコア

ルービンの因果効果を, もし無作為割り当てを行った場合の処置群と対照群の期待値の差とする。 因果効果 = E(y1) – E(y0) となる。[3]

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

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

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

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

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

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

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

Propensity Score Matching

傾向スコアマッチング (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) [5] または 処理群の平均処置効果 (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

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で学ぶ 傾向スコア解析入門 - 無作為割り当てが出来ない時の因果効果推定 -

コメントを残す

必須欄は * がついています