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. > library(survival) > > options(na.action=na.exclude) # preserve missings > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type > aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...) > > # > # Test out all of the routines on a more complex data set > # > temp <- survfit(Surv(time, status) ~ ph.ecog, lung) > summary(temp, times=c(30*1:11, 365*1:3)) Call: survfit(formula = Surv(time, status) ~ ph.ecog, data = lung) 1 observation deleted due to missingness ph.ecog=0 time n.risk n.event survival std.err lower 95% CI upper 95% CI 30 60 3 0.952 0.0268 0.9012 1.000 60 58 2 0.921 0.0341 0.8562 0.990 90 56 2 0.889 0.0396 0.8146 0.970 120 56 0 0.889 0.0396 0.8146 0.970 150 55 1 0.873 0.0419 0.7946 0.959 180 52 2 0.841 0.0461 0.7553 0.936 210 48 2 0.808 0.0498 0.7164 0.912 240 45 0 0.808 0.0498 0.7164 0.912 270 38 2 0.770 0.0543 0.6709 0.884 300 33 2 0.727 0.0591 0.6203 0.853 330 29 2 0.681 0.0637 0.5670 0.818 365 22 6 0.535 0.0728 0.4100 0.699 730 5 11 0.193 0.0707 0.0943 0.396 ph.ecog=1 time n.risk n.event survival std.err lower 95% CI upper 95% CI 30 111 2 0.982 0.0124 0.9583 1.000 60 110 3 0.956 0.0193 0.9186 0.994 90 104 4 0.920 0.0255 0.8718 0.972 120 99 5 0.876 0.0310 0.8174 0.939 150 93 6 0.823 0.0359 0.7556 0.896 180 82 8 0.751 0.0407 0.6756 0.836 210 68 9 0.666 0.0450 0.5831 0.760 240 57 6 0.604 0.0474 0.5176 0.704 270 53 4 0.561 0.0487 0.4729 0.665 300 46 3 0.527 0.0495 0.4384 0.633 330 40 4 0.480 0.0504 0.3903 0.589 365 34 4 0.431 0.0509 0.3417 0.543 730 7 21 0.114 0.0388 0.0582 0.222 ph.ecog=2 time n.risk n.event survival std.err lower 95% CI upper 95% CI 30 46 5 0.9000 0.0424 0.82057 0.987 60 43 2 0.8600 0.0491 0.76900 0.962 90 40 3 0.8000 0.0566 0.69647 0.919 120 34 4 0.7174 0.0641 0.60216 0.855 150 31 3 0.6541 0.0680 0.53342 0.802 180 26 6 0.5275 0.0719 0.40385 0.689 210 21 4 0.4431 0.0717 0.32266 0.608 240 17 3 0.3766 0.0705 0.26100 0.543 270 17 0 0.3766 0.0705 0.26100 0.543 300 13 3 0.3102 0.0677 0.20223 0.476 330 11 2 0.2624 0.0651 0.16135 0.427 365 9 2 0.2147 0.0614 0.12258 0.376 730 1 6 0.0371 0.0345 0.00601 0.229 ph.ecog=3 time n.risk n.event survival std.err lower 95% CI upper 95% CI 30 1 0 1 0 1 1 60 1 0 1 0 1 1 90 1 0 1 0 1 1 > print(temp[2:3]) Call: survfit(formula = Surv(time, status) ~ ph.ecog, data = lung) n events median 0.95LCL 0.95UCL ph.ecog=1 113 82 306 268 429 ph.ecog=2 50 44 199 156 288 > > temp <- survfit(Surv(time, status)~1, lung, type='fleming', + conf.int=.9, conf.type='log-log', error='tsiatis') > summary(temp, times=30 *1:5) Call: survfit(formula = Surv(time, status) ~ 1, data = lung, error = "tsiatis", type = "fleming", conf.int = 0.9, conf.type = "log-log") time n.risk n.event survival std.err lower 90% CI upper 90% CI 30 219 10 0.956 0.0135 0.928 0.974 60 213 7 0.926 0.0173 0.891 0.950 90 201 10 0.882 0.0213 0.842 0.913 120 189 10 0.838 0.0244 0.793 0.874 150 179 10 0.794 0.0268 0.745 0.834 > > temp <- survdiff(Surv(time, status) ~ inst, lung, rho=.5) > print(temp, digits=6) Call: survdiff(formula = Surv(time, status) ~ inst, data = lung, rho = 0.5) n=227, 1 observation deleted due to missingness. N Observed Expected (O-E)^2/E (O-E)^2/V inst=1 36 21.190058 17.455181 0.799149708 1.171232977 inst=2 5 3.173330 1.964395 0.744007932 0.860140808 inst=3 19 10.663476 11.958755 0.140294489 0.200472362 inst=4 4 2.245347 3.559344 0.485085848 0.677874608 inst=5 9 5.010883 4.500982 0.057765161 0.077128402 inst=6 14 8.862602 7.078516 0.449665221 0.582743947 inst=7 8 4.445647 4.416133 0.000197254 0.000253632 inst=10 4 2.901923 2.223283 0.207150016 0.249077097 inst=11 18 7.807867 9.525163 0.309611863 0.422142221 inst=12 23 14.009656 12.216768 0.263117640 0.365712493 inst=13 20 9.140983 11.863298 0.624699853 0.874238212 inst=15 6 3.170744 3.558447 0.042241456 0.057938955 inst=16 16 8.870360 9.992612 0.126038005 0.175170113 inst=21 13 9.263733 4.460746 5.171484268 6.149354145 inst=22 17 8.278566 11.971473 1.139171459 1.645863937 inst=26 6 1.627074 3.542694 1.035821659 1.286365543 inst=32 7 1.792468 2.679904 0.293869782 0.343966668 inst=33 2 0.929177 0.416202 0.632249272 0.676682390 Chisq= 15.1 on 17 degrees of freedom, p= 0.5904 > > # verify that the zph routine does the actual score test > dtime <- lung$time[lung$status==2] > lung2 <- survSplit(Surv(time, status) ~ ., lung, cut=dtime) > > cfit1 <-coxph(Surv(time, status) ~ ph.ecog + ph.karno + pat.karno + wt.loss + + sex + age + strata(inst), lung) > cfit2 <-coxph(Surv(tstart, time, status) ~ ph.ecog + ph.karno + pat.karno + + wt.loss + sex + age + strata(inst), lung2) > all.equal(cfit1$loglik, cfit2$loglik) [1] TRUE > all.equal(coef(cfit1), coef(cfit2)) [1] TRUE > > # the above verifies that the data set is correct > zp1 <- cox.zph(cfit1, transform="log") > zp2 <- cox.zph(cfit2, transform="log") > # everything should match but the call > icall <- match("Call", names(zp1)) > all.equal(unclass(zp2)[-icall], unclass(zp1)[-icall]) [1] TRUE > > # now compute score tests one variable at a time > ncoef <- length(coef(cfit2)) > check <- double(ncoef) > cname <- names(coef(cfit2)) > for (i in 1:ncoef) { + temp <- log(lung2$time) * lung2[[cname[i]]] + # score test for this new variable + tfit <- coxph(Surv(tstart, time, status) ~ ph.ecog + ph.karno + pat.karno + + wt.loss + sex + age + strata(inst) + + temp, lung2, init=c(cfit2$coef, 0), iter=0) + check[i] <- tfit$score + } > aeq(check, zp1$table[1:ncoef,1]) # skip the 'global' test [1] TRUE > > # > # Tests of using "." > # > fit1 <- coxph(Surv(time, status) ~ . - meal.cal - wt.loss - inst, lung) > fit2 <- update(fit1, .~. - ph.karno) > fit3 <- coxph(Surv(time, status) ~ age + sex + ph.ecog + pat.karno, lung) > all.equal(fit2, fit3) [1] TRUE > > proc.time() user system elapsed 2.07 0.14 2.20