[連載]フリーソフトによるデータ解析・マイニング第37回

Rと生存時間分析(2)

1.セミノンパラメトリックモデル

(1)コックス比例ハザードモデル

 イベントに影響を及ぼす複数の因子(共変量)の影響を解析することを前提としたノンパラメトリックモデルをセミノンパラメトリックモデルと呼ぶ。最も広く用いられているモデルはコックス比例ハザードモデル(Cox proportional hazard model)である。
 2つのハザード(瞬間死亡率)関数h1(t)、h2(t)が、すべて可能なt>0に対し、関係式h1(t)=ch2(t)が成り立つとき、2つのハザード関数は比例すると言う。関係式の中のcは時間tと無関係な定数である。
 共変量x=(x1,x2,…,xm)Tを持つハザード関数をh(t|x)とし、xを変数とする関数r(x)と、あるハザード関数h0がすべてのt>0とxについて関係式

関係式


が成り立つとき、この式を比例ハザードモデル、あるいは比例危険度モデルと呼ぶ。
 共変量xは、年齢や血圧のような連続変数、性別や結婚の有無のようなカテゴリカル変数、これらの交差項などを含む変数ベクトルである。
 コックス比例ハザードモデルは、次の式で定義される。

コックス比例ハザードモデル式の定義


 このモデルは、生存時間を目的変数とした回帰モデルである。共変量 がすべて0の場合のハザードをベースラインハザード(baseline hazard)と呼ぶ。コックス比例ハザードモデルは、ベースラインハザードを基準とした分析手法である。
 モデルのパラメータβ=(β1,β2,…,βm) は、次に示す条件付きの確率の部分尤度(partial likelihood)を最大化する方法で推定できる[1][2]。

条件付きの確率の部分尤度


 式の中のj=1,2,…,Mはイベントの数、tjはイベントの時点、R(tj)は各時点のリスクの集合である。

(2)パラメータの推定

 コックス比例ハザードモデルのパラメータの推定方法としては、直接法、Breslowの近似法、Efronの近似法などが提案されている。直接法は、イベントの数が多くなると、部分尤度の分母の項が煩雑になるので、近似法としてBreslowの近似法、Efronの近似法が提案されている。これらの近似法は、直接法と比較して計算が簡単であるが、同時に起こるイベントの数が多くなった場合、推定されるβが0に偏った値を取り、妥当性を失うと言われている。
 パッケージsurvivalには、コックス比例ハザードモデルのパラメータを推定する関数coxphがある。関数 coxphの書き式を次に示す。

coxph(formula, data, method,…)

 引数formulaには、モデルに用いる共変量などを指定する。引数methodには、3種類(efron, breslow, exact)のオプションの中から1つを選択して指定する。デフォルトは"efron"法が指定されている。
 パッケージsurvivalの中に用意されているデータkidneyを用いた使用例を次に示す。データkidneyは、ポータブル透析装置の使用と腎臓患者の生存時間に関して、38ペア(使用と不使用)に対する実験データである。データセットの中の変数を次に示す。

patient  ID
time 時間
status 打ち切りは0、その他は1
age 年齢
sex 男性=1、女性=2
disease 病気の種類(GN,AN,PKD,Other)
frail オリジナル論文からのフレイルティ(frailty)の推定値

 データの中の性別(sex)と病気の種類(disease)を説明変数としたコックス比例ハザードモデルの解析例を次に示す。

>library(survival)
>data(kidney)
>kidney.cox<-coxph( Surv(time, status) ~ sex+disease, data=kidney)
>summary(kidney.cox)
summary(kidney.cox)の結果

 関数coxphは、変数ごとの係数β、exp(β)、係数の標準誤差、z値、p値、信頼区間、仮説H0=0,H1≠0の検定統計量などを返す。
 関数coxphには3種類(尤度比の検定、ワルド検定、スコア検定)の検定が用いられている。これらは、係数が正規分布に従う仮定の下で異なる角度から求めている。
 モデルによる生存時間の当てはめは、次のように関数survfitを用いると便利である。

