Rで生存時間解析　ーコックス比例ハザードモデル, coxph関数ー

今回は、生存時間解析として広く使用されるコックス比例ハザードモデルについて解説します。
以前、生存時間解析として、カプラン・マイヤー曲線やログランク検定などの手法を紹介しました。

Rで生存時間解析　ーカプラン・マイヤー曲線,survfit関数と、ログランク検定,survdiff関数ー

臨床現場で白血病患者の治療をしない群（コントロール群）と治療群の生存時間のデータを取得している。取得データから、各群の白血病患者の生存時間を予測したい。また、コントロール群と治療群の2群のデータがあり、それらの2群の比較も合わせて行...

■目次
1. コックス比例ハザードモデルとは（カプラン・マイヤー曲線との違い）
2. Rでコックス比例ハザードモデル
3. Rによるデータの可視化

1. コックス比例ハザードモデルとは（カプラン・マイヤー曲線との違い）

生存時間解析とは、生物の死や機械システムの故障など、1つの事象（event、イベント）が発生するまでの予想される期間を分析する統計学の一分野である。
引用：Wikipedia

生存時間解析の手法として、コックス比例ハザードモデルは広く使われます。カプラン・マイヤー曲線との違いについて解説します。

カプラン・マイヤー曲線(Kaplan-Meier method)とは？　
・目的: 生存率の推定と視覚化
・概要: 生存データ（例えば、患者の生存時間）に基づき、生存関数を推定するための非パラメトリックな手法。時間の経過に伴う生存確率を示す曲線を描く。
・使用例: 治療後の生存率を示すために使用される。例えば、癌患者の生存期間を、治療方法のちがい（単変量）で比較するなど。
・出力: 時間とともに生存する確率の推移を示すグラフ。

💡カプランマイヤー曲線は、単一の因子（変数）が生存時間に与える影響を評価するのに使われる。

コックス比例ハザードモデル（Cox Proportional-Hazards Model）とは？
・目的: 生存データにおける共変量の影響評価
・概要: 生存時間に対する共変量（例えば、年齢、性別、治療方法など）の影響を調べるための回帰モデル。ハザード比を用いて、各共変量が生存時間に与える影響を定量化する。
・使用例: 複数の要因が生存時間にどのように影響するかを調べる際に使用される。
・出力: 共変量の影響を示すハザード比（例えば、1より大きい場合はリスクが増加、1より小さい場合はリスクが減少することを示す）。

💡コックス比例ハザードモデルは、複数因子（変数）が生存時間に与える影響を評価するのに使われる。

★ログランク検定について
・目的：生存曲線を評価するための検定。
・概要：カプラン・マイヤー曲線で得られる生存曲線の差異を検定し、異なる治療法や条件が生存率に与える影響を確認する。具体的には、2つ以上のグループ間で生存曲線が統計的に有意に異なるかを示すp値を算出する。
コックス比例ハザードモデルの予測と実際の生存曲線の一致を確認するために、ログランク検定が参考にされる。p値が有意水準を下まわる場合、モデルの予測がデータと適合しており、各変数が統計学的有意にモデルに寄与していると判断される。

💡ログランク検定は、カプラン・マイヤー曲線でも、コックス比例ハザードモデルでも評価に使われる。

2. Rでコックス比例ハザードモデル

今回は、肺がん患者のデータセットにコックス比例ハザードモデルを使用します。

survival パッケージに入っているlung（進行肺がん患者の生存率）のデータセットを使用します。
　inst：機関コード、time：生存期間（日）、status: 打ち切りステータス 1=打ち切り, 2=死亡、
　age：年齢（歳）、sex: 性別：男性=1 女性=2、
　ph.ecog： ECOGパフォーマンスステータス（身体機能レべルの指標　医師の診断）0=無症状、1=症状はあるが完全に歩行可能、2=1日の50%未満はベッドにいる、3=1日の50%以上はベッドにいるが寝たきりではない、4=寝たきり、
　ph.karno： KSPスコア（身体機能レベルの指標　医師の診断）（悪い＝0-良い＝100）、
　pat.karno： KSPスコア（身体機能レベルの指標　患者の診断）、meal.cal：食事時の消費カロリー、
　wt.loss：過去6ヶ月間の体重減少（ポンド）

コックス比例ハザードモデルを適用し、年齢（age）、性別（sex）、およびECOGパフォーマンスステータス（ph.ecog）が生存期間（time）に与える影響を解析します。

install.packages("survival")
library(survival)

head(lung) 　#進行肺癌患者の生存率　データセット
#  inst 　time 　status 　age　 sex 　ph.ecog　 ph.karno　 pat.karno 　meal.cal　 wt.loss
#1    3  　306    　  2  　 74 　　  1    　   1    　　   90      　　 100     　　1175   　    NA
#2    3  　455      　2  　 68 　　  1    　   0    　　   90    　　    90  　　    1225    　   15
#3    3 　1010      　1  　56  　　 1    　   0   　　    90    　　    90    　 　   NA     　   15
#4    5  　210      　2  　57　　   1      　 1    　　   90      　　  60   　　    1150           11
#5    1  　883      　2  　60  　　 1      　 0    　　  100     　　   90  　 　     NA       　  0
#6   12   1022      　1  　74 　　  1    　   1   　　    50       　　 80    　　    513            0

