群馬大学 | 医学部 | サイトトップ | Software Tips | RについてのTips

生存時間解析入門

最終更新: 2009年 6月 18日 (木曜日) 20時41分 :作成

このページは,群馬大学公衆衛生学教室において毎週行っている疫学・社会医学勉強会の月曜日枠「R実践活用」で,2008年11月17日から2009年6月15日まで,約半年をかけて読んだ,「観察疫学研究における生存時間解析」というチュートリアル論文(Bull K, Spiegelhalter DJ: Tutorial in biostatistics: Survival analysis in observational studies. Statistics in Medicine, 16: 1041-1074, 1997)の内容を紹介し,Rを使ってそれを実行する方法と,出力結果の解釈について解説する目的で設置する。Rを使った生存時間解析の一助になれば幸いである。

1. 緒言

1.1. 背景

多施設データベースが確立したものになってくるにつれて,リスク因子と時間に臨床的な結果を関連付ける報告が多数出現しつつある。これらの研究は,多様な目的をもっている。ある一群の患者が経験することの記述であるかもしれないし,リスク因子の同定であるかもしれないし,意思決定のための個人についての予測であるかもしれない。そして,増加しつつあるのが,施設間比較,あるいは施術者間比較のための標準化(リスク層別化)という目的である。研究者はまた,代替的な介入あるいは臨床戦略の効果についても結論を引き出したいと思うかもしれない(データベースの分析から治療の利益についての判断を下すことの危険性は,既によく議論されてきてはいるけれども)。

1.2. この論文の構成

2. データ

2.1. バイアスについての一般的問題

2.2. 実例となるデータ

Patientagepresagelastageop1deadsexpaanatadfolfollow-upopfpresunopageunopfprepreopdedhadopded-lyrppage-presx
111274-100011273-1127412730000
221234010111213840380110
32119-11101117-11191171010
43120-10100117-11201170020
5610-110014-11041010
665415194001154091881941880100
7732611041010132541034104110340100
881819-101011811-1181918110000
911696-10001685-16966850000
1013641529011164021629160100
11293127144100130981151441150100
12304234710013931747170100
13355794-101015759-1579457590000
14452926211112471762170110
155468-1110114-168141010
16581849-110011791-1184917911000
1768343-11101275-13432751010
1810932761294000131671185129411850100
1911920720710118888207880110
20121143012301011309212320100
21231308237101177623760110
22258347-1000089-1347890020
2334933553830001300634383340100
2436933512826101129822457282624570101
254375474410000110444140121
2677138348680001306397868970101
2712857209-101015924-1720959240001
28245535553532111111001077353210770101
29516153545353111119319253531920111
30549756395633101114213656331360111

表のままだと欠損値が-1になっているので,Rで扱うにはNAに置き換えねばならない。学生がExcelで入力してくれたファイルをタブ区切りテキスト形式で保存し,-1をNAに置換したファイルを,surv.txtとして置いておく。

3. 時間間隔を固定した場合の結果指標の単純な記述

3.1. 分析の特定

3.1.1. データを解析に含める基準
3.1.2. 結果の指標(イベント,反応,従属変数としても知られる)

3.2. 作業例:来院から1年以内に死亡した割合

Fig.1. and Fig.2

上図を描くためのコードは,plotx.Rである。画像はファイル出力されてしまうので注意されたい。

3.3. 注意事項と推論

4. 基本的な生存時間解析(説明変数なしの場合)

4.1. はじめに

4.2. 分析の特定

4.2.1. 時間の始点
4.2.2. 研究への参加開始
4.2.3. 研究からの脱落(打ち切り[censoring])
さまざまな理由で,観察対象イベントが発生する前に,追跡中なのに対象者が消えてしまうことがある。予定された追跡期間が終わっても観察対象イベントが発生していないこともある。追跡終了時に,これらの対象者がどういう状態にあるかわかっていないことになるので,その扱いには注意が必要である。つまり,この人たちが研究からいなくなった理由がリスクと無関係だと仮定できなくてはいけない(この仮定は,情報量の無い打ち切り[non-informative censoring]として知られている)。もし,リスクの高い人あるいはリスクの低い人が選択的に追跡から脱落していたら,バイアスの入った結果が得られてしまう。
(注)打ち切り例を分析から完全に除去してしまうと,打ち切り例には「予定された追跡期間が終わってもイベントが起こっていない人」が含まれていることから,打ち切りでなかった人よりも生存時間が長い人が多く含まれている可能性があることになり,その結果として,生存時間を低めに見積もってしまうバイアスが生まれる可能性がある。従って,データそのものから打ち切り例を除くのは情報量のあるバイアス[informative bias]の原因となって不適切である。そのため,打ち切り例を含めて適切に分析する方法が特別に必要になる。
4.2.4. 生存時間
4.2.5. 観察期間

