1
生物学や医学分野において、個体（例えば、培養細胞、動物あるいは被験者・患者） に対して化合物や医薬品といった処理を与えることにより惹き起こされる反応は、あ らかじめ定められたある一定の期間のなかで個体が追跡され、反応の有無と反応まで の時間がともに観測されることが一般的である。例えば、がんの臨床試験では、試験 参加者が試験に組み入れられた後に、試験参加者が関心のある事象（例えば、死亡や 再発）を体験･発現する、試験が完了する、あるいはデータが解析されるいずれかの 時点まで追跡され、事象発現の有無とともに事象発現までの時間が記録される。そし て、ある試験参加者については関心のある事象が試験期間内で観測されるが、ある試 験参加者については試験の途中で転院･転居などの理由で追跡不能となったり、また 別の試験参加者については十分な試験期間が設けられていたにも関わらずその事象を 発現せず、試験の終了とともに追跡が打ち切られる。後者の、関心のある事象を観測 できないままに追跡が打ち切られるような観測の打ち切りは、中途打ち切り（特に右 側中途打ち切り）と呼ばれ、実際の事象発現までの時間は未知であるが、観測された （打ち切りまでの）時間は実際の事象発現までの時間よりは短い。試験期間が長期にわ たれば、関心のある事象を発現する試験参加者数は増えるが、その一方で追跡不能と なる試験参加者数が増える可能性が高まる。このとき、関心のある事象の発現率や平 均的な生存時間を計算したり、さらには計算した発現率を治療群間で比較することに おいて、追跡が中途打ち切られた試験参加者のデータを適切に処理せねば、発現率や 事象発現までの時間を過小に（あるいは過大に）評価する恐れがある。第5回では、関 心のある事象と事象発現までの時間の解析において用いられる代表的な手法として、 「Kaplan-Meier法」と「ログランク検定」を概説する。また、事例を通して統計計算ソ フトウェアでのこれらの手法の実行の仕方を概説する。一般に、関心のある事象は死 亡に限らず、疾患の増悪、疾患の再発といったように、死亡以外である場合も考えら れるが、ここでは多くの教科書にならい、このような場面で得られるデータの総称と して生存時間データと呼ぶ。2
Kaplan-Meier法（Kaplan&Meier,1958）は、生存時間データに基づいて所与 の時間t における生存率（生存率は時間t の関数であることから生存時間関数と呼ば■ 連載 第 5 回
朝倉こう子
＊・濱﨑俊光
＊ ＊KokoAsakura,ToshimitsuHamasaki OfficeofBiostatisticsandDataManagement DepartmentofAdvancedMedicalTechnologyDevelopment NationalCerebralandCardiovascularCenter 国立循環器病研究センター先進医療・治験推進部 DM/ 統計室 〒 565―8565 吹田市藤白台 5―7―1E-mail:[email protected]医学データの統計解析の基本
生存時間データの解析
Fundamentalsofstatisticalanalysisinbiomedicalresearch：Analysisofsurvivaldataれることもあるが、本稿では単に生存率と呼ぶ）を求めるための1つの方法である。 Hosmeretal.（2008）に掲載された WHAS100試験の被験者5名の生存時間データ （表1）を用いて、Kaplan-Meier法の計算を例解する（生存時間の昇順に並びかえてい る）。表1において、被験者番号3は、生存時間が中途打ち切り＊1されており、少なく とも21カ月までの生存が確認されている。残りの4名の被験者は死亡が確認され、生 存時間が観測されている。標本サイズを n、昇順の第i番目の生存時間の観測値を t（i）、 時点t における生存率を S（t）とすると、Kaplan-Meier法に基づく生存率の推定値Ŝ（t）は Ŝ（t）=
∏
（
n－i）
δ(i) n－i＋1 t （i）≤t で与えられる。δ（i）は昇順の第i番目の観測値が中途打ち切りであるときにゼロ、非中 途打ち切りのときに1をとる死亡・打ち切りを表す指標関数である。この式に基づけ ば、表2のように Ŝ（t）を求めることができる。 これらの Kaplan-Meier法に基づく生存率を95%信頼区間とともに描けば図1が得 られる。図中の実線は生存率の推定値を、点線は95％信頼区間を表している。図1 から死亡が観測された6カ月、14カ月、44カ月、62カ月で生存率は減少し、それ の以外の時間区間では21カ月で打ち切られたとしても一定であり、階段関数として 与えられることがわかる。Kaplan-Meier法に基づく生存率の推定値Ŝ（t）の標準誤差表1 WHAS 試験の被験者5 名の生存時間データ（Hosmer et al.,2 0 0 8）
被験者番号 生存時間（月） 死亡・中途打ち切り（打ち切り指標） 1 t（1） 6 死亡 δ（1）=1 4 t（2） 1 4 死亡 δ（2）=1 3 t（3） 2 1 中途打ち切り δ（3）=0 2 t（4） 4 4 死亡 δ（4）=1 5 t（5） 6 2 死亡 δ（5）=1 表2 Kaplan―Meier 法による WHAS 試験の生存率の推定 時間区間 生存率Ŝ（t） 0≤t<6 t（i）≤t を満たす観測値はないため Ŝ（t）=1.0とする。 6≤t<14 t（1）がt（i）≤t であるから Ŝ（t）=
（
5 －1）
1 = 4 = 0.8 5 －1 ＋1 5 14≤t<21 t（1）とt（2）がt（i）≤tであるから Ŝ（t）=（
5 －1）
1（
5 －2）
1 =（
4）
（
3）
= 0.6 5 －1 ＋1 5 －2 ＋1 5 4 21≤t<44 t（1），t（2）およびt（3）がt（i）≤t であるから Ŝ（t）=（
5 －1）
1（
5 －2）
1（
5 －3）
0 =（
4）
（
3）
1 = 0.6 5 －1 ＋1 5 －2 ＋1 5 －3 ＋1 5 4 44≤t<62 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）
= 0.3 5 －1 ＋1 5 －2 ＋1 5 －3 ＋1 5 －4 ＋1 5 4 2 62≤t 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）
= 0.0 5 －1 ＋1 5 －2 ＋1 5 －3 ＋1 5 －4 ＋1 5 －5 ＋1 5 4 2 1 ＊1 中 途 打 ち 切 り に は 、 あ ら か じ め 打 ち 切 り の 時 間 が 決 ま っ て い る「 時 間 打 ち 切 り 」（ 試 験 期 間 な ど ）、 関 心 の あ る 事 象 が あ ら か じ め 定 め た 数 ま で 確 認 さ れ た 時 点 で 観 測 を 終 了 す る「 個 数 打 ち 切 り 」、 観 測 期 間 中 の（ 転 院 や 転 居 に よ る ） 「 無 作 為 な 打 ち 切 り 」な ど が あ る が 、 こ れ ら は「 無 情 報 な（ n o n - i n f o r m a t i v e ） 打 ち 切 り 」と し て 、 解 析 の 際 に こ れ ら 3 つ の 中 途 打 ち 切 り を 区 別 す る こ と は な い 。 他 方 、 状 態 悪 化 の た め の 治 療・ 投 薬 の 中 止 に よ る 中 途 打 ち 切 り は「 情 報 （ i n f o r m a t i v e ）打 ち 切 り 」 と 呼 ば れ 、 そ の 取 り 扱 い に は 注 意 が い る 。 こ こ で は 、 中 途 打 ち 切 り は「 無 情 報 な 打 ち 切 り 」で あ る と 仮 定 す る 。SE［Ŝ（t）］は、Greenwood の公式（Greenwood,1926）によれば SE［Ŝ（t）］=Ŝ（t）
∑
δ(i) （n－i）（n－i＋1） t （i）≤t で与えられる。例えば、t=14における Kaplan-Meier法に基づく生存率の推定値の標 準誤差は SE［Ŝ（14）］= 35 （5－1）（5－1＋1） （5－2）1 ＋ （5－2＋1）1 ＝（
35 20 121＋ 1）
＝ 0.283 である。標本サイズが大きくなるにともない、Kaplan-Meier法に基づく生存率の推 定値Ŝ（t）の分布は正規分布に近づく性質を利用して、図1の点線のように生存率に 95%信頼区間の帯を付すことができる（Gross&Clark,1973）。たいていの統計計算 ソフトウェアで採用されている Greenwood の公式は、生存率がゼロあるいは1に近 いときに、実際の分散よりも過小評価することが知られており、その場合には代替の ものを用いることが推奨される（Petoetal.,1977）。 Kaplan-Meier法のほかに、Kaplan-Meier法に比べて計算量の小さい生命表に基づ き生存率を推定することもできる（例えば、Gross&Clark（1976）を参照）。生命表に よる生存率の推定は、保険数理推定値としても知られている。Kaplan-Meier法では 事象が生じた時点で時間区間を構成するが、生命表では、事象を観測するかどうかに かかわらず、最初に観測期間をいくつかの時間区間に分割し、その開始時点で生存率 を推定する。このため、その区間のなかの生存率が過大に推定されることに注意がい る。 2つの生存率を比較する代表的な解析手法として、ログランク検定がある。ログラ ンク検定は、死亡が観測された各時点で、朝倉・濱﨑（2015b）で紹介したような治 療群×生存・死亡の2×2分割表を生成する。ここで、周辺度数が固定されていると 図1 Kaplan-Meier 法に基づく生存率の推定値のプロット 生存率 生存時間 0 6 14 21 44 62 1.0 0.8 0.6 0.3 0.0 死亡 中途打ち切り考えれば4つのセルのうち、1つのセルの頻度に注目すればよく、この頻度は超幾何 分布に従う。この性質に基づき各時点のセル頻度の観測度数と期待度数からカイ2乗 検定統計量を計算し、これを Mantel-Haenszel検定の仕方（例えば、Mantel（1966）を 参照）で統合したものがログランク検定統計量である。ログランク検定の名称は Peto &Peto（1972）による。試験治療（E）と対照治療（C）の2つの生存率S（t）と SE （t）の比C 較において、ログランク検定における帰無仮説と対立仮説は、片側検定であれば
｛
H0：すべての t についてS（t）=SE （t）C H1：ある t についてS（t）>SE （t）C である。ログランク検定の p値があらかじめ定めた有意水準よりも小さければ、試験 治療の生存率は時間を通して対照治療よりも高いと解釈できる。ログランク検定は、 朝倉・濱﨑（2015b）で紹介した割合の比較と同様に、正規近似に基づく手法である。 シミュレーションに基づく標本サイズと検定の実質的サイズ（シミュレーションによ り計算した第1種の過誤確率）の関係を図2に示す（片側検定のもとで有意水準は2.5% とし、シミュレーションの反復回数は100万回とした）。図2から示唆されるように、 1群あたりの標本サイズが50例よりも小さい場合に、第1種の過誤確率はあらかじめ 定めた有意水準2.5% よりもわずかに大きくなり、2つの生存率に違いがないにもか かわらずその違いを検出する可能性が若干高まる。標本サイズが小さい場合のログラ ンク検定の適用には、十分に注意を払う必要がある。 ログランク検定は、各時点のカイ2乗検定の統合において、すべての重みを等し く1としたものであるが、諸種の重みを考慮し統合することができる。ログランク検 定の代替法の1つとして、朝倉・濱﨑（2015a）で紹介した Wilcoxon の順位和検定を 中途打ち切りデータまでとり扱えるように拡張した一般化Wilcoxon検定（あるいは Gehan検定）（Gehan,1965;Breslow,1970）は、各時点直前までの生存数（Numberof“at Risk”）で重みを付したものである。実地でよく観測される生存率の形状を図3に示す。 図2 シミュレーションに基づく1 群あたりの標本サイズ（等例数）とログランク検定の実質的サイズの関係 実線は検定の実質的サイズを、点線は有意水準を表す。片側検定の有意水準は2.5% とし、シミュレーションの 反復回数は100万回とした。また、シミュレーションでは試験の登録期間を2年、追跡期間を3年とした。 検定の実質的サイズ（％） 1 群の標本サイズ 20 40 50 100 150 2.5 2.0 1.0 0図中の2本の曲線は、2つの群の生存率を表している。この2つの生存率を比較すると き、ログランク検定は（a）の示す生存率の違いを検出することには優れているが、す べての重みを等しく1として統合するため、（b）の生存率の違いを検出することには 不向きである。この場合、各時点直前までの生存数で重み付けし最初の生存率の違い を強調できる一般化Wilcoxon検定がログランク検定よりも優れている。（c）はログラ ンク検定と一般化Wilcoxon検定のいずれでも違いを検出することが難しい。
3
先に紹介した Kaplan-Meier法に基づく生存率、ログランク検定と一般化 Wilcoxon検定を統計計算ソフトウェアで実行する。SAS®UniversityEdition、R、エ クセル統計、JMP®の4つの統計計算ソフトウェアを使用し、R の survival パッケー ジ＊2に組み込まれているテストデータ“ovarian”を解析する（図4）。“ovarian”は卵巣が んに対する外科手術後に異なる2つの化学療法を実施した場合の有効性を検討した確 率化臨床試験のデータである。化学療法の種類は“rx”という変数で得られており、「シ 図3 実地でよく観測される生存率の形状 時間 （a） 時間 （b） 時間 （c） 生存率 1 0 ＊2 生 存 時 間 に 関 す る テ ス ト デ ー タ と 解 析 の た め の 関 数 が ま と め ら れ た パ ッ ケ ー ジ で あ り 、 C R A N（ R 本 体 や パ ッ ケ ー ジ を ダ ウ ン ロ ー ド す る た め の w e b サ イ ト ）か ら ダ ウ ン ロ ー ド し イ ン ス ト ー ル す る こ と で 、“ o v a r i a n ”デ ー タ も 得 ら れ る 。 パ ッ ケ ー ジ の イ ン ス ト ー ル の 方 法 は 例 え ば 舟 尾（ 2 0 0 9 ）を 参 照 さ れ た い 。 図4 テストデータ“ovarian”クロホスファミド単独（rx=1）」もしくは「アドリアマイシンとの併用（rx=2）」である。 治療の有効性を検討する指標として、生存時間または打ち切りまでの時間が“futime”、 打ち切りに関する指標が“fustat”（0：打ち切り、1：死亡）として得られている。いま、 Kaplan-Meier法に基づく生存率の推定値を計算、描画し、2つの化学療法の効果の比 較のために、2群の生存率についてログランク検定または一般化Wilcoxon検定にて比 較する。それぞれのソフトウェアについて SAS®Studio3.4、Rversion3.2.0、エク セル統計2015および JMP®Pro10.0.2にて動作確認を行っている。以降ではそれぞ れのソフトウェアにおける実行と出力される結果について個々に述べる。
① SAS® University Edition：生存時間解析を適用する LIFETIME プロシジャを用い
て Kaplan-Meier法に基づく生存率を描画し、化学療法に関する2群比較のための検 定を実行する手順を紹介する。 PLOTS ステートメントで「PLOTS=S」と指定すれば、Kaplan-Meier法に基づく生 存率の推定値が描画される。このとき、CL オプションを指定（「PLOTS=S（CL）」）す れば、生存率の95%信頼区間が同時に描画される。TIME ステートメントで生存時 間を示す変数（時間変数）、打ち切りの指標および打ち切りを示す値を指定し（ここで は「futime＊fustat（0）」となる）、STRATA ステートメントで群を示す変数（群変数、 ここでは“rx”）を指定する。ただし現時点では、SAS®UniversityEdition の表示言語 が日本語になっている場合に STRATA ステートメントを指定すると、エラーにより 結果が適切に表示されないため、そのような場合には表示言語を英語にする。ここで 示す出力結果は、言語を英語にしてプログラムを実行した場合のものである。 以上の手順を実行すれば、各群について生存率の推定値が図5のように描画され、 ログランク検定と一般化Wilcoxon検定の結果が図6のように出力される。図5、図6 図5 LIFETEST プロシジャによる Kaplan―Meier 法に基づく生存率の推定値のプロット
は、すべての出力のうち例示に用いる部分のみを抜粋している。これより、視覚的に はアドリアマイシンとの併用群がシクロホスファミド単独群に比べて生存率が高いよ うに見えるが、ログランク検定による p値は0.3026であり、有意水準0.05のもとで、 2群の生存率に差があるとはいえないと結論づけられる。また、一般化Wilcoxon検定 による p値は0.1665であり、ログランク検定に基づく p値に比べ小さいが、ログラ ンク検定と一般化Wilcoxon検定で導かれる結論は変わらない。 ② R：ここでは survival パッケージを用いた手順を紹介する。Kaplan-Meier法、ロ グランク検定のいずれにおいても survival パッケージを用いるため、事前にパッ ケージを読み込んでおく。まず survfit関数を用い、Kaplan-Meier法に基づく生存率 の推定値を計算する。その際、survfit関数の引数として指定する式は「Surv（時間変 数,打ち切りの指標）」となる。時間変数としてここでは“futime”、打ち切りの指標と して“fustat”を指定する。さらに群変数として“rx”、およびデータ名として“ovarian” を指定する。これにより群ごとに生存率の推定値が計算される。結果がオブジェク トとして保存されるよう名前を付け（ここでは、“km”とする）実行する。次に plot 関数を用い、先ほど保存したオブジェクトを描画する（「plot（km）」）。これにより、
図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（事象発現）でも可能）。
図8 survdiff 関数によるログランク検定の実行結果（R）
図7 survfit 関数と plot 関数による Kaplan―Meier 法に基づく生存率の推定値のプロット（R）
0 200 400 600 800 1000 1200 1.0 0.6 0.8 0.4 0.2 0.0 図9 survdiff 関数による一般化 Wilcoxon 検定の実行結果（R）
出力された結果より、ログランク検定に基づく p値は SAS®による実行結果に一致 するが、一般化Wilcoxon検定については、SAS®において規定値で適用される検定 とはカイ2乗検定の統合における重みが異なるため、結果が若干異なる。しかしここ では導かれる結論は同じである。 ③エクセル統計：｢生存分析・ハザード分析｣のなかの「カプラン = マイヤー法（表形式）｣ または｢カプラン = マイヤー法（データベース形式）｣を用いる。いずれを用いるかは データの形式によるが、ここでは他のソフトウェアで解析する場合と同じ形式（図4 の形式）でデータが入力されているとし、「カプラン = マイヤー法（データベース形 式）」を用いる。出てくるウィンドウの「変数」タブで「時間」に時間変数“futime”を、「状 態」に打ち切りの指標“fustat”を、「グループ変数」に群変数“rx”を指定する。エクセ ル統計では、打ち切りの指標は0（打ち切り）と1（事象発現）で与える。すると Kaplan-Meier法に基づく生存率の推定値とその描画、検定結果が出力される（図10）。出力 結果より、ログランク検定については“Cochran-Mantel-Haenszel”が SAS®や R にて 前記の手順で実行した検定と対応している。一般化Wilcoxon検定については“Gehan-Breslow”が SAS®での前記の手順で実行した検定と対応しており、それらの結果は 一致している。 図1 0 ｢生存分析・ハザード分析｣による Kaplan-Meier 法に基づく生存率の推定値のプロットとログランク検定 および一般化 Wilcoxon 検定の実行結果（エクセル統計）
図1 1 「生存時間分析｣による Kaplan-Meier 法に基づく生存率の推定値のプロットとログランク検定および一般 化 Wilcoxon 検定の実行結果（JMP®） ④ JMP®：｢分析｣の｢信頼性/生存時間分析｣のなかの「生存時間分析」を用い、「Y，事 象発現までの時間」に時間変数“futime”、「グループ変数」に群変数“rx”、「打ち切り」 に打ち切りの指標“fustat”を指定する。また「打ち切りの値」には打ち切りを表す値（こ こでは0）を指定する。すると、Kaplan-Meier法に基づく生存率の推定値の描画や検 定結果などが出力される（図11）。なお、図11は、出力のうち例示に必要な部分のみ
文献 1）Breslow,N.E.（1970）.AgeneralizedKruskal-WallistestforcomparingKsamplessubjecttounequal patternsofcensorship.Biometrika57,579-594. 2）Gehan,E.A.（1965）.AgeneralizedWilcoxontestforcomparingarbitrarilysingly-censoredsamples. Biometrika52,203-223. 3）Greenwood,M.（1926）.Thenaturaldurationofcancer.ReportsonPublicHealthandMedicalSubjects 33,1-26. 4）Gross,A.J.,Glark,V.A.（1976）.SurvivalDistributions:ReliabilityApplicationsintheBiomedical Sciences.NewYork:JohnWiley&Sons. 5）Hosmer,D.W.,Lemeshow,S.,May,S.（2008）.AppliedSurvivalAnalysis:RegressionModelingofTime-to-EventData.Hoboken:JohnWiley&Sons. 6）Kaplan,E.L.,Meier,P.（1958）.Nonparametricestimationfromincompleteobservations.Journalofthe AmericanStatisticalAssociation53,457-481. 7）Mantel,N.（1966）.Evaluationofsurvivaldataandtwonewrankorderstatisticsarisinginits consideration.CancerChemotherapyReports50,163-170. 8）Peto,R.,Peto,J.（1972）.Asymptoticallyefficientrankinvarianttestprocedure.JournaloftheRoyal StatisticalSociety.SeriesA135,185-207. 9）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）.Designandanalysisofrandomizedclinicaltrialsrequiringprolonged observationofeachpatient.PartII:analysisandexamples.BritishJournalofCancer35,1-39. 10）朝倉こう子・濱﨑俊光（2015a）.2 つの平均の比較 .DrugDeliverySystem30（2）,149-158. 11）朝倉こう子・濱﨑俊光（2015b）.2 つの割合の比較 .DrugDeliverySystem30（3）,235-245. 12）舟尾暢男（2009）.TheRTips:データ解析環境 R の基本技・グラフィック活用集 .第2 版 .オーム社 を抜粋したものである。出力結果より、ログランク検定の結果は SAS®、R、エクセ ル統計による実行結果と一致し、一般化Wilcoxon検定の結果は SAS®による実行結 果（すなわちエクセル統計では“Gehan-Breslow”）と一致している。
