連載 第 5 回 医学データの統計解析の基本 生存時間データの解析 Fundamentals of statistical analysis in biomedical research:analysis of survival data 朝倉こう子 * * 濱﨑俊光 1 生物学や医学分野において

Transcription

1 連載 第 5 回 医学データの統計解析の基本 生存時間データの解析 Fundamentals of statistical analysis in biomedical research:analysis of survival data 朝倉こう子 * * 濱﨑俊光 1 生物学や医学分野において 個体 ( 例えば 培養細胞 動物あるいは被験者 患者 ) に対して化合物や医薬品といった処理を与えることにより惹き起こされる反応は あらかじめ定められたある一定の期間のなかで個体が追跡され 反応の有無と反応までの時間がともに観測されることが一般的である 例えば がんの臨床試験では 試験参加者が試験に組み入れられた後に 試験参加者が関心のある事象 ( 例えば 死亡や再発 ) を体験 発現する 試験が完了する あるいはデータが解析されるいずれかの時点まで追跡され 事象発現の有無とともに事象発現までの時間が記録される そして ある試験参加者については関心のある事象が試験期間内で観測されるが ある試験参加者については試験の途中で転院 転居などの理由で追跡不能となったり また別の試験参加者については十分な試験期間が設けられていたにも関わらずその事象を発現せず 試験の終了とともに追跡が打ち切られる 後者の 関心のある事象を観測できないままに追跡が打ち切られるような観測の打ち切りは 中途打ち切り ( 特に右側中途打ち切り ) と呼ばれ 実際の事象発現までの時間は未知であるが 観測された ( 打ち切りまでの ) 時間は実際の事象発現までの時間よりは短い 試験期間が長期にわたれば 関心のある事象を発現する試験参加者数は増えるが その一方で追跡不能となる試験参加者数が増える可能性が高まる このとき 関心のある事象の発現率や平均的な生存時間を計算したり さらには計算した発現率を治療群間で比較することにおいて 追跡が中途打ち切られた試験参加者のデータを適切に処理せねば 発現率や事象発現までの時間を過小に ( あるいは過大に ) 評価する恐れがある 第 5 回では 関心のある事象と事象発現までの時間の解析において用いられる代表的な手法として Kaplan-Meier 法 と ログランク検定 を概説する また 事例を通して統計計算ソフトウェアでのこれらの手法の実行の仕方を概説する 一般に 関心のある事象は死亡に限らず 疾患の増悪 疾患の再発といったように 死亡以外である場合も考えられるが ここでは多くの教科書にならい このような場面で得られるデータの総称として生存時間データと呼ぶ 2 Kaplan-Meier 法 (Kaplan & Meier,1958) は 生存時間データに基づいて所与 の時間 t における生存率 ( 生存率は時間 t の関数であることから生存時間関数と呼ば * Koko Asakura, Toshimitsu Hamasaki Office of Biostatistics and Data Management Department of Advanced Medical Technology Development National Cerebral and Cardiovascular Center 国立循環器病研究センター先進医療 治験推進部 DM/ 統計室 吹田市藤白台 toshi.hamasaki@ncvc.go.jp 474 Drug Delivery System 30 5, 2015

2 *1 中途打ち切りには あらかじめ打ち切りの時間が決まっている 時間打ち切り ( 試験期間など ) 関心のある事象があらかじめ定めた数まで確認された時点で観測を終了する 個数打ち切り 観測期間中の ( 転院や転居による ) 無作為な打ち切り などがあるが これらは 無情報な (non- informative) 打ち切り として 解析の際にこれら 3 つの中途打ち切りを区別することはない 他方 状態悪化のための治療 投薬の中止による中途打ち切りは 情報 (informative) 打ち切り と呼ばれ その取り扱いには注意がいる ここでは 中途打ち切りは 無情報な打ち切り であると仮定する れることもあるが 本稿では単に生存率と呼ぶ ) を求めるための 1 つの方法である Hosmer et al.(2008) に掲載された WHAS100 試験の被験者 5 名の生存時間データ ( 表 1) を用いて Kaplan-Meier 法の計算を例解する ( 生存時間の昇順に並びかえてい る ) 表 1 において 被験者番号 3 は 生存時間が中途打ち切り *1 されており 少なく とも 21 カ月までの生存が確認されている 残りの 4 名の被験者は死亡が確認され 生 存時間が観測されている 標本サイズを n 昇順の第 i 番目の生存時間の観測値を t (i) 時点 t における生存率を S(t) とすると Kaplan-Meier 法に基づく生存率の推定値 Ŝ(t) は Ŝ(t)= ( n - i n - i+1 t(i) t で与えられる δ (i) は昇順の第 i 番目の観測値が中途打ち切りであるときにゼロ 非中 途打ち切りのときに 1 をとる死亡 打ち切りを表す指標関数である この式に基づけ ば 表 2 のように Ŝ(t) を求めることができる 表 1 WHAS 試験の被験者 5 名の生存時間データ (Hosmer et al., ) 被験者番号 生存時間 ( 月 ) 死亡 中途打ち切り ( 打ち切り指標 ) 1 t (1) 6 死亡 δ (1) =1 4 t (2) 14 死亡 δ (2) =1 3 t (3) 21 中途打ち切り δ (3) =0 2 t (4) 44 死亡 δ (4) =1 5 t (5) 62 死亡 δ (5) =1 ) δ (i) 表 2 Kaplan Meier 法による WHAS 試験の生存率の推定 時間区間 0 t<6 6 t<14 14 t<21 21 t<44 44 t<62 62 t t (i) t を満たす観測値はないため Ŝ(t)=1.0 とする t (1) が t (i) t であるから Ŝ(t)=( 5-1 ) 1 = 4 = 生存率 Ŝ(t) t (1) と t (2) が t (i) t であるから Ŝ(t)=( 5-1 ) 1 ( 5-2 ) 1 =( 4 )( 3 ) = t (1),t (2) および t (3) が t (i) t であるから Ŝ(t)=( 5-1 ) 1 ( 5-2 ) 1 ( 5-3 ) 0 =( 4 )( 3 )1 = t (1),t (2),t (3) および t (4) が t (i) t であるから Ŝ(t)=( 5-1 ) 1 ( 5-2 ) 1 ( 5-3 ) 0 ( 5-4 ) 1 =( 4 )( 3 )1( 1 ) = t (1),t (2),t (3),t (4) および t (5) が t (i) t であるから Ŝ(t)=( 5-1 ) 1 ( 5-2 ) 1 ( 5-3 ) 0 ( 5-4 ) 1 ( 5-5 ) 1 =( 4 )( 3 )1( 1 )( 0 ) = これらの Kaplan-Meier 法に基づく生存率を95% 信頼区間とともに描けば図 1が得られる 図中の実線は生存率の推定値を 点線は95% 信頼区間を表している 図 1 から死亡が観測された6カ月 14カ月 44カ月 62カ月で生存率は減少し それの以外の時間区間では21カ月で打ち切られたとしても一定であり 階段関数として与えられることがわかる Kaplan-Meier 法に基づく生存率の推定値 Ŝ(t) の標準誤差 Drug Delivery System 30 5,

