生物科学研究所 井口研究室
Laboratory of Biology, Okaya, Nagano, Japan
R バグ: Cox 比例ハザード Schoenfeld 検定 ggcoxzph グラフ
井口豊（生物科学研究所，長野県岡谷市）
最終更新： 2023 年 3 月 11 日
グラフの信頼区間のバグ
生存時間解析 Cox 回帰モデルで，比例ハザードが成り立つかどうか検定する方法として， Schoenfeld 残差検定がある。統計ソフト R のパッケージ survminer にある ggcoxzph 関数によるグラフが利用されることも多い。しかし，この信頼区間（2 SE で近似）にはバグがあり， GitHub から修正プログラムが出ている。要するに，信頼区間が広すぎる，のである。
日本語で検索しても，このバグは話題になっていないようで，気づいていないようなサイトも見られる。例えば，以下のサイト。
利用したデータやスクリプトが丁寧に書かれてるので，修正プログラムで検証してみた。
#############
library(survival)
library(survminer)
library(ggplot2)
library(tidyverse)
d <- veteran %>%
mutate(age_cat = if_else((age<60), "under60", "over60"),
age_cat = factor(age_cat),
trt = factor(trt, labels=c("standard", "test")),
prior = factor(prior, labels=c("N"," Y")))
fit <- coxph(
Surv(time, status)~trt + celltype + age_cat + prior,
data=d
)
fit_test <- cox.zph(fit)
# バグ修正前
g.unfix<- ggcoxzph(fit_test)
ggpar(g.unfix, subtitle = "Bug unfixed", theme(aspect.ratio = 1))
# バグ修正後
g.fix<- ggcoxzph(fit_test)
ggpar(g.fix, subtitle = "Bug fixed")
################
結果は，以下のように，バグ修正の前（Bug unfixed）と後（Bug fixed）で，それぞれ図 1，図 2 として示してある。
図 1. Schoenfeld の比例ハザード性検定（バグ修正前）
図 2. Schoenfeld の比例ハザード性検定（バグ修正後）
比例ハザード性の仮定が成立しているか検証する[R] の Schoenfeld 検定グラフは，やはり未修正のプログラムを利用しているようであり，前述の図 1 と同じく，ほぼ全てのデータ点が信頼区間の中に入ってしまっている。