4.3. ノンパラメトリックな生存関数

カプラン=マイヤ法による生存曲線の推定法を説明している。打ち切り例は,リスク集合から順次除かれていくことにより,適切に処理される。

4.4. 作業例:来院からの生存(図2に対応)

表3にカプラン=マイヤ法の計算が示されている。表3と図3に掲載されている情報がわかるプログラムコードをweibull.Rに示す。実行すると下図がpng形式のファイルとして出力されてしまうので注意されたい。

Fig.3. KM and Weibull estimates of survival.

ちなみに,後で詳述するが,ワイブル関数をあてはめた結果は以下の通り。

Call:
survreg(formula = eventtime ~ 1, dist = "weibull")
            Value Std. Error     z        p
(Intercept) 8.494      0.600 14.15 1.81e-45
Log(scale)  0.775      0.221  3.51 4.52e-04

Scale= 2.17 

Weibull distribution
Loglik(model)= -128.6   Loglik(intercept only)= -128.6
Number of Newton-Raphson Iterations: 5 
n= 30

カプラン=マイヤ法を数学的に表現すると,時点tkにおいてk番目の死が起こり,その時点でのリスク集合のサイズがrkだったとし,その時点で起こった死亡数がfkだったとすれば,時点tkまでの生存確率の推定値pkは,
pk=(1-f1/r1)(1-f2/r2)...(1-fk/rk)
となる。

4.5. ノンパラメトリックな生存関数:遅れて参加すること(late entry)の影響

左側切捨て[left-truncation]あるいは遅れて参加すること[late entry]の問題。例えば子供時代の生存については,成人になってから研究に参加した患者についての情報では何もいえないが,通常,右側打ち切りと同じく,左側切捨ても情報量が無い[non-informative late entry]と仮定する。即ち,同じ年齢で既に研究に参加していた人と本質的に同等であったと見なす。

4.6. 遅れて参加した場合を含む生存時間の作業例:誕生からの生存(図1に対応)

表4(a)に示すように,誕生からの生存年数は普通にカプラン=マイヤ法で計算できるが,リスク集合の中に,そのイベント発生年齢では研究にまだ参加していなかった人を含めてしまうのは不合理である。そこで,そういう人をリスク集合から除いてカプラン=マイヤ法をやり直すと,表4(b)のように計算でき,下図のようにグラフが描ける。これを実行するRコードは,table4fig4.Rである(これも実行するとカレントディレクトリにpngファイルができてしまうので注意)。この計算はあまりなされないため,survfit()関数に元々は組み込まれておらず,いちいち計算させる必要があって,多少プログラムが複雑になっている。

Fig. 4. KM estimates of survival from birth: time origin birth contrasted with time origin at presentation.

4.7. パラメトリックな生存関数

ノンパラメトリックなカプラン=マイヤ法は,データを使って直接生存曲線を描くもので,どんな生存曲線にも対応する。しかし,データによっては,何らかの仮定をおけば指数分布とかワイブル分布といった単純な分布を生存曲線に当てはめることが可能なことがある。ワイブル曲線を当てはめ,信頼区間付きで描画したのが図3である。パラメトリックな曲線の信頼区間は,ノンパラメトリックな曲線の信頼区間よりも一般に狭いが,この精度向上は,仮定を増やすことを代償として得られたものであり,バイアスを増やしてしまうことにつながるかもしれない。

4.8. ハザード関数

4.9. 計算

4.10. 推論と注意事項

(a) 情報量のある右側打ち切り
(b) 情報量のある遅れての参加

5. ある固定した時間での結果:ある時点での1因子の曝露

固定時間間隔内でのイベントの単純な記述は,可能な説明変数の影響を調べる場合に容易に拡張できて,臨床的に,より役に立つようになる。

5.1. 分析の特定

ランダム化試験では,すべての説明変数(リスク因子,共変量,予測因子あるいは独立変数としても知られている)が,ランダム化の時点において定義されていて,ランダム化という行為によって(偶然の変動を別にすれば)異なる治療群間でバランスがとれていることが保証されている。

5.2. 作業例:来院登録1年以内の死亡

ここで調べるのは,肺動脈の解剖学的所見(paanat)と来院時年齢が1歳未満かどうか(agepresx)の影響である。表5は,これらの死亡イベントへの影響を1つずつ別々にみた場合(左側の単変量解析と書かれている方)と,同時にみた場合(右側の多変量解析と書かれている方)をまとめたものである。

5.3. 計算

Rで実行するコードはtable5.Rである。

5.4. 注意事項と推論

6. 固定的な説明変数1つのある生存(時間の分析)