head(Surv(lung$time, lung$status))　#　コックス比例ハザードモデルでも打ち切りデータを使用
[1] 306 455 1010+ 210 883 1022+
cox_model <- coxph(Surv(time, status)~ age + sex + ph.ecog, data = lung)　#コックス比例ハザードモデル
summary(cox_model)
#Call:
#coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)
#  n= 227, number of events= 164 　　# n=モデルに使用された総サンプル数、number of events=観測されたイベント（死亡）数
#   ( 1 個の観測値が欠損のため削除されました )
#            　 　  coef　 exp(coef) 　 se(coef)      　z 　　　Pr(>|z|)    　# coef：回帰係数、exp(coef)：ハザード比、 se(coef)：回帰係数の標準誤差
#age  　    0.011067   1.011128 　 0.009267  　1.194　  0.232416    
#sex  　   -0.552612   0.575445 　 0.167739　 -3.294 　0.000986 ***
#ph.ecog  0.463728   1.589991  　0.113577  　 4.083 　4.45e-05 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#       　　 exp(coef) 　exp(-coef)　 lower .95　 upper .95　　# ハザード比の95%信頼区間　※１
#age　        1.0111 　    0.9890   　　 0.9929   　 1.0297
#sex  　      0.5754 　    1.7378   　　 0.4142   　 0.7994
#ph.ecog    1.5900    　 0.6289  　　  1.2727    　1.9864
#Concordance= 0.637  (se = 0.025 )　　# モデルの適合度　※２
#Likelihood ratio test= 30.5  on 3 df,   p=1e-06
#Wald test            = 29.93  on 3 df,   p=1e-06
#Score (logrank) test = 30.5  on 3 df,   p=1e-06

・年齢（age）: 回帰係数（coef）は0.011067で、p値は0.232416。この因子は統計学的に有意であるといえない。
・性別（sex）: 回帰係数（coef）は-0.552612で、p値は0.000986。この因子は統計学的に有意であるといえる。ハザード比が約0.575であり、女性（通常、sexが0）に比べて男性（通常、sexが1）のリスクが約0.575倍になることを示す。
・ECOGパフォーマンスステータス（ph.ecog）: 回帰係数（coef）は0.463728で、p値は4.45e-05。この因子も統計学的に有意であるといえる。ハザード比が約1.59であり、ECOGが1単位増加するごとにリスクが約1.59倍になることを示す。

※１　ハザード比の95%信頼区間
真のハザード比がこの範囲にある確率が95％であることを示す。
例えば性別（sex）だと0.4142～0.7994の間に真のハザード比がある確率が95%であることを意味する。

※２　モデルの適合度　
Concordance（一致係数）※３は0.637で、これは予測の精度を示し、1に近いほど良い予測を意味する。
Likelihood ratio test（尤度比検定）、Wald test、Score (logrank) test(ログランク検定)は、回帰係数に対する検定である。
帰無仮説H0：回帰係数βは0である。　対立仮説H1：回帰係数βは0でない　⇒p<0.05で、帰無仮説棄却
今回結果はすべて同じく、p値が1e-06で統計学的に有意であることを示す。これらのテストはモデル全体の有効性を評価し、変数がモデルに有意に寄与していることを示し、モデルは適合していると結論付けられる。

上記の結果より、
予測モデル式から、以下のことが分かりました。
・性別は生存期間に統計学的に有意な影響を与えており、男性のリスクは女性の約0.575倍です。
・ECOGパフォーマンスステータスも生存期間に統計学的に有意な影響を与えており、パフォーマンスステータスが1単位増加するごとにリスクは約1.59倍増加します。
・年齢は生存期間に統計学的に有意な影響を与えていません。
これらの結果から、性別とECOGパフォーマンスステータスは生存期間に影響を与える重要な要因であると結論づけることができます。
さらに、モデルの適合度も評価されました。
・一致係数（Concordance）は0.637で、モデルの予測力を示します。
・Likelihood Ratio Test、Wald Test、およびScore (Logrank) Testのすべてでp値が1e-06未満であり、これらのテストはモデル全体が統計学的に有意であることを示します。
これにより、提案されたコックス比例ハザードモデルはデータに適合しており、信頼できる予測を提供していると評価できます。
以上のことから、性別とECOGパフォーマンスステータスは生存期間に影響を与える重要な要因であり、モデル全体としてもデータに対して適切にフィットしていることが確認されました。

※３Concordance（C-index 一致係数）での判断についても諸説あるので、以下の文献を参考にされたい。
・Harrell FE Jr, Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Stat Med. 1996 Feb 28;15(4):361-87. doi: 10.1002/(SICI)1097-0258(19960229)15:4<361::AID-SIM168>3.0.CO;2-4.
・Hartman N, Kim S, He K, Kalbfleisch JD. Pitfalls of the concordance index for survival outcomes. Stat Med. 2023 Jun 15;42(13):2179-2190. doi: 10.1002/sim.9717.

3. Rによるデータの可視化

library("survminer")
fit<-survfit(cox_model)　# survfitオブジェクトを作成
ggsurvplot(fit, data = lung, ggtheme = theme_bw())

95%信頼区間が色付きで表示されます。