从0开始学习R语言--Day51--PH检验

在用cox回归做分析时，我们一般会得出各种变量在结局的风险影响（HR大于1，就代表变量值增大，对应结局影响的风险就随之增大），但是这里有个坏处是，cox回归得到的是瞬时风险值，我们最多得到一段时间，这个变量的影响。

而PH检验，就是为了验证变量对结局的风险影响是否会随时间变化，如果不会，那么求得的p值就会大于0.05，也就说明回归里的HR值的可解释性高。

以下是一个例子：

# 加载生存分析包
library(survival)
library(survminer) # 用于PH检验和绘图

# 生成数据（100个样本）
set.seed(123)
n <- 100
group <- sample(c("Treatment", "Placebo"), n, replace = TRUE)
time <- ifelse(group == "Treatment", 
               rexp(n, rate = 0.1),  # 治疗组风险更低
               rexp(n, rate = 0.2))   # 安慰剂组风险更高
status <- rbinom(n, 1, 0.8)  # 80%的事件观察到（1=事件发生，0=删失）

# 创建数据框
df <- data.frame(time, status, group)
head(df)

# 拟合Cox模型
cox_model <- coxph(Surv(time, status) ~ group, data = df, ties = "breslow")
summary(cox_model)

# PH检验（全局检验）
ph_test <- cox.zph(cox_model)
print(ph_test)

# 可视化检验结果（ Schoenfeld残差 vs 时间）
ggcoxzph(ph_test)

输出：

Call:
coxph(formula = Surv(time, status) ~ group, data = df, ties = "breslow")

  n= 100, number of events= 77 

                  coef exp(coef) se(coef)      z Pr(>|z|)
groupTreatment -0.2469    0.7813   0.2440 -1.012    0.312

               exp(coef) exp(-coef) lower .95 upper .95
groupTreatment    0.7813       1.28    0.4843      1.26

Concordance= 0.52  (se = 0.032 )
Likelihood ratio test= 1  on 1 df,   p=0.3
Wald test            = 1.02  on 1 df,   p=0.3
Score (logrank) test = 1.03  on 1 df,   p=0.3

       chisq df    p
group   0.34  1 0.56
GLOBAL  0.34  1 0.56

结果表明，虽然group在cox回归的p值较低，HR也小于1，但是group和global的PH是大于0.05的，这可能意味着方向是对的，只是特征要处理一下，或者是要剖析一下结果。