K-M(カプラン=マイヤ法)で推定された生存関数は,患者の2つのカテゴリについて容易に計算できることは明らかである。

6.1. 作業例:リスクの異なるグループにおける生存

表6と図6に肺動脈の解剖学的所見が異なる2群(paanat=0とpaanat=1)別々にK-Mで生存関数を推定した結果を示した。Rで書いたコード[table6figure6.R]を実行すると,表6を作るために必要な情報と下図が得られる(直接ファイル出力されるので注意)。

図6

6.2. ハザードとその標準誤差の計算

表6にはイベントが起こった瞬間のハザードと1年当たりのハザード及びその標準誤差(Cox and Oakesによる)が,2群別々に,1年未満と1〜10年までの2区間に分けて計算して示されている。例えばpaanat=0群について来院1年未満での死亡イベントが発生した時点でのリスク集合は20, 19, 16, 14人で,各1名ずつの死亡が起こっているので,

f <- c(1,1,1,1)
r <- c(20,19,16,14)
sum(f/r)
sqrt(sum(f/(r*(r-f))))

6.3. 2群間での生存の比較

選んだ時間の生存,例えば1年間の生存を比較することは,観察された生存の差とその推定された標準誤差に基づいて近似的なp値を計算することによって可能である。多くの他の目的のために,このような2群間での生存経験全体の比較は望ましいことである。そのためには,ログランク検定あるいはコックス=マンテル検定を要する。

6.4. 作業例:異なる肺動脈解剖所見をもつ患者群間での生存の比較のためのログランク検定

表7にログランク検定の結果とハザード比が示されている。この検定を実行するためのRのコードはlogrank.Rだが,出力は以下のようになる。

Call:
survdiff(formula = Surv(followup, dead) ~ paanat, data = dat)

          N Observed Expected (O-E)^2/E (O-E)^2/V
paanat=0 20        7    10.48      1.15      3.87
paanat=1 10        8     4.52      2.68      3.87

 Chisq= 3.9  on 1 degrees of freedom, p= 0.0493 

[1] 2.648137

表7と微妙に異なる結果なのは丸め誤差であり,もちろんRの結果の方が正しい。

7. ある固定した時間における結果:1つより多い説明因子がある場合

7.1. ロジスティック回帰分析を用いて調整されたオッズ比

複数の説明変数の目的変数への影響を,互いに他の変数の影響を調整しながら求めるための多変量解析として,ロジスティック回帰分析が行われる。モデルは,
p/(1-p)=d0×d1×…×dIである。

7.2. 計算

モデルの両辺の対数をとって右辺を線型結合に直せば,係数は比較的容易に求められる。求めた係数は,対数オッズ比なので,指数をとってオッズ比に直す。結果は表8に示されている。

この結果を得るためのRのコードは,既に5.3.で示した通り,table5.Rである。ロジスティック回帰分析に該当する部分の出力は下記の通り。

(Intercept)      paanat    agepresx 
  0.3382121   6.9484700   0.3450070 

                 2.5 %     97.5 %
(Intercept) 0.09415143  0.9765344
paanat      1.13747110 62.0705248
agepresx    0.02601249  2.9421095

Call:
glm(formula = dedlyrpp ~ paanat + agepresx, family = binomial, 
    data = dat2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5550  -0.7633  -0.7633   0.8421   1.6586  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -1.0841     0.5800  -1.869   0.0616 .
paanat        1.9385     0.9864   1.965   0.0494 *
agepresx     -1.0642     1.1712  -0.909   0.3636  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 35.594  on 26  degrees of freedom
Residual deviance: 31.115  on 24  degrees of freedom
AIC: 37.115

Number of Fisher Scoring iterations: 4

表8と若干数値が違うのは,ほぼ丸め誤差によると思われる。

7.3. 注意事項

8. 生存−固定因子がたくさんあるとき

8.1. すべての生存経験を用いたコックス回帰

比例ハザード性を仮定すれば,時点tにおけるハザード比は,
h(t) = h0(t)×h1×…×hiと表現できる。分析実行のためのRコードはtable9.Rであり,出力をまとめ直せば表9ができる。

8.2. 計算

8.3. 注意事項

9. 時間に依存する因子をもつ生存

(読者への警告:もしこれまでの節が難しく感じられたら,たぶん,いまが第11節の考察に飛ぶべきところだ!)とあったので,この勉強会で読むのはやめた。

9.1. 時間とともに変わる因子

9.2. 避けるべきエラー

9.3. 注意事項

10. 複雑な生存分析

10.1. 遅れて参加した事例を含むコックスモデリングを使った,固定効果をもつ因子と時間に依存した効果

10.2. 計算

10.3. 注意事項

11. 考察


リンクと引用について

Correspondence to: minato@phi.med.gunma-u.ac.jp.