ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rでsurvival】Kaplan-Meier曲線推定とログランク検定

 今回は、Rのsurvivalパッケージを用いて、Kaplan-Meier曲線の推定ログランク検定を行う方法について説明します。生存時間解析は、たとえばある集団において「ある時点まで生存している確率」を求めたり、異なるグループ間で生存状況に差があるかどうかを検定する手法です。

2つグループの生存関数の比較イメージ

1. ちょっとだけ基本

Kaplan-Meier曲線(Kaplan-Meier Curve)

  • 目的:
    ある集団で、時間の経過とともに死亡などのイベントが発生していく状況を考えます。Kaplan-Meier曲線は、時間 \displaystyle{
t
} まで生存している確率、すなわち生存関数 \displaystyle{
S(t)
} を推定する方法です。

  • 推定方法:
    観察された各イベント発生時刻 \displaystyle{
t _ j
} について、その時点直前に生存していた人数 \displaystyle{
n _ j
} と、その時刻で実際にイベント (死亡など)が起こった人数 \displaystyle{
d _ j
} を使い、以下の式で生存確率を求めます。

\displaystyle{
S(t) = \prod_{t_j \leq t} \left(1 - \frac{d_j}{n_j}\right)
}

 この計算方法により得られた生存曲線は、時間の経過とともに段階的(階段状)に下がっていきます。

 以下の記事も参考にしてください。

chaosmemo.com

ログランク検定(Log-Rank Test)

  • 目的:
    2つ以上のグループの生存曲線に、統計的に有意な差があるかどうかを調べます。

  • 検定の考え方:
    各時刻ごとに、各グループで実際に観察されたイベント数と、もし全グループの生存曲線が同じであったとした場合に期待されるイベント数との差を計算し、それを累積して検定統計量を求めます。帰無仮説は「グループ間で生存曲線に差がない」というもので、統計量はカイ二乗分布に従います。

2. Rスクリプトの例

以下は、シミュレーションデータ(コンピュータで生成した架空の生存時間データ)を用いてKaplan-Meier曲線の推定とログランク検定を行うRスクリプトです。ここでは、サンプルデータとして生存時間(指数分布に従う)、イベントの有無(打ち切りかイベント発生か)、および2つのグループ情報を使っています。

# 1. パッケージの読み込み
# survivalパッケージがインストールされていない場合は自動でインストールし、読み込みます。
if (!require("survival", quietly = TRUE)) {
  install.packages("survival", dependencies = TRUE)
  library(survival)
}

# 2. シミュレーションデータの作成
# 生存時間(time)は指数分布、イベントの有無(status: 1=イベント発生, 0=打ち切り)はランダム、グループ(group)は2群に分けています。
n <- 100
time <- rexp(n, rate = 0.1)                # 生存時間(指数分布)
status <- sample(0:1, n, replace = TRUE)   # イベントの有無
group <- sample(c("Group A", "Group B"), n, replace = TRUE)  # グループ分け

# データフレームにまとめる
data <- data.frame(time = time, status = status, group = group)

# 3. Kaplan-Meier曲線の推定
# Surv()関数で生存オブジェクトを作成し、survfit()関数でグループごとに生存関数を推定します。
km_fit <- survfit(Surv(time, status) ~ group, data = data)

# 推定結果の表示
print(km_fit)

# Kaplan-Meier曲線のプロット
plot(km_fit, col = c("red", "blue"), lwd = 2,
     xlab = "Time", ylab = "Survival Probability",
     main = "Kaplan-Meier Curve by Group")
legend("topright", legend = levels(factor(group)), col = c("red", "blue"), lty = 1, lwd = 2)

# 4. ログランク検定の実施
# survdiff()関数を使ってグループ間の生存曲線の差を検定します。
logrank_test <- survdiff(Surv(time, status) ~ group, data = data)
print(logrank_test)

3. Rスクリプトの解説

3.1 パッケージの読み込みとデータ準備

  • パッケージ読み込み:
    library(survival)を実行することで、サバイバル解析に必要な関数(Survsurvfitsurvdiffなど)が使えるようになります。
    このスクリプトでは、もしパッケージがインストールされていなければ、自動でインストールする仕組みも取り入れています。

  • シミュレーションデータ:
    コンピュータで作った架空のデータを使っています。

    • time: 生存時間(ここでは指数分布に従う乱数で生成)。
    • status: イベントが起こったかどうか(1ならイベント発生、0なら打ち切り)。
    • group: 2つのグループに分けた情報。