>kidney.fit<-survfit(kidney.cox)
>summary(kidney.fit)
summary(kidney.fit)の結果
<後略>

 関数survfitで推定された生存曲線および信頼区間は、関数plotを用いて図示することができる。その例を次に示す。

>plot(kidney.fit)


図1 推測されたkidney.fitの生存曲線

図1 推測されたkidney.fitの生存曲線

(3)残差分析

 生存分析のデータは、打ち切りデータが含まれているので、残差分析は一般線形回帰モデルの残差分析ほど単純ではない。そのため、マルチンゲール残差(martingale residuals)、シェーンフィールド残差(Schoenfeld residuals)、スコア残差(score residuals)、デヴィアンス残差(deviance residuals)など幾つかの定義が提案されている。それぞれ特徴を持っているのでどれが優れているかに関する評価はしがたいが、より多く用いられているのはマルチンゲール残差である。
 coxphオブジェクトの残差は、関数residuals.coxph(略してresiduals)で求めることができる。デフォルトには、マルチンゲール残差が指定されている。関数residualsを使用する際には、短縮した文字列residを用いることも可能である。既に求めたcoxphのオブジェクトkidney.coxの残差のプロットを、関数scatter.smoothを用いて作成する例を次に示す。

>scatter.smooth(residuals(kidney.cox))
>abline(h=0,lty=3,col=2)


図2 マルチンゲール残差プロット

図2 マルチンゲール残差プロット

 残差プロットから分かるように、マルチンゲール残差の最大値は1で、最小値は下限がない歪んだ分布である。この歪んだ分布を標準化して扱いやすくするために考案されたのがデヴィアンス残差である。デヴィアンス(尤離度)残差は、関数residualsに引数type="deviance"を用いることで、求めることができる。その例を次に示す。

>scatter.smooth(residuals(kidney.cox, type="deviance"))
> abline(h=0,lty=3,col=2)


図3 デヴィアンス残差プロット

図3 デヴィアンス残差プロット

 シェーンフィールド残差、スコア残差は関数residualsにそれぞれ引数type=" schoenfeld "、type=" score "を指定することで求めることができる。

(4)ハザードの比例性の分析

 コックス比例ハザードモデルの前提は、ハザード比が時間によらず一定であることである。よって、比例ハザードモデルを用いる場合は、この仮定を吟味する必要がある。
 パッケージsurvivalには、比例性を分析する関数cox.zphが用意されている。関数cox.zphは、シェーンフィールド残差を用いて、β(t)=β+θg(t)における仮説H0:θ=0の検定に関する計算結果を返す。式の中のg(t)は滑らかな関数である。また、関数cox.zphはg(ti)と標準化されたシェーンフィールド残差との相関係数をも返す。
 関数cox.zphには引数transformがあり、式β(t)=β+θg(t)の中のg(t)を指定する。デフォルトにはKaplan-Meier推定量1-S(t) が指定されている。用意されているオプションから選択することも、各自が関数を指定することもできる。
 引数transformの関数g(t)をデフォルト値にし、kidney.coxを用いた例を次に示す。

>kidney.zph<- cox.zph(kidney.cox)
>kidney.zph
kidney.zphの結果

 仮説が棄却されると比例ハザードの仮定が充たされていない可能性があることを示唆する。ただし、仮説検定の結果はg(t)に依存することに留意してほしい。
 関数plot.cox.zph (略してplot)を用いると変数ごとのシェーンフィールド残差プロットに、スプライン平滑(関数はsmooth.spline)化された曲線と標準誤差の2倍の信頼区間が示されるグラフを作成することができる。スプライン平滑化の自由度は引数df="n"で指定できる。デフォルトは、n=4になっている。
 オブジェクトkidney.coxを用いた例を次に示し、その結果を図4に示す。図4の各図はシェーンフィールド残差の各成分のプロットで、横軸は時間である。ここでは、スプライン平滑化曲線の傾向を観察しやすくするため自由度df=2にした。

>par(mfrow=c(2,2))
>plot(kidney.zph,df=2)


図4 比例性の診断プロット

