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) > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y)) > # > # These results can be found in Miller > # > fit <- coxph(Surv(aml$time, aml$status) ~ aml$x, method='breslow') > fit Call: coxph(formula = Surv(aml$time, aml$status) ~ aml$x, method = "breslow") coef exp(coef) se(coef) z p aml$xNonmaintained 0.9042 2.4700 0.5122 1.765 0.0775 Likelihood ratio test=3.3 on 1 df, p=0.06945 n= 23, number of events= 18 > resid(fit, type='mart') 1 2 3 4 5 6 0.86225539 0.79200985 -0.20799015 0.74818869 0.65652976 -0.39796610 7 8 9 10 11 12 0.45424957 0.25475051 -1.05400917 -0.55400917 -1.55400917 0.87844483 13 14 15 16 17 18 0.87844483 0.74006941 0.74006941 0.57677292 -0.51373647 0.15162716 19 20 21 22 23 0.01702219 -0.14897252 -0.56448258 -1.15185244 -1.60340676 > resid(fit, type='score') 1 2 3 4 5 6 -0.546856248 -0.492501830 0.141063944 -0.479907930 -0.447416819 0.268453990 7 8 9 10 11 12 -0.235908976 -0.072655945 0.640826596 0.640826596 0.640826596 0.237767767 13 14 15 16 17 18 0.237767767 0.232585063 0.232585063 0.203878910 -0.165307985 0.044923326 19 20 21 22 23 0.007079721 -0.039651990 -0.181184547 -0.395076175 -0.472116894 > resid(fit, type='scho') 5 5 8 8 9 12 13 0.2706690 0.2706690 0.3081229 0.3081229 -0.6423931 0.3360212 -0.6335658 18 23 23 27 30 31 33 -0.6494307 -0.6791937 0.3208063 0.3269751 0.3360212 -0.5970995 0.3505693 34 43 45 48 -0.5525731 0.3778334 0.5484457 0.0000000 > > # Test the drop of an itercept: should have no effect > fit2 <- coxph(Surv(time, status) ~ x -1, method='breslow', + data=aml) > aeq(fit$loglik, fit2$loglik) [1] TRUE > aeq(coef(fit), coef(fit2)) [1] TRUE > aeq(fit$var, fit2$var) [1] TRUE > > fit <- survfit(Surv(aml$time, aml$status) ~ aml$x) > fit Call: survfit(formula = Surv(aml$time, aml$status) ~ aml$x) n events median 0.95LCL 0.95UCL aml$x=Maintained 11 7 31 18 NA aml$x=Nonmaintained 12 11 23 8 NA > summary(fit) Call: survfit(formula = Surv(aml$time, aml$status) ~ aml$x) aml$x=Maintained time n.risk n.event survival std.err lower 95% CI upper 95% CI 9 11 1 0.909 0.0867 0.7541 1.000 13 10 1 0.818 0.1163 0.6192 1.000 18 8 1 0.716 0.1397 0.4884 1.000 23 7 1 0.614 0.1526 0.3769 0.999 31 5 1 0.491 0.1642 0.2549 0.946 34 4 1 0.368 0.1627 0.1549 0.875 48 2 1 0.184 0.1535 0.0359 0.944 aml$x=Nonmaintained time n.risk n.event survival std.err lower 95% CI upper 95% CI 5 12 2 0.8333 0.1076 0.6470 1.000 8 10 2 0.6667 0.1361 0.4468 0.995 12 8 1 0.5833 0.1423 0.3616 0.941 23 6 1 0.4861 0.1481 0.2675 0.883 27 5 1 0.3889 0.1470 0.1854 0.816 30 4 1 0.2917 0.1387 0.1148 0.741 33 3 1 0.1944 0.1219 0.0569 0.664 43 2 1 0.0972 0.0919 0.0153 0.620 45 1 1 0.0000 NaN NA NA > survdiff(Surv(aml$time, aml$status)~ aml$x) Call: survdiff(formula = Surv(aml$time, aml$status) ~ aml$x) N Observed Expected (O-E)^2/E (O-E)^2/V aml$x=Maintained 11 7 10.69 1.27 3.4 aml$x=Nonmaintained 12 11 7.31 1.86 3.4 Chisq= 3.4 on 1 degrees of freedom, p= 0.07 > > # > # Test out the weighted K-M > # > # First, equal case weights- shouldn't change the survival, but will > # halve the variance > temp2 <-survfit(Surv(aml$time, aml$status)~1, weight=rep(2,23)) > temp <-survfit(Surv(time, status)~1, aml) > aeq(temp$surv, temp2$surv) [1] TRUE > aeq(temp$std.err^2, 2*temp2$std.err^2) [1] TRUE > > # Risk weights-- use a null Cox model > tfit <- coxph(Surv(aml$time, aml$status) ~ offset(log(1:23))) > sfit <- survfit(tfit, stype=2, ctype=1, censor=FALSE) > > # Now compute it by hand. The survfit program will produce a curve > # corresponding to the mean offset. > # Ties are a nuisance, the line above forced the Nelson rather than Efron > # to make it easier > rscore <- exp(log(1:23) - mean(log(1:23)))[order(aml$time)] > atime <- sort(aml$time) > denom <- rev(cumsum(rev(rscore))) > denom <- denom[match(unique(atime), atime)] > deaths <- tapply(aml$status, aml$time, sum) > chaz <- cumsum(deaths/denom) > all.equal(sfit$surv, as.vector(exp(-chaz[deaths>0]))) [1] TRUE > > # And the Efron result > summary(survfit(tfit)) Call: survfit(formula = tfit) time n.risk n.event survival std.err lower 95% CI upper 95% CI 5 23 2 0.932 0.0461 0.8463 1.000 8 21 2 0.863 0.0637 0.7467 0.997 9 19 1 0.827 0.0704 0.6999 0.977 12 18 1 0.793 0.0755 0.6576 0.955 13 17 1 0.757 0.0801 0.6152 0.931 18 14 1 0.719 0.0846 0.5709 0.905 23 13 2 0.645 0.0907 0.4893 0.849 27 11 1 0.607 0.0929 0.4496 0.819 30 9 1 0.565 0.0955 0.4054 0.787 31 8 1 0.519 0.0982 0.3579 0.752 33 7 1 0.474 0.0994 0.3140 0.715 34 6 1 0.423 0.1009 0.2649 0.675 43 5 1 0.373 0.1006 0.2198 0.633 45 4 1 0.312 0.1009 0.1657 0.588 48 2 1 0.199 0.1102 0.0674 0.589 > > # Lots of ties, so its a good test case > x1 <- coxph(Surv(time, status)~x, aml, method='efron') > x1 Call: coxph(formula = Surv(time, status) ~ x, data = aml, method = "efron") coef exp(coef) se(coef) z p xNonmaintained 0.9155 2.4981 0.5119 1.788 0.0737 Likelihood ratio test=3.38 on 1 df, p=0.06581 n= 23, number of events= 18 > x2 <- coxph(Surv(rep(0,23),time, status) ~x, aml, method='efron') > aeq(x1$coef, x2$coef) [1] TRUE > > > proc.time() user system elapsed 0.82 0.15 0.98