こんにちは、リブセンスでデータサイエンティストをしている北原です。今回はCovariate Balancing Propensity Score(CBPS)の紹介をします。また、Rのmomentfitパッケージを利用したCBPSの実装も扱います。
CBPSは共変量のバランスも考慮して傾向スコアを算出する方法です。以前の記事で紹介したGeneralized Method of Moments(GMM、一般化モーメント法)を利用しているところも特徴の一つになっています。CBPS論文著者らが開発したRパッケージがあるので簡単に使うことができますが、CBPSは拡張性の高い手法なのでモーメント条件をカスタマイズしたいこともあるでしょう。そこで今回はGMMを実行できるmomentfitパッケージを使ってCBPSの実装を行います。momentfitパッケージを使ったGMMの実行方法については下記の記事をご参照ください。
傾向スコア
観察データを使って因果推論を行う場合、ランダム割付と同等の状況を作り出して比較を行う必要があります。例えば、マッチングやInverse Probability Weighting（IPW、逆確率重み付け法）では、処置群と対照群の共変量のバランスをとることでランダム割付と同等の状況を作り出します。複数の共変量のままでは扱いにくいときにこれらの情報を１変数に集約したいときに使われるのが傾向スコアです。
傾向スコアは共変量で条件づけた処置の確率です。$i$番目の観測データの処置を$T_i \in \{0, 1\}$、共変量を${\mathbf X}_i$で表し、${\mathbf X}_i$以外に割付に影響を与える情報はないとき
を傾向スコアと呼びます。このとき条件付き独立
となります。$Y_i(1)$は処置があったとき$Y_i(0)$は処置がなかったときの潜在的な結果です。要するに割付の情報を全て集約したものが傾向スコアです。傾向スコアが同じならばランダム割付と同じ状況になっていると考えられます。
傾向スコアに該当する情報を入手できない場合は推定する必要があります。推定する方法はいろいろありますが、ドメイン知識によるモデルの調整がしやすいロジスティック回帰が使われることが多いようです。傾向スコアを利用するときに算出に用いたモデルの解釈性が必要になることはないのでノンパラメトリックな機械学習手法が使われることもあります。
ランダム割付と同等の状況を作り出すのが目的であることをふまえると、傾向スコアの推定では処置群と対照群の共変量の分布を等しくできる方法を用いることが望ましいです。処置確率の予測精度が高いほど共変量のバランスがよくなるとは限らないので、予測性能の高い手法を使えばよいわけではありません。
Covariate Balancing Propensity Score
Covariate Balancing Propensity Score(CPBS)は、共変量のバランスをとりつつ処置確率を推定する方法です。理論上は傾向スコアで条件付ければ共変量のバランスはとれるのですが、傾向スコアの推定を適切にできるとは限らないため使う共変量やモデルによっては処置確率を推定しただけでは共変量のバランスがとりにくいことがあります。CBPSは共変量のバランスを推定時の条件に入れることで、この問題に対処しています。
CBPSの特徴は誤ったモデルを使ったときでも通常のロジスティック回帰と比較して安定した性能が得られることです。モデルが誤っていると尤度を最大化しても共変量のバランスがとりにくくなるため、尤度を多少犠牲にしても共変量のバランスを改善するCBPSが相対的に有利になると考えられます。このことはシミュレーション実験とLaLondeのデータを使った検証で示されています。そのため現象を適切にモデリングしにくいときに適した手法と考えられます。ただし、どのようなモデルを使うかで推定結果が大きく変わることも示されており、CBPSを使ってもバイアスの小さい推定ができるとは限らないことには注意が必要です。
実際に傾向スコアを使う場合はこちらにもあるように地道に共変量のバランスをとる必要があるので大変ですが、共変量のバランスも考慮されているCBPSを使うことで負担が減ることがあります。Google ScholarなどでCBPS論文を引用している文献を調べると実際の分析で使っているケースがいくつか見つかりますし、CBPSの利用を推奨しているチュートリアルもあります。国内でCBPSが使われている例としてはこちらがあります。また、2016年のAtlantic Causal Inference Conferenceのコンペでも使われていて、上位ではないもののパラメトリックな手法の中では好成績であることがわかります。
CBPSでは、バランスをとりたい共変量の関数について傾向スコアの逆数で重み付けしたものをモーメント条件としてGMMでパラメータ推定を行います。パラメータ数よりバランスをとりたい関数の数のほうが多いことがあるのでGMMが使われています。処置確率の推定と共変量のバランシングを同時に行う方法と考えることもできます。処置確率の推定あるいは共変量のバランシングの一方だけでパラメータ数と同数のモーメント条件を使うので、双方を同時に行うときはGMMが必要になります。なお、GMMについては過去の記事で紹介したのでそちらをご覧ください。
パラメータを${\mathbf \beta} \in {\mathbb R}^L$、パラメトリックな傾向スコアモデルを$\pi_{\mathbf \beta}({\mathbf X}_i)$とするとき
を仮定します。このとき共変量バランシングのモーメント条件
を最小化することで傾向スコアを算出する手法がCBPSです。ここで$\tilde{\mathbf X}_{i} = f({\mathbf X}_i)$は釣り合わせたい共変量の関数の長さ$M$のベクトルです。シンプルには$\tilde{\mathbf X}_{i} = {\mathbf X}_i$と指定すれば、共変量の一次モーメントをバランシングさせることになります。論文のシミュレーション実験やCBPSパッケージでは一次モーメントのバランシングが使われています。どのような関数をバランシングさせればよいかはまだ明らかにされていないので、この部分をカスタマイズしたい場合は試行錯誤が必要です。
通常は、共変量の一次モーメントのバランシングに加えて、傾向スコアの推定を行うスコアのモーメント条件も使います。傾向スコアモデルにロジスティック回帰モデルを使う場合は
となります。データ数を$N$として最尤法を使うと
でパラメータを推定できます。${\mathbf \beta}$で微分し整理すると
となります。($\ref{s_moment}$)は標本モーメントですし、($\ref{cbps_s}$)は($\ref{cbps_moment}$)においてモデルの１階微分$\pi^{'}_{\mathbf \beta}({\mathbf X}_i)=\partial \pi({\mathbf X}_i)/\partial {\mathbf \beta}$を釣り合わせるのと同じです。つまり、ロジスティック回帰モデルは$\tilde{\mathbf X}_{i} = \pi^{'}_{\mathbf \beta}({\mathbf X}_i)$とした共変量バランシングと解釈できます。なお、ロジスティック回帰のように($\ref{cbps_moment}$)の形で表現できなくとも、モーメント条件がわかりGMMで確率を推定できるモデルならば傾向スコアのモデルとして利用できると考えられます。
共変量バランシングのモーメント条件($\ref{cbps_moment}$)を少し変形して標本モーメントの形で書くと
となります。なお、これはAverage Treatment Effect（ATE）のケースです。
パラメータ${\mathbf \beta}$は($\ref{s_moment}$)と($\ref{w_moment}$)をモーメント条件の標本モーメントとしてGMMを実行することで推定できます。
重み行列には、効率的な重みと呼ばれる$g_{\beta}\left(T_{i}, {\mathbf X}_{i}\right)$の分散共分散行列の逆行列ではなく
を使っています。期待値は処置$T_i$についてとっていることに注意し($\ref{ps}$)を条件付き確率として計算すると
となります。この重みでは処置の判別が難しいときはスコアのモーメント条件の重みが相対的に大きくなり、逆に簡単なときは共変量バランシングの重みが相対的に大きくなっていることがわかります。
通常のGMMと同様にJ検定も使えます。J検定ではモーメント条件を満たせないと帰無仮説が棄却されます。傾向スコアのモデルが適切ではない場合にはモーメント条件を満たせなくなるため帰無仮説が棄却されます。CBPSはモデルが誤っていても安定した結果が得られることが特徴なので、帰無仮説が棄却されるようなスコアモデルがデータに適合しないケースでも単純に推定失敗とはみなしません。
標準化と特異値分解
モーメント条件と重み行列がわかれば実装はできるものの、推定が不安定になりやすいという問題があります。論文には書いていないのですが、Rパッケージのソースコードを見ると標準化と特異値分解による前処理が行われており、これらが推定の安定化に必要と考えられます。そこで本節では標準化と特異値分解を使った前処理と後処理について紹介します。
標準化と特異値分解のいずれも推定を安定化させたり予測精度を向上させる前処理としてよく知られています。CBPSのスコアは前処理後のデータから推定されたパラメータを使えば得られるのですが、推定されたパラメータを知りたい場合は前処理後のデータから推定されたパラメータを前処理前のデータから得られるパラメータに逆変換する必要があります。
処理の流れとしては以下のようになります。
- データを標準化
- 標準化したデータを特異値分解
- 分解したデータを使ってGMM
- GMMで得られたパラメータを特異値分解前のデータでのパラメータに逆変換
- 特異値分解前のデータでのパラメータを標準化前のデータでのパラメータに逆変換
前処理前のデータから得られるパラメータへの逆変換はあまり知られていないと思うので、ここでは逆変換を中心に説明します。標準化と特異値分解それぞれについて説明し、最後にそれらを組み合わせたものを紹介します。
ロジスティック回帰モデルは非線形モデルなので標準化前のパラメータへの逆変換は自明ではなく、いくつかの方法があるようです。CBPSパッケージでは、こちらで議論されているAgresti methodが使われています。
$j$番目のデータ$x_j$について、平均$\mu_j$と標準偏差$\sigma_j$を使って
と標準化したデータ${\mathbf x}^{\ast}$から得られたロジスティック回帰のパラメータを${\mathbf \beta}^{\ast}_j$とすると、標準化前のパラメータ${\mathbf \beta}_j$は次のようにして得られます。
なお$\beta_0$は定数項のパラメータです。定数項は標準化しないので、定数項のパラメータだけ計算方法が異なります。
サンプルデータで確認してみましょう。次のようなサンプルデータを使います。
library(tidyverse) # サンプルデータ set.seed(1) n <- 1000 x <- tibble( x1 = rnorm(n, mean = 1, sd = 5), x2 = rnorm(n, mean = 1, sd = 5), x3 = rnorm(n, mean = 1, sd = 5)) %>% mutate( y = rbernoulli(n, p = 1 / (1 + exp(-(x1 + 2 * x2 + 3 * x3)))) %>% as.integer() ) # 標準化前のデータから得られるパラメータ fit_glm <- glm(y ~ x1 + x2 + x3, x, family = binomial) coef_glm <- coefficients(fit_glm)
標準化後のデータでロジスティック回帰を行い得られたパラメータを逆変換して、標準化前のパラメータと比較します。
# 標準化 x_std <- x %>% dplyr::select(x1, x2, x3) %>% scale() x_sig <- attributes(x_std)[["scaled:scale"]] x_mean <- attributes(x_std)[["scaled:center"]] x_std <- cbind(x$y, x_std) %>% data.frame() colnames(x_std)[1] <- 'y' # 標準化後のデータでロジスティック回帰 fit_glm_std <- glm(y ~ x1 + x2 + x3, x_std, family = binomial) summary(fit_glm_std) coef_std <- coefficients(fit_glm_std) # 逆変換で元のパラメータに戻す rev_coef_std <- coef_std[-1] / x_sig rev_coef_std0 <- coef_std[1] - sum(rev_coef_std * x_mean) rev_coef_std <- c(rev_coef_std0, rev_coef_std) print( data.frame( 元のパラメータ = coef_glm, 標準化後のパラメータ = coef_std, 元に戻したパラメータ = rev_coef_std ) )
元のパラメータ 標準化後のパラメータ 元に戻したパラメータ (Intercept) 0.2855891 6.300383 0.2855891 x1 1.0299750 5.329687 1.0299750 x2 1.9463096 10.120630 1.9463096 x3 3.0251862 15.596453 3.0251862
元のパラメータと同じになっており、標準化の逆変換ができていることがわかります。
次に特異値分解の逆変換を紹介します。データの行列$X$を
のように特異値分解し、$U$から得られたロジスティック回帰のパラメータを${\mathbf \beta}^{\ast}$とすると、特異値分解前のパラメータは
となります。対数尤度を微分して0とおいた式にこれらを代入すると成り立つことが確認できます。
特異値分解についてもデータで確認してみましょう。元のパラメータと同じになっており逆変換ができていることがわかります。
# 特異値分解 x_svd <- x %>% mutate(intercept = 1) %>% dplyr::select(intercept, x1:x3) %>% svd() # 特異値分解後のデータでロジスティック回帰 fit_glm_svd <- glm(x$y ~ . -1, data.frame(x_svd$u), family = binomial) summary(fit_glm_svd) coef_svd <- coefficients(fit_glm_svd) # 逆変換で元のパラメータに戻す rev_coef_svd <- c(x_svd$v %*% diag(1 / x_svd$d) %*% coef_svd) print( data.frame( 元のパラメータ = coef_glm, 特異値分解後のパラメータ = coef_svd, 元に戻したパラメータ = rev_coef_svd ) )
元のパラメータ 特異値分解後のパラメータ 元に戻したパラメータ (Intercept) 0.2855891 -621.917035 0.2855891 x1 1.0299750 -81.327997 1.0299750 x2 1.9463096 -189.719793 1.9463096 x3 3.0251862 -2.581429 3.0251862
最後に標準化と特異値分解の組み合わせについて見てみましょう。
# 逆変換の関数 rev_trans_coef <- function(coef, v, d, mu, sig) { coef_ <- c(v %*% diag(1 / d) %*% coef) rev_coef <- coef_[-1] / sig rev_coef0 <- coef_[1] - sum(rev_coef * mu) return(c(rev_coef0, rev_coef)) } # 標準化 x_std <- x %>% dplyr::select(x1, x2, x3) %>% scale() x_sig <- attributes(x_std)[["scaled:scale"]] x_mean <- attributes(x_std)[["scaled:center"]] # 特異値分解 x_svd <- cbind(1, x_std) %>% svd() # ロジスティック回帰 fit_glm_std_svd <- glm(x$y ~ . - 1, data.frame(x_svd$u), family = binomial) coef_std_svd <- coef(fit_glm_std_svd) # 逆変換 rev_coef <- rev_trans_coef(coef_std_svd, x_svd$v, x_svd$d, x_mean, x_sig) print( data.frame( 元のパラメータ = coef_glm, 前処理後のパラメータ = coef_std_svd, 元に戻したパラメータ = rev_coef ) )
元のパラメータ 前処理後のパラメータ 元に戻したパラメータ (Intercept) 0.2855891 573.0577 0.2855891 x1 1.0299750 -199.2356 1.0299750 x2 1.9463096 -183.3706 1.9463096 x3 3.0251862 -166.3748 3.0251862
CBPSの実装
CBPSの実装にはRのmomentfitパッケージを利用します。momentfitパッケージを使ったGMMの実行方法については過去の記事をご参照ください。なお、論文での検証やRパッケージと同じく、共変量バランシングの関数には１次モーメントを使います。
momentfitパッケージを使ってGMMを実行するにはモーメント条件の関数を作成する必要があります。($\ref{gmm_moment}$)をそのまま実装すればよいので次のようになります。$s_{\mathbf \beta}\left(T_{i}, {\mathbf X}_{i}\right)$はロジスティック回帰のモーメント条件と同じです。
# モーメント条件の関数 g <- function(theta, X_) { tr <- X_[, 1] X <- X_[, -1] p <- c(1 / (1 + exp(-X %*% theta))) # s w <- (tr - p) / (p * (1 - p)) # w return( cbind( (tr - p) * X, w * X )) }
GMMでは重み行列に($\ref{weight}$)を使うので、重み行列の作成も実装する必要があります。重み行列は次のように作成できます。
# 重み行列を作成する関数 cbps_weight <- function(x_, theta) { p <- c(1 / (1 + exp(-x_ %*% theta))) x_p <- sqrt(p * (1 - p)) * x_ # 行列演算時にp * (1 - p)になるようsqrt() x_rp <- x_ / sqrt(p * (1 - p)) # 行列演算時に1/(p * (1 - p)になるようsqrt() S_11 <- t(x_p) %*% x_p S_12 <- t(x_) %*% x_ S_21 <- S_12 S_22 <- t(x_rp) %*% x_rp return( rbind( cbind(S_11, S_12), cbind(S_21, S_22) ) / nrow(x_) ) }
ではサンプルデータを使って実際にパラメータを推定してみましょう。次のようなサンプルデータを使います。t1が処置の有無、yが結果です。処置、結果のいずれも説明変数の非線形関数で生成されるようにしています。CBPSはモデルの誤りに影響を受けにくい推定ができるのが特徴なので、非線形関数で生成されたデータに対して、傾向スコアと結果に線形モデルを指定して推定したときにどうなるかを見てみましょう。
# データ作成 set.seed(1) n <- 1000 x <- tibble( x1 = rnorm(n, mean = 1, sd = 5), x2 = rnorm(n, mean = 1, sd = 5), x3 = rnorm(n, mean = 1, sd = 5), x4 = rnorm(n, mean = 1, sd = 5)) %>% mutate( t1 = rbernoulli(n, p = 1 / (1 + exp(5 * x2^2 - 2 * x3^3 - 3 * exp(x4)))) %>% as.integer() ) %>% mutate(y = rnorm(n, mean = 1 + 2 * x1^2 + 3 * x2 + 4 * x3 + 5 * t1, sd = 10)) # 処置ごとの平均 x %>% group_by(t1) %>% summarise_all(mean) %>% print()
バイアスのあるデータになっていることがわかります。
# A tibble: 2 x 6
t1 x1 x2 x3 x4 y
<int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 0.738 1.07 -1.97 -1.61 53.5
2 1 1.11 0.794 3.54 3.26 75.9
まずは一般的なロジスティック回帰による傾向スコアの結果を見てみましょう。ここではweightItパッケージのweightit()関数を利用して共変量のバランス確認と効果推定を行います。共変量のバランス確認にはcobaltパッケージのtab.tab()関数、推定にはlmtestパッケージのcoeftest()関数、推定時に頑健標準誤差を計算するのにsandwitchパッケージのvcovHC()関数を利用します。処置の生成にはx1は寄与していませんが、誤ったモデルを指定したときの結果を見るためx1も使っています。
# 共変量のバランス確認と効果推定で利用するパッケージ library(WeightIt) library(cobalt) library(lmtest) library(sandwich) # ロジスティック回帰 fit_glm <- glm(t1 ~ x1 + x2 + x3 + x4, x, family = binomial) print(fit_glm$coefficients) ps_glm <- fit_glm$fitted.values w.out_glm <- weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), ps = ps_glm, data = x) # 共変量のバランス確認 print(bal.tab(w.out_glm)) # ATEの推定 fit_outcome_glm <- lm(y ~ t1, x, weights = w.out_glm$weights) print(coeftest(fit_outcome_glm, vcov. = vcovHC))
(Intercept) x1 x2 x3 x4 -0.435056936 -0.008405366 -0.043157239 0.472980076 0.414328802
Call
weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), data = x,
ps = ps_glm)
Balance Measures
Type Diff.Adj
prop.score Distance -0.1987
x1 Contin. -0.1285
x2 Contin. 0.2263
x3 Contin. 0.3100
x4 Contin. -0.5140
Effective sample sizes
Control Treated
Unadjusted 447. 553.
Adjusted 50.91 34.17
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.6785 9.3987 5.3921 8.695e-08 ***
t1 13.7539 11.3649 1.2102 0.2265
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
共変量のバランスはbal.tab()の出力のDiff.Adjで確認します。バランスがとれているとDiff.Adjの絶対値が小さくなります。Diff.Adjの絶対値が0.1を下回るとバランスがとれていると判断することが多いです。この結果では、いずれの共変量でも0.1を超えているのでバランスはよくないと判断できます。
平均処置効果の推定値と標準誤差はcoeftest()の出力のt1のところで確認します。推定された平均処置効果は13.75とデータ作成時に指定した5よりも大きめなことがわかります。なお、共変量のバランス確認方法や推定方法については後日別記事で紹介予定です。
次にCBPSでスコアを算出してみましょう。前節で説明したように標準化と特異値分解の前処理を行ってからGMMを実行し、逆変換でパラメータを元に戻します。重み行列を指定するのでevalWeights()関数で重み行列のオブジェクトを作成し、gmmFit()関数呼び出し時に渡します。初期値にはロジスティック回帰の結果を使います。CBPSでも誤ったモデルを指定したときの結果を見るためx1も使っています。
library(momentfit) # 標準化 x_std <- x %>% dplyr::select(x1, x2, x3, x4) %>% scale() x_sig <- attributes(x_std)[["scaled:scale"]] x_mean <- attributes(x_std)[["scaled:center"]] # 特異値分解 x_svd <- cbind(1, x_std) %>% svd() # 重み行列 Sig <- cbps_weight(x_svd$u, fit_glm$coefficients) # GMM fit_glm <- glm(x$t1 ~ . -1, data.frame(x_svd$u), family = binomial) model <- momentModel( g, cbind(x$t1, x_svd$u), theta0 = fit_glm$coefficients) wObj <- evalWeights(model, w = Sig) res <- gmmFit(model, weights = wObj) print(summary(res)) # 係数の逆変換 cbps_coefs <- rev_trans_coef(coef(res), x_svd$v, x_svd$d, x_mean, x_sig) print(cbps_coefs)
GMMについては次のような推定結果が得られます。J検定の帰無仮説が棄却されているのでモデル指定が誤っている可能性があることがわかります。
Model based on moment conditions
*********************************
Moment type: function
Covariance matrix: iid
Number of regressors: 5
Number of moment conditions: 10
Number of Endogenous Variables: 0
Sample size: 1000
Estimation: One-Step GMM with fixed weights
Convergence Optim: 1
Sandwich vcov: TRUE
coefficients:
Estimate Std. Error t value Pr(>|t|)
X1 73.3472 6.2428 11.7492 < 2.2e-16 ***
X2 -19.3409 5.1822 -3.7322 0.0001898 ***
X3 -3.7984 4.3216 -0.8789 0.3794391
X4 -23.2500 4.5142 -5.1504 2.599e-07 ***
X5 72.1747 8.1686 8.8356 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
J-Test
Statistics df pvalue
Test E(g)=0: 12.93 5 0.024041
逆変換で得られたパラメータは次のとおりです。ロジスティック回帰とは少し異なる値になっています。
x1 x2 x3 x4 -0.31403992 -0.02975958 -0.02459620 0.51202782 0.38985264
共変量のバランスと推定結果は以下のとおりです。
# CBPSの出力 cbps_scores <- c(1 / (1 + exp(-as.matrix(cbind(1, x %>% dplyr::select(x1:x4))) %*% cbps_coefs))) w.out_cbps <- weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), ps = cbps_scores, data = x) # 共変量のバランス確認 print(bal.tab(w.out_cbps)) # ATEの推定 fit_outcome_cbps <- lm(y ~ t1, x, weights = w.out_cbps$weights) print(coeftest(fit_outcome_cbps, vcov. = vcovHC))
Call
weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), data = x,
ps = cbps_scores)
Balance Measures
Type Diff.Adj
prop.score Distance -0.1244
x1 Contin. -0.0019
x2 Contin. 0.0011
x3 Contin. 0.0006
x4 Contin. -0.0017
Effective sample sizes
Control Treated
Unadjusted 447. 553.
Adjusted 48.61 72.24
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.2452 8.7709 5.7286 1.34e-08 ***
t1 8.9609 11.4917 0.7798 0.4357
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
傾向スコアのDiff.Adjの絶対値は0.1をわずかに上回っているものの、いずれの共変量もDiff.Adjの絶対値は0.1を下回っており、推定値もロジスティック回帰と比較するとデータ作成時に設定した5に近い値になっていることがわかります。
momentfitパッケージには複数の推定アルゴリズムが実装されているので使ってみましょう。ここでは推定を収束するまで繰り返すIterated GMMと、重み行列とパラメータの推定を一度に行うContinuously Updated GMMを使います。Iterated GMMは次のようになります。
# Iterated GMM res <- gmmFit(model, type = "iter") cbps_itr_coefs <- rev_trans_coef(coef(res), x_svd$v, x_svd$d, x_mean, x_sig) cbps_itr_scores <- c(1 / (1 + exp(-as.matrix(cbind(1, x %>% dplyr::select(x1:x4))) %*% cbps_itr_coefs))) w.out_cbps_itr <- weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), ps = cbps_itr_scores, data = x) # 共変量のバランス確認 print(bal.tab(w.out_cbps_itr)) # ATEの推定 fit_outcome_itr <- lm(y ~ t1, x, weights = w.out_cbps_itr$weights) print(coeftest(fit_outcome_itr, vcov. = vcovHC))
Call
weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), data = x,
ps = cbps_itr_scores)
Balance Measures
Type Diff.Adj
prop.score Distance -0.4702
x1 Contin. -0.0937
x2 Contin. -0.1338
x3 Contin. -0.0339
x4 Contin. -0.3820
Effective sample sizes
Control Treated
Unadjusted 447. 553.
Adjusted 47.2 35.36
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.2506 8.2712 6.1963 8.442e-10 ***
t1 5.1796 11.9148 0.4347 0.6639
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
推定値はデータ生成時に指定した5に近い結果が得られていますが、Diff.Adjの絶対値は0.1を超えており、ロジスティック回帰ほどではないですが共変量バランスはよくありません。偶然よい推定値が得られたと考えたほうがよいでしょう。Continuously Updated GMMは次のようになります。
# Continuously Updated GMM res <- gmmFit(model, type = "cue") cbps_cue_coefs <- rev_trans_coef(coef(res), x_svd$v, x_svd$d, x_mean, x_sig) cbps_cue_scores <- c(1 / (1 + exp(-as.matrix(cbind(1, x %>% dplyr::select(x1:x4))) %*% cbps_cue_coefs))) w.out_cbps_cue <- weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), ps = cbps_cue_scores, data = x) # 共変量のバランス確認 print(bal.tab(w.out_cbps_cue)) # ATEの推定 fit_outcome_cue <- lm(y ~ t1, x, weights = w.out_cbps_cue$weights) print(coeftest(fit_outcome_cue, vcov. = vcovHC))
Call
weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), data = x,
ps = cbps_cue_scores)
Balance Measures
Type Diff.Adj
prop.score Distance -0.4537
x1 Contin. -0.1190
x2 Contin. -0.0136
x3 Contin. 0.1107
x4 Contin. -0.5623
Effective sample sizes
Control Treated
Unadjusted 447. 553.
Adjusted 47.98 29.13
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.9046 8.6993 5.8516 6.597e-09 ***
t1 8.7252 11.7225 0.7443 0.4569
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Diff.Adjの絶対値は0.1を超えており、ロジスティック回帰よりはよいものの、CBPSの重み行列を使ったときほど共変量のバランスはよくありません。この結果だけで判断はできないもののCBPSの重み行列は使ったほうがよいようです。
念のためCBPSパッケージを使ったときの結果も確認しましょう。
# CBPSパッケージ library(CBPS) fit_cbps_pkg <- CBPS(t1 ~ x1 + x2 + x3 + x4, x, ATT = 0) cbps_pkg_scores <- fit_cbps_pkg$fitted.values w.out_cbps_pkg <- weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), ps = cbps_pkg_scores, data = x) # 共変量のバランス確認 print(bal.tab(w.out_cbps_pkg)) # ATEの推定 fit_outcome_pkg <- lm(y ~ t1, x, weights = w.out_cbps_pkg$weights) print(coeftest(fit_outcome_pkg, vcov. = vcovHC))
Call
weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), data = x,
ps = cbps_pkg_scores)
Balance Measures
Type Diff.Adj
prop.score Distance 0.0178
x1 Contin. -0.0931
x2 Contin. 0.1682
x3 Contin. 0.3106
x4 Contin. -0.2647
Effective sample sizes
Control Treated
Unadjusted 447. 553.
Adjusted 62.31 51.25
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.3962 8.3765 6.1357 1.221e-09 ***
t1 12.8225 10.2224 1.2543 0.21
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
数値計算上の問題なのか、傾向スコアのバランスはよいものの、共変量のバランスが悪い結果が得られました。基本的には同じアルゴリズムを使っているはずですが、細かな実装は異なっているために違いが生じたのかもしれません。例えば最適化で使うBFGSの実行時に、CBPSパッケージは微分関数を指定していますが、momentfitパッケージではデフォルトで数値微分を使うなどといった違いがあります。なお、共変量バランシングのモーメント条件($\ref{cbps_w}$)のみを使うと丁度識別のCBPSを計算することもできます。今回momentfitパッケージを使って推定した結果はCBPSパッケージで丁度識別のCBPSを推定した結果に近く、ロジスティック回帰のモーメント条件($\ref{cbps_s}$)が十分機能していない可能性があります。
最後に標準化と特異値分解の前処理をしない場合の結果を見てみましょう。
# 標準化と特異値分解をしない場合 fit_glm <- glm(t1 ~ x1 + x2 + x3 + x4, x, family = binomial) x_ <- x %>% mutate(intercept = 1) %>% dplyr::select(intercept, x1:x4) %>% as.matrix() Sig <- cbps_weight(x_, fit_glm$coefficients) model <- momentModel( g, cbind(x$t1, x_), theta0 = fit_glm$coefficients) wObj <- evalWeights(model, w = Sig) res <- gmmFit(model, weights = wObj) print(summary(res)) cbps_nopreproc_scores <- c(1 / (1 + exp(-as.matrix(cbind(1, x %>% dplyr::select(x1:x4))) %*% coef(res)))) w.out_cbps_nopreproc <- weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), ps = cbps_nopreproc_scores, data = x) # 共変量のバランス確認 print(bal.tab(w.out_cbps_nopreproc)) # ATEの推定 fit_outcome_nopreproc <- lm(y ~ t1, x, weights = w.out_cbps_nopreproc$weights) print(coeftest(fit_outcome_nopreproc, vcov. = vcovHC))
Model based on moment conditions
*********************************
Moment type: function
Covariance matrix: iid
Number of regressors: 5
Number of moment conditions: 10
Number of Endogenous Variables: 0
Sample size: 1000
Estimation: One-Step GMM with fixed weights
Convergence Optim: 0
Sandwich vcov: TRUE
coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.312492 0.174999 -1.7857 0.07415 .
x1 -0.030253 0.034534 -0.8760 0.38101
x2 -0.024407 0.028927 -0.8438 0.39881
x3 0.512560 0.044645 11.4809 < 2e-16 ***
x4 0.389714 0.037976 10.2622 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
J-Test
Statistics df pvalue
Test E(g)=0: 12.967 5 0.023689
Call
weightit(formula = formula(t1 ~ x1 + x2 + x3 + x4), data = x,
ps = cbps_nopreproc_scores)
Balance Measures
Type Diff.Adj
prop.score Distance -0.1242
x1 Contin. -0.0002
x2 Contin. 0.0000
x3 Contin. -0.0002
x4 Contin. -0.0003
Effective sample sizes
Control Treated
Unadjusted 447. 553.
Adjusted 48.61 72.37
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.2236 8.7718 5.7256 1.363e-08 ***
t1 8.9178 11.5120 0.7747 0.4387
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
今回のサンプルデータでは共変量のバランスもよく前処理をしたときとほぼ同じ推定結果になっており問題ないように見えます。しかし、重み行列を確認すると極端に大きな値が含まれており、最適化計算時に不安定になりやすいことがわかります。
# 重み行列の要素を大きい順に5個表示 print(c(Sig) %>% sort(decreasing = TRUE) %>% head(10))
[1] 37609.36 35418.97 32109.23 32109.23 10106.04
まとめ
今回はCBPSという共変量のバランスも考慮して最適化された傾向スコアについて紹介しました。共変量のバランシングに使う関数を分析者が選ぶことができる拡張性の高い手法なので、Rのmomentfitパッケージを利用した実装についても紹介しました。共変量のバランシングに使う関数にどのようなものを使えばよいかは知られていませんが、例えば結果変数への寄与の大きい共変量のバランシングには強い条件を設定するなどの拡張を行うことで、よりバイアスの小さい推定が可能になると考えられます。