3 S E[ Ŝ(t)] は Greenwood の公式 (Greenwood,1926) によれば δ (i) SE[Ŝ(t)]=Ŝ(t) t (i) t(n -i)( n - i+1) で与えられる 例えば t=14における Kaplan-Meier 法に基づく生存率の推定値の標準誤差は SE[Ŝ(14)]= (5-1)( 5-1+1) (5-2)( 5-2+1) = 3 ( ) = である 標本サイズが大きくなるにともない Kaplan-Meier 法に基づく生存率の推定値 Ŝ(t) の分布は正規分布に近づく性質を利用して 図 1の点線のように生存率に 95% 信頼区間の帯を付すことができる (Gross & Clark,1973) たいていの統計計算ソフトウェアで採用されている Greenwood の公式は 生存率がゼロあるいは1に近いときに 実際の分散よりも過小評価することが知られており その場合には代替のものを用いることが推奨される (Peto et al.,1977) 生存率 死亡中途打ち切り 生存時間 図 1 Kaplan-Meier 法に基づく生存率の推定値のプロット Kaplan-Meier 法のほかに Kaplan-Meier 法に比べて計算量の小さい生命表に基づき生存率を推定することもできる ( 例えば Gross & Clark(1976) を参照 ) 生命表による生存率の推定は 保険数理推定値としても知られている Kaplan-Meier 法では事象が生じた時点で時間区間を構成するが 生命表では 事象を観測するかどうかにかかわらず 最初に観測期間をいくつかの時間区間に分割し その開始時点で生存率を推定する このため その区間のなかの生存率が過大に推定されることに注意がいる 2つの生存率を比較する代表的な解析手法として ログランク検定がある ログランク検定は 死亡が観測された各時点で 朝倉 濱﨑 (2015b) で紹介したような治療群 生存 死亡の2 2 分割表を生成する ここで 周辺度数が固定されていると 476 Drug Delivery System 30 5, 2015

