R言語

【R言語】生存時間解析survival カプランマイヤー曲線・Cox比例ハザードモデル・ログランク検定

  1. HOME >
  2. R言語 >

【R言語】生存時間解析survival カプランマイヤー曲線・Cox比例ハザードモデル・ログランク検定

スポンサーリンク

R言語で生存時間解析を行う方法を紹介していきます。

生存率の推定方法の一種であるカプランマイヤー曲線のプロットやCoxの比例ハザードモデル、ログランク検定

を解説していきます。

今回扱うプログラミングコードは以下からダウンロードできます。

Rでデータ解析を始めるならコレ
¥3,300 (2022/06/01 17:20時点 | Amazon調べ)
\楽天ポイント5倍セール!/
楽天市場
RStudioで解析したい方へ
¥3,278 (2022/06/02 22:52時点 | Amazon調べ)
\楽天ポイント5倍セール!/
楽天市場

生存時間解析

R言語で生存時間解析を行うためにはパッケージsurvivalをインストールする必要があります。

重回帰分析などのように標準で関数がインストールされているわけではないので注意が必要です。

早速、以下のコードを実行しsurvivalをインストールします。

パッケージsurvivalにはcancerという癌に関するデータセットが含めれています。時間(time)イベント発生の有無(status)などの変数があり、まさに生存時間解析の練習用のデータセットといった感じです。

このデータセットは既にデータ整形されているため、少し加工し実際のデータに近い形にします。(生存時間や生存と死亡のフラグが既にnumeric型として扱えますが、実際はcharacter型だったり、日付のままだったりするので解説のために加工していきます。)

以下を実行すると、cancerに観測開始日時と終了日時が追加されたデータセットを作成することができます。

以下のように、生存時間ではなく開始と終了に関する日時が追加されているのが分かります。また、生存・死亡のフラグであるstatusも0, 1から"survived", "dead"に変更しました。

多分、実務でデータ入力する際に開始日と終了日の差を考慮しないことがほとんどなので、上のようなデータ日付データが一般的だと思います。

生存時間のプロット

まず、生存時間のグラフの描き方についてみていきます。

生存時間解析の本でよく見かける観測開始時点と終了時点に関するグラフを描いていきます。

生存時間を取得する際に用いるDateクラスは次の記事で紹介しています。基本的な日付・時間に関する操作をまとめています。

R言語 日付・時間 DateクラスとPOSIXクラス

R言語の日付と時間に関してよく使うテクニックをまとめました。 Rに標準でインストールされているパッケージbaseの中のDateクラスとPOSIXクラスの操作方法について解説していきます。 日付や時間か ...

続きを見る

まず、パッケージggplot2をインストールします。Rのbaseに標準で入っているplotを使うとx軸とy軸を反転できないため、ggplot2を使って描きます。

観測開始時点と終了時点に関するグラフは次のようにしてプロットすることができます。変数statusによって色を分けるように指定しているので、死亡("dead")が青打ち切り("survived")が赤に色分けされています。

生存時間プロット

また、グラフのy軸のstartを0に設定すると次のような左端を固定したグラフを描くこともできます。

生存時間プロット(左端固定)

カプランマイヤー曲線

カプランマイヤー法を用いて生存率を推定していきます。

カプランマイヤー法を実行するにはパッケージsurviveの関数survfitを用います。

関数survfit

関数survfitとその引数を以下にまとめます。

survfit(formula, data, weights, subset, na.action, newdata, individual=F, conf.int=.95, se.fit=T, type=c("kaplan-meier","fleming-harrington", "fh2"), error=c("greenwood","tsiatis"), conf.type=c("log","log-log","plain","none"), conf.lower=c("usual", "peto", "modified"))
関数survfitの引数
引数説明
formulaSurvオブジェクトまたはcoxphオブジェクト。他のformulaと同様に目的変数を~の左側に、説明変数を右側に指定する。
data引数formula、subset、weightsで列名を指定する際にどの列なのかを解釈するためのデータフレーム
weights非負の重み。正で与えることを推奨している。
subsetモデルを適合する際にデータのどの行を用いるかを指定する。
na.action欠損値に対して実行する関数。デフォルトでoptions()$na.actionを行う。
newdataformulaで指定した変数と同じ変数を持つデータフレーム。formulaがcoxphオブジェクトである場合に適用できる。
individuallogical型。データの行がぞれぞれ異なる個々の時間を表す場合TRUE、個々が複数の時間を持つ場合はFALSE。
conf.int信頼区間の信頼水準を指定する。
se.fitlogical型。標準誤差を計算するかどうか。
typecharacter型。どの手法の生存曲線を求めるかどうかを指定する。"kaplan-meier"でカプランマイヤー法、"fleming-harrington"でfleming-harrington法、"fh2"でfh2法。
error誤差の計算方法。"greenwood"でGreenwoodの式、"tsiatis"でTsiatisの式によって計算する。
conf.type信頼区間の計算方法。"none"で信頼区間を計算しない。デフォルトで"log"。
conf.lower信頼区間の下限の設定。

