[連載] フリーソフトによるデータ解析・マイニング 第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)

Call:
coxph(formula = Surv(time, status) ~ sex + disease, data = kidney)
n= 76, number of events= 58

coef exp(coef) se(coef) z Pr(>|z|)
sex
-1.4774 0.2282 0.3569 -4.140 3.48e-05 ***
diseaseGN
0.1392 1.1494 0.3635 0.383 0.7017
diseaseAN
0.4132 1.5116 0.3360 1.230 0.2188
diseasePKD
-1.3671 0.2549 0.5889 -2.321 0.0203 *

---
Signif. codes: 0 ‘***’  0.001 ‘**’  0.01 ‘*’  0.05 ‘.’  0.1 ‘ ’  1

exp(coef) exp(-coef) lower.95 upper.95
sex
0.2282 4.3815 0.11339 0.4594
diseaseGN
1.1494 0.8700 0.56368 2.3437
diseaseAN
1.5116 0.6616 0.78245 2.9202
diseasePKD
0.2549 3.9238 0.08035 0.8084

Concordance= 0.696 (se = 0.045 )
Rsquare= 0.206 (max possible = 0.993 )
Likelihood ratio test = 17.56 on 4 df, p=0.001501
Wald test        = 19.77 on 4 df, p=0.0005533
Score (logrank) test = 19.97 on 4 df, p=0.0005069

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

> kidney.fit<-survfit(kidney.cox)
> summary(kidney.fit)

Call: survfit(formula = kidney.cox)

time n.risk n.event survival std.err lower 95% CI upper 95% CI
2 76 1 0.99018 0.00985 0.971059 1.000
7 71 2 0.96856 0.01831 0.933335 1.000
8 69 2 0.94641 0.02422 0.900102 0.995
9 65 1 0.93499 0.02685 0.883820 0.989
12 64 2 0.91106 0.03177 0.850869 0.976
13 62 1 0.89844 0.03411 0.834007 0.968
<後略>

  関数 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 ) における仮説 H 0 : θ = 0 の検定に関する計算結果を返す。式の中の g (t ) は滑らかな関数である。また、関数 cox.zph は g (t i ) と標準化されたシェーンフィールド残差との相関係数をも返す。
  関数 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

rho chisq p
sex
0.18839 2.60676 0.106
diseaseGN
-0.01392 0.01123 0.916
diseaseAN
0.06162 0.20036 0.654
diseasePKD
0.00701 0.00438 0.947
GLOBAL
NA 4.20781 0.379

  仮説が棄却されると比例ハザードの仮定が充たされていない可能性があることを示唆する。ただし、仮説検定の結果は 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

Call:
coxph(formula = Surv(time, status) ~ (sex + age + disease)^2, data = kidney)

coef exp(coef) se(coef) z p
sex
-2.61e+00 7.33e-02 1.21e+00 -2.16 0.031
age
2.11e-02 1.02e+00 7.53e-02 0.28 0.780
diseaseGN
-7.85e-01 4.56e-01 3.04e+00 -0.26 0.796
diseaseAN
-1.88e+00 1.52e-01 2.58e+00 -0.73 0.466
diseasePKD
-1.44e+01 5.76e-07 5.59e+00 -2.57 0.010
sex:age
-4.40e-03 9.96e-01 4.04e-02 -0.11 0.913
sex:diseaseGN
9.95e-01 2.71e+00 1.34e+00 0.74 0.459
sex:diseaseAN
1.30e+00 3.67e+00 1.22e+00 1.06 0.288
sex:diseasePKD
3.71e+00 4.08e+01 1.69e+00 2.19 0.028
age:diseaseGN
-2.31e-02 9.77e-01 2.88e-02 -0.80 0.422
age:diseaseAN
-3.38e-03 9.97e-01 3.00e-02 -0.11 0.910
age:diseasePKD
1.40e-01 1.15e+00 1.11e-01 1.25 0.210

Likelihood ratio test=32.9 on 12 df, p=0.00102 n= 76, number of events= 58

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

> library(MASS); stepAIC(kidney.cox2)

Start: AIC=366.95
Surv(time, status) ~ (sex + age + disease)^2

Df AIC
- age:disease
3 362.99
- sex:age
1 364.96
<none>
366.95
- sex:disease
3 367.26

Step: AIC=362.99
Surv(time, status) ~ sex + age + disease + sex:age + sex:disease

Df AIC
- sex:age
1 361.36
<none>
362.99
- sex:disease
3 368.67
<中略>

Step: AIC=359.54
Surv(time, status) ~ sex + disease + sex:disease

Df AIC
<none>
359.54
- sex:disease
3 366.24

Call:
coxph(formula = Surv(time, status) ~ sex + disease + sex:disease, data = kidney)

coef exp(coef) se(coef) z p
sex
-2.525278 0.080036 0.565962 -4.46 8.1e-06
diseaseGN
-1.341579 0.261432 1.303772 -1.03 0.30348
diseaseAN
-1.418077 0.242179 1.484055 -0.96 0.33930
diseasePKD
-7.381915 0.000622 1.949378 -3.79 0.00015
sex:diseaseGN
0.831292 2.296284 0.760444 1.09 0.27432
sex:diseaseAN
1.075400 2.931166 0.815635 1.32 0.18734
sex:diseasePKD
4.214473 67.658476 1.150109 3.66 0.00025

Likelihood ratio test=30.3 on 7 df, p=8.5e-05
n= 76, number of events= 58

  このように機械的に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")

Call:
survreg(formula = Surv(time, status) ~ sex + disease, data = kidney, dist = "lognormal")
Coefficients:

(Intercept) sex diseaseGN diseaseAN diseasePKD
1.7923643 1.5062960 -0.3334601 -0.5321264 0.6810495

Scale= 1.129394
Loglik(model)= -329.1 Loglik(intercept only)= -340
Chisq= 21.8 on 4 degrees of freedom, p= 0.00022 n= 76

  どの確率分布を用いた方がよいかに関しては、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による統計解析 、シュプリンガー・フェアラーク東京