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())