実行例

関数survfitを使ってカプランマイヤー曲線をプロットするコードを紹介します。

まず、さきほど作ったdatasetに生存時間の変数を追加します。観察開始日と終了日の変数は日付を表すcharacter型となっていますが、as.Dateを使ってDate型に変換し、さらにas.numericによってnumeric型に変換することで数値に変換することができます。

datasetの一番右の列にtimeが追加されているのが確認できます。

関数survfitを使ってstatus別の生存曲線を描く前に、変数statussexnumeric型に変換する必要があります。

注意ポイント

関数survfitの引数formulaにSurv(time, event)~sexのようにSurvオブジェクトを指定しますが、time、event、sexに関する変数がnumeric型ではない場合、正しく生存曲線が描画されないので注意が必要です。生存曲線ではなく死亡曲線が描かれたりします。

変数statussexの整形が完了したら、survfit実行してみましょう。次のように記述することでカプランマイヤー法で生存率を推定することができます。

他の解析法と同様にfitには生存率やその信頼区間などの変数が含まれており、fitの後ろに$を付けることで参照することができます。

データフレームにこれらの値を代入することで次の画像のように、csvファイルに生存率とその信頼区間を出力することができます。

カプランマイヤー推定

カプランマイヤー生存曲線をプロットする際は、survivalパッケージに含まれている関数plot.survfitを使います。実際に使用する際にはplot()のように使うため、一見、Rのbaseに標準で実装されている関数plotだと思ってしまいますが、こちらの関数は第一引数にSurvオブジェクトを指定するため全く別物となります。

下記のコードを実行することで、通常のカプランマイヤー曲線とsex別のカプランマイヤー曲線が描けます。

1行目のように引数formulaをSurv(time, status) ~ 1と指定することで、グループ別に分けることなく生存曲線を描くことが可能です。factor型の変数のlevelごとに曲線を描きたいときはSurv(time, status) ~ sexのようにします。ここでは性別ごとの生存曲線を描いています。

プロットに打ち切りの線を引きたい場合はmark.T = TRUEとし、信頼区間を描きたい場合はconf.int = TRUEとします。(画像参照)

カプランマイヤー生存曲線
カプランマイヤー生存曲線(Sex別)
¥3,509 (2022/06/11 02:16時点 | Amazon調べ)
\楽天ポイント5倍セール!/
楽天市場
¥4,036 (2022/06/11 02:41時点 | Amazon調べ)
\楽天ポイント5倍セール!/
楽天市場

Coxの比例ハザードモデル

続いてCoxの比例ハザードモデルの実行方法について解説します。

Coxの比例ハザードモデルを実行するには関数coxphを用います。

関数coxph

以下、関数coxphとその引数です。

coxph(formula, data=, weights, subset, na.action, init, control, ties=c("efron","breslow","exact"), singular.ok=TRUE, robust, model=FALSE, x=FALSE, y=TRUE, tt, method=ties, id, cluster, istate, statedata, nocenter=c(-1, 0, 1), ...)
関数coxphの引数
引数説明
formulaSurvオブジェクト。他のformulaと同様に目的変数を~の左側に、説明変数を右側に指定する。
data引数formula、subset、weightsで列名を指定する際にどの列なのかを解釈するためのデータフレーム。
weightsnumeric型のベクトル。モデリングの際の重み。
subsetモデルを適合する際にデータのどの行を用いるかを指定する。
na.action欠損値に対して実行する関数。デフォルトでoptions()$na.actionを行う。
initnumeric型のベクトル。回帰係数を求める際の反復法の初期値。
controlcox.controlオブジェクト。反復法の際の極限やそのほかのオプションを指定する。
tie生存時間にtieが存在するときにどのように部分尤度を計算するかどうか。デフォルトはEfronの近似。 "Breslow"でBreslowの近似。"exact"で正確な部分尤度。
singular.oklogical型。計画行列が共線性を持つ場合にどうするか。TRUEの場合、行列の特定の行が前の行の線形結合であるとき、自動的に無視する。
robustlogical型。ロバスト分散を計算するかどうか。
idどの行が同じ対象であるかを指定する。1つの対象が2つ以上のデータを持つ場合指定する。
cluster観測値をクラスタリングするためのオプション値。
istate状態の初期値。
statedata多状態モデルを表現する際の別のデータ。
modellogical型。TRUEの場合、モデルを要素ごとに返す。
xlogical型。TRUEの場合、計画行列Xを要素ごとに返す。
ylogical型。TRUEの場合、応答変数Yを要素ごとに返す。
tt時間を変換する関数のリスト。
metod引数tieの別の名前。
nocenter行列Xの列が指定したベクトルの範囲に存在する場合、センタリングしない。

