Rで統計解析するための備忘録的なやつです。

コックス比例ハザード回帰分析の式

\log_e\frac{h(t)}{h_0(t)} = a_1x_1 + a_2x_2 + a_3x_3 + ...\\
\log_e\frac{h(t)}{h_0(t)}：ハザード対数比\\
x_i：説明変数\\
a_i：偏回帰係数

必要なパッケージ

survivalとsurvminerパッケージをインストールし、読み込み。

# ライブラリの読み込み
library(survival)
library(survminer)

使用するデータ

今回使うデータは、これ。

kidney

組み込みデータです。
人工透析装置の利用が、各腎疾患の患者の生存時間にどう影響するかを調べたデータらしいです。

まずは読み込んで確認します。

# データの読み込みと確認
DF <- kidney
summary(DF)

実行結果は、

> #データの読み込みと確認
> DF <- kidney
> summary(DF)
       id            time           status            age            sex         disease       frail      
 Min.   : 1.0   Min.   :  2.0   Min.   :0.0000   Min.   :10.0   Min.   :1.000   Other:26   Min.   :0.200  
 1st Qu.:10.0   1st Qu.: 16.0   1st Qu.:1.0000   1st Qu.:34.0   1st Qu.:1.000   GN   :18   1st Qu.:0.600  
 Median :19.5   Median : 39.5   Median :1.0000   Median :45.5   Median :2.000   AN   :24   Median :1.100  
 Mean   :19.5   Mean   :101.6   Mean   :0.7632   Mean   :43.7   Mean   :1.737   PKD  : 8   Mean   :1.184  
 3rd Qu.:29.0   3rd Qu.:149.8   3rd Qu.:1.0000   3rd Qu.:54.0   3rd Qu.:2.000              3rd Qu.:1.500  
 Max.   :38.0   Max.   :562.0   Max.   :1.0000   Max.   :69.0   Max.   :2.000              Max.   :3.000

それぞれのパラメーターの説明

id：患者ID
time：生存期間
status：期間内に患者が死亡したかどうか〈0：生存（打ち切り）、1：死亡〉
age：年齢
sex：性別（1：男、2：女）
disease：腎疾患の種類
frail：虚弱性（frailty）の推定値

コックス比例ハザード回帰分析

性別と腎疾患の種類による生存率への影響を検討するモデルを構築。
**coxph()**の引数をSurv(time, status) ~ sex + disease, data = DFとします。
目的変数は、timeとstatusから計算されるハザード比の対数。
説明変数は、sexとdisease。

# コックス比例ハザード回帰分析
DFcox <- coxph(Surv(time, status) ~ sex + disease, data = DF)
summary(DFcox) #分析結果の確認

実行結果は、

> #コックス比例ハザード回帰分析
> DFcox <- coxph(Surv(time, status) ~ sex + disease, data = DF)
> summary(DFcox) #分析結果の確認
Call:
coxph(formula = Surv(time, status) ~ sex + disease, data = DF)

  n= 76, number of events= 58 

              coef exp(coef) se(coef)      z Pr(>|z|)    
sex        -1.4774    0.2282   0.3569 -4.140 3.48e-05 ***
diseaseGN   0.1392    1.1494   0.3635  0.383   0.7017    
diseaseAN   0.4132    1.5116   0.3360  1.230   0.2188    
diseasePKD -1.3671    0.2549   0.5889 -2.321   0.0203 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

           exp(coef) exp(-coef) lower .95 upper .95
sex           0.2282     4.3815   0.11339    0.4594
diseaseGN     1.1494     0.8700   0.56368    2.3437
diseaseAN     1.5116     0.6616   0.78245    2.9202
diseasePKD    0.2549     3.9238   0.08035    0.8084

Concordance= 0.696  (se = 0.039 )
Likelihood ratio test= 17.56  on 4 df,   p=0.002
Wald test            = 19.77  on 4 df,   p=6e-04
Score (logrank) test = 19.97  on 4 df,   p=5e-04

coefは変数ごとの係数なので、回帰モデルはこんな感じになります。性別の効果は有意で、値が大きいほどハザードが起きる確率が低い、つまり（値が2である）女性ほど死ぬ確率が低いということです。

\log_e\frac{h(t)}{h_0(t)} = -1.474 × sex + 0.1392 × diseaseGN + 0.4132 × diseaseAN -1.3671 × diseasePKD\\

カプラン-マイヤー曲線

腎疾患ごとのカプラン-マイヤー曲線を描いてみます。
survivalパッケージのsurvfit()関数とsurvminerパッケージの**ggsurvplot()**関数を使います。

DFfit <- survfit(Surv(time, status) ~ disease, data = DF)
ggsurvplot(fit = DFfit, data = DF)

こんな感じの描画

見てくれをいい感じに。

# カプラン-マイヤー曲線の完成形（病気の種類ごと）
ggsurvplot(fit = DFfit, data = DF,
           size = 1, #size:線の太さ
           conf.int = FALSE, #conf.int:信頼区間
           pval = TRUE, pval.method = TRUE, #pval:p値
           risk.table = TRUE, #risk.table:リスクテーブル
           risk.table.height = 0.3, #risk.table.height:リスクテーブルの高さ
           risk.table.col = "strata", #risk.table.col:リスクテーブルの色
           ggtheme = theme_bw())

最終的なスクリプト

# コックス比例ハザード回帰分析

# Rをきれいにする
rm(list = ls())

# ライブラリの読み込み
library(survival)
library(survminer)

# データの読み込みと確認
DF <- kidney
summary(DF)

# コックス比例ハザード回帰分析
DFcox <- coxph(Surv(time, status) ~ sex + disease, data = DF)
summary(DFcox) #分析結果の確認

# カプラン-マイヤー曲線（病気の種類ごと）
DFfit <- survfit(Surv(time, status) ~ disease, data = DF)
ggsurvplot(fit = DFfit, data = DF)

# カプラン-マイヤー曲線の完成形（病気の種類ごと）
ggsurvplot(fit = DFfit, data = DF,
           size = 1, #size:線の太さ
           conf.int = FALSE, #conf.int:信頼区間
           pval = TRUE, pval.method = TRUE, #pval:p値
           risk.table = TRUE, #risk.table:リスクテーブル
           risk.table.height = 0.3, #risk.table.height:リスクテーブルの高さ
           risk.table.col = "strata", #risk.table.col:リスクテーブルの色
           ggtheme = theme_bw())