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'), + show.signif.stars=FALSE) #ensure constrast type > library(survival) > > # > # Fit the models found in Wei et. al. > # > wfit <- coxph(Surv(stop, event) ~ (rx + size + number)* strata(enum), + cluster=id, bladder, ties='breslow') > wfit Call: coxph(formula = Surv(stop, event) ~ (rx + size + number) * strata(enum), data = bladder, ties = "breslow", cluster = id) coef exp(coef) se(coef) robust se z p rx -0.51762 0.59594 0.31576 0.30750 -1.683 0.09231 size 0.06789 1.07025 0.10125 0.08529 0.796 0.42602 number 0.23599 1.26617 0.07608 0.07208 3.274 0.00106 rx:strata(enum)enum=2 -0.10182 0.90319 0.50427 0.32654 -0.312 0.75518 rx:strata(enum)enum=3 -0.18226 0.83339 0.55790 0.39162 -0.465 0.64165 rx:strata(enum)enum=4 -0.13317 0.87531 0.65813 0.49680 -0.268 0.78865 size:strata(enum)enum=2 -0.14401 0.86588 0.16800 0.11190 -1.287 0.19812 size:strata(enum)enum=3 -0.27920 0.75639 0.20862 0.15115 -1.847 0.06472 size:strata(enum)enum=4 -0.27106 0.76257 0.25146 0.18563 -1.460 0.14422 number:strata(enum)enum=2 -0.09843 0.90625 0.11931 0.11439 -0.861 0.38951 number:strata(enum)enum=3 -0.06616 0.93598 0.12983 0.11672 -0.567 0.57083 number:strata(enum)enum=4 0.09280 1.09724 0.14657 0.11754 0.790 0.42982 Likelihood ratio test=29.39 on 12 df, p=0.003443 n= 340, number of events= 112 > > # Check the rx coefs versus Wei, et al, JASA 1989 > rx <- c(1,4,5,6) # the treatment coefs above > cmat <- diag(4); cmat[1,] <- 1; #contrast matrix > wfit$coef[rx] %*% cmat # the coefs in their paper (table 5) [,1] [,2] [,3] [,4] [1,] -0.5176209 -0.6194404 -0.6998771 -0.6507935 > t(cmat) %*% wfit$var[rx,rx] %*% cmat # var matrix (eqn 3.2) [,1] [,2] [,3] [,4] [1,] 0.09455501 0.06017669 0.05677331 0.0437777 [2,] 0.06017669 0.13242834 0.13011557 0.1160420 [3,] 0.05677331 0.13011557 0.17235879 0.1590865 [4,] 0.04377770 0.11604200 0.15908650 0.2398112 > > # Anderson-Gill fit > fita <- coxph(Surv(start, stop, event) ~ rx + size + number, cluster=id, + bladder2, ties='breslow') > summary(fita) Call: coxph(formula = Surv(start, stop, event) ~ rx + size + number, data = bladder2, ties = "breslow", cluster = id) n= 178, number of events= 112 coef exp(coef) se(coef) robust se z Pr(>|z|) rx -0.45979 0.63142 0.19996 0.25801 -1.782 0.07474 size -0.04256 0.95833 0.06903 0.07555 -0.563 0.57317 number 0.17164 1.18726 0.04733 0.06131 2.799 0.00512 exp(coef) exp(-coef) lower .95 upper .95 rx 0.6314 1.5837 0.3808 1.047 size 0.9583 1.0435 0.8264 1.111 number 1.1873 0.8423 1.0528 1.339 Concordance= 0.634 (se = 0.032 ) Likelihood ratio test= 16.77 on 3 df, p=8e-04 Wald test = 11.76 on 3 df, p=0.008 Score (logrank) test = 18.57 on 3 df, p=3e-04, Robust = 11.44 p=0.01 (Note: the likelihood ratio and score tests assume independence of observations within a cluster, the Wald and robust score tests do not). > > # Prentice fits. Their model 1 a and b are the same > fit1p <- coxph(Surv(stop, event) ~ rx + size + number, bladder2, + subset=(enum==1), ties='breslow') > fit2pa <- coxph(Surv(stop, event) ~ rx + size + number, bladder2, + subset=(enum==2), ties='breslow') > fit2pb <- coxph(Surv(stop-start, event) ~ rx + size + number, bladder2, + subset=(enum==2), ties='breslow') > fit3pa <- coxph(Surv(stop, event) ~ rx + size + number, bladder2, + subset=(enum==3), ties='breslow') > #and etc. > fit1p Call: coxph(formula = Surv(stop, event) ~ rx + size + number, data = bladder2, subset = (enum == 1), ties = "breslow") coef exp(coef) se(coef) z p rx -0.51762 0.59594 0.31576 -1.639 0.10115 size 0.06789 1.07025 0.10125 0.671 0.50253 number 0.23599 1.26617 0.07608 3.102 0.00192 Likelihood ratio test=9.66 on 3 df, p=0.02164 n= 85, number of events= 47 > fit2pa Call: coxph(formula = Surv(stop, event) ~ rx + size + number, data = bladder2, subset = (enum == 2), ties = "breslow") coef exp(coef) se(coef) z p rx -0.424214 0.654284 0.402200 -1.055 0.292 size -0.125033 0.882467 0.117088 -1.068 0.286 number 0.001987 1.001989 0.093760 0.021 0.983 Likelihood ratio test=2.02 on 3 df, p=0.5688 n= 46, number of events= 29 > fit2pb Call: coxph(formula = Surv(stop - start, event) ~ rx + size + number, data = bladder2, subset = (enum == 2), ties = "breslow") coef exp(coef) se(coef) z p rx -0.259112 0.771737 0.405110 -0.640 0.522 size -0.116363 0.890152 0.119236 -0.976 0.329 number -0.005709 0.994307 0.096671 -0.059 0.953 Likelihood ratio test=1.27 on 3 df, p=0.7353 n= 46, number of events= 29 > fit3pa Call: coxph(formula = Surv(stop, event) ~ rx + size + number, data = bladder2, subset = (enum == 3), ties = "breslow") coef exp(coef) se(coef) z p rx -0.89854 0.40716 0.55352 -1.623 0.105 size 0.08497 1.08869 0.20861 0.407 0.684 number -0.01716 0.98298 0.12802 -0.134 0.893 Likelihood ratio test=4.16 on 3 df, p=0.2452 n= 27, number of events= 22 > rm(rx, cmat, wfit, fita, fit1p, fit2pa, fit2pb, fit3pa) > > proc.time() user system elapsed 0.89 0.07 0.96