3.2 Kaplan-Meier曲線の推定

  • 生存オブジェクトの作成:
    Surv(time, status)によって、生存時間とイベント情報をひとまとめにしたオブジェクトを作ります。

  • 生存関数の推定:
    survfit(Surv(time, status) ~ group, data = data)により、各グループごとに生存関数(生存確率の推移)を推定します。

  • プロット:
    推定された各グループの生存曲線を色分けして描画し、どのグループの曲線かを分かりやすく示します。

3.3 ログランク検定の実施

  • 検定:
    survdiff(Surv(time, status) ~ group, data = data)を使い、各グループで観察されたイベント数と期待されるイベント数を比較して、統計的に有意な差があるかどうかを調べます。
    結果には検定統計量(カイ二乗値)やp値が表示され、p値が小さい場合は「グループ間に有意な差がある」と判断されます。

  • 注意点:
    このサンプルデータでは、2つのグループ間に本来の違いは設定していません。そのため、真の生存曲線はどちらも同じであり、見かけ上の差はランダムなばらつきによるものです。

4. survfit関数とsurvdiff関数の使い方・オプションについて

survfit関数の使い方

survfit関数は、Kaplan-Meier曲線などの生存曲線を推定するために使います。使い方のポイントは以下の通りです。

  • 基本の構文:

    • 単一群の場合:
      r fit <- survfit(Surv(time, status) ~ 1, data = mydata) ここで~ 1はグループ分けをせずに全体の生存曲線を求めることを意味します。

    • グループごとに解析する場合:
      r fit_group <- survfit(Surv(time, status) ~ group, data = mydata) 右辺にグループ変数を指定することで、各グループごとの生存曲線が得られます。

  • 主なオプション:

    • data: 使用するデータフレームを指定します。
    • subset: 解析対象の部分集合を指定できます(例:subset = (group == "A"))。
    • na.action: 欠損値の扱い方を指定します(例:na.action = na.exclude)。
    • conf.int: 信頼区間の信頼水準を指定します(デフォルトは0.95ですが、conf.int = 0.90なども可能)。
    • conf.type: 信頼区間の計算方法を指定します。たとえば、"log"(対数変換、デフォルト)、"plain"(変換なし)、"log-log"(対数-対数変換)などがあります。
  • 出力内容:
    実行結果は"survfit"クラスのオブジェクトで、以下の情報が含まれます。

    • time: イベントが発生した時刻。
    • n.risk: 各時刻のリスクにさらされた人数。
    • n.event: 各時刻のイベント発生数。
    • n.censor: 打ち切りとなった人数。
    • surv: 推定された生存確率。
    • lower, upper: 信頼区間の下限・上限。
    • 複数グループの場合、strataとして各グループの情報がまとめられます。
      summary()関数を使えば、より詳細な各時点の情報を確認できます。

survdiff関数の使い方

survdiff関数は、複数グループ間で生存曲線に有意な差があるかどうかを検定するためのものです。

  • 基本の構文:
    r result <- survdiff(Surv(time, status) ~ group, data = mydata) ここで、Surv(time, status)は生存時間とイベント情報、~ groupはグループ分けを意味します。

  • 主なオプション:

    • data: 使用するデータフレームを指定します。
    • subset: 解析対象のデータを絞る場合に使用します。
    • na.action: 欠損値の扱いを指定します。
    • rho: 重み付けパラメータで、デフォルトはrho = 0(標準的なログランク検定)です。rhoの値を変えると、検定が特定の時点(例:早期のイベント)に重点を置くようになります。
  • 出力内容:
    survdiffの結果は、以下の情報を含むオブジェクトとして返されます。

    • chisq: 検定統計量(カイ二乗値)。
    • n: 各グループの症例数。
    • obs: 各グループで実際に観察されたイベント数。
    • exp: 各グループで期待されるイベント数(帰無仮説に基づく)。
    • var: 検定統計量の分散。
      これらの値からp値を計算し、通常は有意水準(例えば0.05)と比較してグループ間に差があるかどうかを判断します。

ログランク検定の結果の読み方

 ログランク検定の結果は、グループ間で生存曲線に有意な差があるかどうかを判断するための統計的な情報を提供します。分析結果の読み方とその解釈について簡単に説明します。

1. 出力結果の主な要素

