今年もこの季節がやってきました。
この記事は、R Advent Calendar 2013 21日目の記事です。
(遅れてしまって申し訳ない。。)
年に一度しかブログを書かない末席医療系研究者の@Hiro_macchanことゲス野郎です。
私は、統計学の専門家でない上数学の素養もない、さらに、エンジニアとしてのスキルに長けているわけではないデータサイエンティストの対局に位置する人ですのでまさかりは優しく投げていただけると幸いです。
（まさかり自体はありがたいです。）
今年の記事は、最近業務で必要となった、モデルの感度分析の勉強過程の話をしたいと思います。
医学研究の領域で感度というと、ある検査（判定）に於いて「陽性（正例）と判定されるべきものを正しく陽性（正例）と判定する確率」をがイメージされますが、別物のお話です。
今回扱う、感度分析とは測定できなかった共変量を明示的にモデルに組み込み、 その影響力を変化させて、推計した因果効果がどれくらい変化するかを確認する方法とです（星野：調査観察データの統計科学，P133を参考）。
もう少し具体的に考えます。
観察データに基づいてある治療法（\( Z \)）の病気に対する治療効果（\( Y \)）を検証したいとします。
\( Y \)に影響する因子は\( Z \)の他に\( X \)と\( U \)があり、治療法\( Z \)を受けた群、受けなかった群でその存在が偏っています。
ここで、治療法\( Z \)の\( Y \)に対する効果をモデル化するとすると \[ Y=a+bZ+cX+dU+\eta \] を推計し、bの大きさを見てあげれば良いと言うことになります。 このモデル式を真のモデル式と呼称しましょう。
さて、ここでよくある話ですが、Uが測定できなかったとします。
すると、推計可能なモデル式は以下の様になります。
\[
Y=a+b'Z+c'X+\eta
\]
この推計結果\( b' \)は\( b \)とは異なり、Biasされている事になります。
このモデル式を誤ったモデル式と呼称しましょう。
解析者「この結果を発表しても、各方面から袋だたき。死ぬしかないのか？」
いえ、そんなことはないですね。
いくつか戦略が考えられるでしょうが、今回は感度解析を行いましょう。
実は、この\( b' \)の大きさは以下の様になっています（Angrist:「「ほとんど無害な」計量経済学」,p61より）。 \[ b'=b+d\delta \] ここで\( d \)を真のモデル式でのZの回帰係数,\( δ \)を\( U \)を\( Z \)に回帰した際の回帰係数とします。
よって、\( U \)の情報を外挿しながら論理を組み立てて、Biasの度合いを報告して 上げれば良いと言えます。致命的な場合はどうしようもないですが。。。
さて、次に実際の観察データの解析を考えましょう。
以下は、HIV罹患者に対する住宅支援サービスによってHIV罹患者の予後が改善するかどうかを検証した論文です。
Schwarcz, S. K., Hsu, L. C., Vittinghoff, E., Vu, A., Bamberger, J. D., & Katz, M. H. (2009). Impact of housing on the survival of persons with AIDS. BMC public health, 9, 220. doi:10.1186/1471-2458-9-220
この論文では、Coxの比例ハザードモデルを使ってHIV罹患者に対する住宅支援サービスがHIV罹患者の死亡ハザードに与える効果を調べています。
モデル式は
\[
λ(t|X)=λ0(t)exp(bZ+\beta X)
\]
の線型モデルになっています。
報告されているモデルでのハザード比は0.2（80％のリスク減）、なので、\( b \)の値は\( log(0.2)=-1.609438 \)となります。
Zの存在比率は\( P=70/676 \) なので分散は\( P(1-P)=0.0928 \)
ここで、まず観測できなかったConfounder があるものとしたモデルを考えます。
先ほどの以下の式を思い出しましょう。
\[
b'=b+d\delta
\]
もし、\( b \)の値がが本当は0である、つまり治療法\( Z \)には効果がないとしましょう。
とすると、 \[ b'=d\delta \] となります。 よってこういった状況が生じうる観測できなかった Coufounder の効果の大きさ\( d \)は以下の式で表す事が出来ます。
\[ d=b'/\delta \]
回帰式の係数は以下の様に表せるので \[ \delta = \frac{Cov(Z,U)}{Var(U)} \]
\[ d=\frac{b'\sqrt{Var(Z)}}{r\sqrt{Var(U)}} \] ただし、\( r \)は観測できなかったConfounderと介入\( Z \)の相関係数で \[ r=\frac{Cov(Z,U)}{\sqrt{Var(Z)}\sqrt{Var(U)}} \]
モデルで推計された\( Z \)の係数\( b' \)と\( Z \)の分散は観察データからわかるので、 \( Z \)とConfounder の関連（相関係数）とUの分散を仮定して、\( d \) がどれくらい の幅で存在するのかを見てみましょう。
beta <- log(0.2)
p_2 = 70/676
Var_Z <- p_2 * (1 - p_2)
r <- seq(-0.9, 0.9, length = 50)
p <- seq(0.1, 0.9, length = 50)
f <- function(r, p) {
(beta/r) * (sqrt(Var_Z)/sqrt(p * (1 - p)))
}
d <- outer(r, p, f)
persp(r, p, d, theta = 30, phi = 30, expand = 1, col = rainbow(50), border = NA)
論文中で彼らは、観測出来なかったConfounderの影響を、Confounder と住宅提供との関連が強く（相関係数0.5）、存在率が70%のConfounder で住宅提供を受けなかった症例の死亡ハザードを9倍に上げる効果のあるものと見積もっています。
（下の計算では、8.5倍でした。論文に出てない数字から出しているのか、私の計算が根本的に間違っているのか。。。。どっちだろ？）
1/exp(f(0.5, 0.7))
## [1] 8.5
全症例の7割で観察出来て、そこまでリスクに影響する因子は現場勘として存在しないよねという結論で議論をしめています。
以上が、感度分析の流れになります。
多分、分析屋として使う事はあまりない気がしますが、ドメイン知識に基づいたモデルの検証をするためにはこう言う形のアプローチもあるんじゃないかなと思っています。
ドメインの知識・経験超大事。
あと、あまり今回Rのコード出てきませんでした。すいません。
数式R Markdownで書いたんですが、楽に書けますね。お勧めです。
今年は色々な事がありました。しんどかったこと、うれしかったこと様々です。
来年は、どんな年になるかわかりませんが、取りあえず頑張ります。
来年はF#でもチャレンジしてみるかなぁ。
では、皆さん、メリークリスマス、良いお年を。