実行例

生存率をセンサーデータstatusその他の説明変数でモデリングしていきます。

  • センサーデータ:status
  • 説明変数:sex, age, ph.ecog, meal.cal, wt.loss

次のコードを実行することで、上記の変数のCoxの比例ハザードモデルを行うことができます。

fitを参照すると次のように回帰係数p値を見ることができます。Likelihood ratio test=9.33 on 5 df, p=0.09643と書いてあることから、timeはsex, age, ph.ecog, meal.cal, wt.lossで変化するとは言えないことが分かりました。また、各説明変数の検定結果からph.ecogのp値が0.0506<0.1であり、timeに影響を及ぼす変数であることがうかがえます。

Coxの比例ハザードモデルの結果を抽出するときは、関数summaryを使うのをお勧めします。

以下のようにsummaryを用いることで、先ほどの回帰係数やp値、各種検定結果等を抽出することが可能です。

モデルや検定結果などをデータフレームに格納し、csvファイルに出力するには以下のように実行しましょう。

上を実行すると、次の画像のcsvファイルが作業ディレクトリに出力されます。

Coxの比例ハザードモデル

ログランク検定

最後にログランク検定についてみていきます。

パッケージsurvivalの関数survdiffを用いることで、ログランク検定を実行することができます。

関数survdiff

以下、関数survdiffとその引数となります。

survdiff(formula, data, subset, na.action, rho=0, timefix=TRUE)
関数survdiffの引数
引数説明
formulaSurvオブジェクト。Surv(time, status) ~ predictorsのように記述する。
data引数formula、subset、weightsで列名を指定する際にどの列なのかを解釈するためのデータフレーム
subsetnumeric型のベクトル。モデルを適合する際にデータのどの行を用いるかを指定する。マイナスの場合そのインデックスのデータを除外する。
na.action欠損値に対して実行する関数。デフォルトでoptions()$na.actionを行う。
rhonumeric型。検定の種類を指定する。
timefixlogical型。関数aeqSurvが丸め誤差を排除する際のい時間を有限にするかどうか。

実行例

ログランク検定を用いて性別(sex)によって生存時間(time)が異なるかどうかを検定してみます。

次のように関数survdiff第一引数にSurv(time, status) ~ sexを代入することで、生存時間が性別によって異なるかというログランク検定を実行することができます。

注意ポイント

注意として、Survオブジェクトeventに関するデータをnumeric型にするのを忘れないようにしましょう。

testResultを参照すると、観測度数や期待度数、p値などのカイ2乗検定に関する値がコンソールに表示されます。また、Chisq= 2.2 on 1 degrees of freedom, p= 0.1 より有意水準0.05で生存時間(time)は性別(sex)によって変化するとは言えないことが分かりました。

testResultの後ろに$を付けることで、これらの検定結果の値を取得することができます。

Coxの比例ハザードモデルでもやったように、データフレームに代入することで検定結果をcsvファイルに出力することができます。

ログランク検定

まとめ

R言語で生存時間解析を解説しました。

この記事ではカプランマイヤー曲線や、Coxの比例ハザードモデル、ログランク検定など有名な解析を紹介しました。

今回、関数の引数などの詳細はあまり触れませんでしたが、詳細についてはRdocumentationで調べるのをお勧めします。

スポンサーリンク

  • この記事を書いた人
  • 最新記事

usagi-san

統計学とゲームとかをメインに解説していくよ。 数式とかプログラミングコードにミスがあったり質問があったりする場合はコメントで受け付けます。すぐに対応します。

-R言語
-, , ,

© 2022 ウサギさんの統計学サロン Powered by AFFINGER5