R を使った生存時間分析:
Kaplan-Meier の作り方とログランク検定、論文への記載方法など

UB3/statistics/survival/survival_r

このページの最終更新日: 2021/07/10

  1. R を使った Log-rank 検定 (自作データセットを使用)
  2. R を使った生存曲線の作成 (自作データセットを使用)
  3. 組み込みデータセットを使った Log-rank 検定

広告

R を使った Log-rank test

統計パッケージ R を使い、実際に Log-rank test を行う方法を解説する。データフレームの作り方については、テキストファイルからのデータ読み込み を参照のこと。csv で以下のデータテーブルのような形で保存し、read.csv でそのファイルを読み込むのが手軽だと思う。

R で Log-rank test を行うには、少なくとも以下のデータが dataframe に含まれる必要がある。

  • 生存時間
  • 打ち切りの有無
  • 群 (上に書いたように 2 群が基本)

簡略化のため、3 個体 x 2 群でデータフレームの例を示しておく。念のため、作り方は次の通り。

R 生存時間分析

データテーブルは以下の通り。

group time censor
1 A 10 1
2 A 11 1
3 A 10 1
4 B 19 1
5 B 25 1
6 B 23 0 0 が打ち切りを示す

グループ B は A の約 2 倍長生きで、6 番の個体のみ打ち切りであったというデータである。このデータの並べ方は重要。この他、性別や年齢など関連するデータがあるなら、列を追加してもよい (2)。

このようなデータフレーム test がある状態で、以下のように survival パッケージを読み込み、ログランク検定をかける。

生存時間分析では、まず Surv 関数を使い、これを survdiff 関数で解析するという方法をとる。

library(survival)
survdiff(Surv(time, censor) ~ group, data = test)

以下の例では、いったん survdiff の結果を test_stat に格納し、test_stat と入力することで結果を表示させている。

R 生存時間分析

なお、上の結果からわかるように、得られる検定統計量は chi square value である。統計検定の英語表現 のページにまとめたように、上記の Log-rank 検定の結果を paper に記載するときの正式な表記法は、Log-rank test, chi-square = 5.2, p = 0.02 のようになる。

また、coxph にすると Cox ハザードモデルでの解析が行える。

library(survival)
coxph(Surv(time, censor) ~ group, data = test)

R を使った生存曲線の作成

生存曲線

survdiff でなく survfit という関数を使ってこれをプロットすれば、Kaplan-Meier カーブを作ることができる。

test_fit = survfit(Surv(time,censor)~group, data=test) plot(test_fit)

描画データは test_fit$surv というカラムに保存されているので、これを表示すれば最終的な生存率が何%であったかを知ることができる。


広告

組み込みデータセットを使った Log-rank 検定

このページ によると、R で生存時間分析に適したデータセットは以下の 3 つである。

  • MASS パッケージに含まれている gehan
  • survival パッケージに含まれている leukemia
  • survival パッケージに含まれている veteran

ここでは、veteran パッケージを使った例を示しておく。これは肺がん治療のデータを示したデータセットで、以下のような項目を含んでいる (2)。本当は 137 行あるが、1 から 10 行目のみ記載する。

Rのveteranデータセット
  • trt: 治療の種類。1 が standard (1-69 行)、2 が test (70-137 行)。
  • celltype: がん細胞の種類。
  • time: 生存時間。
  • status: 1 が打ち切りなし、0 が打ち切り。
  • karno: Karnofsky performance score
  • diagtime: 診断されてから clinical trial までの時間
  • age: 年齢
  • prior: Clinical trial 前の治療の有無。0 が治療なし、10 が治療あり。

trt, time, status を使って解析する。

library(survival)
survdiff(Surv(time, status) ~ trt, data = veteran)

Rのveteranデータセット Log-rank検定結果

p = 0.9 なので、有意差なしという結果になる。グラフは以下の通り。線の色を変更するなど、グラフの見た目については plot 関数のページ を参照のこと。

test_fit = survfit(Surv(time, status) ~ trt, data = veteran) plot(test_fit)


Rのveteranデータセット 生存曲線
広告

References

  1. 佐藤弘樹、市川度. 2013.

生存時間解析 について平易に書いた数少ない解説書。

統計のなかでも、生存時間解析はそれだけで 1 冊の本になるほど複雑なわりに、ANOVAや t 検定などと違い使用頻度が低いため、とっつきにくい検定である。

この本では、とくに Kalpan-Meier 生存曲線、Log-rank 検定、Cox 比例ハザードモデルを重点的に解説しているが、prospective study と retrospective study, 選択バイアス、プラセボなど、臨床統計実験で重要な概念についても詳しい説明がある。臨床でない、基礎生物学の実験ではあまり意識しない重要な点であるので押さえておきたい。

  1. R で生存時間分析を行う. Link: Last access 2018/09/20.
  2. Veterans' Administration Lung Cancer study. Link: Last access 2021/07/05.

コメント欄

各ページのコメント欄を復活させました。スパム対策のため、以下の禁止ワードが含まれるコメントは表示されないように設定しています。レイアウトなどは引き続き改善していきます。「管理人への質問」「フォーラム」へのバナーも引き続きご利用下さい。

禁止ワード: http, the, м (ロシア語のフォントです)

このページにコメント

Name:


Comment:



これまでに投稿されたコメント

Date Name Comment

↑