傾向スコアの基礎
UPDATE: 2022-12-19 20:26:53
はじめに
このノートではRを使って因果推論に関する基礎的な理論から分析を実行するための方法をまとめている。ここでは傾向スコアについてまとめておく。
傾向スコアの数理
傾向スコアの数理については下記が詳しいので、ここでは扱わない。
傾向スコア
傾向スコアとバランシング
傾向スコアとIPW
傾向スコアの層別解析
傾向スコアとマッチング
準備
必要なライブラリを読み込んでおく。
library(tidyverse)
library(broom)
library(MatchIt)
library(WeightIt)
library(tableone)
library(cobalt)
library(lmtest)
library(sandwich)ここで使用するサンプルデータは、できる!傾向スコア分析SPSS・Stata・Rを用いた必勝マニュアルで利用されているサンプルデータをお借りする。ここでは、このデータをもとに、書籍では紹介されていないMatchItやWeightItなどのパッケージを使って傾向スコアを用いた分析手法をまとめておく。
# id: 患者ID
# sequela_y:退院時後遺症(0/1)
# adl_y: 退院時ADLスコア(0-100)
# treat: 治療薬の処置(0/1)
# sex: 男性=1
# age: 年齢
# ht: 高血圧(0/1)
# dm: 糖尿病(0/1)
# stroke: 脳卒中の既往(0/1)
# mi: 心筋梗塞の既往(0/1)
df <- read_csv('~/Desktop/PSbook_data.csv', show_col_types = FALSE, locale = locale(encoding = "shift-jis")) %>%
select(id, sequela_y = sequela, adl_y = ADL_disc, treat = TreatmentX, sex, age = Age, ht = HT, dm = DM, stroke = Stroke, mi = MI) %>%
mutate(sex = if_else(sex == '男', 1, 0))
df## # A tibble: 15,000 × 10
## id sequela_y adl_y treat sex age ht dm stroke mi
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 45 0 1 74 0 0 0 0
## 2 2 1 40 0 0 73 1 0 0 1
## 3 3 0 40 0 1 75 0 0 0 0
## 4 4 0 40 0 1 76 0 0 0 0
## 5 5 0 40 0 0 80 0 0 0 0
## 6 6 0 40 0 1 76 0 0 0 0
## 7 7 0 40 0 0 74 0 0 0 0
## 8 8 0 40 0 0 78 0 0 0 0
## 9 9 0 45 0 0 72 0 0 0 0
## 10 10 0 45 1 0 75 0 0 0 0
## # … with 14,990 more rows基本的な統計量を確認しておく。tableoneパッケージのCreateTableOne()関数が便利。これを見る限り、age
、ht、dm、stroke、miなどに差があることがわかる。この状態で比較しても、効果が交絡している状態となるので、本来知りたい処置の効果を推定できない。
CreateTableOne(
data = df,
vars = c("sequela_y", "adl_y", "sex", "age", "ht", "dm", "stroke", "mi"),
factorVars = c("sequela_y", "sex", "ht", "dm", "stroke", "mi"),
strata = "treat",
test = FALSE
)## Stratified by treat
## 0 1
## n 9653 5347
## sequela_y = 1 (%) 839 ( 8.7) 403 ( 7.5)
## adl_y (mean (SD)) 43.83 (4.97) 38.37 (6.14)
## sex = 1 (%) 5825 (60.3) 3215 (60.1)
## age (mean (SD)) 72.59 (5.18) 79.35 (4.82)
## ht = 1 (%) 1170 (12.1) 2522 (47.2)
## dm = 1 (%) 487 ( 5.0) 1554 (29.1)
## stroke = 1 (%) 110 ( 1.1) 877 (16.4)
## mi = 1 (%) 298 ( 3.1) 1140 (21.3)何も考慮せずに処置効果を計算してみると、処置を行うことで、退院時ADLスコアadl_yが低下するという結果になっており、
tidy(lm(adl_y ~ treat, data = df))[2,2]## # A tibble: 1 × 1
## estimate
## <dbl>
## 1 -5.46退院時後遺症sequela_yは起こりにくいもものオッズ比が1に近い結果となっている。
exp(tidy(glm(formula = sequela_y ~ treat, family = binomial(link = "logit"), data = df))[2,2])## # A tibble: 1 × 1
## estimate
## <dbl>
## 1 0.856まずは傾向スコアを計算する。MatchItやWeightItパッケージを使えば、傾向スコアを計算する必要はないが、ここでは勉強のために計算し、可視化しておく。
傾向スコアの分布が2群で重なっていることがわかる。この重なりがない状態だと共有サポート(common surpport)を受けられず、場合に応じてはATTやATEが計算できない。共有サポートがないのであれば、比較する群に似ている個体がない状態なので、比較できず、因果推論もできない。
formula_ps <- formula("treat ~ age + sex + ht + dm + stroke + mi" )
fit_ps <- glm(
formula = formula_ps,
data = df,
family = binomial(link = "logit")
)
df$ps <- fit_ps$fitted.values
ggplot(df, aes(x = ps, col = factor(treat), fill = factor(treat))) +
geom_histogram(binwidth = 0.1, alpha = .5, position = "identity") +
labs(x = "傾向スコア", y = "カウント", fill = "処置") +
guides(col = "none") +
theme_classic()
処置ごとにわけて傾向スコアを可視化するとこのようになる。対照群(=0)では傾向スコアが大きい個体があまりいないことがわかる。
ggplot(df, aes(x = ps)) +
geom_histogram(color = "black", binwidth = 0.1, position = "identity") +
facet_grid(rows = vars(treat), scales = "free_y") +
labs(x = "傾向スコア", y = "カウント") +
theme_classic()
傾向スコアと重み付け法
傾向スコアを利用した重み付けで、因果効果を推定する方法をまとめる。推定対象(estimand)がATEかATEなのかによって重みの計算式が異なるので注意。ATTを推定するために、傾向スコアは下記の通り。
\[ w_{i}^{ATT} = Z_{i} + (1 - Z_{i}) \frac{e_{i}(X)}{1 - e_{i}(X)} \] この重みは処置群(Z=1) の個体の重みは1となるので、処置群の個体には、傾向スコアに関係なく等しい重みがつけられる。一方で、処置群(Z=1)の個体の重みは傾向スコアが高いほど、大きな重みがつけられるようになる。イメージとして、傾向スコアが高く処置群に入る確率が高いにも拘らず対照群に入った個体は、処置群の比較対象として重みが大きなり、傾向スコアが低い個体は、処置群の比較対象としてあまり重要ではないので、重みが小さくなる。
ATEを推定するための重みはIPW(inverse probability weighting)という方法を利用する。処置群と対照群について、その群に割り付けられる傾向スコアの逆数で重みが決まるため、それぞれの群で珍しい個体は重みが大きくなる。
\[ w_{i}^{ATE} = \frac{Z_{i}}{e_{i}(X)} + \frac{(1 - Z_{i}) }{1 - e_{i}(X)} \]
ATTとATEの重みを計算する場合は下記の通り。処置と対照に合わせて傾向スコアを利用する。
df_weight <- df %>%
mutate(
weight_att = if_else(treat == 1, 1 , ps/(1-ps)),
weight_ate = if_else(treat == 1, 1/ps, 1/(1-ps))
)
df_weight %>%
select(treat, ps, weight_att, weight_ate)## # A tibble: 15,000 × 4
## treat ps weight_att weight_ate
## <dbl> <dbl> <dbl> <dbl>
## 1 0 0.220 0.283 1.28
## 2 0 0.266 0.362 1.36
## 3 0 0.263 0.358 1.36
## 4 0 0.312 0.453 1.45
## 5 0 0.537 1.16 2.16
## 6 0 0.312 0.453 1.45
## 7 0 0.220 0.282 1.28
## 8 0 0.420 0.723 1.72
## 9 0 0.150 0.176 1.18
## 10 1 0.263 1 3.80
## # … with 14,990 more rows処置群に対する平均処置効果ATTの推定
さきほどのように手計算しなくても、WeightItパッケージのweightit()関数で、傾向スコアと重みの計算が一気通貫でできる。また、cobaltパッケージのbal.plot()関数を合わせると、重みによる調整前後の傾向スコア分布も可視化できる。対照群の分布が調整されたことがわかる。
weight_att <- weightit(
formula_ps,
data = df,
method = "ps",
estimand = "ATT"
)
bal.plot(
x = weight_att,
var.name = "prop.score",
which = "both",
type = "histogram",
mirror = TRUE,
sample.names = c("pre", "post")
) +
scale_fill_discrete(name = "処置") +
labs(x = "傾向スコア", title = "傾向スコア分布(ATT)")
各共変量がバランシングしたかどうかは、cobaltパッケージのlove.plot()関数を利用する。ageの以外の共変量は重みよって処置群と対照群のバランスが改善し、標準化平均差(standardized
mean
difference)の絶対値が0.1未満に収まっていることがわかる。基準としては、\(|d|<0.1\)や\(|d|<0.25\)が提案されている。
標準化平均差は下記の通り計算される。
\[ d = \frac{\bar X_{Z=1} - \bar X_{Z=0}}{\sqrt{\frac{Var(\bar X_{Z=1}) + Var(\bar X_{Z=0})}{2}}} \]
love.plot(weight_att,
threshold = 0.1,
abs = TRUE,
grid = TRUE,
sample.names = c("pre", "post"), title = "バランシング(ATT)") +
labs(x = "標準化平均差の絶対値")
バランシングしていなければ、因果効果は計算しても交絡の可能性が疑われるので、傾向スコアのモデルを改善するなどを行う必要がある。ここではパッケージの利用方法をまとめているので、このまま進める。この重みを利用して、退院時ADLスコアadl_yのATTを計算する。
# E[Y(1)|Z=1]
y1_att <- df_weight %>% filter(treat == 1) %>% pull(adl_y) %>% mean()
# E[Y(0)|Z=0]を重みをつけてE[Y(0)|Z=1]とする
y0_att <- df_weight %>% filter(treat == 0) %>% with(weighted.mean(adl_y, w = weight_att))
list(
y1_att = y1_att,
y0_att = y0_att,
ATE = y1_att - y0_att
)## $y1_att
## [1] 38.37011
##
## $y0_att
## [1] 35.6895
##
## $ATE
## [1] 2.680611この結果は、lm()関数のweight = weight_attとした重みづけ回帰分析でも同じ結果が得られる。
tidy(lm(adl_y ~ treat, data = df_weight, weight = weight_att))## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 35.7 0.0750 476. 0
## 2 treat 2.68 0.109 24.7 9.58e-132ロバスト分散を使用した効果を推定するのであれば、lmtestとsandwithパッケージを利用する。処置効果ATTとして、退院時ADLスコアは2.68[1.39-3.96]ポイント改善していることがわかる。
fit_gaussian_att <- glm(
formula = adl_y ~ treat,
family = gaussian(link = "identity"),
data = df_weight,
weights = weight_att
)
robust_fit_gaussian_att <- coeftest(fit_gaussian_att, vcov. = sandwich)
res_robust_fit_gaussian_att <- c(
robust_fit_gaussian_att[2],
robust_fit_gaussian_att[2] - 1.96*robust_fit_gaussian_att[4],
robust_fit_gaussian_att[2] + 1.96*robust_fit_gaussian_att[4],
robust_fit_gaussian_att[8]
)
names(res_robust_fit_gaussian_att) <- c("treat", "lower95CI", "upper95CI", "pvalue")
# マッチングした個体のペアをクラスターとして反映させたクラスターに頑丈な標準誤差
# coeftest(fit_gaussian_att, vcov. = vcovCL, weight = ~weight_att)
# coefci( fit_gaussian_att, vcov. = vcovCL, weight = ~weight_att, leverl = 0.95)
res_robust_fit_gaussian_att## treat lower95CI upper95CI pvalue
## 2.680611e+00 1.392828e+00 3.968393e+00 4.505909e-05退院時後遺症sequela_yが改善したかどうかを確認するためには、ロジスティック回帰分析を行えば良い。処置効果ATTとして、退院時後遺症のオッズ比は0.63[0.48-0.84]と改善していることがわかる。
fit_logit_att <- glm(
formula = sequela_y ~ treat,
family = binomial(link = "logit"),
data = df_weight,
weights = weight_att
)
robust_fit_logit_att <- coeftest(fit_logit_att, vcov. = sandwich)
res_robust_fit_logit_att <- c(
exp(c(
robust_fit_logit_att[2],
robust_fit_logit_att[2] - 1.96*robust_fit_logit_att[4],
robust_fit_logit_att[2] + 1.96*robust_fit_logit_att[4])),
robust_fit_logit_att[8]
)
names(res_robust_fit_logit_att) <- c("treat", "lower95CI", "upper95CI", "pvalue")
res_robust_fit_logit_att## treat lower95CI upper95CI pvalue
## 0.637919702 0.482283159 0.843781374 0.001630531平均処置効果ATEの推定
ATEも同じように計算しておく。処置群にあわせて対照群が調整されていることがわかる。
weight_ate <- weightit(
formula_ps,
data = df,
method = "ps",
estimand = "ATE"
)
bal.plot(
x = weight_ate,
var.name = "prop.score",
which = "both",
type = "histogram",
mirror = TRUE,
sample.names = c("pre", "post")
) +
scale_fill_discrete(name = "処置") +
labs(x = "傾向スコア", title = "傾向スコア分布(ATE)")
各共変量がバランシングしたかどうかも確認しておく。すべての共変量は重みよって処置群と対照群のバランスが改善し、標準化平均差の絶対値が0.1未満に収まっていることがわかる。
love.plot(weight_ate,
threshold = 0.1,
abs = TRUE,
grid = TRUE,
sample.names = c("pre", "post"), title = "バランシング(ATE)") +
labs(x = "標準化平均差の絶対値")
この重みを利用して、退院時ADLスコアadl_yのATEを計算する。ATEはプラスの値になっており、処置を行うことで、退院時ADLスコアが改善されたことがわかる。
# E[Y(1)]
y1_ate <- df_weight %>% filter(treat == 1) %>% with(weighted.mean(adl_y, w = weight_ate))
# E[Y(0)]
y0_ate <- df_weight %>% filter(treat == 0) %>% with(weighted.mean(adl_y, w = weight_ate))
list(
y1_ate = y1_ate,
y1_at0 = y0_ate,
ATE = y1_ate - y0_ate
)## $y1_ate
## [1] 42.26167
##
## $y1_at0
## [1] 40.7518
##
## $ATE
## [1] 1.509876この結果は、lm()関数のweight = weight_ateとした重みづけ回帰分析でも同じ結果が得られる。
tidy(lm(adl_y ~ treat, data = df_weight, weight = weight_ate))## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 40.8 0.0744 548. 0
## 2 treat 1.51 0.108 14.0 2.40e-44処置効果ATEとして、退院時ADLスコアは1.5[0.78-2.23]ポイント改善していることがわかる。ロバスト分散を使用した効果を推定する方法はこちら。
fit_gaussian_ate <- glm(
formula = adl_y ~ treat,
family = gaussian(link = "identity"),
data = df_weight,
weights = weight_ate
)
robust_fit_gaussian_ate <- coeftest(fit_gaussian_ate, vcov = sandwich)
res_robust_fit_gaussian_ate <- c(
robust_fit_gaussian_ate[2],
robust_fit_gaussian_ate[2] - 1.96*robust_fit_gaussian_ate[4],
robust_fit_gaussian_ate[2] + 1.96*robust_fit_gaussian_ate[4],
robust_fit_gaussian_ate[8]
)
names(res_robust_fit_gaussian_ate) <- c("treat", "lower95CI", "upper95CI", "pvalue")
res_robust_fit_gaussian_ate## treat lower95CI upper95CI pvalue
## 1.509876e+00 7.815383e-01 2.238213e+00 4.841141e-05退院時後遺症sequela_yが改善したかどうかを確認するためには、ロジスティック回帰分析を行えば良い。処置効果ATEとして、退院時後遺症のオッズ比は0.67[0.55-0.81]と改善していることがわかる。
fit_logit_ate <- glm(
formula = sequela_y ~ treat,
family = binomial(link = "logit"),
data = df_weight,
weights = weight_ate
)
robust_fit_logit_ate <- coeftest(fit_logit_ate, vcov = sandwich)
res_robust_fit_logit_ate <- c(
exp(c(
robust_fit_logit_ate[2],
robust_fit_logit_ate[2] - 1.96*robust_fit_logit_ate[4],
robust_fit_logit_ate[2] + 1.96*robust_fit_logit_ate[4])),
robust_fit_logit_ate[8]
)
names(res_robust_fit_logit_ate) <- c("treat", "lower95CI", "upper95CI", "pvalue")
res_robust_fit_logit_ate## treat lower95CI upper95CI pvalue
## 6.711194e-01 5.521882e-01 8.156663e-01 6.140711e-05傾向スコアとマッチング
ここではMatchItパッケージの利用方法を中心に傾向スコアマッチングの方法をまとめておく。
傾向スコアマッチングの推定対象(estimand)はATTであり、ATEではない。これは処置群の個体について、対照群からマッチングする候補を傾向スコアをもとに選んでくるため、マッチング後の処置群の個体を中心に再構成されるため。層別解析や重み付け法を使えばATEもATTも計算可能になる。
距離とマッチング方法
MatchItパッケージのmatchit関数では、distance引数で傾向スコアを計算する方法を指定できる。デフォルトはglmのロジスティック回帰。
- 一般家線形モデル(
glm) - 一般化加法モデル(
gam) - 決定木(
rpart) - ランダムフォレスト(
randomforest) - ニューラルネットワーク(
nnet) - 共変量バランシング傾向スコア(
cbps) - ベイズ加法回帰木(
bart) - マハラノビス距離(
mahalanobis)
MatchItパッケージのmatchit関数では、method引数で傾向スコアのマッチングする方法を指定できる。デフォルトはnearestの最近隣法マッチング。最近隣法マッチングはわかりやすく、傾向スコアの値が近い個体をマッチングすることで個体の距離を最小化する。個別の個体距離の最小化を目指し、全体での距離最小化は目指さない。
- 最近隣法マッチング(
nearest) - 最適マッチング(
optimal) - 遺伝的マッチング(
genetic) - 厳密マッチング(
exact) - 単純厳密マッチング(
cem) - 層化マッチング(
subclass) - 最適フルマッチング(
full)
処置群に対する平均処置効果ATTの推定
MatchItパッケージのmatchit関数のformula引数にモデル式を渡す。ここでは、replace = TRUEで復元抽出を選択している。match.data関数は、マッチングしたデータフレームを作成してくれる関数。
fit_match <- matchit(
formula = formula_ps,
data = df,
replace = TRUE, # 復元
distance = "glm",
method = "nearest",
estimand = "ATT"
)
fit_match## A matchit object
## - method: 1:1 nearest neighbor matching with replacement
## - distance: Propensity score
## - estimated with logistic regression
## - number of obs.: 15000 (original), 5628 (matched)
## - target estimand: ATT
## - covariates: age, sex, ht, dm, stroke, mibal.plot(
x = fit_match,
var.name = "distance",
which = "both",
type = "histogram",
mirror = TRUE,
sample.names = c("pre", "post")
) +
scale_fill_discrete(name = "処置") +
labs(x = "傾向スコア", title = "傾向スコア分布(ATT)")
すべての共変量はマッチングによって処置群と対照群のバランスが改善し、標準化平均差の絶対値が0.1未満に収まっていることがわかる。
love.plot(fit_match,
threshold = 0.1,
abs = TRUE,
grid = TRUE,
sample.names = c("pre", "post"), title = "バランシング(ATT)") +
labs(x = "標準化平均差の絶対値")
マッチングしたデータを使うのであれば単純に差を計算すれば良い。ここでは、sequela_yはオッズ比ではなく単純な比率。
df_matched <- match.data(fit_match)
df_matched %>%
group_by(treat) %>%
summarise(
mean_adl_y = mean(adl_y),
mean_sequela_y = mean(sequela_y)
) %>%
pivot_longer(.,
cols = c(mean_adl_y, mean_sequela_y),
names_to = 'outcome',
values_to = 'value') %>%
pivot_wider(names_from = treat, names_prefix = 'Z_') %>%
mutate(diff = abs(Z_1 - Z_0))## # A tibble: 2 × 4
## outcome Z_0 Z_1 diff
## <chr> <dbl> <dbl> <dbl>
## 1 mean_adl_y 36.4 38.4 2.00
## 2 mean_sequela_y 0.117 0.0754 0.0421ロバスト分散を使用した効果を推定する。処置効果ATTとして、退院時ADLスコアは1.90[0.64-3.16]ポイント改善していることがわかる。
fit_gaussian_att <- glm(
formula = adl_y ~ treat,
family = gaussian(link = "identity"),
data = df_matched,
weights = weights
)
robust_fit_gaussian_att <- coeftest(fit_gaussian_att, vcov. = sandwich)
res_robust_fit_gaussian_att <- c(
robust_fit_gaussian_att[2],
robust_fit_gaussian_att[2] - 1.96*robust_fit_gaussian_att[4],
robust_fit_gaussian_att[2] + 1.96*robust_fit_gaussian_att[4],
robust_fit_gaussian_att[8]
)
names(res_robust_fit_gaussian_att) <- c("treat", "lower95CI", "upper95CI", "pvalue")
# マッチングした個体のペアをクラスターとして反映させたクラスターに頑丈な標準誤差
# coeftest(fit_gaussian_att, vcov. = vcovCL, weight = ~weight_att)
# coefci( fit_gaussian_att, vcov. = vcovCL, weight = ~weight_att, leverl = 0.95)
res_robust_fit_gaussian_att## treat lower95CI upper95CI pvalue
## 1.903871330 0.643290607 3.164452052 0.003074261退院時後遺症sequela_yが改善したかどうかを確認するためには、ロジスティック回帰分析を行えば良い。処置効果ATTとして、退院時後遺症のオッズ比は0.50[0.28-0.88]と改善していることがわかる。
fit_logit_att <- glm(
formula = sequela_y ~ treat,
family = binomial(link = "logit"),
data = df_matched,
weights = weights
)
robust_fit_logit_att <- coeftest(fit_logit_att, vcov. = sandwich)
res_robust_fit_logit_att <- c(
exp(c(
robust_fit_logit_att[2],
robust_fit_logit_att[2] - 1.96*robust_fit_logit_att[4],
robust_fit_logit_att[2] + 1.96*robust_fit_logit_att[4])),
robust_fit_logit_att[8]
)
names(res_robust_fit_logit_att) <- c("treat", "lower95CI", "upper95CI", "pvalue")
res_robust_fit_logit_att## treat lower95CI upper95CI pvalue
## 0.50430657 0.28695122 0.88630087 0.01733403傾向スコアと層化解析
傾向スコアを使った層化解析の方法をまとめておく。傾向スコアを使った層化解析では、傾向スコアが似た個体を層に分類する。下記のスライドのp33からが分かりやすい。
その層の中では処置以外の共変量は似たような個体が集まるため、各層ごとに平均値を計算して、加重平均を使って平均処置効果を推定する。層ごとに統合してATEを計算する方法は下記が詳しい。
層化解析もMatchItパッケージのmatchit関数のmethod引数にsubclassを渡すことで計算可能。バランシングはsummary()関数でも可能だが、可視化したほうが判断しやすい。各セクションの意味は下記の通り。
Summary of Balance for All Data：マッチング前のバランシングSummary of Balance for Matched Data：マッチング後のバランシングSample Sizes：matched cohortに採用された症例数
各カラムの値は下記の通り。
Means Treated：処置群の平均Means Control：対照群の平均Std. Mean Diff：群間の標準化平均差(SMD)。0に近いほどバランシングされている。Var. Ratio：群間の分散比。連続変数に対してのみ計算される。バランシングされていると1に近く。eCDF Mean：群間でeCDFの平均乖離度合い(経験的累積分布関数 empirical cumulative distribution function: eCDF)eCDF Max：群間でeCDFが最大乖離度合い
n_subclass <- 10
fit_subclass <- matchit(
formula = formula_ps,
data = df,
method = "subclass",
subclass = n_subclass,
estimand = "ATE"
)
summary(fit_subclass)##
## Call:
## matchit(formula = formula_ps, data = df, method = "subclass",
## estimand = "ATE", subclass = n_subclass)
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5707 0.2378 1.4105 1.7385 0.2552
## age 79.3540 72.5853 1.3526 0.8675 0.1471
## sex 0.6013 0.6034 -0.0044 . 0.0022
## ht 0.4717 0.1212 0.8310 . 0.3505
## dm 0.2906 0.0505 0.6739 . 0.2402
## stroke 0.1640 0.0114 0.5603 . 0.1526
## mi 0.2132 0.0309 0.5800 . 0.1823
## eCDF Max
## distance 0.5151
## age 0.5007
## sex 0.0022
## ht 0.3505
## dm 0.2402
## stroke 0.1526
## mi 0.1823
##
## Summary of Balance Across Subclasses
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.3600 0.3523 0.0325 1.0402 0.0093
## age 75.1611 74.9064 0.0509 0.9822 0.0059
## sex 0.5522 0.6119 -0.1219 . 0.0597
## ht 0.2573 0.2444 0.0306 . 0.0129
## dm 0.1411 0.1232 0.0502 . 0.0179
## stroke 0.0694 0.0463 0.0848 . 0.0231
## mi 0.1004 0.0841 0.0521 . 0.0164
## eCDF Max
## distance 0.0229
## age 0.0215
## sex 0.0597
## ht 0.0129
## dm 0.0179
## stroke 0.0231
## mi 0.0164
##
## Sample Sizes:
## Control Treated
## All 9653 5347.
## Matched (ESS) 6300 1434.75
## Matched 9653 5347.
## Unmatched 0 0.
## Discarded 0 0.共変量のバランシングはstrokeのバランシングが少し悪いが今回はこのまま進める。
# plot(fit_subclass, type = "density")
# plot(fit_subclass, type = "jitter")
love.plot(
fit_subclass,
stats = "m", # 平均差
binary = "std", # 標準化
abs = TRUE, # 絶対値
disp.subclass = TRUE,
threshold = 0.1,
grid = TRUE,
title = "共変量のバランス") +
labs(x = "標準化平均差の絶対値")
傾向スコアを使った層化解析からATEを計算するためには少し工夫がいる。ここでは、下記の書籍を参考にしている。
coef_subclass <- NULL
samplesize_subclass <- NULL
robust_var <- NULL
df_subclass <- match.data(fit_subclass)
for (i in 1:n_subclass) {
# 各層のデータを抽出
dataps <- df_subclass[df_subclass$subclass == i, ]
# 各層のデータでモデリング
fit <- lm(adl_y ~ treat + age + sex + ht + dm + stroke + mi, data = dataps)
# 回帰モデルの処置効果の偏回帰係数を取得
coef_subclass[i] <- summary(fit)$coefficients[2,1]
# 各層のデータサイズを取得
samplesize_subclass[i] <- nrow(dataps)
# クラスタに頑丈な標準誤差を計算
robust_var[i] <- coeftest(fit, vcov. = vcovCL, cluster = ~weights)[2,2]
}
list(
coef_subclass = coef_subclass,
samplesize_subclass = samplesize_subclass,
robust_var = robust_var
)## $coef_subclass
## [1] 1.600960 1.854417 1.471620 1.187132 1.485429 1.907515 1.240757 1.381447
## [9] 1.643934 1.912480
##
## $samplesize_subclass
## [1] 1339 1628 1503 1419 1208 1781 1507 1604 1488 1523
##
## $robust_var
## [1] 0.0009859430 0.0201471762 0.0273109711 0.0146781969 0.0093261449
## [6] 0.0093296867 0.0009972078 0.0095593519 0.0199511628 0.0083378377処置効果ATEとして、退院時ADLスコアは1.57[1.51-1.64]ポイント改善していることがわかる。
# 層化する前のサンプルサイズを取得
samplesize <- nrow(df)
# 偏回帰係数を加重平均で計算(p174 12.1 12.2)
tau_hat <- sum((samplesize_subclass/samplesize) * coef_subclass)
# ロバスト分散を加重平均で計算(p174 12.1 12.2)
robust_var_tau <- sum((samplesize_subclass/samplesize)^2 * robust_var)
# ロバスト標準誤差を計算
robust_se_tau <- sqrt(robust_var_tau)
# ここでのモデルの推定パラメタは8。t分布を使って95％信頼区間を計算
tval <- qt(0.975, samplesize - 8)
upper95CI <- tau_hat + tval * robust_se_tau
lower95CI <- tau_hat - tval * robust_se_tau
list(
tau_hat = tau_hat,
lower95CI = lower95CI,
upper95CI = upper95CI,
robust_se_tau = robust_se_tau
)## $tau_hat
## [1] 1.579686
##
## $lower95CI
## [1] 1.51054
##
## $upper95CI
## [1] 1.648832
##
## $robust_se_tau
## [1] 0.03527625