4 考えれば 4 つのセルのうち 1 つのセルの頻度に注目すればよく この頻度は超幾何 分布に従う この性質に基づき各時点のセル頻度の観測度数と期待度数からカイ 2 乗 検定統計量を計算し これを Mantel-Haenszel 検定の仕方 ( 例えば Mantel(1966) を 参照 ) で統合したものがログランク検定統計量である ログランク検定の名称は Peto & Peto(1972) による 試験治療 (E) と対照治療 (C) の 2 つの生存率 S E (t) と S C (t) の比 較において ログランク検定における帰無仮説と対立仮説は 片側検定であれば { H 0: すべての t について S E (t)= S C (t) H 1 : ある t について S E (t)> S C (t) である ログランク検定の p 値があらかじめ定めた有意水準よりも小さければ 試験 治療の生存率は時間を通して対照治療よりも高いと解釈できる ログランク検定は 朝倉 濱﨑 (2015b) で紹介した割合の比較と同様に 正規近似に基づく手法である シミュレーションに基づく標本サイズと検定の実質的サイズ ( シミュレーションによ り計算した第 1 種の過誤確率 ) の関係を図 2 に示す ( 片側検定のもとで有意水準は 2.5% とし シミュレーションの反復回数は 100 万回とした ) 図 2 から示唆されるように 1 群あたりの標本サイズが 50 例よりも小さい場合に 第 1 種の過誤確率はあらかじめ 定めた有意水準 2.5% よりもわずかに大きくなり 2 つの生存率に違いがないにもか かわらずその違いを検出する可能性が若干高まる 標本サイズが小さい場合のログラ ンク検定の適用には 十分に注意を払う必要がある 検定の実質的サイズ (%) 群の標本サイズ 図 2 シミュレーションに基づく 1 群あたりの標本サイズ ( 等例数 ) とログランク検定の実質的サイズの関係 実線は検定の実質的サイズを 点線は有意水準を表す 片側検定の有意水準は 2.5% とし シミュレーションの反復回数は 万回とした また シミュレーションでは試験の登録期間を 2 年 追跡期間を 3 年とした ログランク検定は 各時点のカイ2 乗検定の統合において すべての重みを等しく1としたものであるが 諸種の重みを考慮し統合することができる ログランク検定の代替法の1つとして 朝倉 濱﨑 (2015a) で紹介した Wilcoxon の順位和検定を中途打ち切りデータまでとり扱えるように拡張した一般化 Wilcoxon 検定 ( あるいは Gehan 検定 )(Gehan,1965;Breslow,1970) は 各時点直前までの生存数 (Number of at Risk ) で重みを付したものである 実地でよく観測される生存率の形状を図 3に示す Drug Delivery System 30 5,

5 1 生存率 0 時間 (a) 時間 (b) 時間 (c) 図 3 実地でよく観測される生存率の形状 図中の2 本の曲線は 2つの群の生存率を表している この2つの生存率を比較するとき ログランク検定は (a) の示す生存率の違いを検出することには優れているが すべての重みを等しく1として統合するため (b) の生存率の違いを検出することには不向きである この場合 各時点直前までの生存数で重み付けし最初の生存率の違いを強調できる一般化 Wilcoxon 検定がログランク検定よりも優れている (c) はログランク検定と一般化 Wilcoxon 検定のいずれでも違いを検出することが難しい *2 生存時間に関するテストデータと解析のための関数がまとめられたパッケージであり CRAN(R 本体やパッケージをダウンロードするための web サイト ) からダウンロードしインストールすることで ovarian データも得られる パッケージのインストールの方法は例えば舟尾 (2009) を参照されたい 3 先に紹介した Kaplan-Meier 法に基づく生存率 ログランク検定と一般化 Wilcoxon 検定を統計計算ソフトウェアで実行する SAS University Edition R エクセル統計 JMP の4つの統計計算ソフトウェアを使用し R の survival パッケージ *2 に組み込まれているテストデータ ovarian を解析する ( 図 4) ovarian は卵巣がんに対する外科手術後に異なる2つの化学療法を実施した場合の有効性を検討した確率化臨床試験のデータである 化学療法の種類は rx という変数で得られており シ 図 4 テストデータ ovarian 478 Drug Delivery System 30 5, 2015

6 クロホスファミド単独 (rx=1) もしくは アドリアマイシンとの併用(rx=2) である 治療の有効性を検討する指標として 生存時間または打ち切りまでの時間が futime 打ち切りに関する指標が fustat (0: 打ち切り 1: 死亡 ) として得られている いま Kaplan-Meier 法に基づく生存率の推定値を計算 描画し 2つの化学療法の効果の比較のために 2 群の生存率についてログランク検定または一般化 Wilcoxon 検定にて比較する それぞれのソフトウェアについて SAS Studio 3.4 R version エクセル統計 2015および JMP Pro にて動作確認を行っている 以降ではそれぞれのソフトウェアにおける実行と出力される結果について個々に述べる 1SAS University Edition: 生存時間解析を適用する LIFETIME プロシジャを用いて Kaplan-Meier 法に基づく生存率を描画し 化学療法に関する2 群比較のための検定を実行する手順を紹介する PLOTS ステートメントで PLOTS=S と指定すれば Kaplan-Meier 法に基づく生存率の推定値が描画される このとき CL オプションを指定 ( PLOTS=S(CL) ) すれば 生存率の95% 信頼区間が同時に描画される TIME ステートメントで生存時間を示す変数 ( 時間変数 ) 打ち切りの指標および打ち切りを示す値を指定し( ここでは futime*fustat(0) となる) STRATA ステートメントで群を示す変数 ( 群変数 ここでは rx ) を指定する ただし現時点では SAS University Edition の表示言語が日本語になっている場合に STRATA ステートメントを指定すると エラーにより結果が適切に表示されないため そのような場合には表示言語を英語にする ここで示す出力結果は 言語を英語にしてプログラムを実行した場合のものである 以上の手順を実行すれば 各群について生存率の推定値が図 5のように描画され ログランク検定と一般化 Wilcoxon 検定の結果が図 6のように出力される 図 5 図 6 図 5 LIFETEST プロシジャによる Kaplan Meier 法に基づく生存率の推定値のプロット (SAS University Edition) Drug Delivery System 30 5,

7 図 6 LIFETEST プロシジャによるログランク検定と一般化 Wilcoxon 検定の実行結果 (SAS University Edition) は すべての出力のうち例示に用いる部分のみを抜粋している これより 視覚的にはアドリアマイシンとの併用群がシクロホスファミド単独群に比べて生存率が高いように見えるが ログランク検定による p 値は0.3026であり 有意水準 0.05のもとで 2 群の生存率に差があるとはいえないと結論づけられる また 一般化 Wilcoxon 検定による p 値は0.1665であり ログランク検定に基づく p 値に比べ小さいが ログランク検定と一般化 Wilcoxon 検定で導かれる結論は変わらない 2R: ここでは survival パッケージを用いた手順を紹介する Kaplan-Meier 法 ログランク検定のいずれにおいても survival パッケージを用いるため 事前にパッケージを読み込んでおく まず survfit 関数を用い Kaplan-Meier 法に基づく生存率の推定値を計算する その際 survfit 関数の引数として指定する式は Surv( 時間変数, 打ち切りの指標 ) となる 時間変数としてここでは futime 打ち切りの指標として fustat を指定する さらに群変数として rx およびデータ名として ovarian を指定する これにより群ごとに生存率の推定値が計算される 結果がオブジェクトとして保存されるよう名前を付け ( ここでは km とする ) 実行する 次に plot 関数を用い 先ほど保存したオブジェクトを描画する ( plot(km) ) これにより 480 Drug Delivery System 30 5, 2015

8 図 7のように Kaplan-Meier 法に基づく生存率の推定値が描画される この際 conf.int=true と指定すれば 95% 信頼区間があわせて描画される また検定を行う場合は survdiff 関数を用い survfit 関数の場合と同様に式 群変数およびデータを指定し実行すれば 治療群間での生存率の比較のための検定が実行される 検定の種類を特に指定しなければ 規定値でログランク検定の結果が出力される ( 図 8) また survdiff 関数の実行の際に rho=1 と指定すれば 一般化 Wilcoxon 検定の結果が出力される ( 図 9) ただし R で適用される一般化 Wilcoxon 検定は Peto-Peto の修正版 (Peto & Peto,1972) である なお R では打ち切りの指標を通常は0( 打ち切り ) と1( 事象発現 ) で与える必要がある ( 他に FALSE( 打ち切り )/TRUE( 事象発現 ) 1( 打ち切り )/2( 事象発現 ) でも可能 ) 図 7 survfit 関数と plot 関数による Kaplan Meier 法に基づく生存率の推定値のプロット (R) 図 8 survdiff 関数によるログランク検定の実行結果 (R) 図 9 survdiff 関数による一般化 Wilcoxon 検定の実行結果 (R) Drug Delivery System 30 5,

9 出力された結果より ログランク検定に基づく p 値は SAS による実行結果に一致するが 一般化 Wilcoxon 検定については SAS において規定値で適用される検定とはカイ2 乗検定の統合における重みが異なるため 結果が若干異なる しかしここでは導かれる結論は同じである 3エクセル統計 : 生存分析 ハザード分析 のなかの カプラン = マイヤー法 ( 表形式 ) または カプラン = マイヤー法 ( データベース形式 ) を用いる いずれを用いるかはデータの形式によるが ここでは他のソフトウェアで解析する場合と同じ形式 ( 図 4 の形式 ) でデータが入力されているとし カプラン = マイヤー法 ( データベース形式 ) を用いる 出てくるウィンドウの 変数 タブで 時間 に時間変数 futime を 状態 に打ち切りの指標 fustat を グループ変数 に群変数 rx を指定する エクセル統計では 打ち切りの指標は0( 打ち切り ) と1( 事象発現 ) で与える すると Kaplan- Meier 法に基づく生存率の推定値とその描画 検定結果が出力される ( 図 10) 出力結果より ログランク検定については Cochran-Mantel-Haenszel が SAS やRにて前記の手順で実行した検定と対応している 一般化 Wilcoxon 検定については Gehan- Breslow が SAS での前記の手順で実行した検定と対応しており それらの結果は一致している 図 10 生存分析 ハザード分析 による Kaplan-Meier 法に基づく生存率の推定値のプロットとログランク検定および一般化 Wilcoxon 検定の実行結果 ( エクセル統計 ) 482 Drug Delivery System 30 5, 2015

10 4JMP : 分析 の 信頼性/ 生存時間分析 のなかの 生存時間分析 を用い Y, 事象発現までの時間 に時間変数 futime グループ変数 に群変数 rx 打ち切り に打ち切りの指標 fustat を指定する また 打ち切りの値 には打ち切りを表す値 ( ここでは0) を指定する すると Kaplan-Meier 法に基づく生存率の推定値の描画や検定結果などが出力される ( 図 11) なお 図 11は 出力のうち例示に必要な部分のみ 図 11 生存時間分析 による Kaplan-Meier 法に基づく生存率の推定値のプロットとログランク検定および一般化 Wilcoxon 検定の実行結果 (JMP ) Drug Delivery System 30 5,

11 を抜粋したものである 出力結果より ログランク検定の結果は SAS R エクセ ル統計による実行結果と一致し 一般化 Wilcoxon 検定の結果は SAS による実行結 果 ( すなわちエクセル統計では Gehan-Breslow ) と一致している 文献 1)Breslow, N.E. (1970). A generalized Kruskal-Wallis test for comparing K samples subject to unequal patterns of censorship. Biometrika 57, )Gehan, E.A. (1965). A generalized Wilcoxon test for comparing arbitrarily singly-censored samples. Biometrika 52, )Greenwood, M. (1926). The natural duration of cancer. Reports on Public Health and Medical Subjects 33, )Gross, A.J., Glark, V.A. (1976). Survival Distributions: Reliability Applications in the Biomedical Sciences. New York: John Wiley & Sons. 5)Hosmer, D.W., Lemeshow, S., May, S. (2008). Applied Survival Analysis: Regression Modeling of Timeto-Event Data. Hoboken: John Wiley & Sons. 6)Kaplan, E.L., Meier, P. ( ). Nonparametric estimation from incomplete observations. Journal of the American Statistical Association 53, )Mantel, N. (1966). Evaluation of survival data and two new rank order statistics arising in its consideration. Cancer Chemotherapy Reports 50, )Peto, R., Peto, J. (1972). Asymptotically efficient rank invariant test procedure. Journal of the Royal Statistical Society. Series A135, )Peto, R., Pike, M.C., Armitage, P., Breslow, N.E., Cox, D.R., Howard, S.V., Mantel, N., McPherson, K., Peto, J., Smith, P.G. (1977). Design and analysis of randomized clinical trials requiring prolonged observation of each patient. Part II: analysis and examples. British Journal of Cancer 35, ) 朝倉こう子 濱﨑俊光 ( a). 2 つの平均の比較. Drug Delivery System 30(2), ) 朝倉こう子 濱﨑俊光 ( b). 2 つの割合の比較. Drug Delivery System 30(3), ) 舟尾暢男 (2009). The R Tips: データ解析環境 R の基本技 グラフィック活用集. 第 2 版. オーム社 X 線還元反応 (Radiolytic Reduction) DDS 用語解説 No.300 X 線照射は極めて高いエネルギーを持つ外部刺激であり 物質のイオン化 ラジカル形成さらに高エネルギー励起活性種の形成を通じて 化学結合の切断および形成に関わる 生体を構成する重要な成分である水に X 線を照射すると イオン化や励起が引き起こされ 数段階を経て 水素原子 水和電子 ヒドロキシルラジカルが主要な分解生成物として そして反応活性種として形成される このうち 水和電子は非常に高い還元反応特性を示し 水中に溶解している溶質を還元することが知られている 筆者らはこれまでに インドールキノン基 アジドメチル基 ジスルフィド結合などの官能基を持つ化合物の水溶液に X 線を照射すると 水和電子による還元反応が進行し 化学結合開裂が生じることを見出している これら置換基を活用して X 線照射によって活性化されるプロドラッグの構築を進めている 田邉一仁 (Kazuhito Tanabe) 青山学院大学理工学部化学 生命科学科 484 Drug Delivery System 30 5, 2015