Cox回归HR分组

Question

我想对以下问题进行Cox回归:一组患者是否接受“药物”治疗(0 / 1)。我的时间变量“时间”告诉我，病人被观察了多少天，以及如果病人存活或死亡(死亡= 1，存活= 0)的“状态”。

library(survival)

set.seed(123)

df <- data.frame(time = round(runif(100, min = 1, max = 70)), 
                 status = round(runif(100, min = 0, max = 1)),
                 drug = round(runif(100, min = 0, max = 1)),
                 age40 = round(runif(100, min = 0, max = 1)), 
                 stringsAsFactors = FALSE)

object <- Surv(df$time, df$status)
model <- coxph(object ~ drug, data = df)
summary(model)

这对我来说很好，并且告诉我，HR是0.89，所以这种药物可以防止病人死亡。

现在我想做一些子群分析，，f.e。如果病人是<= 40岁或> 40岁，HR如何变化(age40: 0 vs 1)。

要将变量"age40“包含到coxph中，我所要做的就够了吗？

object2 <- Surv(df$time, df$status)
model2 <- coxph(object2 ~ drug + age40, data = df)
summary(model2)

如果这样做，我在drug1摘要中的HR略有更改为0.86，而age40 (1.12)则会得到另一个。

现在我的问题是:如果病人是<= 40岁或40岁以上，那么在治疗中死亡的危险比率是多少(药物= 1)。

编辑:另一个问题是用图形显示药物对森林地块，f.e状态影响的不同HRs。就像这样：https://rpkgs.datanovia.com/survminer/reference/ggforest-2.png。而不是“性”、"rx“、”坚持“等等。我想显示Age40 =0vs.1和其他变量的HRs，比如高血压=0比1，吸烟者=0比1。

谢谢!

Answer 1

您需要在predict上使用的函数是model2，它需要提供一个新数据参数，其中包括您想要考虑的所有情况：

exp( predict(model2, newdata=expand.grid(drug=c(0,1), age40=c(0,1))) )
#        1         2         3         4 
#1.0000000 0.8564951 1.1268713 0.9651598 

现在您有所有4例可能的药物和age40的组合。基本案例有一个统一的值，因为您正在估算风险比率，形成了{drug=0，age40=0}的基线情况，您可以看到其他风险比率与什么相关。

expand.grid(drug=c(0,1), age40=c(0,1))
  drug age40
1    0     0
2    1     0
3    0     1
4    1     1

请注意，对于单独考虑的每个年龄类别，drug=0与drug=1的比值是相同的。如果你想知道药物在这两种年龄类别中的作用是否不同，你可以使用一个交互作用模型：

model3 <- coxph(object2 ~ drug * age40, data = df)
summary(model3)
#----------------
Call:
coxph(formula = object2 ~ drug * age40, data = df)

  n= 100, number of events= 50 

               coef exp(coef) se(coef)      z Pr(>|z|)
drug       -0.18524   0.83091  0.45415 -0.408    0.683
age40       0.09611   1.10089  0.39560  0.243    0.808
drug:age40  0.05679   1.05843  0.63094  0.090    0.928

           exp(coef) exp(-coef) lower .95 upper .95
drug          0.8309     1.2035    0.3412     2.024
age40         1.1009     0.9084    0.5070     2.390
drug:age40    1.0584     0.9448    0.3073     3.645

Concordance= 0.528  (se = 0.042 )
Likelihood ratio test= 0.34  on 3 df,   p=1
Wald test            = 0.33  on 3 df,   p=1
Score (logrank) test = 0.33  on 3 df,   p=1

现在的影响估计有点不同：

 exp( predict(model3, newdata=expand.grid(drug=c(0,1), age40=c(0,1))) )
#        1         2         3         4 
#1.0000000 0.8309089 1.1008850 0.9681861 

Answer 2

使用参数strata。

coxph(object ~ drug + strata(age40), data = df)