検定統計量(chisq)

  • 意味:
    これは、各グループで観測されたイベント数と、もしグループ間に差がなかった場合に期待されるイベント数との差を、各時点ごとに累積して計算した統計量です。
  • 解釈:
    検定統計量が大きいほど、グループ間の差が大きいことを示唆します。たとえば、グループ間で実際のイベント数と期待されるイベント数のズレが大きければ、chisq値も大きくなります。

自由度(df)

  • 意味:
    自由度は「グループの数 - 1」で決まります。たとえば、2群の場合は自由度は1、3群の場合は自由度は2となります。
  • 解釈:
    この自由度をもとに、カイ二乗分布を参照してp値を計算します。

p値

  • 意味:
    検定統計量(chisq値)が、帰無仮説「グループ間に生存曲線の差がない」のもとで観察される確率です。
  • 解釈:
    • p値があらかじめ設定した有意水準(例えば0.05)より小さい場合、グループ間の生存曲線に有意な差があると判断します。
    • 逆に、p値が有意水準以上の場合、統計的には「差があるとは言えない」と結論づけます。

各グループのObserved(obs)とExpected(exp)の値

  • Observed(obs):
    各グループで実際に観察されたイベント数です。
  • Expected(exp):
    仮にすべてのグループで生存曲線が同じだった場合に、各グループで予測されるイベント数です。
  • 解釈:
    各グループのobsとexpの差が大きいほど、帰無仮説(グループ間に差がない)が成立しにくいと考えられます。この差の大きさは、chisq値の計算にも反映されています。

2. 結果の読み方

 例えば、次のような出力があったとします。

Call: survdiff(formula = Surv(time, status) ~ group, data = data)

             N Observed Expected (O-E)^2/E
group=A     50       12     10.5      0.214
group=B     50        8      9.5      0.211

Chisq= 0.43  on 1 degrees of freedom, p= 0.512

この出力から読み取れることは:

  • グループごとのデータ:

    • グループA:50人中実際のイベント数は12、期待値は10.5
    • グループB:50人中実際のイベント数は8、期待値は9.5
      それぞれのグループで観測値と期待値の差は小さく、(O-E)2/Eの値も低いです。
  • 検定統計量と自由度:

    • 検定統計量(chisq)は0.43で、自由度は1(2群の場合)。
  • p値:

    • p値は0.512となっており、これは一般的な有意水準(0.05)よりはるかに大きいため、「グループ間で有意な差は認められない」という結論になります。

3. 結果の解釈

  • 有意な差がある場合:
    有意水準を0.05とするのであれば、p値が0.05未満であれば、グループ間の生存曲線に統計的に有意な差があると判断します。たとえば、p値が0.02なら、帰無仮説(「グループ間に差がない」)を棄却し、グループ間で生存のパターンに明らかな違いがあると結論づけます。

  • 有意な差がない場合:
    有意水準を0.05としたとき、p値が0.05以上の場合、グループ間の生存曲線に有意な差がない、またはデータのばらつきが偶然の範囲内であると判断します。たとえば、上の例ではp値が0.512であり、これは「差があるとは言えない」と判断します。

ログランク検定の結果を読むときは、以下の点に注意してください。

  • 検定統計量(chisq): 差の大きさを示す数値。大きいほどグループ間の差がある可能性が高い。
  • 自由度: グループ数に依存し、検定統計量を評価するために必要。
  • p値: 統計的に有意かどうかの判断基準。有意水準(通常0.05)と比較して判断する。
  • ObservedとExpectedの値: 各グループでの実際のイベント数と、仮に差がなかった場合の期待値を比較することで、どの程度の差があるかを把握する。

このようにして、ログランク検定の結果から、グループ間の生存パターンの違いが偶然か否かを判断します。

5. まとめ

  • Kaplan-Meier曲線:
    生存時間データをもとに、ある時点までの生存確率 \displaystyle{
S(t)
} を推定する方法です。イベントが起こるたびに生存率が階段状に下がっていきます。

  • ログランク検定:
    異なるグループ間で生存曲線に差があるかどうかを、各時点での実際のイベント数と期待されるイベント数の差から検定する方法です。

  • Rでの実際:
    シミュレーションデータを用いて、survivalパッケージのsurvfitsurvdiffを使うことで、生存曲線の推定とグループ間の差の検定が簡単に実施できます。
    また、survfitsurvdiffにはさまざまなオプションがあり、データの特性や解析目的に合わせて柔軟に設定できます。

このように、Rのsurvivalパッケージを利用することで、実際の臨床データや観察データに基づく生存解析を簡単かつ詳細に行うことができます。