R Under development (unstable) (2024-12-09 r87433 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 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) > > expect <- survexp(futime ~ 1, + rmap = list(age=(accept.dt - birth.dt), sex=1, + year=accept.dt, race='white'), jasa, cohort=F, + ratetable=survexp.usr) > > survdiff(Surv(jasa$futime, jasa$fustat) ~ offset(expect)) Call: survdiff(formula = Surv(jasa$futime, jasa$fustat) ~ offset(expect)) Observed Expected Z p 75.000 0.587 -97.119 0.000 > # Now fit the 6 models found in Kalbfleisch and Prentice, p139 > sfit.1 <- coxph(Surv(start, stop, event)~ (age + surgery)*transplant, + jasa1, method='breslow') > sfit.2 <- coxph(Surv(start, stop, event)~ year*transplant, + jasa1, method='breslow') > sfit.3 <- coxph(Surv(start, stop, event)~ (age + year)*transplant, + jasa1, method='breslow') > sfit.4 <- coxph(Surv(start, stop, event)~ (year +surgery) *transplant, + jasa1, method='breslow') > sfit.5 <- coxph(Surv(start, stop, event)~ (age + surgery)*transplant + year , + jasa1, method='breslow') > sfit.6 <- coxph(Surv(start, stop, event)~ age*transplant + surgery + year, + jasa1, method='breslow') > > summary(sfit.1) Call: coxph(formula = Surv(start, stop, event) ~ (age + surgery) * transplant, data = jasa1, method = "breslow") n= 170, number of events= 75 coef exp(coef) se(coef) z Pr(>|z|) age 0.01386 1.01395 0.01813 0.765 0.445 surgery -0.54652 0.57896 0.61091 -0.895 0.371 transplant 0.11572 1.12268 0.32729 0.354 0.724 age:transplant 0.03473 1.03534 0.02725 1.274 0.202 surgery:transplant -0.29037 0.74799 0.75819 -0.383 0.702 exp(coef) exp(-coef) lower .95 upper .95 age 1.014 0.9862 0.9786 1.051 surgery 0.579 1.7272 0.1748 1.917 transplant 1.123 0.8907 0.5911 2.132 age:transplant 1.035 0.9659 0.9815 1.092 surgery:transplant 0.748 1.3369 0.1692 3.306 Concordance= 0.595 (se = 0.037 ) Likelihood ratio test= 12.45 on 5 df, p=0.03 Wald test = 11.62 on 5 df, p=0.04 Score (logrank) test = 12.02 on 5 df, p=0.03 > sfit.2 Call: coxph(formula = Surv(start, stop, event) ~ year * transplant, data = jasa1, method = "breslow") coef exp(coef) se(coef) z p year -0.2652 0.7671 0.1051 -2.522 0.0117 transplant -0.2871 0.7504 0.5137 -0.559 0.5762 year:transplant 0.1371 1.1469 0.1409 0.973 0.3306 Likelihood ratio test=8.61 on 3 df, p=0.03488 n= 170, number of events= 75 > summary(sfit.3) Call: coxph(formula = Surv(start, stop, event) ~ (age + year) * transplant, data = jasa1, method = "breslow") n= 170, number of events= 75 coef exp(coef) se(coef) z Pr(>|z|) age 0.01558 1.01571 0.01734 0.899 0.36887 year -0.27413 0.76023 0.10588 -2.589 0.00962 ** transplant -0.59388 0.55218 0.54222 -1.095 0.27339 age:transplant 0.03380 1.03438 0.02795 1.209 0.22653 year:transplant 0.20228 1.22419 0.14247 1.420 0.15566 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 exp(coef) exp(-coef) lower .95 upper .95 age 1.0157 0.9845 0.9818 1.0508 year 0.7602 1.3154 0.6178 0.9356 transplant 0.5522 1.8110 0.1908 1.5981 age:transplant 1.0344 0.9668 0.9792 1.0926 year:transplant 1.2242 0.8169 0.9259 1.6185 Concordance= 0.63 (se = 0.036 ) Likelihood ratio test= 14.85 on 5 df, p=0.01 Wald test = 13.77 on 5 df, p=0.02 Score (logrank) test = 14.06 on 5 df, p=0.02 > sfit.4 Call: coxph(formula = Surv(start, stop, event) ~ (year + surgery) * transplant, data = jasa1, method = "breslow") coef exp(coef) se(coef) z p year -0.2541 0.7756 0.1077 -2.360 0.0183 surgery -0.2365 0.7894 0.6282 -0.377 0.7065 transplant -0.2969 0.7431 0.5054 -0.588 0.5568 year:transplant 0.1653 1.1797 0.1416 1.167 0.2431 surgery:transplant -0.5499 0.5770 0.7759 -0.709 0.4785 Likelihood ratio test=12.36 on 5 df, p=0.03021 n= 170, number of events= 75 > sfit.5 Call: coxph(formula = Surv(start, stop, event) ~ (age + surgery) * transplant + year, data = jasa1, method = "breslow") coef exp(coef) se(coef) z p age 0.01504 1.01516 0.01760 0.855 0.3928 surgery -0.42025 0.65689 0.61564 -0.683 0.4948 transplant 0.07412 1.07694 0.33113 0.224 0.8229 year -0.13631 0.87258 0.07097 -1.921 0.0548 age:transplant 0.02691 1.02728 0.02712 0.992 0.3210 surgery:transplant -0.29656 0.74337 0.75800 -0.391 0.6956 Likelihood ratio test=16.2 on 6 df, p=0.0127 n= 170, number of events= 75 > sfit.6 Call: coxph(formula = Surv(start, stop, event) ~ age * transplant + surgery + year, data = jasa1, method = "breslow") coef exp(coef) se(coef) z p age 0.01527 1.01539 0.01750 0.873 0.3828 transplant 0.04457 1.04558 0.32171 0.139 0.8898 surgery -0.62113 0.53734 0.36786 -1.688 0.0913 year -0.13608 0.87278 0.07090 -1.919 0.0550 age:transplant 0.02703 1.02740 0.02714 0.996 0.3193 Likelihood ratio test=16.06 on 5 df, p=0.006687 n= 170, number of events= 75 > > # Survival curve for an "average" subject, > # done once as overall, once via individual method > surv1 <- survfit(sfit.1, newdata=list(age=-2, surgery=0, transplant=0)) > newdata <- data.frame(start=c(0,50,100), stop=c(50,100, max(jasa1$stop)), + event=c(1,1,1), age=rep(-2,3), surgery=rep(0,3), + transplant=rep(0,3), name=c("Smith", "Smith", "Smith")) > surv2 <- survfit(sfit.1, newdata, id=name) > # Have to use unclass to avoid [.survfit trying to pick curves, > # remove the final element "call" because it won't match, nor will newdata > ii <- match(c("newdata", "call"), names(surv1)) > all.equal(unclass(surv1)[-ii], + unclass(surv2)[-ii]) [1] TRUE > > > # Survival curve for a subject of age 50, with prior surgery, tx at 6 months > # Remember that 'age' in jasa 1 was centered at 48 > data <- data.frame(start=c(0,183), stop=c(183,3*365), event=c(1,1), + age=c(2,2), surgery=c(1,1), transplant=c(0,1), id=c(1,1)) > # This output changed in version 3.8-0; the drop in std(surv) at 183 was > # incorrect > summary(survfit(sfit.1, data, id=id)) Call: survfit(formula = sfit.1, newdata = data, id = id) time n.risk n.event survival std.err lower 95% CI upper 95% CI 0.5 103 1 0.994 0.00722 0.980 1.000 1.0 102 3 0.975 0.01860 0.939 1.000 2.0 99 3 0.956 0.02914 0.900 1.000 4.0 96 2 0.943 0.03605 0.875 1.000 5.0 94 2 0.930 0.04286 0.849 1.000 7.0 92 1 0.923 0.04623 0.837 1.000 8.0 91 1 0.917 0.04959 0.824 1.000 11.0 89 1 0.910 0.05294 0.812 1.000 15.0 88 3 0.890 0.06278 0.775 1.000 16.0 85 1 0.883 0.06608 0.763 1.000 17.0 84 1 0.877 0.06928 0.751 1.000 20.0 83 2 0.864 0.07538 0.728 1.000 27.0 81 1 0.857 0.07849 0.716 1.000 29.0 80 1 0.850 0.08160 0.705 1.000 31.0 78 1 0.844 0.08473 0.693 1.000 34.0 77 1 0.837 0.08786 0.681 1.000 35.0 76 1 0.830 0.09098 0.669 1.000 36.0 75 1 0.823 0.09412 0.658 1.000 38.0 74 1 0.816 0.09727 0.646 1.000 39.0 72 2 0.802 0.10349 0.623 1.000 42.0 70 1 0.795 0.10664 0.611 1.000 44.0 69 1 0.788 0.10982 0.600 1.000 49.0 68 1 0.781 0.11300 0.588 1.000 50.0 67 1 0.774 0.11614 0.577 1.000 52.0 66 1 0.767 0.11925 0.565 1.000 57.0 65 1 0.760 0.12238 0.554 1.000 60.0 64 1 0.752 0.12552 0.542 1.000 65.0 63 1 0.745 0.12866 0.531 1.000 67.0 62 2 0.730 0.13494 0.508 1.000 68.0 60 1 0.722 0.13809 0.497 1.000 71.0 59 2 0.707 0.14420 0.474 1.000 76.0 57 1 0.699 0.14729 0.463 1.000 77.0 56 1 0.691 0.15043 0.451 1.000 79.0 55 1 0.683 0.15362 0.439 1.000 80.0 54 1 0.674 0.15680 0.428 1.000 84.0 53 1 0.666 0.16005 0.416 1.000 89.0 52 1 0.657 0.16326 0.404 1.000 95.0 51 1 0.648 0.16648 0.392 1.000 99.0 50 1 0.639 0.16972 0.380 1.000 101.0 49 1 0.630 0.17293 0.368 1.000 109.0 47 1 0.621 0.17611 0.356 1.000 148.0 45 1 0.611 0.17927 0.344 1.000 152.0 44 1 0.601 0.18236 0.332 1.000 164.0 43 1 0.592 0.18551 0.320 1.000 185.0 41 1 0.583 0.18292 0.315 1.000 187.0 40 1 0.574 0.18040 0.310 1.000 206.0 39 1 0.565 0.17790 0.305 1.000 218.0 38 1 0.556 0.17543 0.299 1.000 262.0 37 1 0.546 0.17298 0.294 1.000 284.0 35 2 0.527 0.16813 0.282 0.985 307.0 33 1 0.517 0.16575 0.276 0.969 333.0 32 1 0.507 0.16341 0.270 0.954 339.0 31 1 0.497 0.16110 0.263 0.938 342.0 29 1 0.486 0.15876 0.256 0.922 583.0 21 1 0.471 0.15581 0.246 0.900 674.0 17 1 0.452 0.15262 0.233 0.876 732.0 16 1 0.433 0.14944 0.220 0.852 851.0 14 1 0.410 0.14600 0.204 0.824 979.0 11 1 0.383 0.14229 0.185 0.793 995.0 10 1 0.356 0.13883 0.166 0.765 1031.0 9 1 0.330 0.13567 0.147 0.739 > > # These should all give the same answer > # When there are offsets, the default curve is always for someone with > # the mean offset. > j.age <- jasa$age -48 > fit1 <- coxph(Surv(futime, fustat) ~ j.age, data=jasa) > fit2 <- coxph(Surv(futime, fustat) ~ j.age, jasa, init=fit1$coef, iter=0) > fit3 <- coxph(Surv(start, stop, event) ~ age, jasa1) > fit4 <- coxph(Surv(start, stop, event) ~ offset(age*fit1$coef), jasa1) > > s1 <- survfit(fit1, list(j.age=fit3$means), censor=FALSE) > s2 <- survfit(fit2, list(j.age=fit3$means), censor=FALSE) > s3 <- survfit(fit3, censor=FALSE) > s4 <- survfit(fit4, censor=FALSE) > > all.equal(s1$surv, s2$surv) [1] TRUE > all.equal(s1$surv, s3$surv) [1] TRUE > all.equal(s1$surv, s4$surv) [1] TRUE > > # Still the same answer, fit multiple strata at once > # Strata 1 has independent coefs of strata 2, so putting in > # the other data should not affect it > ll <- nrow(jasa1) > ss <- rep(0:1, c(ll,ll)) > tdata <- with(jasa1, data.frame(start=rep(start,2), stop=rep(stop,2), + event=rep(event,2), ss=ss, age=rep(age,2), + age2 = (rep(age,2))^2 * ss)) > fit <- coxph(Surv(start, stop, event) ~ age*strata(ss) + age2, tdata) > # Above replaced these 2 lines, which kill Splus5 as of 8/98 > # Something with data frames, I expect. > #fit <- coxph(Surv(rep(start,2), rep(stop,2), rep(event,2)) ~ > # rep(age,2)*strata(ss) + I(rep(age,2)^2*ss) ) > all.equal(fit$coef[1], fit3$coef) [1] TRUE > s5 <- survfit(fit, data.frame(age=fit3$means, age2=0, ss=0), censor=FALSE) > all.equal(s5$surv[1:(s5$strata[1])], s3$surv) [1] TRUE > > > > proc.time() user system elapsed 0.89 0.15 0.98