ここでは、ベースラインのHbA1cが高い患者ほど強化療法を受けやすく、その層では自然経過でも検査値が下がりやすい、という状況を置いたランダムに生成したデータセットを使います。傾向スコアを推定するところまでは多くの観察研究で共通しますが、その先でマッチングと重み付けのどちらを選ぶかによって、論文に載る数値そのものが変わってしまいます。
この選択について、両方を当てはめて近い方を採用する進め方を勧めないのは、二つの方法が別の集団の平均効果を推定しているためです。以下では、重み付けが1人を何人分として数え直すことなのかを定義から追い、マッチングとの違いを4,000例で比べてみます。
それにより、ATE重み、ATT重み、重複重みのどれを使うかという選択が、解析手法の好みではなく対象集団の指定だという点に気づかれると思います。
重み付けとマッチングは前提を共有します
傾向スコアは、共変量から見て治療を受ける確率として定義されます。
$$e(\mathbf{x}) = \Pr(T=1 \mid \mathbf{X}=\mathbf{x}) \tag{1}$$
式(1)は、共変量の値が$\mathbf{x}$の患者が治療群に入る確率を表しています。$T_i$は患者$i$が治療を受けたかどうかの指標で、$Y_i$はアウトカム、$\mathbf{x}_i$は調整に使う共変量のベクトルです。
この定義についてRosenbaum PR, Rubin DB（1983）が示したのは、この1つの確率で条件づければ共変量の分布が治療群と対照群で揃うという性質で、多数の共変量を1次元へ要約したまま調整に使える理由もここにあります。
ただし、この性質が成り立つには、治療の割り付けが観測した共変量だけで説明できること（条件付き交換可能性）と、どの共変量の値でも両方の治療を受ける確率が0と1の間にあること（正値性）が必要です。マッチングも重み付けもこの前提を共有し、どちらでも未測定交絡は残りますが、マッチングで誰が解析から消えるかは傾向スコアマッチング：誰が解析から消えたのかで詳しく扱っています。
もっとも、違いが出るのは前提ではなく、前提を満たしたあとで何を推定するかの方です。
1人を何人分として数え直すかが重みです
重み付けでは、各患者に自分が受けた治療を受ける確率の逆数を重みとして与えます。
$$w_i^{\mathrm{ATE}} = \frac{T_i}{e(\mathbf{x}_i)} + \frac{1-T_i}{1-e(\mathbf{x}_i)} \tag{2}$$
式(2)は、治療群の患者には傾向スコアの逆数を、対照群の患者には1から傾向スコアを引いた値の逆数を割り当てるという意味です。$w_i^{\mathrm{ATE}}$は患者$i$に与える重みを表します。
実際、数を入れると効き方が見えてきて、傾向スコアが0.1の患者が治療を受けていれば、治療を受けにくい層を代表する1例として10人分に数え直されます。今回のデータセットでATE重みの最大値は28.1で、1人が28人分として数えられていました。
この重みを使った平均の差は、次の形で計算します。
$$\hat{\Delta} = \frac{\sum_i w_i T_i Y_i}{\sum_i w_i T_i} – \frac{\sum_i w_i (1-T_i) Y_i}{\sum_i w_i (1-T_i)} \tag{3}$$
式(3)は、治療群と対照群のそれぞれで重み付き平均を計算し、その差を取る操作です。分母を重みの合計にしたこの形はHajek型と呼ばれ、治療の指標を説明変数、アウトカムを目的変数として各患者に重みを与えた重み付き単回帰の傾きと一致します。
そのため得られるのは、全員が治療を受けた場合と全員が受けなかった場合の平均アウトカムの差、つまり集団全体の平均治療効果（ATE）です。推定対象の違いが結果をどう動かすかは、生存時間解析における因果推論：回帰調整とIPTWは違う量を推定するでも取り上げました。
重みを変えると対象集団が変わります
治療群に合わせた効果を知りたい場合は、治療群の重みを1に固定し、対照群だけを治療群の構成へ寄せます。
$$w_i^{\mathrm{ATT}} = T_i + (1-T_i)\frac{e(\mathbf{x}_i)}{1-e(\mathbf{x}_i)} \tag{4}$$
式(4)は、治療群をそのまま残し、対照群の患者を治療を受ける確率と受けない確率の比、つまりオッズに応じて数え直すという意味です。これで推定されるのは、実際に治療を受けた人たちにおける平均治療効果（ATT）になります。
もう一点、別の選び方としてLi F, Morgan KL, Zaslavsky AM（2018）が提案した重複重みがあります。
$$w_i^{\mathrm{ATO}} = T_i\{1-e(\mathbf{x}_i)\} + (1-T_i)e(\mathbf{x}_i) \tag{5}$$
式(5)は、どちらの治療にもなり得た患者、つまり傾向スコアが0.5付近の患者を重く、どちらか一方にほぼ決まっていた患者を軽く数える重みです。ロジスティック回帰で傾向スコアを推定した場合、モデルに入れた共変量の重み付き平均が2群で厳密に一致するという性質を持ちます。
重みの選び方の違いは真値の側にも表れており、今回のデータセットでは真のATEが-0.34%ポイント、真のATTが-0.62%ポイントでした。ベースラインHbA1cが高い患者ほど改善が大きくなる効果修飾を置いてあり、しかも治療群にはそうした患者が集まっているため、同じデータから計算しても両者は倍近く離れます。
この差はどちらかが偏っているということではなく、見ている集団が別だということです。
マッチングで外れた549例が推定対象を動かします
ロジット尺度で0.2標準偏差のキャリパーを置き、復元なしの貪欲法で1対1最近傍マッチングを当てはめると、1,657例の治療群のうち1,108例だけが相手を見つけました。処理順に依存しないよう治療群の並びは無作為にしてありますが、それでも残る549例は解析から外れ、対照群2,343例のうち1,235例も使われないまま残ります。
実際、図1の両端が二つの方法で扱いの最も分かれる領域です。重み付けはこの少数例に大きな重みを与えて残し、キャリパーつきマッチングは相手が見つからないものとして落とします。
そこで問題になるのが推定対象の記述で、マッチングはATTを推定する方法として紹介されがちですが、キャリパーで治療群の3割が落ちた時点で、対象は治療を受けた人全体から相手が見つかった治療群へ移っています。
実際、今回のマッチング推定値は-0.42%ポイントにとどまり、真のATT（-0.62%ポイント）から0.20%ポイント離れました。解析に入るのは2,216例で、4,000例全員を残す重み付けとは使う人数も変わります。
なお、この現象をマッチングの偏りと読むのは原因の取り違えです。脱落が減るほど推定対象は治療群全体へ戻りますが、キャリパーを緩めれば相手との傾向スコアの隔たりは広がるため、均衡と推定対象は逆方向に動きます。
均衡が揃っても推定値は揃いません
共変量の均衡は標準化差（SMD）で確認するのが一般的で、Austin PC（2011）は絶対値0.1を目安として置いています。調整前のこのデータセットでは、ベースラインHbA1cのSMDが1.08と大きく離れていました。
この調整後の最大絶対SMDは、ATE重みで0.042、マッチングで0.019、重複重みでは0になります。一方でATT重みだけはベースラインHbA1cで-0.118と目安を外れており、この方法については均衡が取り切れていないと読むべきです。
この点は実務で迷いやすいところです。
もっとも、均衡が目安を満たしたATE重み、重複重み、マッチングの間でも、推定値は-0.35から-0.42まで割れます。均衡は調整がうまくいったかを測る指標であって、どの集団の効果かを教える指標ではないため、査読で重みの選択を問われた時にSMDの表を示すだけでは答えになりません。
有効標本サイズが重み付けの代償を映します
重み付けで全例を残したとしても、情報量まで4,000例分あるわけではありません。重みのばらつきが大きいほど、実質的に使える人数は減っていきます。
$$\mathrm{ESS} = \frac{(\sum_i w_i)^2}{\sum_i w_i^2} \tag{6}$$
式(6)は、重みの合計の2乗を重みの2乗和で割った値で、Kishの有効標本サイズと呼ばれます。全員の重みが等しければ実際の人数と一致し、一部に重みが集中するほど小さくなります。
この指標を当てはめると、ATE重みで2,260.7、ATT重みで1,499.6、重複重みで2,900.5となりました。ATE重みが4,000例を使いながら実質2,261例分にとどまるのは上位1%の患者が重みの合計の7.2%を占めるためで、ATT重みがさらに小さいのは対照群の一部だけが治療群の構成を代表しているためです。
逆に重複重みが最も大きいのは、重みが0から1に収まり、極端な値が生じない構造のためです。ただし、この重みが答えるのは重なりの良い部分集団の効果で、研究計画で定めた対象集団と一致するとは限らないため、有効標本サイズの大きさだけを理由に重複重みを選ぶことは適切ではありません。
報告では重みの定式化と分散推定が問われます
重み付けを使った研究でまず求められるのは、どの重みを使ったかの明記です。Austin PC, Stuart EA（2015）は、重みの定式化、重み付き標準化差による均衡確認、重みの分布の要約、分散推定の方法を報告事項として整理しています。
そのなかで見落とされやすいのが分散の扱いで、重み付けを無視した通常の標準誤差は区間を狭く見せます。重みを既知として扱う頑健分散はその問題を避けますが、傾向スコア自体の推定誤差は織り込まずATEに対して保守的になりやすいため、今回の区間は傾向スコアの当てはめごと500回再標本化して求めました。
さらに、極端な重みにはCole SR, Hernán MA（2008）の安定化重みがあり、重みの平均が1付近に収まるかを診断として使えます。もっとも、安定化は分散を抑える工夫であって、正値性が崩れた領域の情報不足を補うものではありません。
実務では、私なら重みを決める前に図1のような分布を見て、重なりが広ければATE重みで計画通りの対象集団を推定し、片方の群しかいない領域が目立つ時は、Desai RJ, Franklin JM（2019）の手引きにあるように対象集団を治療群へ限定するか、重みの上限を決めたうえでその判断を論文に書きます。時間とともに治療が変わる設定では、周辺構造モデル：検査値で調整すると治療効果が半分になるで扱った時間依存の重みにつながります。
Rで重みと推定値を確かめます
以下では、重み付けの中心部分にあたるデータの設定、傾向スコア、3種類の重み、Hajek型の推定値を追います。base Rとstatsだけを使い、追加パッケージはありません。
set.seed(42)
n <- 4000
age <- rnorm(n, mean = 62, sd = 9)
hba1c0 <- rnorm(n, mean = 8.4, sd = 1.1)
duration <- pmax(rnorm(n, mean = 9, sd = 4), 0.5)
ps_true <- 1 / (1 + exp(-(-0.40 + 0.045 * (age - 62) + 1.15 * (hba1c0 - 8.4) + 0.05 * (duration - 9))))
treat <- rbinom(n, size = 1, prob = ps_true)
tau <- -0.35 - 0.45 * (hba1c0 - 8.4) # 個人ごとの真の効果(効果修飾)
change <- 0.10 - 0.30 * (hba1c0 - 8.4) + 0.010 * (age - 62) +
tau * treat + rnorm(n, mean = 0, sd = 0.60)ここで用意したのは、年齢、ベースラインHbA1c、罹病期間を共変量とする4,000例で、強化療法が1,657例、標準療法が2,343例です。アウトカムは24週時点のHbA1c変化量（%ポイント、負が改善）とし、ベースラインHbA1cが1%ポイント高いほど真の効果が0.45%ポイント大きい改善になる効果修飾を置きました。
fit_ps <- glm(treat ~ age + hba1c0 + duration, family = binomial(),
control = glm.control(epsilon = 1e-12, maxit = 100))
ps <- fitted(fit_ps)
w_ate <- ifelse(treat == 1, 1 / ps, 1 / (1 - ps))
w_att <- ifelse(treat == 1, 1, ps / (1 - ps))
w_ato <- ifelse(treat == 1, 1 - ps, ps)
hajek <- function(w) {
sum(w[treat == 1] * change[treat == 1]) / sum(w[treat == 1]) -
sum(w[treat == 0] * change[treat == 0]) / sum(w[treat == 0])
}
round(c(unadjusted = mean(change[treat == 1]) - mean(change[treat == 0]),
ate = hajek(w_ate), att = hajek(w_att), ato = hajek(w_ato)), 3)#> unadjusted ate att ato
#> -0.914 -0.352 -0.584 -0.405出力の読み方はこうなります。未調整の差-0.914は真のATE（-0.342）から0.57%ポイント離れており、交絡の分だけ改善を大きく見せています。
一方、重み付けの3行はそれぞれの真値に近い位置へ落ち着いており、ATE重みの-0.352は真のATEと、ATT重みの-0.584は真のATT（-0.621）と、重複重みの-0.405は重なり集団での真値（-0.407）と対応しています。
読み違えやすいのは、-0.352と-0.584の差を精度の問題として受け取ってしまうことです。両者は別の集団の平均を推定しているため、どちらの集団について語りたい時かを先に決めることになります。
ess <- function(w) sum(w)^2 / sum(w^2)
round(c(ate = ess(w_ate), att = ess(w_att), ato = ess(w_ato)), 1)#> ate att ato
#> 2260.7 1499.6 2900.5重みによってこれだけ開くのが、同じ4,000例から得られる情報量の違いで、信頼区間の幅も重複重みが最も狭く、ATT重みが最も広いという同じ順序になりました。
png("fig1-ps-overlap.png", width = 1200, height = 800, res = 150)
breaks_ps <- seq(0, 1, by = 0.05)
hist_treated <- hist(ps[treat == 1], breaks = breaks_ps, plot = FALSE)
hist_control <- hist(ps[treat == 0], breaks = breaks_ps, plot = FALSE)
plot(hist_control, col = rgb(0.2, 0.4, 0.8, 0.5), border = "white",
ylim = c(0, max(hist_treated$counts, hist_control$counts)),
main = "Propensity score distribution by group",
xlab = "Estimated propensity score", ylab = "Number of patients")
plot(hist_treated, col = rgb(0.9, 0.4, 0.1, 0.5), border = "white", add = TRUE)
dev.off()なお、マッチングと信頼区間のコードは分量の都合で省きますが、算出条件はロジット尺度で0.2標準偏差のキャリパーを置いた復元なしの1対1最近傍、区間はブートストラップ500回、SMDはAustin PC（2011）の定義です。5つの推定値を並べたのが図2になります。
また、図2で目を引くのは、マッチングの点が真のATTではなく真のATEの近くに来ていることです。マッチングがATTを推定するという説明を、キャリパーで対象を絞った後もそのまま当てはめると、推定対象の記述を誤ります。
選択は対象集団から決まります
重み付けとマッチングのどちらを選ぶかは、推定したい集団を先に決めれば自動的に絞られます。集団全体について語りたい時はATE重み、実際に治療を受けた患者について語りたい時はATT重みかマッチングで、後者では誰が外れたかまで報告します。
次に見るのが傾向スコアの重なりで、片方の群しかいない領域が広ければATE重みは少数例に頼った推定になり、重複重みへ移るなら推定対象が重なりの良い部分集団へ変わることを本文と要旨の両方に書くのが条件です。
ここまで読んで、自分の研究ではどの重みを選べばよいのか、と思われた方もいるかもしれませんが、実際のデータでは共変量の重なり、外れる患者の特徴、報告先の雑誌が求める記載事項まで含めて判断することになります。Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
- Rosenbaum PR, Rubin DB. The central role of the propensity score in observational studies for causal effects. Biometrika. 1983;70(1):41-55.
- Austin PC. An introduction to propensity score methods for reducing the effects of confounding in observational studies. Multivariate Behavioral Research. 2011;46(3):399-424.
- Austin PC, Stuart EA. Moving towards best practice when using inverse probability of treatment weighting (IPTW) using the propensity score to estimate causal treatment effects in observational studies. Statistics in Medicine. 2015;34(28):3661-3679.
- Cole SR, Hernán MA. Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology. 2008;168(6):656-664.
- Li F, Morgan KL, Zaslavsky AM. Balancing covariates via propensity score weighting. Journal of the American Statistical Association. 2018;113(521):390-400.
- Desai RJ, Franklin JM. Alternative approaches for confounding adjustment in observational studies using weighting based on the propensity score: a primer for practitioners. BMJ. 2019;367:l5657.
