File: C:/Program Files/R/R-4.4.0/library/survival/tests/doaml.Rout.save
R Under development (unstable) (2019-05-15 r76504) -- "Unsuffered Consequences"
Copyright (C) 2019 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
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.820 0.028 0.848