治療の効果を比べる際には、アウトカムに影響する患者特性がそろった集団同士を比較しなければならない。治療割り付けとアウトカムの両方に関連を持った因子は交絡因子となって治療効果の推定を歪めるため、何らかの方法で調整する必要がある。
傾向スコア（propensity score, PS）とは「観察された因子のもとで、ある治療を受ける確率」のことで、患者背景のバランスをとって治療群間の比較可能性を高めるために用いられる。
傾向スコアがどのように治療効果の適切な推定につながるかを理解するためには、潜在アウトカム（potential outcome）といった因果推論のフレームワークを勉強する必要があるが、平たく言えば、
「（もし未測定交絡がなければ）傾向スコアが同じ人では（平均的に見て）治療群間で共変量の偏りがなくなり、擬似的に治療をランダムに割り付けた状態を仮想できる」
ということ（あまり平たくなってないかも）。
傾向スコアを使った解析の流れ
傾向スコアを使った解析の流れは、大きく分けて次の3つのステップから成る。
- 傾向スコアを推定する
- 推定された傾向スコアを評価する
- 傾向スコアを使って治療効果を推定する
この記事で扱うのはステップ1〜2の一部だけ。
傾向スコアの推定方法
まず第1のステップとして、患者背景因子など、治療が割り付けられる前の情報を使って、その患者が「ある治療を受ける確率」を推定する。最も一般的に用いられる方法は、ロジスティック回帰モデルを使った方法である。
ロジスティック回帰モデルでは、共変量Xが与えられた下である治療Z=1へ割り付けられる確率 は次の式で表される（詳細は割愛）。
ここではロジスティック回帰モデルを用いて推定する方法を紹介するが、決定木をベースにした方法で推定することも可能。
傾向スコアの推定に用いる変数
傾向スコアを推定するためのモデル（PSモデル）の従属変数は治療割り付け変数だが、説明変数にはどのような変数を含めるべきか。
まず、全ての交絡因子はPSモデルの説明変数に含まれていなければならない（傾向スコアを使っても未測定交絡に対処できない）。
次に、強い予後予測因子も含めた方がよい（検出力が上がるから）。
逆に、治療の割り付けにしか関係しない因子を含めると、検出力を下げてしまうので含めてはいけない。
治療の割り付けの下流にある因子（中間因子）を含めてはいけない。
推定したPSがイマイチだったら、交互作用項（）や高次項（
など）を加えて推定してもよい。PSモデルでは治療とアウトカムの関連を見ている訳ではないので、何回繰り返してもP-fishing（=有意差が出るまで解析を繰り返すこと）には当たらない。
パッケージとデータの準備
PSAgraphicsパッケージのlindnerデータセットを使う。
lindnerデータは、Lindnerセンターで経皮的冠動脈インターベンション（percutaneous cardiac intervention, PCI）を受けた996例において、abciximabという抗血小板薬を追加したときの効果をみた研究のデータ。*1
下のコードでPCI後6か月時点の予後をdeathとして追加した。
library(tidyverse) library(PSAgraphics) data("lindner") d <- lindner %>% mutate(death = if_else(lifepres==0,1,0))
> head(d) lifepres cardbill abcix stent height female diabetic acutemi ejecfrac ves1proc death 1 0.0 14301 1 0 163 1 1 0 56 1 1 2 11.6 3563 1 0 168 0 0 0 56 1 0 3 11.6 4694 1 0 188 0 0 0 50 1 0 4 11.6 7366 1 0 175 0 1 0 50 1 0 5 11.6 8247 1 0 168 1 0 0 55 1 0 6 11.6 8319 1 0 178 0 0 0 50 1 0
このうち、使用する変数は以下のとおり。
death：PCI後6か月時点の予後（0=生存, 1=死亡）abcix：abciximab使用有無（0=PCI単独, 1=PCI+abciximab）stent：冠動脈ステント留置（0=なし, 1=あり）height：身長（cm）female：性別（0=男性, 1=女性）diabetic：糖尿病有無（0=なし, 1=あり）acutemi：7日前以内の急性心筋梗塞有無（0=なし, 1=あり）ejecfrac：左室収縮率（%）ves1proc：病変血管数
元データで群間のバランスを見てみる
元のデータで治療群間（abcximabなし vs. あり）の患者背景を比較してみる。
library(tableone) table_crude <- CreateTableOne(vars = c("stent", "height", "female", "diabetic", "acutemi", "ejecfrac", "ves1proc","death"), factorVars = c("stent", "female", "diabetic", "acutemi","death"), strata= "abcix", data = d, test = FALSE, smd=TRUE) print(table_crude, smd=TRUE)
Stratified by abcix
0 1 SMD
n 298 698
stent = 1 (%) 174 (58.4) 492 (70.5) 0.255
height (mean (SD)) 171.45 (10.59) 171.44 (10.69) <0.001
female = 1 (%) 115 (38.6) 231 (33.1) 0.115
diabetic = 1 (%) 80 (26.8) 143 (20.5) 0.150
acutemi = 1 (%) 18 ( 6.0) 125 (17.9) 0.372
ejecfrac (mean (SD)) 52.29 (10.30) 50.40 (10.42) 0.182
ves1proc (mean (SD)) 1.20 (0.48) 1.46 (0.71) 0.427
death = 1 (%) 15 ( 5.0) 11 ( 1.6) 0.194
ここでSMDとはstandardized mean difference（標準化平均差）のこと。カテゴリー変数については、各群の割合を とすると、
と定義される。連続変数については、各群の平均を 、標準偏差を
とすると、次のようになる。
SMD<0.1であれば、その変数は群間でバランスが取れていると考えてよい。
マッチング用の関数を集めたMatchItパッケージなら、バランスに関する他の指標も一括して表示できて便利そう（別の記事で説明予定 まとめ済み）。
ロジスティック回帰モデルを使って傾向スコアを推定する
glm()を使ってロジスティック回帰モデルに基づいたPSを計算し、元のデータセットに付け加える。
ps_model <- glm(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data = d, family = binomial(link="logit")) d <- d %>% mutate(ps_logit = predict(ps_model, type = "link"), ps = predict(ps_model, type = "response"))
> d %>% select(abcix, death, ps_logit, ps) %>% head() abcix death ps_logit ps 1 1 1 -0.37258319 0.4079170 2 1 0 0.31645572 0.5784602 3 1 0 0.09786545 0.5244469 4 1 0 -0.10918387 0.4727311 5 1 0 -0.02781550 0.4930466 6 1 0 0.25152728 0.5625524
ここで、psは傾向スコアそのもの、ps_logitは傾向スコアをロジット変換したもので に相当する。
推定された傾向スコアを確認する
推定された傾向スコアが妥当なものかを確認する方法はいくつか提案されているが、
- 治療割り付け予測因子としての性能を評価する方法
- 共変量の分布を均衡化しているかを評価する方法
に分けられる。
傾向スコアの分布
まずは、それぞれの治療群における傾向スコアの分布を見てみる。
ggplot(aes(x = ps, fill = factor(abcix)),data=d) + geom_density(alpha=.25) + theme_bw()
当たり前だが、abciximab使用あり群の方が高い。別記事で説明するが、ある程度の重なりが必要。
ROC曲線
予測因子としての識別能をROC曲線で評価する。
library(pROC) roc_ps <- roc(abcix ~ ps_logit, data=d, ci=TRUE) ggroc(roc_ps, size=1, legacy.axes = TRUE) + geom_abline(color="dark grey", size=0.5)
ROC曲線下面積（area under the ROC, AUROC）は下のコードで計算できる。
> roc_ps Call: roc.formula(formula = abcix ~ ps_logit, data = d, ci = TRUE) Data: ps_logit in 298 controls (abcix 0) < 698 cases (abcix 1). Area under the curve: 0.6807 95% CI: 0.6457-0.7156 (DeLong)
この例ではAUCは0.68で、予測能としてはそれほど高くない。「AUR>0.8が必要」と書かれていることもあるが、AUCは高すぎてもダメ。AUCが高い場合、同じくらいの傾向スコアはどちらか一方の治療に割り付けられやすくなり、比べる相手がいなくなってしまうので、0.6-0.8くらいが丁度いいみたい。
※ ROC曲線を描く方法については以前まとめました（2022-07-17 追加）
Hosmer-Lemeshow検定
Hosmer-Lemeshow検定は、予測モデルの当てはまり具合を評価する一般的な方法。検定が棄却されなければ、「予測値と観測値は違うとは言えない」から当てはまっていると言っていいんじゃない？という論理。傾向スコアでは、あんまりチェックされてないような気がする。
ResourceSelectionパッケージのhoslem.test( )が簡単そうだった。
library(ResourceSelection) hoslem.test(d$abcix,fitted(ps_model),g=10)
Hosmer and Lemeshow goodness of fit (GOF) test data: d$abcix, fitted(ps_model) X-squared = 8.7059, df = 8, p-value = 0.3677
共変量のバランスを確認する
傾向スコアの役割は、治療の割り付けを予測することではなく。治療群間で共変量のバランスを取ることである。傾向スコアを導入することで、共変量がどれくらい均衡化されたかを評価することが重要。
前述のSMDのほか、分散比、経験的累積密度関数（empirical cumulative density function, eCDF）の差などで評価する。量が多そうなので別の機会に（MatchItパッケージで一括表示可能）。
※ MatchItパッケージについての記事を書きました（2022-07-17 追加）
おわりに
- 解析のやり方ばっかり頭に残ると、解析の前提条件とか限界を忘れがちなので注意。
- ロジスティック回帰以外の方法での傾向スコア推定方法についてはそのうちまとめたい。
参考資料
潜在アウトカムなどの因果推論の解説を読みたい方はこちらがオススメ。 www.krsk-phs.com
傾向スコアを使った解析の詳細な概要（どっちも）。 www.krsk-phs.com
星野先生、岡田先生の総説。傾向スコア解析の全体像がつかみやすいのに、細かいところまで解説されています。 ci.nii.ac.jp
*1:Kereiakes DJ, Obenchain RL, Barber BL, Smith A, McDonald M, Broderick TM, Runyon JP, Shimshak TM, Schneider JF, Hattemer CR, Roth EM, Whang DD, Cocks D, Abbottsmith CW. Abciximab provides cost-effective survival advantage in high-volume interventional practice. Am Heart J. 2000 Oct;140(4):603-10.
