Cox 比例风险模型-交互答案

【问题标题】：Cox proportional hazard model-interactionCox 比例风险模型-交互
【发布时间】：2020-04-01 14:07:16
【问题描述】：

我想使用俄亥俄州立大学骨髓移植研究数据的主效应和交互作用项来检验移植类型和疾病类型之间的交互作用（对于 Cox 比例风险模型）。

这是数据的使用代码：

time_Allo_NHL<- c(28,32,49,84,357,933,1078,1183,1560,2114,2144)
censor_Allo_NHL<- c(rep(1,5), rep(0,6))
df_Allo_NHL <- data.frame(group = "Allo NHL", 
                          time = time_Allo_NHL,
                          censor = censor_Allo_NHL,
                          Z1 = c(90,30,40,60,70,90,100,90,80,80,90),
                          Z2 = c(24,7,8,10,42,9,16,16,20,27,5))

time_Auto_NHL<- c(42,53,57,63,81,140,176,210,252,476,524,1037)
censor_Auto_NHL<- c(rep(1,7), rep(0,1), rep(1,1), rep(0,1), rep(1,1), rep(0,1))
df_Auto_NHL <- data.frame(group = "Auto NHL", 
                          time = time_Auto_NHL, 
                          censor = censor_Auto_NHL,
                          Z1 = c(80,90,30,60,50,100,80,90,90,90,90,90),
                          Z2 = c(19,17,9,13,12,11,38,16,21,24,39,84))

time_Allo_HOD<- c(2,4,72,77,79)
censor_Allo_HOD<- c(rep(1,5))
df_Allo_HOD <- data.frame(group = "Allo HOD", 
                          time = time_Allo_HOD, 
                          censor = censor_Allo_HOD,
                          Z1 = c(20,50,80,60,70),
                          Z2 = c(34,28,59,102,71))

time_Auto_HOD<- c(30,36,41,52,62,108,132,180,307,406,446,484,748,1290,1345)
censor_Auto_HOD<- c(rep(1,7), rep(0,8))
df_Auto_HOD <- data.frame(group = "Auto HOD", 
                          time = time_Auto_HOD, 
                          censor = censor_Auto_HOD,
                          Z1 = c(90,80,70,60,90,70,60,100,100,100,100,90,90,90,80),
                          Z2 = c(73,61,34,18,40,65,17,61,24,48,52,84,171,20,98))

myData <- Reduce(rbind, list(df_Allo_NHL, df_Auto_NHL, df_Allo_HOD, df_Auto_HOD))

这里是交互的代码，但是我不确定应该在下面写什么（myData$(here?) 从下面的代码中可以运行它。

n<-length(myData$time)
n

for (i in 1:n){
  if (myData$(here?)[i]==2)
    myData$W1[i] <-1
  else myData$W1[i]<-0
}

for (i in 1:n){
  if (myData$(here?)[i]==2)
    myData$W2[i] <-1
  else myData$W2[i]<-0
}

myData

Coxfit.W<-coxph(Surv(time,censor)~W1+W2+W1*W2, data = myData)
summary(Coxfit.W)

【问题讨论】：

标签： r cox-regression

【解决方案1】：

一个简单的方法是使用 tidyr 包中的 separate 函数将四个组变量分开。

library(tidyr)

myData <- separate(myData, col=group, into=c("disease","transpl"))
head(myData)
  disease transpl time censor Z1 Z2
1    Allo     NHL   28      1 90 24
2    Allo     NHL   32      1 30  7
3    Allo     NHL   49      1 40  8
4    Allo     NHL   84      1 60 10
5    Allo     NHL  357      1 70 42
6    Allo     NHL  933      0 90  9

然后您可以将这两个新变量（disease 和 transpl）放入 Cox 模型中，并带有交互项。

Coxfit.W<-coxph(Surv(time,censor)~transpl*disease, data = myData)
summary(Coxfit.W)

Call:
coxph(formula = Surv(time, censor) ~ transpl * disease, data = myData)

  n= 43, number of events= 26 

                          coef exp(coef) se(coef)      z Pr(>|z|)   
transplNHL             -1.8212    0.1618   0.6747 -2.699  0.00695 **
diseaseAuto            -1.6628    0.1896   0.6188 -2.687  0.00721 **
transplNHL:diseaseAuto  2.3050   10.0244   0.8494  2.714  0.00665 **

                       exp(coef) exp(-coef) lower .95 upper .95
transplNHL                0.1618    6.17946   0.04312    0.6073
diseaseAuto               0.1896    5.27387   0.05638    0.6377
transplNHL:diseaseAuto   10.0244    0.09976   1.89700   52.9720

【讨论】：