CODEWITHSUNDEEP
rcox-regression
Dana
Dana

Reputation: 75

Cox proportional hazard model-interaction

I want to test for an interaction ( for Cox proportional hazard model) between type of transplant and disease type using main effects and interaction terms on the data bone marrow transplant study at Ohio State University.

Here is the used code for the data:

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

Here is the code for interaction, but I'm not sure what it should be written in here (myData$(here?) from the following code to be able to run it.

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)

Upvotes: 0

Views: 227

Answers (1)

Edward
Edward

Reputation: 18683

An easy way is to separate the four groups variables using the separate function from the tidyr package.

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

Then you can put these two new variables (disease and transpl) into the Cox model, with interaction term.

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

Upvotes: 1

Related Questions