Propensity score

MikuHatsune2014-05-11

Propensity Scoreというものを知った。近年Propensity scoreを使った報告がめっちゃ増えている。

医学研究において最強の比較試験と言われているのが、RCT(Randomize Control Trial)と呼ばれる研究である。
これは、研究計画のプロトコールを満足する患者が登録されたときに、患者がどんな背景を持っていてもそれは考慮せずに、ランダムに新薬や新しい治療法で治療が行われる群と、プラセボや既存の治療法で治療される群に割りつける。そうすることで、介入する/しない 以外の条件は、2群でおよそ均等に分布する(はず)ので、治療によって予後が良くなった/ならない が、介入によるものだということがかなり直接的に推定できる。
 
しかしながら、治療をする/しない が完全にランダムに割りつけられるかというと、いつもできるとは限らない。というのも、実際に病気の人に「治療をしない」ということはかなりチャレンジングで倫理的にアレである。なので実際のところ、プロトコールに合致する患者がいても、現場の医師の判断で当時の医療水準にあった治療が行われることが多い。これは観察研究に多くあるパターンで、こうなると集団に偏りができて、治療した/しない で予後が良くなっても、「治療される人たちはそもそも予後が悪い人に治療がされるのだから、治療して予後が悪いのはそもそも予後が悪い人のせい」といちゃもんが付けれられる。
観察研究に見られる、集団内の偏りをなんとか是正しようとするのが、Propensity scoreらしい。
ランダム化臨床試験でなくとも N Engl J Med に掲載される
Propensity Score | 大阪大学腎臓内科
Propensity Score 図解まとめ “MIT白熱教室 これからの因果推論を考えよう”
概念としては1983年に出たらしいが、臨床研究に出たのは1996年の、右心カテでICU予後がどうなるか調べた論文(JAMA 276: 889-897, 1996)らしい。
これ、実は読んでないので、最近読んだやつを見る(Anesthesiology. 2014 May;120(5):1098-108.)。
心臓の手術で、たいてい心臓がへタるから、inotropeという変力作用のある薬剤、つまり心筋の収縮力を上げて脈拍・血圧を上昇される処置が行われる。
しかし、inotropeの投与には、「心臓をがんばらさなければならないほど、そもそも心臓がヘタっている」という状況が既にあり、「inotropeが投与される患者はそもそも心機能が低い」という偏りがある。なので、「inotrope投与で予後が悪い」とは昔からそれとなくみんな思っていたが、「inotropeで予後が悪いのはそもそも心機能が低いから」と言われ続け、かといってRCTで確かめようにも、「心機能が悪いのにinotropeを使わないなんてそれなんて無理ゲー」ということでそんなにできなかったのが実情。
 