図4 比例性の診断プロット

 スプライン平滑化曲線の傾向に、時間に伴う明らかな変化パターンが見られなければ比例ハザードの仮定には問題がないと言われている。個別の変数に比例ハザードの仮定が充たされていない場合は、変数の交差項をモデルに導入するような試みが必要であろう。

(5)交互作用と変数の選択

 一般回帰分析の場合と同じく、コックス比例ハザードモデルの場合でも、説明変数の交互作用効果を取り入れたモデルの構築が可能である。次にその例を示す。

>kidney.cox2<-coxph( Surv(time, status) ~ (sex+age+ disease)^2, data=kidney)
>kidney.cox2
kidney.cox2の結果

 返されたp値から分かるように、年齢(age)はモデルへの寄与度が低い。このようなモデルの作成にあまり役に立たない変数を取り除く(役に立つ変数を選択する)方法として、AIC情報量を用いる関数stepAICを用いることが可能である。関数stepAICは、パッケージMASSに含まれている。関数stepAICの応用例として、求めたオブジェクトkidney.cox2を用いた例とその結果を次に示す。

>library(MASS);stepAIC(kidney.cox2)
library(MASS);stepAIC(kidney.cox2)の結果
<中略>
library(MASS);stepAIC(kidney.cox2)の結果

 このように機械的にAIC情報量に基づいて最適と判断された結果が返される。

2.パラメトリックモデル

 生存時間が確率分布に従う仮定の下でモデルを構築したモデルをパラメトリックモデルと呼ぶ。パラメトリックモデルは、生存時間の確率分布を仮定する制約条件があるので、応用範囲が限定される。パラメトリックモデルは、コックス比例ハザードモデルと比べると計算速度が速いのが特徴であるが、今日のコンピュータでは、もはや大きい問題として取り上げる必要がないであろう。
生存時間モデルに多く使われる確率分布は、指数分布、ワイブル分布、対数正規分布、対数ロジスティック分布などである。これらの分布の確率密度関数、生存関数、ハザード関数を表1に示す。

表1 主な確率分布の生存関数とハザード関数
分布名称 関数
指数(exponential) 指数(exponential)
ワイブル(weibull) ワイブル(weibull)
対数正規(log-normal) 対数正規(log-normal)
対数ロジスティック(log-logistic) 対数ロジスティック(log-logistic)

 指数分布はワイブル分布の ときの特殊なケースである。
 パッケージsurvivalには、パラメトリック生存モデルを当てはめる関数survregがある。その書き式を次に示す。

survreg(formula=formula(data),dist="weibull", ...)

 引数distで確率分布を指定する。デフォルトにはワイブル分布が指定されている。オプションとしては、指数分布(exponential)、ガウス分布(gaussian) 、対数正規分布 (log-normal) 、ロジスティック分布(logistic)、対数ロジスティック分布(log-logistic)が用意されている。
 次に対数正規分布を用いた例を示す。

>survreg(Surv(time, status) ~sex+disease, kidney,dist="lognormal")
survreg(Surv(time, status) ~sex+disease, kidney,dist=

 どの確率分布を用いた方がよいかに関しては、AICを用いて評価することができる。関数survregは、直接AICを返さないが、対数尤度Loglike(model)を返すので、AIC求めることは困難ではない。
 パッケージsurvivalには、生存分析に関する豊富なオブジェクトが用意されている。日本語に訳された一覧表は次のサイトからパッケージ名前survivalで検索できる。
http://www.okada.jp.org/RWiki/index.php?cmd=search
 パッケージsurvivalに関する邦文の参考文献としては[3][4]などがある。

参考文献
[1] 中村 剛(2004): Cox比例ハザードモデル、朝倉書店
[2] 丹後俊朗(2004): 統計モデル入門、朝倉書店
[3] 岡田 昌史 編(2004): The R Book―データ解析環境Rの活用事例集 (単行本) 、九天社
[4] W.N. ヴェナブルズ ・その他著(2001): S‐PLUSによる統計解析 、シュプリンガー・フェアラーク東京