R Under development (unstable) (2023-08-12 r84939 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > options(na.action=na.exclude) # preserve missings > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type > library(survival) > > # > # Test the coxph program on the Ovarian data > summary(survfit(Surv(futime, fustat)~1, data=ovarian), censor=TRUE) Call: survfit(formula = Surv(futime, fustat) ~ 1, data = ovarian) time n.risk n.event survival std.err lower 95% CI upper 95% CI 59 26 1 0.962 0.0377 0.890 1.000 115 25 1 0.923 0.0523 0.826 1.000 156 24 1 0.885 0.0627 0.770 1.000 268 23 1 0.846 0.0708 0.718 0.997 329 22 1 0.808 0.0773 0.670 0.974 353 21 1 0.769 0.0826 0.623 0.949 365 20 1 0.731 0.0870 0.579 0.923 377 19 0 0.731 0.0870 0.579 0.923 421 18 0 0.731 0.0870 0.579 0.923 431 17 1 0.688 0.0919 0.529 0.894 448 16 0 0.688 0.0919 0.529 0.894 464 15 1 0.642 0.0965 0.478 0.862 475 14 1 0.596 0.0999 0.429 0.828 477 13 0 0.596 0.0999 0.429 0.828 563 12 1 0.546 0.1032 0.377 0.791 638 11 1 0.497 0.1051 0.328 0.752 744 10 0 0.497 0.1051 0.328 0.752 769 9 0 0.497 0.1051 0.328 0.752 770 8 0 0.497 0.1051 0.328 0.752 803 7 0 0.497 0.1051 0.328 0.752 855 6 0 0.497 0.1051 0.328 0.752 1040 5 0 0.497 0.1051 0.328 0.752 1106 4 0 0.497 0.1051 0.328 0.752 1129 3 0 0.497 0.1051 0.328 0.752 1206 2 0 0.497 0.1051 0.328 0.752 1227 1 0 0.497 0.1051 0.328 0.752 > > # Various models > coxph(Surv(futime, fustat)~ age, data=ovarian) Call: coxph(formula = Surv(futime, fustat) ~ age, data = ovarian) coef exp(coef) se(coef) z p age 0.16162 1.17541 0.04974 3.249 0.00116 Likelihood ratio test=14.29 on 1 df, p=0.0001564 n= 26, number of events= 12 > coxph(Surv(futime, fustat)~ resid.ds, data=ovarian) Call: coxph(formula = Surv(futime, fustat) ~ resid.ds, data = ovarian) coef exp(coef) se(coef) z p resid.ds 1.2092 3.3507 0.6724 1.798 0.0721 Likelihood ratio test=3.76 on 1 df, p=0.05251 n= 26, number of events= 12 > coxph(Surv(futime, fustat)~ rx, data=ovarian) Call: coxph(formula = Surv(futime, fustat) ~ rx, data = ovarian) coef exp(coef) se(coef) z p rx -0.5964 0.5508 0.5870 -1.016 0.31 Likelihood ratio test=1.05 on 1 df, p=0.3052 n= 26, number of events= 12 > coxph(Surv(futime, fustat)~ ecog.ps, data=ovarian) Call: coxph(formula = Surv(futime, fustat) ~ ecog.ps, data = ovarian) coef exp(coef) se(coef) z p ecog.ps 0.3984 1.4894 0.5864 0.679 0.497 Likelihood ratio test=0.47 on 1 df, p=0.4935 n= 26, number of events= 12 > > coxph(Surv(futime, fustat)~ resid.ds + rx + ecog.ps, data=ovarian) Call: coxph(formula = Surv(futime, fustat) ~ resid.ds + rx + ecog.ps, data = ovarian) coef exp(coef) se(coef) z p resid.ds 1.3466 3.8444 0.6803 1.980 0.0478 rx -0.7493 0.4727 0.5949 -1.260 0.2078 ecog.ps 0.4528 1.5727 0.5903 0.767 0.4431 Likelihood ratio test=6.03 on 3 df, p=0.1101 n= 26, number of events= 12 > coxph(Surv(futime, fustat)~ age + rx + ecog.ps, data=ovarian) Call: coxph(formula = Surv(futime, fustat) ~ age + rx + ecog.ps, data = ovarian) coef exp(coef) se(coef) z p age 0.1470 1.1583 0.0463 3.175 0.0015 rx -0.8146 0.4428 0.6342 -1.285 0.1990 ecog.ps 0.1032 1.1087 0.6064 0.170 0.8649 Likelihood ratio test=15.92 on 3 df, p=0.00118 n= 26, number of events= 12 > coxph(Surv(futime, fustat)~ age + resid.ds + ecog.ps, data=ovarian) Call: coxph(formula = Surv(futime, fustat) ~ age + resid.ds + ecog.ps, data = ovarian) coef exp(coef) se(coef) z p age 0.14240 1.15303 0.05196 2.740 0.00614 resid.ds 0.66261 1.93984 0.75045 0.883 0.37727 ecog.ps 0.16634 1.18097 0.61486 0.271 0.78675 Likelihood ratio test=15.1 on 3 df, p=0.001732 n= 26, number of events= 12 > coxph(Surv(futime, fustat)~ age + resid.ds + rx, data=ovarian) Call: coxph(formula = Surv(futime, fustat) ~ age + resid.ds + rx, data = ovarian) coef exp(coef) se(coef) z p age 0.1285 1.1372 0.0473 2.718 0.00657 resid.ds 0.6964 2.0065 0.7585 0.918 0.35858 rx -0.8489 0.4279 0.6392 -1.328 0.18416 Likelihood ratio test=16.77 on 3 df, p=0.0007889 n= 26, number of events= 12 > > # Residuals > fit <- coxph(Surv(futime, fustat)~ age + resid.ds + rx + ecog.ps, data=ovarian) > resid(fit) 1 2 3 4 5 6 0.84103277 0.54424388 0.59670824 -0.11281376 0.75111588 -0.32609026 7 8 9 10 11 12 0.59998927 0.29570718 -2.15325805 0.76243469 0.06474272 -0.11680752 13 14 15 16 17 18 -1.22562781 -0.63474839 -0.07535824 -0.17058905 -0.22986038 -0.14654862 19 20 21 22 23 24 -0.18762920 -0.12771548 -0.53373114 -0.65480022 0.95866131 0.82111675 25 26 0.55136554 -0.09154014 > resid(fit, 'dev') 1 2 3 4 5 6 1.41281595 0.69505907 0.78916003 -0.47500266 1.13106322 -0.80757694 7 8 9 10 11 12 0.79532966 0.33122166 -2.07521471 1.16179002 0.06619519 -0.48333740 13 14 15 16 17 18 -1.56564862 -1.12671948 -0.38822221 -0.58410453 -0.67802711 -0.54138455 19 20 21 22 23 24 -0.61258338 -0.50540178 -1.03318066 -0.54976346 2.11059000 1.34157009 25 26 0.70736314 -0.42787881 > resid(fit, 'scor') age resid.ds rx ecog.ps 1 2.26503249 0.05686357 -0.10565379 -0.42661688 2 3.02525428 0.04641312 -0.08623662 -0.34821275 3 -0.06851355 0.07131430 -0.13250357 0.06167527 4 0.94597623 -0.02541510 -0.06423496 0.05971729 5 -5.41507168 0.21605962 -0.32258092 -0.39333909 6 1.48999552 0.24899474 0.14035143 -0.15380664 7 -0.68612431 0.13740891 0.28392482 0.29196506 8 0.93116906 0.08428957 0.16040160 0.18430641 9 -8.20092595 -0.51356176 0.95647608 1.11337112 10 0.95287510 -0.31078224 0.21463992 0.17363388 11 2.85526159 0.09417730 -0.14186603 -0.07586086 12 0.92721107 0.07495002 -0.05400751 0.07061578 13 -1.93962967 -0.43919871 -0.56668535 -0.48467672 14 0.63185387 -0.22745949 -0.29348437 0.38373600 15 1.41495195 0.04835392 0.04051535 0.04555769 16 2.54591188 0.10945916 0.09171493 -0.06745975 17 4.40282381 -0.08236953 0.12358137 -0.09089870 18 1.97071836 0.09403352 0.07878991 0.08859570 19 0.77692371 0.12039304 -0.08675286 0.11343089 20 1.76784279 -0.04576632 -0.05905095 0.07721016 21 -0.82272526 0.34247077 -0.24677770 -0.21106494 22 -3.48057998 -0.03965965 0.07368852 -0.26669335 23 -14.86623758 0.28137017 -0.52279208 -0.43881151 24 3.96084273 -0.56566921 0.34648950 0.44907410 25 4.30025715 0.15241262 0.22417527 -0.20390438 26 0.31490641 0.07091764 -0.05212198 0.04845623 > resid(fit, 'scho') age resid.ds rx ecog.ps 59 2.69315603 0.06761160 -0.1256239 -0.5072536 115 5.36390105 0.08039116 -0.1493686 -0.6031318 156 -0.89877512 0.10683985 -0.1985108 0.1984379 268 6.95664326 0.12857949 -0.2389036 0.2388157 329 -15.73656605 0.28889883 -0.5367805 -0.4634169 353 4.06104389 -0.70587654 0.4535120 0.5282024 365 5.50035833 0.25348264 0.4796230 -0.4413864 431 -8.06809505 0.27490176 -0.4297023 -0.5248323 464 -2.15471559 0.23158421 0.5066040 0.4814387 475 0.57065051 0.25226659 0.5518479 0.5244351 563 0.06487219 -0.47274522 0.3319974 0.2747028 638 1.64752655 -0.50593437 -0.6446947 0.2939883 > > fit <- coxph(Surv(futime, fustat) ~ age + ecog.ps + strata(rx), data=ovarian) > summary(fit) Call: coxph(formula = Surv(futime, fustat) ~ age + ecog.ps + strata(rx), data = ovarian) n= 26, number of events= 12 coef exp(coef) se(coef) z Pr(>|z|) age 0.13853 1.14858 0.04801 2.885 0.00391 ** ecog.ps -0.09670 0.90783 0.62994 -0.154 0.87800 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 age 1.1486 0.8706 1.0454 1.262 ecog.ps 0.9078 1.1015 0.2641 3.120 Concordance= 0.819 (se = 0.058 ) Likelihood ratio test= 12.71 on 2 df, p=0.002 Wald test = 8.43 on 2 df, p=0.01 Score (logrank) test = 12.24 on 2 df, p=0.002 > summary(survfit(fit)) Call: survfit(formula = fit) rx=1 time n.risk n.event survival std.err lower 95% CI upper 95% CI 59 13 1 0.978 0.0266 0.9275 1 115 12 1 0.951 0.0478 0.8620 1 156 11 1 0.910 0.0760 0.7722 1 268 10 1 0.862 0.1055 0.6776 1 329 9 1 0.737 0.1525 0.4909 1 431 8 1 0.627 0.1704 0.3680 1 638 5 1 0.333 0.2296 0.0865 1 rx=2 time n.risk n.event survival std.err lower 95% CI upper 95% CI 353 13 1 0.943 0.0560 0.839 1.000 365 12 1 0.880 0.0812 0.735 1.000 464 9 1 0.789 0.1143 0.594 1.000 475 8 1 0.697 0.1349 0.477 1.000 563 7 1 0.597 0.1494 0.366 0.975 > sfit <- survfit(fit, list(age=c(30,70), ecog.ps=c(2,3))) #two columns > sfit Call: survfit(formula = fit, newdata = list(age = c(30, 70), ecog.ps = c(2, 3))) n events median 0.95LCL 0.95UCL rx=1, 1 13 7 NA NA NA rx=2, 1 13 5 NA NA NA rx=1, 2 13 7 268 115 NA rx=2, 2 13 5 365 353 NA > summary(sfit) Call: survfit(formula = fit, newdata = list(age = c(30, 70), ecog.ps = c(2, 3))) rx=1 time n.risk n.event survival1 survival2 59 13 1 0.999 0.87905 115 12 1 0.999 0.74575 156 11 1 0.998 0.57398 268 10 1 0.996 0.41764 329 9 1 0.992 0.16673 431 8 1 0.988 0.06489 638 5 1 0.973 0.00161 rx=2 time n.risk n.event survival1 survival2 353 13 1 0.999 0.7092 365 12 1 0.997 0.4738 464 9 1 0.994 0.2494 475 8 1 0.991 0.1207 563 7 1 0.987 0.0489 > > > # Check of offset + surv, added 7/2000 > fit1 <- coxph(Surv(futime, fustat) ~ age + rx, ovarian, + control=coxph.control(eps=1e-8)) > fit2 <- coxph(Surv(futime, fustat) ~ age + offset(rx*fit1$coef[2]), ovarian, + control=coxph.control(eps=1e-8)) > all.equal(fit1$coef[1], fit2$coef[1]) [1] TRUE > > fit <- coxph(Surv(futime, fustat) ~ age + offset(rx), ovarian) > survfit(fit, censor=FALSE)$surv^exp(-1.5) [1] 0.9977751 0.9951975 0.9917927 0.9881504 0.9825769 0.9770280 0.9704304 [8] 0.9603196 0.9499085 0.9385539 0.9217097 0.9031334 > > # Check it by hand -- there are no tied times > # Remember that offsets from survfit are centered, which is 1.5 for > # this data set. > eta <- fit$coef*(ovarian$age - fit$mean) + (ovarian$rx - 1.5) > ord <- order(ovarian$futime) > risk <- exp(eta[ord]) > rsum <- rev(cumsum(rev(risk))) # cumulative risk at each time point > dead <- (ovarian$fustat[ord]==1) > baseline <- cumsum(1/rsum[dead]) > all.equal(survfit(fit, censor=FALSE)$surv, exp(-baseline)) [1] TRUE > > rm(fit, fit1, fit2, ord, eta, risk, rsum, dead, baseline, sfit) > > proc.time() user system elapsed 0.92 0.14 1.06