Load survival package and look at AML data.
library(survival)
aml
## time status x
## 1 9 1 Maintained
## 2 13 1 Maintained
## 3 13 0 Maintained
## 4 18 1 Maintained
## 5 23 1 Maintained
## 6 28 0 Maintained
## 7 31 1 Maintained
## 8 34 1 Maintained
## 9 45 0 Maintained
## 10 48 1 Maintained
## 11 161 0 Maintained
## 12 5 1 Nonmaintained
## 13 5 1 Nonmaintained
## 14 8 1 Nonmaintained
## 15 8 1 Nonmaintained
## 16 12 1 Nonmaintained
## 17 16 0 Nonmaintained
## 18 23 1 Nonmaintained
## 19 27 1 Nonmaintained
## 20 30 1 Nonmaintained
## 21 33 1 Nonmaintained
## 22 43 1 Nonmaintained
## 23 45 1 Nonmaintained
From the Maintained/Nonmaintained factor, create a continuous variable
set.seed(2016-03-25)
aml$xc <- sapply(as.numeric(aml$x), function(x) rnorm(1, mean=x*5))
aml
## time status x xc
## 1 9 1 Maintained 5.135729
## 2 13 1 Maintained 4.661955
## 3 13 0 Maintained 3.961140
## 4 18 1 Maintained 4.515426
## 5 23 1 Maintained 6.257935
## 6 28 0 Maintained 5.939816
## 7 31 1 Maintained 5.939539
## 8 34 1 Maintained 6.683239
## 9 45 0 Maintained 2.072196
## 10 48 1 Maintained 4.889319
## 11 161 0 Maintained 4.760556
## 12 5 1 Nonmaintained 9.997055
## 13 5 1 Nonmaintained 10.715022
## 14 8 1 Nonmaintained 10.469295
## 15 8 1 Nonmaintained 11.157180
## 16 12 1 Nonmaintained 9.339657
## 17 16 0 Nonmaintained 10.829025
## 18 23 1 Nonmaintained 9.686948
## 19 27 1 Nonmaintained 10.582921
## 20 30 1 Nonmaintained 8.568784
## 21 33 1 Nonmaintained 10.657252
## 22 43 1 Nonmaintained 9.302998
## 23 45 1 Nonmaintained 9.996449
See what that looks like:
boxplot(xc~x, data=aml, ylab="xc", xlab="x")
Create survival object.
s <- Surv(aml$time, aml$status)
s
## [1] 9 13 13+ 18 23 28+ 31 34 45+ 48 161+ 5 5 8
## [15] 8 12 16+ 23 27 30 33 43 45
Fit the Cox PH model to the factor variable:
fitx <- coxph(s ~ aml$x)
summary(fitx)
## Call:
## coxph(formula = s ~ aml$x)
##
## n= 23, number of events= 18
##
## coef exp(coef) se(coef) z Pr(>|z|)
## aml$xNonmaintained 0.9155 2.4981 0.5119 1.788 0.0737
##
## exp(coef) exp(-coef) lower .95 upper .95
## aml$xNonmaintained 2.498 0.4003 0.9159 6.813
##
## Concordance= 0.619 (se = 0.073 )
## Rsquare= 0.137 (max possible= 0.976 )
## Likelihood ratio test= 3.38 on 1 df, p=0.06581
## Wald test = 3.2 on 1 df, p=0.07371
## Score (logrank) test = 3.42 on 1 df, p=0.06454
To the continuous variable you made up:
fitxc <- coxph(s ~ aml$xc)
summary(fitxc)
## Call:
## coxph(formula = s ~ aml$xc)
##
## n= 23, number of events= 18
##
## coef exp(coef) se(coef) z Pr(>|z|)
## aml$xc 0.2184 1.2441 0.1016 2.15 0.0316
##
## exp(coef) exp(-coef) lower .95 upper .95
## aml$xc 1.244 0.8038 1.019 1.518
##
## Concordance= 0.657 (se = 0.085 )
## Rsquare= 0.204 (max possible= 0.976 )
## Likelihood ratio test= 5.25 on 1 df, p=0.02197
## Wald test = 4.62 on 1 df, p=0.03159
## Score (logrank) test = 4.99 on 1 df, p=0.0255
If you have broom installed the models are nicer.
library(broom)
rbind(tidy(fitx), tidy(fitxc))
## term estimate std.error statistic p.value conf.low
## 1 aml$xNonmaintained 0.9155326 0.5119343 1.788379 0.07371486 -0.08784017
## 2 aml$xc 0.2184150 0.1016089 2.149566 0.03158956 0.01926525
## conf.high
## 1 1.9189053
## 2 0.4175648