今回は6005人の観察研究で、2097人(35%)にinotropeが投与され、3098人(65%)はinotropeを投与されなかった。
ここからPropensity score(PS)の算出にとりかかるわけだが、要はPSは[0 - 1]のロジスティック回帰の結果である。inotrope投与をする/しない の判断に関わりそうな臨床的な項目30弱を選んで、これらを回帰して「この患者にinotropeが投与される確率」を計算する。ロジスティック回帰は興味のある事象、例えば予後で死んだ/生きた Y=\{\begin{matrix}1&death\\0&alive\end{matrix}を、自分がデータとして持っている検査値などの各種パラメータX_iで回帰して、各パラメータの重要度である係数\beta_{i}を決めにかかる作業
Y=\beta_0+\beta_{1}X_{1}+\beta_{2}X_{2}+\dots+\beta_{n}X_{n}
で、ロジット変換とか使えば、[0 - 1]のロジスティック曲線に乗って、確率となる。これがPSである。
すると、治療をした群と治療をしていない群で、それぞれの患者のPSの分布が取れることになる。
PSを計算しただけの集団の概念図。Inotropeが投与される人(右側)の分布は、Inotropeが投与されない人(左側)の分布に比べて、PS(縦軸)が大きい、つまり治療されやすい方に偏っている。

# 適当なPS分布の作成
s1 <- seq(0, 10, length=100)
ino1 <- dnorm(s1, 6, 1.2) # treat
ino2 <- dnorm(s1, 4, 1.4) # w/o treat

xl <- c(-max(ino2), max(ino1))
yl <- c(0, 1)
par(mar=c(5, 4.5, 2, 2))
plot(0, type="n", xlim=xl, ylim=yl, xlab="", ylab="Propensity score (PS)", cex.lab=1.6, cex.axis=1.6, las=1)
abline(v=0, lty=2)
y <- seq(0, 1, length=length(s1))
points(ino1, y, type="l", col="red", lwd=2)
points(-ino2, y, type="l", col="blue", lwd=2)
text(mean(c(0, par()$usr[1])), par()$usr[3], "No Inotrope", xpd=TRUE, adj=c(NA, 4), cex=1.5)
text(mean(c(0, par()$usr[2])), par()$usr[3], "Inotrope", xpd=TRUE, adj=c(NA, 4), cex=1.5)


ここからPS マッチングを行うが、PSが等しい人を治療群と非治療群から一人ずつ選びだす。すると、非常に事後解析な見方だが、「PS(=治療を行う確率)が等しい患者を、治療した場合と治療しなかった場合」という、いわば「たられば」的状況を作り出すことが可能となる。自分的解釈では、RCTでは自然に割りつけたら最終的に患者背景は2群で均等に分布するが、PSではほぼ同一の背景を持つ、対応のある2人を均等に分けていく作業だと思う。なので、検定をするときにも、対応のある検定としてMcNemar検定や、対応のある生存曲線を考える必要がある。
しかも、「観察研究でありながらRCT化できそう」という印象があるが、ロジスティック回帰に持ち込むときに、パラメータが抜けていたらそのパラメータがバイアスとなりうる。なので、感度分析が必要になるが、例えばこの論文と同じ号で、GNASのバリアント次第ではCABG後にβブロッカーを投与すると予後が悪くなる(Anesthesiology. 2014 May;120(5):1109-17.)というような、強力な遺伝素因があったりするとバイアスになる。ものすごいテーラーメイド医療っぽい。
また、PSのマッチング処理としては、今回は似たような人を探してマッチング、つまり、シミュレーションなのでPSがほぼ同一の人をマッチングさせたが、実際には数千人規模でしかも2群で差があれば数が減ってしまう。なので、5等分くらいに分割したカテゴリーでマッチングさせているやつが多いし、また「PSが0.95なのに治療しない」「PSが0.01なのに治療しちゃう」という状況は、まあ非常に稀だし、2群の分布の山に重なりがなければそもそもマッチングしなくなる。というかそういう場合はそもそも研究デザインがアレだと思うが。
というわけで、解析人数を減らさないためにも、分布に重み付けしてふくらませるような手法があったりするっぽいけど、今回はマッチングだけ。なので、PS補正した集団では、灰色で示される各々1170人の集団に減ってしまっている。

minPS <- apply(cbind(ino1, ino2), 1, min)
par(mar=c(5, 4.5, 2, 2))
plot(0, type="n", xlim=xl, ylim=yl, xlab="", ylab="Propensity score (PS)", cex.lab=1.6, cex.axis=1.6, las=1)
polygon(c(minPS, rev(-minPS)), c(y, rev(y)), border=NA, col=grey(0.85))
abline(v=0, lty=2)
y <- seq(0, 1, length=length(s1))
points(ino1, y, type="l", col="red", lwd=1.5, lty=3)
points(-ino2, y, type="l", col="blue", lwd=1.5,lty=3)
text(mean(c(0, par()$usr[1])), par()$usr[3], "No Inotrope", xpd=TRUE, adj=c(NA, 4), cex=1.5)
text(mean(c(0, par()$usr[2])), par()$usr[3], "Inotrope", xpd=TRUE, adj=c(NA, 4), cex=1.5)
points(minPS, y, type="l", col="red", lwd=2)
points(-minPS, y, type="l", col="blue", lwd=2)
idx <- 60
segments(-10^10, y[idx], minPS[idx], lwd=2, lty=3)
segments(-minPS[idx], y[idx], minPS[idx], lwd=3)


 
PS マッチングを行うと、2群で偏りがあったものが、なんとか是正されるっぽい。たぶん。
特に、LVEF(左室駆出率。左室内の収縮によってどれくらい血液を前進に供給できるかの指標) < 30% や、CPB(人工心肺を稼動している時間) > 120分 などは心機能にかなり影響する因子なので、2群で偏っていたっぽいが、これらが均一になる。
Anesthesiology. 2014 May;120(5):1098-108.

 
PS マッチングを行った補正集団での解析結果。元の集団では、各種臨床パラメータが回帰分析の結果、ほとんどが統計学的有意な結果になっている。が、これはもともと2群で偏っていたからで当たり前である。一方、補正集団では統計学的有意なパラメータは消えている。しかし、LVEF > 30% のようなきつい偏りがあったパラメータは、補正してもちょっと偏りがあるっぽい気がしないこともない。
Anesthesiology. 2014 May;120(5):1098-108.

 
今回の論文では、候補に挙げたパラメータが十分だったかを感度分析で確かめるようなことはしていないっぽいが、それでもまあ以前から言われていたとおり、inotropeを投与したらやっぱり予後が悪いのではないかと言ってる。しかし、これは統計学的になんやかんやしたp-hacking的側面は否めないし、結局はRCTやろうぜ、みたいなことがググっても出てくるので面白いけど適応は慎重に、といったところ?
Rではいくつかパッケージが出てくるので、余力があればやる。