Running 'aareg.R' [2s]
  Comparing 'aareg.Rout' to 'aareg.Rout.save' ... OK
  Running 'anova.R' [149s]
  Running 'bladder.R' [172s]
  Running 'book1.R' [170s]
  Running 'book2.R' [147s]
  Running 'book3.R' [169s]
  Running 'book4.R' [2s]
  Comparing 'book4.Rout' to 'book4.Rout.save' ... OK
  Running 'book5.R' [2s]
  Comparing 'book5.Rout' to 'book5.Rout.save' ... OK
  Running 'book6.R' [2s]
  Comparing 'book6.Rout' to 'book6.Rout.save' ... OK
  Running 'book7.R' [2s]
  Comparing 'book7.Rout' to 'book7.Rout.save' ... OK
  Running 'brier.R' [195s]
  Running 'cancer.R' [3s]
  Comparing 'cancer.Rout' to 'cancer.Rout.save' ... OK
  Running 'checkSurv2.R' [176s]
  Running 'clogit.R' [2s]
  Comparing 'clogit.Rout' to 'clogit.Rout.save' ... OK
  Running 'concordance.R' [155s]
  Running 'concordance2.R' [160s]
  Running 'concordance3.R' [167s]
  Running 'counting.R' [1s]
  Comparing 'counting.Rout' to 'counting.Rout.save' ... OK
  Running 'coxsurv.R' [170s]
  Running 'coxsurv2.R' [160s]
  Running 'coxsurv3.R' [1s]
  Comparing 'coxsurv3.Rout' to 'coxsurv3.Rout.save' ... OK
  Running 'coxsurv4.R' [173s]
  Running 'coxsurv5.R' [2s]
  Comparing 'coxsurv5.Rout' to 'coxsurv5.Rout.save' ... OK
  Running 'coxsurv6.R' [3s]
  Comparing 'coxsurv6.Rout' to 'coxsurv6.Rout.save' ... OK
  Running 'detail.R' [153s]
  Running 'difftest.R' [154s]
  Running 'doaml.R' [171s]
  Running 'doweight.R' [2s]
  Comparing 'doweight.Rout' to 'doweight.Rout.save' ... OK
  Running 'dropspecial.R' [1s]
  Comparing 'dropspecial.Rout' to 'dropspecial.Rout.save' ... OK
  Running 'expected.R' [161s]
  Running 'expected2.R' [1s]
  Comparing 'expected2.Rout' to 'expected2.Rout.save' ... OK
  Running 'factor.R' [1s]
  Comparing 'factor.Rout' to 'factor.Rout.save' ... OK
  Running 'factor2.R' [142s]
  Running 'finegray.R' [2s]
  Comparing 'finegray.Rout' to 'finegray.Rout.save' ... OK
  Running 'fr_cancer.R' [167s]
  Running 'fr_kidney.R' [2s]
  Comparing 'fr_kidney.Rout' to 'fr_kidney.Rout.save' ... OK
  Running 'fr_lung.R' [175s]
  Running 'fr_ovarian.R' [167s]
  Running 'fr_rat1.R' [1s]
  Comparing 'fr_rat1.Rout' to 'fr_rat1.Rout.save' ... OK
  Running 'fr_resid.R' [1s]
  Comparing 'fr_resid.Rout' to 'fr_resid.Rout.save' ... OK
  Running 'fr_simple.R' [171s]
  Running 'frailty.R' [1s]
  Comparing 'frailty.Rout' to 'frailty.Rout.save' ... OK
  Running 'frank.R' [1s]
  Comparing 'frank.Rout' to 'frank.Rout.save' ... OK
  Running 'infcox.R' [1s]
  Comparing 'infcox.Rout' to 'infcox.Rout.save' ... OK
  Running 'jasa.R' [2s]
  Comparing 'jasa.Rout' to 'jasa.Rout.save' ... OK
  Running 'model.matrix.R' [1s]
  Comparing 'model.matrix.Rout' to 'model.matrix.Rout.save' ... OK
  Running 'mstate.R' [164s]
  Running 'mstrata.R' [2s]
  Comparing 'mstrata.Rout' to 'mstrata.Rout.save' ... OK
  Running 'multi2.R' [2s]
  Comparing 'multi2.Rout' to 'multi2.Rout.save' ... OK
  Running 'multi3.R' [165s]
  Running 'multistate.R' [173s]
  Running 'neardate.R' [164s]
  Running 'nested.R' [2s]
  Comparing 'nested.Rout' to 'nested.Rout.save' ... OK
  Running 'nsk.R' [167s]
  Running 'ovarian.R' [1s]
  Comparing 'ovarian.Rout' to 'ovarian.Rout.save' ... OK
  Running 'overlap.R' [174s]
  Running 'prednew.R' [2s]
  Comparing 'prednew.Rout' to 'prednew.Rout.save' ... OK
  Running 'pseudo.R' [2s]
  Comparing 'pseudo.Rout' to 'pseudo.Rout.save' ... OK
  Running 'pspline.R' [1s]
  Comparing 'pspline.Rout' to 'pspline.Rout.save' ... OK
  Running 'pyear.R' [1s]
  Comparing 'pyear.Rout' to 'pyear.Rout.save' ... OK
  Running 'quantile.R' [2s]
  Comparing 'quantile.Rout' to 'quantile.Rout.save' ... OK
  Running 'r_lung.R' [2s]
  Comparing 'r_lung.Rout' to 'r_lung.Rout.save' ... OK
  Running 'r_resid.R' [2s]
  Comparing 'r_resid.Rout' to 'r_resid.Rout.save' ... OK
  Running 'r_sas.R' [162s]
  Running 'r_scale.R' [2s]
  Comparing 'r_scale.Rout' to 'r_scale.Rout.save' ... OK
  Running 'r_stanford.R' [1s]
  Comparing 'r_stanford.Rout' to 'r_stanford.Rout.save' ... OK
  Running 'r_strata.R' [1s]
  Comparing 'r_strata.Rout' to 'r_strata.Rout.save' ... OK
  Running 'r_tdist.R' [163s]
  Running 'r_user.R' [1s]
  Comparing 'r_user.Rout' to 'r_user.Rout.save' ... OK
  Running 'ratetable.R' [1s]
  Comparing 'ratetable.Rout' to 'ratetable.Rout.save' ... OK
  Running 'residsf.R' [2s]
  Comparing 'residsf.Rout' to 'residsf.Rout.save' ... OK
  Running 'royston.R' [1s]
  Comparing 'royston.Rout' to 'royston.Rout.save' ... OK
  Running 'rttright.R' [150s]
  Running 'singtest.R' [2s]
  Comparing 'singtest.Rout' to 'singtest.Rout.save' ... OK
  Running 'strata2.R' [1s]
  Comparing 'strata2.Rout' to 'strata2.Rout.save' ... OK
  Running 'stratatest.R' [1s]
  Comparing 'stratatest.Rout' to 'stratatest.Rout.save' ... OK
  Running 'summary_survfit.R' [2s]
  Comparing 'summary_survfit.Rout' to 'summary_survfit.Rout.save' ... OK
  Running 'surv.R' [145s]
  Running 'survSplit.R' [150s]
  Running 'survcheck.R' [1s]
  Comparing 'survcheck.Rout' to 'survcheck.Rout.save' ... OK
  Running 'survfit1.R' [154s]
  Running 'survfit2.R' [143s]
  Running 'survreg1.R' [1s]
  Comparing 'survreg1.Rout' to 'survreg1.Rout.save' ... OK
  Running 'survreg2.R' [2s]
  Comparing 'survreg2.Rout' to 'survreg2.Rout.save' ... OK
  Running 'survtest.R' [127s]
  Running 'testci.R' [1s]
  Comparing 'testci.Rout' to 'testci.Rout.save' ... OK
  Running 'testci2.R' [1s]
  Comparing 'testci2.Rout' to 'testci2.Rout.save' ... OK
  Running 'testnull.R' [1s]
  Comparing 'testnull.Rout' to 'testnull.Rout.save' ... OK
  Running 'testreg.R' [1s]
  Comparing 'testreg.Rout' to 'testreg.Rout.save' ... OK
  Running 'tiedtime.R' [97s]
  Running 'tmerge.R' [1s]
  Comparing 'tmerge.Rout' to 'tmerge.Rout.save' ... OK
  Running 'tmerge2.R' [78s]
  Running 'tmerge3.R' [62s]
  Running 'tt.R' [1s]
  Comparing 'tt.Rout' to 'tt.Rout.save' ... OK
  Running 'tt2.R' [63s]
  Running 'turnbull.R' [1s]
  Comparing 'turnbull.Rout' to 'turnbull.Rout.save' ... OK
  Running 'update.R' [61s]
  Running 'yates0.R' [59s]
  Running 'yates1.R' [58s]
  Running 'yates2.R' [60s]
  Running 'zph.R' [83s]
Running the tests in 'tests/anova.R' failed.
Complete output:
  > #
  > # Test out anova, with strata terms
  > #
  > options(na.action=na.omit)
  > library(survival)
  > 
  > fit1 <- coxph(Surv(time, status) ~ ph.ecog + wt.loss + strata(sex) + 
  +               poly(age,3), lung)
  > ztemp <- anova(fit1)
  > 
  > tdata <- na.omit(lung[, c('time', 'status', 'ph.ecog', 'wt.loss', 'sex', 'age')])
  > fit2 <- coxph(Surv(time, status)~ ph.ecog + wt.loss + poly(age,3) + strata(sex),
  +               data=tdata)
  > ztemp2 <- anova(fit2)
  > all.equal(ztemp, ztemp2)
  [1] TRUE
  > 
  > 
  > fit2 <-  coxph(Surv(time, status) ~ ph.ecog + wt.loss + strata(sex), tdata)
  > fit3 <-  coxph(Surv(time, status) ~ ph.ecog + strata(sex), tdata)
  > 
  > all.equal(ztemp$loglik, c(fit1$loglik[1], fit3$loglik[2], fit2$loglik[2],
  +                          fit1$loglik[2]))
  [1] TRUE
  > all.equal(ztemp$Chisq[-1], 2* diff(ztemp$loglik))
  [1] TRUE
  > all.equal(ztemp$Df[-1], c(1,1,3))
  [1] TRUE
  > 
  > ztemp2 <- anova(fit3, fit2, fit1)
  > all.equal(ztemp2$loglik, ztemp$loglik[-1])
  [1] TRUE
  > all.equal(ztemp2$Chisq[2:3], ztemp$Chisq[3:4])
  [1] TRUE
  > # Change from ztemp2$P; it's a data frame and in R 3.0.2 abbreviated names
  > #   give a warning
  > all.equal(ztemp2[[4]][2:3], ztemp[[4]][3:4])
  [1] TRUE
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.17    0.14    1.29 
Running the tests in 'tests/bladder.R' failed.
Complete output:
  > 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 
     1.15    0.23    1.35 
Running the tests in 'tests/book1.R' failed.
Complete output:
  > 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))
  > 
  > #
  > # Tests from the appendix of Therneau and Grambsch
  > #  a. Data set 1 and Breslow estimate
  > #  The data below is not in time order, to also test sorting, and has 1 NA
  > #
  > test1 <- data.frame(time=  c(9, 3,1,1,6,6,8),
  +                     status=c(1,NA,1,0,1,1,0),
  +                     x=     c(0, 2,1,1,1,0,0))
  > 
  > # Nelson-Aalen influence
  > s1 <- survfit(Surv(time, status) ~1, test1, id=1:7, influence=TRUE)
  > inf1 <- matrix(c(10, rep(-2,5), 10, -2, 7,7, -11, -11)/72,
  +                ncol=2)
  > indx <- order(test1$time[!is.na(test1$status)])
  > aeq(s1$influence.chaz[indx,], inf1[,c(1,2,2,2)])
  [1] TRUE
  > 
  > # KM influence
  > inf2 <- matrix(c(-20, rep(4,5), -10, 2, -13, -13, 17, 17,
  +                  rep(0,6))/144, ncol=3)
  > aeq(s1$influence.surv[indx,], inf2[, c(1,2,2,3)])
  [1] TRUE
  > 
  > # Fleming-Harrington influence
  > s2 <- survfit(Surv(time, status) ~ 1, test1, id=1:7, ctype=2, influence=2)
  > inf3 <- matrix(c( rep(c(5, -1), c(1, 5))/36, c(5,-1)/36, 
  +                  c(21,21,-29, -29)/144), ncol=2)
  > aeq(s2$influence.chaz[indx,], inf3[,c(1,2,2,2)])
  [1] TRUE
  > 
  > 
  > # Breslow estimate
  > byhand1 <- function(beta, newx=0) {
  +     r <- exp(beta)
  +     loglik <- 2*beta - (log(3*r+3) + 2*log(r+3))
  +     u <- (6 + 3*r - r^2) / ((r+1)*(r+3))
  +     imat <- r/(r+1)^2 + 6*r/(r+3)^2
  + 
  +     x <- c(1,1,1,0,0,0)
  +     status <- c(1,0,1,1,0,1)
  +     xbar <- c(r/(r+1), r/(r+3), 0, 0)  # at times 1, 6, 8 and 9
  +     haz <- c(1/(3*r+3), 2/(r+3), 0, 1 )
  +     ties <- c(1,1,2,2,3,4)
  +     wt <- c(r,r,r,1,1,1)
  +     mart <- c(1,0,1,1,0,1) -  wt* (cumsum(haz))[ties]  #martingale residual
  + 
  +     a <- 3*(r+1)^2; b<- (r+3)^2
  +     score <- c((2*r+3)/a, -r/a, -r/a + 3*(3-r)/b,  r/a - r*(r+1)/b,
  +                r/a + 2*r/b, r/a + 2*r/b)
  + 
  +     # Schoenfeld residual
  +     scho <- c(1/(r+1), 1- (r/(3+r)), 0-(r/(3+r)) , 0)
  + 
  +     surv  <- exp(-cumsum(haz)* exp(beta*newx))
  +     varhaz.g <- cumsum(c(1/(3*r+3)^2, 2/(r+3)^2, 0, 1 ))
  + 
  +     varhaz.d <- cumsum((newx-xbar) * haz)
  + 
  +     varhaz <- (varhaz.g + varhaz.d^2/ imat) * exp(2*beta*newx)
  + 
  +     names(xbar) <- names(haz) <- 1:4
  +     names(surv) <- names(varhaz) <- 1:4
  +     list(loglik=loglik, u=u, imat=imat, xbar=xbar, haz=haz,
  + 	     mart=mart,  score=score,
  + 		scho=scho, surv=surv, var=varhaz,
  + 		varhaz.g=varhaz.g, varhaz.d=varhaz.d)
  +     }
  > 
  > 
  > 
  > fit0 <-coxph(Surv(time, status) ~x, test1, iter=0, method='breslow')
  > truth0 <- byhand1(0,0)
  > aeq(truth0$loglik, fit0$loglik[1])
  [1] TRUE
  > aeq(1/truth0$imat, fit0$var)
  [1] TRUE
  > aeq(truth0$mart, fit0$resid[c(2:6,1)])
  [1] TRUE
  > aeq(truth0$scho, resid(fit0, 'schoen'))
  [1] TRUE
  > aeq(truth0$score, resid(fit0, 'score')[c(3:7,1)])
  [1] TRUE
  > sfit <- survfit(fit0, list(x=0))
  > aeq(sfit$cumhaz, cumsum(truth0$haz))
  [1] TRUE
  > aeq(sfit$surv, exp(-cumsum(truth0$haz)))
  [1] TRUE
  > aeq(sfit$std.err^2, c(7/180, 2/9, 2/9, 11/9))
  [1] TRUE
  > aeq(resid(fit0, 'score'), c(5/24, NA, 5/12, -1/12, 7/24, -1/24, 5/24))
  [1] TRUE
  > 
  > fit1 <- coxph(Surv(time, status) ~x, test1, iter=1, method='breslow')
  > aeq(fit1$coef, 8/5)
  [1] TRUE
  > 
  > # This next gives an ignorable warning message
  > fit2 <- coxph(Surv(time, status) ~x, test1, method='breslow', iter=2)
  Warning message:
  In coxph.fit(X, Y, istrat, offset, init, control, weights = weights,  :
    Ran out of iterations and did not converge
  > aeq(round(fit2$coef, 6), 1.472724)
  [1] TRUE
  > 
  > fit <- coxph(Surv(time, status) ~x, test1, method='breslow', eps=1e-8,
  +              nocenter=NULL)
  > aeq(fit$coef, log(1.5 + sqrt(33)/2))  # the true solution
  [1] TRUE
  > truth <- byhand1(fit$coef, 0)
  > aeq(truth$loglik, fit$loglik[2])
  [1] TRUE
  > aeq(1/truth$imat, fit$var)
  [1] TRUE
  > aeq(truth$mart, fit$resid[c(2:6,1)])
  [1] TRUE
  > aeq(truth$scho, resid(fit, 'schoen'))
  [1] TRUE
  > aeq(truth$score, resid(fit, 'score')[c(3:7,1)])
  [1] TRUE
  > expect <- predict(fit, type='expected', newdata=test1) #force recalc
  > aeq(test1$status[-2] -fit$resid, expect[-2]) #tests the predict function
  [1] TRUE
  > 
  > sfit <- survfit(fit, list(x=0), censor=FALSE)
  > aeq(sfit$std.err^2, truth$var[c(1,2,4)]) # sfit skips time 8 (no events there)
  [1] TRUE
  > aeq(-log(sfit$surv), (cumsum(truth$haz))[c(1,2,4)])
  [1] TRUE
  > sfit <- survfit(fit, list(x=0), censor=TRUE)
  > aeq(sfit$std.err^2, truth$var) 
  [1] TRUE
  > aeq(-log(sfit$surv), (cumsum(truth$haz)))
  [1] TRUE
  > 
  > # 
  > # Done with the formal test, now print out lots of bits
  > #
  > resid(fit)
           1          2          3          4          5          6          7 
  -0.3333333         NA  0.7287136 -0.2712864 -0.4574271  0.6666667 -0.3333333 
  > resid(fit, 'scor')
            1           2           3           4           5           6 
   0.21138938          NA  0.13564322 -0.05049744 -0.12624360 -0.38168095 
            7 
   0.21138938 
  > resid(fit, 'scho')
           1          6          6          9 
   0.1861407  0.4069297 -0.5930703  0.0000000 
  > 
  > predict(fit, type='lp', se.fit=T)
  $fit
           1          2          3          4          5          6          7 
  -0.7376425         NA  0.7376425  0.7376425  0.7376425 -0.7376425 -0.7376425 
  
  $se.fit
          1         2         3         4         5         6         7 
  0.6278672        NA 0.6278672 0.6278672 0.6278672 0.6278672 0.6278672 
  
  > predict(fit, type='risk', se.fit=T)
  $fit
          1         2         3         4         5         6         7 
  0.4782401        NA 2.0910001 2.0910001 2.0910001 0.4782401 0.4782401 
  
  $se.fit
          1         2         3         4         5         6         7 
  0.4342009        NA 0.9079142 0.9079142 0.9079142 0.4342009 0.4342009 
  
  > predict(fit, type='expected', se.fit=T)
  $fit
          1         2         3         4         5         6         7 
  1.3333333        NA 0.2712864 0.2712864 1.4574271 0.3333333 0.3333333 
  
  $se.fit
  [1] 1.0540926        NA 0.2785989 0.2785989 1.1069433 0.3333333 0.3333333
  
  > predict(fit, type='terms', se.fit=T)
  $fit
             x
  1 -0.7376425
  2         NA
  3  0.7376425
  4  0.7376425
  5  0.7376425
  6 -0.7376425
  7 -0.7376425
  
  $se.fit
            x
  1 0.6278672
  2        NA
  3 0.6278672
  4 0.6278672
  5 0.6278672
  6 0.6278672
  7 0.6278672
  
  > 
  > summary(survfit(fit, list(x=2)))
  Call: survfit(formula = fit, newdata = list(x = 2))
  
   time n.risk n.event survival  std.err lower 95% CI upper 95% CI
      1      6       1 3.05e-01 6.50e-01     4.72e-03            1
      6      4       2 1.71e-03 1.98e-02     2.33e-13            1
      9      1       1 8.52e-12 5.29e-10     1.22e-64            1
  > 
  > proc.time()
     user  system elapsed 
     1.50    0.15    1.64 
Running the tests in 'tests/book2.R' failed.
Complete output:
  > library(survival)
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > 
  > #
  > # Tests from the appendix of Therneau and Grambsch
  > #  b. Data set 1 and Efron estimate
  > #
  > test1 <- data.frame(time=  c(9, 3,1,1,6,6,8),
  +                     status=c(1,NA,1,0,1,1,0),
  +                     x=     c(0, 2,1,1,1,0,0))
  > 
  > byhand <- function(beta, newx=0) {
  +     r <- exp(beta)
  +     loglik <- 2*beta - (log(3*r +3) + log((r+5)/2) + log(r+3))
  +     u <- (30 + 23*r - r^3)/ ((r+1)*(r+3)*(r+5))
  +     tfun <- function(x)  x - x^2
  +     imat <- tfun(r/(r+1)) + tfun(r/(r+5)) + tfun(r/(r+3))
  + 
  +     # The matrix of weights, one row per obs, one col per time
  +     #  Time of 1, 6, 6+0 (second death), and 9  
  +     wtmat <- matrix(c(1,1,1,1,1,1,
  +                       0,0,1,1,1,1,
  +                       0,0,.5, .5, 1,1,
  +                       0,0,0,0,0,1), ncol=4)
  +     wtmat <- diag(c(r,r,r,1,1,1)) %*% wtmat
  + 
  +     x <- c(1,1,1,0,0,0)
  +     status <- c(1,0,1,1,0,1)
  +     xbar <- colSums(wtmat*x)/ colSums(wtmat)
  +     haz <-  1/ colSums(wtmat)  # one death at each of the times
  + 
  +     hazmat <- wtmat %*% diag(haz) #each subject's hazard over time
  +     mart <- status - rowSums(hazmat)
  + 
  +     a <- r+1; b<- r+3; d<- r+5  # 'c' in the book, 'd' here
  +     score <- c((2*r + 3)/ (3*a^2),
  +                -r/ (3*a^2),
  +                (675+ r*(1305 +r*(756 + r*(-4 +r*(-79 -13*r)))))/(3*(a*b*d)^2),
  +                r*(1/(3*a^2) - a/(2*b^2) - b/(2*d^2)), 
  +                2*r*(177 + r*(282 +r*(182 + r*(50 + 5*r)))) /(3*(a*b*d)^2), 
  +                2*r*(177 + r*(282 +r*(182 + r*(50 + 5*r)))) /(3*(a*b*d)^2))
  + 
  +     # Schoenfeld residual
  +     d <- mean(xbar[2:3])
  +     scho <- c(1/(r+1), 1- d, 0- d , 0)
  + 
  +     surv  <- exp(-cumsum(haz)* exp(beta*newx))[c(1,3,4)]
  +     varhaz.g <- cumsum(haz^2)  # since all numerators are 1
  + 
  +     varhaz.d <- cumsum((newx-xbar) * haz)
  + 
  +     varhaz <- (varhaz.g + varhaz.d^2/ imat) * exp(2*beta*newx)
  + 
  +     list(loglik=loglik, u=u, imat=imat, xbar=xbar, haz=haz,
  + 	     mart=mart,  score=score, var.g=varhaz.g, var.d=varhaz.d,
  + 		scho=scho, surv=surv, var=varhaz[c(1,3,4)])
  +     }
  > 
  > 
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > fit0 <-coxph(Surv(time, status) ~x, test1, iter=0)
  > truth0 <- byhand(0,0)
  > aeq(truth0$loglik, fit0$loglik[1])
  [1] TRUE
  > aeq(1/truth0$imat, fit0$var)
  [1] TRUE
  > aeq(truth0$mart, fit0$resid[c(2:6,1)])
  [1] TRUE
  > aeq(resid(fit0), c(-3/4, NA, 5/6, -1/6, 5/12, 5/12, -3/4))
  [1] TRUE
  > aeq(truth0$scho, resid(fit0, 'schoen'))
  [1] TRUE
  > aeq(truth0$score, resid(fit0, 'score')[c(3:7,1)])
  [1] TRUE
  > sfit <- survfit(fit0, list(x=0), censor=FALSE)
  > aeq(sfit$std.err^2, truth0$var)
  [1] TRUE
  > aeq(sfit$surv, truth0$surv)
  [1] TRUE
  > 
  > fit <- coxph(Surv(time, status) ~x, test1, eps=1e-8, nocenter=NULL)
  > aeq(round(fit$coef,6), 1.676857)
  [1] TRUE
  > truebeta <- log(cos(acos((45/23)*sqrt(3/23))/3) * 2* sqrt(23/3))
  > truth <- byhand(truebeta, 0)
  > aeq(truth$loglik, fit$loglik[2])
  [1] TRUE
  > aeq(1/truth$imat, fit$var)
  [1] TRUE
  > aeq(truth$mart, fit$resid[c(2:6,1)])
  [1] TRUE
  > aeq(truth$scho, resid(fit, 'schoen'))
  [1] TRUE
  > aeq(truth$score, resid(fit, 'score')[c(3:7,1)])
  [1] TRUE
  > 
  > # Per comments in the source code, the below is expected to fail for Efron 
  > #  at the tied death times.  (When predicting for new data, predict
  > #  treats a time in the new data set that exactly matches one in the original
  > #  as being just after the original, i.e., experiences the full hazard
  > #  jump there, in the same way that censors do.)
  > expect <- predict(fit, type='expected', newdata=test1) #force recalc
  > use <- !(test1$time==6 | is.na(test1$status))
  > aeq(test1$status[use] - resid(fit)[use], expect[use])  
  [1] TRUE
  > 
  > sfit <- survfit(fit, list(x=0), censor=FALSE)
  > aeq(sfit$surv, truth$surv)
  [1] TRUE
  > aeq(sfit$std.err^2, truth$var)
  [1] TRUE
  > 
  > # 
  > # Done with the formal test, now print out lots of bits
  > #
  > resid(fit)
           1          2          3          4          5          6          7 
  -0.3655434         NA  0.7191707 -0.2808293 -0.4383414  0.7310869 -0.3655434 
  > resid(fit, 'scor')
           1          2          3          4          5          6          7 
   0.2208584         NA  0.1132780 -0.0442340 -0.1029199 -0.4078409  0.2208584 
  > resid(fit, 'scho')
          1         6         6         9 
   0.157512  0.421244 -0.578756  0.000000 
  > 
  > predict(fit, type='lp')
  [1] -0.8384287         NA  0.8384287  0.8384287  0.8384287 -0.8384287 -0.8384287
  > predict(fit, type='risk')
  [1] 0.4323894        NA 2.3127302 2.3127302 2.3127302 0.4323894 0.4323894
  > predict(fit, type='expected')
          1         2         3         4         5         6         7 
  1.3655434        NA 0.2808293 0.2808293 1.4383414 0.2689131 0.3655434 
  > predict(fit, type='terms')
             x
  1 -0.8384287
  2         NA
  3  0.8384287
  4  0.8384287
  5  0.8384287
  6 -0.8384287
  7 -0.8384287
  > predict(fit, type='lp', se.fit=T)
  $fit
           1          2          3          4          5          6          7 
  -0.8384287         NA  0.8384287  0.8384287  0.8384287 -0.8384287 -0.8384287 
  
  $se.fit
          1         2         3         4         5         6         7 
  0.6388078        NA 0.6388078 0.6388078 0.6388078 0.6388078 0.6388078 
  
  > predict(fit, type='risk', se.fit=T)
  $fit
          1         2         3         4         5         6         7 
  0.4323894        NA 2.3127302 2.3127302 2.3127302 0.4323894 0.4323894 
  
  $se.fit
          1         2         3         4         5         6         7 
  0.4200565        NA 0.9714774 0.9714774 0.9714774 0.4200565 0.4200565 
  
  > predict(fit, type='expected', se.fit=T)
  $fit
          1         2         3         4         5         6         7 
  1.3655434        NA 0.2808293 0.2808293 1.4383414 0.2689131 0.3655434 
  
  $se.fit
  [1] 1.0649293        NA 0.2864593 0.2864593 1.5922983 0.3661617 0.3661617
  
  > predict(fit, type='terms', se.fit=T)
  $fit
             x
  1 -0.8384287
  2         NA
  3  0.8384287
  4  0.8384287
  5  0.8384287
  6 -0.8384287
  7 -0.8384287
  
  $se.fit
            x
  1 0.6388078
  2        NA
  3 0.6388078
  4 0.6388078
  5 0.6388078
  6 0.6388078
  7 0.6388078
  
  > 
  > summary(survfit(fit))
  Call: survfit(formula = fit)
  
   time n.risk n.event survival std.err lower 95% CI upper 95% CI
      1      6       1   0.8857   0.117     0.683036            1
      6      4       2   0.4294   0.237     0.145743            1
      9      1       1   0.0425   0.116     0.000198            1
  > summary(survfit(fit, list(x=2)))
  Call: survfit(formula = fit, newdata = list(x = 2))
  
   time n.risk n.event survival  std.err lower 95% CI upper 95% CI
      1      6       1 2.23e-01 5.97e-01     1.16e-03            1
      6      4       2 2.87e-05 5.69e-04     3.96e-22            1
      9      1       1 1.08e-17 1.04e-15     1.07e-99            1
  > 
  > proc.time()
     user  system elapsed 
     1.37    0.21    1.57 
Running the tests in 'tests/book3.R' failed.
Complete output:
  > library(survival)
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > 
  > #
  > # Tests from the appendix of Therneau and Grambsch
  > #  c. Data set 2 and Breslow estimate
  > #
  > test2 <- data.frame(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
  +                     stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
  +                     event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
  +                     x    =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0))
  > 
  > byhand <- function(beta, newx=0) {
  +     r <- exp(beta)
  +     loglik <- 4*beta - log(r+1) - log(r+2) - 3*log(3*r+2) - 2*log(3*r+1)
  +     u <- 1/(r+1) +  1/(3*r+1) + 4/(3*r+2) -
  +                  ( r/(r+2) +3*r/(3*r+2) + 3*r/(3*r+1))
  +     imat <- r/(r+1)^2 + 2*r/(r+2)^2 + 6*r/(3*r+2)^2 +
  +             3*r/(3*r+1)^2 + 3*r/(3*r+1)^2 + 12*r/(3*r+2)^2
  + 
  +     hazard <-c( 1/(r+1), 1/(r+2), 1/(3*r+2), 1/(3*r+1), 1/(3*r+1), 2/(3*r+2) )
  +     xbar <- c(r/(r+1), r/(r+2), 3*r/(3*r+2), 3*r/(3*r+1), 3*r/(3*r+1),
  +                 3*r/(3*r+2))
  + 
  +     # The matrix of weights, one row per obs, one col per time
  +     #   deaths at 2,3,6,7,8,9
  +     wtmat <- matrix(c(1,0,0,0,1,0,0,0,0,0,
  +                       0,1,0,1,1,0,0,0,0,0,
  +                       0,0,1,1,1,0,1,1,0,0,
  +                       0,0,0,1,1,0,1,1,0,0,
  +                       0,0,0,0,1,1,1,1,0,0,
  +                       0,0,0,0,0,1,1,1,1,1), ncol=6)
  +     wtmat <- diag(c(r,1,1,r,1,r,r,r,1,1)) %*% wtmat
  + 
  +     x      <- c(1,0,0,1,0,1,1,1,0,0)
  +     status <- c(1,1,1,1,1,1,1,0,0,0)
  +     xbar <- colSums(wtmat*x)/ colSums(wtmat)
  +     n <- length(x)
  + 
  +    # Table of sums for score and Schoenfeld resids
  +     hazmat <- wtmat %*% diag(hazard) #each subject's hazard over time
  +     dM <- -hazmat  #Expected part
  +     for (i in 1:6) dM[i,i] <- dM[i,i] +1  #observed
  +     dM[7,6] <- dM[7,6] +1  # observed
  +     mart <- rowSums(dM)
  + 
  +     # Table of sums for score and Schoenfeld resids
  +     #  Looks like the last table of appendix E.2.1 of the book
  +     resid <- dM * outer(x, xbar, '-')
  +     score <- rowSums(resid)
  +     scho <- colSums(resid)
  +     # We need to split the two tied times up, to match coxph
  +     scho <- c(scho[1:5], scho[6]/2, scho[6]/2)
  +     var.g <- cumsum(hazard*hazard /c(1,1,1,1,1,2))
  +     var.d <- cumsum( (xbar-newx)*hazard)
  + 
  +     surv <- exp(-cumsum(hazard) * exp(beta*newx))
  +     varhaz <- (var.g + var.d^2/imat)* exp(2*beta*newx)
  + 
  +     list(loglik=loglik, u=u, imat=imat, xbar=xbar, haz=hazard,
  + 	     mart=mart,  score=score, rmat=resid,
  + 		scho=scho, surv=surv, var=varhaz)
  +     }
  > 
  > 
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > fit0 <-coxph(Surv(start, stop, event) ~x, test2, iter=0, method='breslow')
  > truth0 <- byhand(0,0)
  > aeq(truth0$loglik, fit0$loglik[1])
  [1] TRUE
  > aeq(1/truth0$imat, fit0$var)
  [1] TRUE
  > aeq(truth0$mart, fit0$resid)
  [1] TRUE
  > aeq(truth0$scho, resid(fit0, 'schoen'))
  [1] TRUE
  > aeq(truth0$score, resid(fit0, 'score')) 
  [1] TRUE
  > sfit <- survfit(fit0, list(x=0), censor=FALSE)
  > aeq(sfit$std.err^2, truth0$var)
  [1] TRUE
  > aeq(sfit$surv, truth0$surv)
  [1] TRUE
  > aeq(fit0$score, truth0$u^2/truth0$imat)
  [1] TRUE
  > 
  > beta1 <- truth0$u/truth0$imat
  > fit1 <- coxph(Surv(start, stop, event) ~x, test2, iter=1, ties="breslow")
  > aeq(beta1, coef(fit1))
  [1] TRUE
  > 
  > truth <- byhand(-0.084526081, 0)
  > fit <- coxph(Surv(start, stop, event) ~x, test2, eps=1e-8, method='breslow',
  +              nocenter= NULL)
  > aeq(truth$loglik, fit$loglik[2])
  [1] TRUE
  > aeq(1/truth$imat, fit$var)
  [1] TRUE
  > aeq(truth$mart, fit$resid)
  [1] TRUE
  > aeq(truth$scho, resid(fit, 'schoen'))
  [1] TRUE
  > aeq(truth$score, resid(fit, 'score'))
  [1] TRUE
  > expect <- predict(fit, type='expected', newdata=test2) #force recalc
  > aeq(test2$event -fit$resid, expect) #tests the predict function
  [1] TRUE
  > 
  > sfit <- survfit(fit, list(x=0), censor=FALSE)
  > aeq(sfit$std.err^2, truth$var) 
  [1] TRUE
  > aeq(-log(sfit$surv), (cumsum(truth$haz)))
  [1] TRUE
  > 
  > # Reprise the test, with strata
  > #  offseting the times ensures that we will get the wrong risk sets
  > #  if strata were not kept separate
  > test2b <- rbind(test2, test2, test2)
  > test2b$group <- rep(1:3, each= nrow(test2))
  > test2b$start <- test2b$start + test2b$group
  > test2b$stop  <- test2b$stop  + test2b$group
  > fit0 <- coxph(Surv(start, stop, event) ~ x + strata(group), test2b, 
  +               iter=0, method="breslow")
  > aeq(3*truth0$loglik, fit0$loglik[1])
  [1] TRUE
  > aeq(3*truth0$imat, 1/fit0$var)
  [1] TRUE
  > aeq(rep(truth0$mart,3), fit0$resid)
  [1] TRUE
  > aeq(rep(truth0$scho,3),  resid(fit0, 'schoen'))
  [1] TRUE
  > aeq(rep(truth0$score,3), resid(fit0, 'score')) 
  [1] TRUE
  > 
  > fit1 <- coxph(Surv(start, stop, event) ~ x + strata(group), test2b, 
  +               iter=1, method="breslow")
  > aeq(fit1$coef, beta1)
  [1] TRUE
  > 
  > fit3 <- coxph(Surv(start, stop, event) ~x + strata(group),
  +              test2b, eps=1e-8, method='breslow')
  > aeq(3*truth$loglik, fit3$loglik[2])
  [1] TRUE
  > aeq(3*truth$imat, 1/fit3$var)
  [1] TRUE
  > aeq(rep(truth$mart,3), fit3$resid)
  [1] TRUE
  > aeq(rep(truth$scho,3), resid(fit3, 'schoen'))
  [1] TRUE
  > aeq(rep(truth$score,3), resid(fit3, 'score'))
  [1] TRUE
  > 
  > # 
  > # Done with the formal test, now print out lots of bits
  > #
  > resid(fit)
            1           2           3           4           5           6 
   0.52111895  0.65741078  0.78977654  0.24738772 -0.60629349  0.36902492 
            7           8           9          10 
  -0.06876579 -1.06876579 -0.42044692 -0.42044692 
  > resid(fit, 'scor')
            1           2           3           4           5           6 
   0.27156496 -0.20696709 -0.45771743 -0.09586133  0.13608234  0.19288983 
            7           8           9          10 
   0.04655651 -0.37389040  0.24367131  0.24367131 
  > resid(fit, 'scho')
           2          3          6          7          8          9          9 
   0.5211189 -0.3148216 -0.5795531  0.2661809 -0.7338191  0.4204469  0.4204469 
  > 
  > predict(fit, type='lp')
   [1] -0.04226304  0.04226304  0.04226304 -0.04226304  0.04226304 -0.04226304
   [7] -0.04226304 -0.04226304  0.04226304  0.04226304
  > predict(fit, type='risk')
   [1] 0.9586176 1.0431688 1.0431688 0.9586176 1.0431688 0.9586176 0.9586176
   [8] 0.9586176 1.0431688 1.0431688
  > predict(fit, type='expected')
          1         2         3         4         5         6         7         8 
  0.4788811 0.3425892 0.2102235 0.7526123 1.6062935 0.6309751 1.0687658 1.0687658 
          9        10 
  0.4204469 0.4204469 
  > predict(fit, type='terms')
               x
  1  -0.04226304
  2   0.04226304
  3   0.04226304
  4  -0.04226304
  5   0.04226304
  6  -0.04226304
  7  -0.04226304
  8  -0.04226304
  9   0.04226304
  10  0.04226304
  attr(,"constant")
  [1] -0.04226304
  > predict(fit, type='lp', se.fit=T)
  $fit
            1           2           3           4           5           6 
  -0.04226304  0.04226304  0.04226304 -0.04226304  0.04226304 -0.04226304 
            7           8           9          10 
  -0.04226304 -0.04226304  0.04226304  0.04226304 
  
  $se.fit
          1         2         3         4         5         6         7         8 
  0.3969086 0.3969086 0.3969086 0.3969086 0.3969086 0.3969086 0.3969086 0.3969086 
          9        10 
  0.3969086 0.3969086 
  
  > predict(fit, type='risk', se.fit=T)
  $fit
          1         2         3         4         5         6         7         8 
  0.9586176 1.0431688 1.0431688 0.9586176 1.0431688 0.9586176 0.9586176 0.9586176 
          9        10 
  1.0431688 1.0431688 
  
  $se.fit
          1         2         3         4         5         6         7         8 
  0.3886094 0.4053852 0.4053852 0.3886094 0.4053852 0.3886094 0.3886094 0.3886094 
          9        10 
  0.4053852 0.4053852 
  
  > predict(fit, type='expected', se.fit=T)
  $fit
          1         2         3         4         5         6         7         8 
  0.4788811 0.3425892 0.2102235 0.7526123 1.6062935 0.6309751 1.0687658 1.0687658 
          9        10 
  0.4204469 0.4204469 
  
  $se.fit
   [1] 0.5182381 0.3982700 0.3292830 0.6266797 1.0255146 0.5852364 0.7341340
   [8] 0.7341340 0.6268550 0.6268550
  
  > predict(fit, type='terms', se.fit=T)
  $fit
               x
  1  -0.04226304
  2   0.04226304
  3   0.04226304
  4  -0.04226304
  5   0.04226304
  6  -0.04226304
  7  -0.04226304
  8  -0.04226304
  9   0.04226304
  10  0.04226304
  attr(,"constant")
  [1] -0.04226304
  
  $se.fit
             x
  1  0.3969086
  2  0.3969086
  3  0.3969086
  4  0.3969086
  5  0.3969086
  6  0.3969086
  7  0.3969086
  8  0.3969086
  9  0.3969086
  10 0.3969086
  
  > 
  > summary(survfit(fit))
  Call: survfit(formula = fit)
  
   time n.risk n.event survival std.err lower 95% CI upper 95% CI
      2      2       1    0.607   0.303       0.2279        1.000
      3      3       1    0.437   0.262       0.1347        1.000
      6      5       1    0.357   0.226       0.1034        1.000
      7      4       1    0.277   0.188       0.0729        1.000
      8      4       1    0.214   0.156       0.0514        0.894
      9      5       2    0.143   0.112       0.0308        0.667
  > summary(survfit(fit, list(x=2)))
  Call: survfit(formula = fit, newdata = list(x = 2))
  
   time n.risk n.event survival std.err lower 95% CI upper 95% CI
      2      2       1    0.644   0.444      0.16657            1
      3      3       1    0.482   0.511      0.06055            1
      6      5       1    0.404   0.504      0.03491            1
      7      4       1    0.322   0.475      0.01801            1
      8      4       1    0.258   0.437      0.00928            1
      9      5       2    0.181   0.377      0.00302            1
  > 
  > proc.time()
     user  system elapsed 
     1.64    0.09    1.70 
Running the tests in 'tests/brier.R' failed.
Complete output:
  > # Tests of the Brier score.
  > #  Start with the example in the vignette
  > library(survival)
  > 
  > rott2 <- rotterdam
  > ignore <- with(rott2, recur ==0 & death==1 & rtime < dtime)
  > rott2$rfs <- with(rott2, ifelse(recur==1 | ignore, recur, death))
  > rott2$rfstime <- with(rott2, ifelse(recur==1 | ignore, rtime, dtime))/365.25
  > 
  > rsurv <- survfit(Surv(rfstime, rfs) ~1, rott2)  #KM
  > rfit <- coxph(Surv(rfstime, rfs) ~ pspline(age) + meno + size + pmin(nodes,12),
  +               rott2)
  > 
  > tau <- c(2,4,6, 8)  # four tau values
  > bfit <- brier(rfit, times=tau)
  > 
  > # Now by hand
  > wtmat <- rttright(Surv(rfstime, rfs) ~ 1, rott2, times=tau)
  > psurv <- survfit(rfit, newdata= rott2) # one curve per subject
  > yhat  <- 1- summary(psurv, times=tau)$surv
  > ybar  <- 1- summary(rsurv, times=tau)$surv
  > 
  > y <- with(rott2, cbind(rfstime <=tau[1] & rfs==1,
  +                        rfstime <=tau[2] & rfs==1,
  +                        rfstime <=tau[3] & rfs==1,
  +                        rfstime <=tau[4] & rfs==1)) * 1L
  > ss1 <- colSums(wtmat * (y - t(yhat))^2)
  > ss2 <- colSums(wtmat * (y - rep(ybar, each=nrow(y)))^2)
  > 
  > all.equal(unname(1- ss1/ss2), bfit$rsquare)
  [1] TRUE
  > all.equal(unname(ss1), bfit$brier)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
    13.98    8.00   21.95 
Running the tests in 'tests/checkSurv2.R' failed.
Complete output:
  > library(survival)
  > #
  > # check of the Surv2 function
  > #
  > # Build a flat form of the mgus2 data set.  Mix up the data set order, to test
  > #  out that part of the underlying code.
  > set.seed(1953)
  > m2 <- mgus2[sample(1:nrow(mgus2), nrow(mgus2),replace=FALSE),]
  > 
  > temp1 <- data.frame(m2[,1:7], ftime=0)
  > temp2 <- with(subset(m2, pstat==1), 
  +               data.frame(id=id, ftime=ptime, event="progression"))
  > 
  > # competing risks: use only the first of death and progression
  > temp3 <- with(subset(m2, pstat==0),
  +               data.frame(id=id, ftime=futime, 
  +                                 event=ifelse(death==0, "censor", "death")))
  > mflat <- merge(temp1, rbind(temp2, temp3), all=TRUE)
  > mflat$event <- factor(mflat$event, c("censor", "progression", "death"))
  > 
  > sfit1 <- survfit(Surv2(ftime, event) ~ sex, mflat, id=id)
  > 
  > # now compare it to the usual way
  > etime <- with(mgus2, ifelse(pstat==1, ptime, futime))
  > estat <- with(mgus2, ifelse(pstat==1, 1, 2*death))
  > estat <- factor(estat, 0:2, c("censor", "progression", "death"))
  > sfit2 <- survfit(Surv(etime, estat) ~ sex, mgus2)
  > 
  > all.equal(sfit1$pstate, sfit2$pstate)
  [1] TRUE
  > 
  > # Cox model
  > cfit1 <- coxph(Surv2(ftime, event) ~ sex + age, data=mflat, id=id)
  > cfit2 <- coxph(Surv(etime, estat)  ~ sex + age, data=mgus2, id=id)
  > all.equal(cfit1[c("coefficients", "var", "loglik", "score")],
  +           cfit2[c("coefficients", "var", "loglik", "score")])
  [1] TRUE
  > 
  > # And using the explicit call to build a data set
  > sdata <- Surv2data(Surv2(ftime, event) ~ ., data=mflat, id=id)
  > cfit3 <- coxph(Surv2.y ~ sex + age, data=sdata, id=id)
  > all.equal(cfit1[c("coefficients", "var", "loglik", "score")],
  +           cfit3[c("coefficients", "var", "loglik", "score")])
  [1] TRUE
  > 
  > # Create a data set with error = two events on the same day
  > #  A model with this data will generate an error.
  > temp4 <- with(m2, 
  +               data.frame(id=id, ftime=futime, 
  +                                 event=ifelse(death==0, "censor", "death")))
  > mflat2 <- merge(temp1, rbind(temp2, temp4), all=TRUE)
  > mflat2$event <- factor(mflat2$event, c("censor", "prog", "death"))
  > stemp <- survcheck(Surv2(ftime, event) ~ sex, data=mflat2, id=id)
  > all.equal(stemp$duplicate$row, which(duplicated(mflat2[,c("id", "ftime")])))
  [1] TRUE
  > 
  > # Full 3 state model.  We need to make progressions that are tied with
  > #  deaths be just a bit sooner.
  > temp2b <- with(subset(m2, pstat==1),
  +                data.frame(id=id, ftime= ifelse(ptime==futime & death==1, ptime-.1, ptime),
  +                event="progression"))
  > temp3b <- with(m2,
  +                data.frame(id=id, ftime=futime, event=ifelse(death==0, "censor", "death")))
  > 
  > mflat3 <- merge(temp1, rbind(temp2b, temp3b), all=TRUE)
  > mflat3$event <- factor(mflat3$event, c("censor", "progression", "death"))
  > 
  > cfit4 <- coxph(Surv2(ftime, event) ~ sex + age + mspike, mflat3, id=id)
  > 
  > # For a standard start-stop data set use tmerge
  > m3 <- tmerge(m2[,1:7], subset(m2,,c(id, futime, death)), id=id, 
  +                                   event= event(futime, 2*death))
  > m3 <- tmerge(m3, temp2b, id=id, event= event(ftime))
  > m3$event <- factor(m3$event, 0:2, c("censor", "progression", "death"))
  > cfit5 <- coxph(Surv(tstart, tstop, event) ~ sex + age + mspike, m3, id=id)
  > 
  > all.equal(cfit4[c("coefficients", "var", "loglik", "score")],
  +           cfit5[c("coefficients", "var", "loglik", "score")])
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.90    0.31    2.20 
Running the tests in 'tests/concordance.R' failed.
Complete output:
  > library(survival)
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > 
  > #
  > # Simple tests of concordance
  > #
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > grank <- function(x, time, grp, wt) 
  +     unlist(tapply(x, grp, rank))
  > grank2 <- function(x, time, grp, wt) {  #for case weights
  +     if (length(wt)==0) wt <- rep(1, length(x))
  +     z <- double(length(x))
  +     for (i in unique(grp)) {
  +         indx <- which(grp==i)
  +         temp <- tapply(wt[indx], x[indx], sum)
  +         temp <- temp/2  + c(0, cumsum(temp)[-length(temp)])
  +         z[indx] <- temp[match(x[indx], names(temp))]
  +     }
  +     z
  + }
  > 
  > # Concordance by brute force.  O(n^2) algorithm, but ok for n<500 or so
  > allpair <- function(x, time, status, wt, all=FALSE) {
  +     if (missing(wt)) wt <- rep(1, length(x))
  +     count <- sapply(which(status==1), function(i) {
  +         atrisk <- (time > time[i]) | (time==time[i] & status==0)
  +         temp <- tapply(wt[atrisk], factor(sign(x[i] -x[atrisk]), c(1, -1, 0)),
  +                        sum)
  +         wt[i]* c(ifelse(is.na(temp), 0, temp),
  +                  (sum(wt[time==time[i] & status==1]) - wt[i])/2) 
  +     })
  +     rownames(count) <- c("concordant", "discordant", "tied.x", "tied.y")
  +     if (all) {
  +         colnames(count) <- time[status==1]
  +         t(count)
  +     }
  +     else rowSums(count)
  + }
  > 
  > 
  > # The std of C = std(numerator)/(number of comparable pairs)
  > # The information matrix of a Cox model is = to the var(C-D)
  > cfun <- function(fit) fit$cvar * sum(fit$count[1:3])^2
  >                 
  > tdata <- aml[aml$x=='Maintained', c("time", "status")]
  > tdata$x <- c(1,6,2,7,3,7,3,8,4,4,5)
  > tdata$wt <- c(1,2,3,2,1,2,3,4,3,2,1)
  > fit <- concordance(Surv(time, status) ~x, tdata)
  > 
  > aeq(fit$count[1:4], c(24,14,2,0))
  [1] TRUE
  > cfit <- coxph(Surv(time, status) ~ tt(x), tdata, tt=grank, method='breslow',
  +               iter=0, x=T)
  > cdt <- coxph.detail(cfit)
  > aeq(sum(cdt$imat), cfun(fit))
  [1] TRUE
  > aeq(sum(2*cdt$score), diff(fit$count[1:2]))
  [1] TRUE
  > aeq(with(tdata, allpair(x, time, status)), c(14,24,2,0))
  [1] TRUE
  > 
  > # Lots of ties
  > tempy <- Surv(c(1,2,2,2,3,4,4,4,5,2), c(1,0,1,0,1,0,1,1,0,1))
  > tempx <- c(5,5,4,4,3,3,7,6,5,4)
  > fit2 <- concordance(tempy ~ tempx)
  > addxy <- function(x) c(x[1:3], sum(x[4:5]))
  > aeq(addxy(fit2$count), allpair(tempx, tempy[,1], tempy[,2]))
  [1] TRUE
  > cfit2 <-  coxph(tempy ~ tt(tempx), tt=grank, method='breslow', iter=0)
  > aeq(cfit2$var, 1/cfun(fit2))
  [1] TRUE
  > 
  > # Direct call
  > fit2b <- concordancefit(tempy, tempx)
  > fit2c <- concordancefit(tempy, tempx, std.err=FALSE)
  > all.equal(fit2[1:5], fit2b)
  [1] TRUE
  > all.equal(fit2b[1:3], fit2c)
  [1] TRUE
  > 
  > # Bigger data
  > fit3 <- concordance(Surv(time, status) ~ age, lung, reverse=TRUE)
  > aeq(addxy(fit3$count), allpair(lung$age, lung$time, lung$status-1))
  [1] TRUE
  > cfit3 <- coxph(Surv(time, status) ~ tt(age), lung, 
  +                iter=0, method='breslow', tt=grank, x=T)
  > cdt <- coxph.detail(cfit3)
  > aeq(sum(cdt$imat), cfun(fit3)) 
  [1] TRUE
  > aeq(2*sum(cdt$score), diff(fit3$count[2:1]))
  [1] TRUE
  > 
  > 
  > # More ties
  > fit4 <- concordance(Surv(time, status) ~ ph.ecog, lung, reverse=TRUE)
  > aeq(addxy(fit4$count), allpair(lung$ph.ecog, lung$time, lung$status-1))
  [1] TRUE
  > aeq(fit4$count[1:5], c(8392, 4258, 7137, 21, 7))
  [1] TRUE
  > cfit4 <- coxph(Surv(time, status) ~ tt(ph.ecog), lung, 
  +                iter=0, method='breslow', tt=grank)
  > aeq(1/cfit4$var, cfun(fit4))
  [1] TRUE
  > 
  > # Case weights
  > fit5 <- concordance(Surv(time, status) ~ x, tdata, weight=wt, reverse=TRUE)
  > fit6 <- concordance(Surv(time, status) ~x, tdata[rep(1:11,tdata$wt),])
  > aeq(addxy(fit5$count), with(tdata, allpair(x, time, status, wt)))
  [1] TRUE
  > aeq(fit5$count[1:4], c(70, 91, 7, 0))  # checked by hand
  [1] TRUE
  > aeq(fit5$count[1:3], fit6$count[c(2,1,3)])  #spurious "tied on time" values, ignore
  [1] TRUE
  > aeq(fit5$std, fit6$std)
  [1] TRUE
  > cfit5 <- coxph(Surv(time, status) ~ tt(x), tdata, weight=wt, 
  +                iter=0, method='breslow', tt=grank2)
  > cfit6 <- coxph(Surv(time, status) ~ tt(x), tdata[rep(1:11,tdata$wt),], 
  +                iter=0, method='breslow', tt=grank)
  > aeq(1/cfit6$var, cfun(fit6))
  [1] TRUE
  > aeq(cfit5$var, cfit6$var)
  [1] TRUE
  > 
  > # Start, stop simplest cases
  > fit7 <- concordance(Surv(rep(0,11), time, status) ~ x, tdata)
  > aeq(fit7$count, fit$count)
  [1] TRUE
  > aeq(fit7$std.err, fit$std.err)
  [1] TRUE
  > fit7 <- concordance(Surv(rep(0,11), time, status) ~ x, tdata, weight=wt)
  > aeq(fit5$count, fit7$count[c(2,1,3:5)])  #one reversed, one not
  [1] TRUE
  > 
  > # Multiple intervals for some, but same risk sets as tdata
  > tdata2 <- data.frame(time1=c(0,3, 5,  6,7,   0,  4,17,  7,  0,16,  2,  0, 
  +                              0,9, 5),
  +                      time2=c(3,9, 13, 7,13, 18, 17,23, 28, 16,31, 34, 45, 
  +                              9,48, 60),
  +                      status=c(0,1, 1, 0,0,  1,  0,1, 0, 0,1, 1, 0, 0,1, 0),
  +                      x = c(1,1, 6, 2,2, 7, 3,3, 7, 3,3, 8, 4, 4,4, 5),
  +                      wt= c(1,1, 2, 3,3, 2, 1,1, 2, 3,3, 4, 3, 2,2, 1))
  > fit8 <- concordance(Surv(time1, time2, status) ~x, tdata2, weight=wt,
  +                     reverse=TRUE)
  > aeq(fit5$count, fit8$count)
  [1] TRUE
  > aeq(fit5$std.err, fit8$std.err)
  [1] TRUE
  > cfit8 <- coxph(Surv(time1, time2, status) ~ tt(x), tdata2, weight=wt, 
  +                iter=0, method='breslow', tt=grank2)
  > aeq(1/cfit8$var, cfun(fit8))
  [1] TRUE
  > 
  > # Stratified
  > tdata3 <- data.frame(time1=c(tdata2$time1, rep(0, nrow(lung))),
  +                      time2=c(tdata2$time2, lung$time),
  +                      status = c(tdata2$status, lung$status -1),
  +                      x = c(tdata2$x, lung$ph.ecog),
  +                      wt= c(tdata2$wt, rep(1, nrow(lung))),
  +                      grp=rep(1:2, c(nrow(tdata2), nrow(lung))))
  > fit9 <- concordance(Surv(time1, time2, status) ~x + strata(grp),
  +                         data=tdata3, weight=wt, reverse=TRUE)
  > aeq(fit9$count[1,], fit5$count)
  [1] TRUE
  > aeq(fit9$count[2,], fit4$count)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.62    0.20    1.81 
Running the tests in 'tests/concordance2.R' failed.
Complete output:
  > library(survival)
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > 
  > #
  > # Tests for the condordance function. 
  > #
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > grank <- function(x, time, grp, wt) {
  +     if (all(wt==1)) unlist(tapply(x, grp, rank))
  +     else unlist(tapply(1:length(x), grp, function(i) {
  +         xx <- x[i]  # x and wts for this subset of the data
  +         ww <- wt[i]
  +         temp <- outer(xx, xx, function(a, b) sign(b-a))
  +         colSums(ww*temp)/2
  +     }))
  + }
  > 
  > # a Cox model using iter=0, ties='breslow' and the above function has a score 
  > #  statistic which is U=(C-D)/2 and score test U^2/H, where H is the Cox model 
  > #  information matrix, with fit$var=1/H.  The concordance is U+ 1/2. 
  > # Pull out the Somers' d and its variance
  > phget <- function(fit) {
  +     c(d =  2*sqrt(fit$score/fit$var), v= 4/fit$var)
  + }
  > fscale <- function(fit) {
  +     if (is.matrix(fit$count)) temp <- colSums(fit$count) else temp <- fit$count
  +     npair <- sum(temp[1:3])
  +     c(d = abs(temp[1]-temp[2]), v=4*fit$cvar*npair^2)
  + }
  > 
  > # Concordance by brute force.  O(n^2) algorithm, but ok for n<500 or so
  > allpair <- function(time, status, x, wt, all=FALSE) {
  +     n <- length(time)
  +     if (missing(wt)) wt <- rep(1, length(x))
  +     count <- sapply(which(status==1), function(i) {
  +         atrisk <- (time > time[i]) | (time==time[i] & status==0)
  +         temp <- tapply(wt[atrisk], factor(sign(x[i] -x[atrisk]), c(-1, 1, 0)),
  +                        sum)
  +         tiedtime <- (time==time[i] & status ==1 & (1:n)>i)
  +         ties <- tapply(wt[tiedtime], factor(x[tiedtime]==x[i], 
  +                                             c(FALSE, TRUE)),sum)
  +         wt[i]* c(ifelse(is.na(temp), 0, temp), ifelse(is.na(ties), 0, ties))
  +     })
  +     rownames(count) <- c("concordant", "discordant", "tied.x", "tied.y",
  +                          "tied.xy")
  +     if (all) {
  +         colnames(count) <- time[status==1]
  +         t(count)
  +     }
  +     else rowSums(count)
  + }
  > 
  > # leverage by brute force
  > leverage <- function(time, status, x, wt, eps=1e-5) {
  +     if (missing(wt)) wt <- rep(1, length(x))
  +     toss <- is.na(time + status + x +wt)
  +     if (any(toss)) {
  +         time <- time[!toss]
  +         status <- status[!toss]
  +         x <- x[!toss]
  +         wt <- wt[!toss]
  +     }
  +     n <- length(time)
  +     influence <- matrix(0, n, 5)
  +     t2 <- time + eps*(status==0)
  +     for (i in 1:n) {
  +         if (status[i] ==0) comparable <- (time<=time[i] & status==1)
  +         else comparable <- ifelse(status==0, time >= time[i], time!=time[i])
  +         temp <- sign((x[i]-x[comparable])*(t2[i] - t2[comparable]))
  +         influence[i,1:3] <-tapply(wt[comparable],factor(temp, c(1,-1,0)), sum)
  +         if (status[i]==1) {
  +             tied <- (time==time[i] & status==1 & (1:n)!= i)
  +             if (any(tied)) {
  +                 itemp<- tapply(wt[tied], factor(x[tied]==x[i], 
  +                                                 c(FALSE, TRUE)), sum)
  +                 influence[i,4:5] <- itemp
  +             }
  +         }
  +     }
  +     dimnames(influence) <- list(as.character(Surv(time, status)), 
  +                                 c("concord", "discord", "tie.x", "tie.y",
  +                                   "tie.xy"))
  +     ifelse(is.na(influence), 0, influence)
  + }
  > 
  > tdata <- aml[aml$x=='Maintained', c("time", "status")]
  > tdata$x <- c(1,6,2,7,3,7,3,8,4,4,5)
  > tdata$wt <- c(1,2,3,2,1,2,3,4,3,2,1)
  > 
  > fit <- concordance(Surv(time, status) ~x, tdata, influence=2)
  > aeq(fit$count, with(tdata, allpair(time, status, x)))
  [1] TRUE
  > aeq(fit$influence, with(tdata, leverage(time, status, x)))
  [1] TRUE
  > 
  > cfit <- coxph(Surv(time, status) ~ tt(x), tdata, tt=grank, ties='breslow',
  +               iter=0, x=T)
  > aeq(phget(cfit), fscale(fit))  # agree with Cox model
  [1] TRUE
  > 
  > # Test 2: Lots of ties
  > tempy <- Surv(c(1,2,2,2,3,4,4,4,5,2), c(1,0,1,0,1,0,1,1,0,1))
  > tempx <- c(5,5,4,4,3,3,7,6,5,4)
  > fit2 <- concordance(tempy ~ tempx, influence=2)
  > aeq(fit2$count, allpair(tempy[,1], tempy[,2], tempx))
  [1] TRUE
  > aeq(fit2$influence, leverage(tempy[,1], tempy[,2], tempx))
  [1] TRUE
  > cfit2 <- coxph(tempy ~ tt(tempx), tt=grank, ties="breslow", iter=0)
  > aeq(phget(cfit2), fscale(fit2))  # agree with Cox model
  [1] TRUE
  > 
  > # Bigger data
  > cox3 <- coxph(Surv(time, status) ~ age + sex + ph.ecog, lung)
  > fit3 <- concordance(Surv(time, status) ~ predict(cox3), lung, influence=2)
  > aeq(fit3$count, allpair(lung$time, lung$status-1,predict(cox3)))
  [1] TRUE
  > aeq(fit3$influence, leverage(lung$time, lung$status-1,predict(cox3)))
  [1] TRUE
  > cfit3 <- coxph(Surv(time, status) ~ tt(predict(cox3)), tt=grank,
  +                ties="breslow", iter=0, data=lung)
  > aeq(phget(cfit3), fscale(fit3))  # agree with Cox model
  [1] TRUE
  > 
  > # More ties
  > fit4  <- concordance(Surv(time, status) ~ ph.ecog, lung, influence=2)
  > fit4b <- concordance(Surv(time, status) ~ ph.ecog, lung, reverse=TRUE)
  > aeq(fit4$count, allpair(lung$time, lung$status-1, lung$ph.ecog))
  [1] TRUE
  > aeq(fit4b$count, c(8392, 4258, 7137, 21, 7))
  [1] TRUE
  > cfit4 <- coxph(Surv(time, status) ~ tt(ph.ecog), lung, 
  +                iter=0, method='breslow', tt=grank)
  > aeq(phget(cfit4), fscale(fit4))  # agree with Cox model
  [1] TRUE
  > 
  > # Case weights
  > fit5 <- concordance(Surv(time, status) ~ x, tdata, weight=wt, influence=2)
  > fit6 <- concordance(Surv(time, status) ~x, tdata[rep(1:11,tdata$wt),])
  > aeq(fit5$count, with(tdata, allpair(time, status, x, wt)))
  [1] TRUE
  > aeq(fit5$count, c(91, 70, 7, 0, 0))  # checked by hand
  [1] TRUE
  > aeq(fit5$count[1:3], fit6$count[1:3])  #spurious "tied.xy" values, ignore
  [1] TRUE
  > aeq(fit5$var[2], fit6$var[2])
  [1] TRUE
  > aeq(fit5$influence, with(tdata, leverage(time, status, x, wt)))
  [1] TRUE
  > cfit5 <- coxph(Surv(time, status) ~ tt(x), tdata, weight=wt, 
  +                iter=0, method='breslow', tt=grank)
  > aeq(phget(cfit5), fscale(fit5))  # agree with Cox model
  [1] TRUE
  > 
  > # Start, stop simplest cases
  > fit6 <- concordance(Surv(rep(0,11), time, status) ~ x, tdata)
  > aeq(fit6$count, fit$count)
  [1] TRUE
  > aeq(fit6$var, fit$var)
  [1] TRUE
  > fit7 <- concordance(Surv(rep(0,11), time, status) ~ x, tdata, weight=wt)
  > aeq(fit7$count, fit5$count)
  [1] TRUE
  > aeq(fit7$var, fit5$var)
  [1] TRUE
  > 
  > # Multiple intervals for some, but same risk sets as tdata
  > tdata2 <- data.frame(time1=c(0,3, 5,  6,7,   0,  4,17,  7,  0,16,  2,  0, 
  +                              0,9, 5),
  +                      time2=c(3,9, 13, 7,13, 18, 17,23, 28, 16,31, 34, 45, 
  +                              9,48, 60),
  +                      status=c(0,1, 1, 0,0,  1,  0,1, 0, 0,1, 1, 0, 0,1, 0),
  +                      x = c(1,1, 6, 2,2, 7, 3,3, 7, 3,3, 8, 4, 4,4, 5),
  +                      wt= c(1,1, 2, 3,3, 2, 1,1, 2, 3,3, 4, 3, 2,2, 1),
  +                      id= c(1,1, 2, 3,3, 4, 5,5, 6, 7,7, 8, 9, 10,10, 11))
  > fit8 <- concordance(Surv(time1, time2, status) ~x, cluster=id, tdata2, 
  +                     weight=wt, influence=2)
  > aeq(fit5$count, fit8$count)
  [1] TRUE
  > # influence has one row per obs, so the next line is false: mismatched lengths
  > # aeq(fit5$influence, fit8$influence) 
  > aeq(fit5$var, fit8$var)
  [1] TRUE
  > cfit8 <- coxph(Surv(time1, time2, status) ~ tt(x), tdata2, weight=wt, 
  +                iter=0, method='breslow', tt=grank)
  > aeq(phget(cfit8), fscale(fit8))  # agree with Cox model
  [1] TRUE
  > 
  > # Stratified
  > tdata3 <- data.frame(time1=c(tdata2$time1, rep(0, nrow(lung))),
  +                      time2=c(tdata2$time2, lung$time),
  +                      status = c(tdata2$status, lung$status -1),
  +                      x = c(tdata2$x, lung$ph.ecog),
  +                      wt= c(tdata2$wt, rep(1, nrow(lung))),
  +                      grp=rep(1:2, c(nrow(tdata2), nrow(lung))),
  +                      id = c(tdata2$id, 100+ 1:nrow(lung)))
  > fit9 <- concordance(Surv(time1, time2, status) ~x + strata(grp), cluster=id,
  +                         data=tdata3, weight=wt, influence=2)
  > aeq(fit9$count, rbind(fit8$count, fit4$count))
  [1] TRUE
  > 
  > # check out case weights, strata, and grouped jackknife;
  > #   force several ties in x, y, and xy (with missing values too for good measure).
  > tdata <- subset(lung, select=-c(meal.cal, wt.loss, sex, age))
  > tdata$wt <- rep(1:25, length=nrow(tdata))/10
  > tdata$time <- ceiling(tdata$time/30)  # force ties in y
  > tfit <- coxph(Surv(time, status) ~ ph.ecog + pat.karno + strata(inst)
  +               + cluster(inst), tdata, weight=wt)
  > tdata$tpred <- predict(tfit)
  > cm4 <- concordance(tfit, influence=3, keepstrata=TRUE)
  > cm5 <- concordance(Surv(time, status) ~ tpred + strata(inst) + cluster(inst),
  +                    data=tdata, weight=wt, reverse=TRUE, influence=3,
  +                    keepstrata=TRUE)
  > all.equal(cm4[1:6], cm5[1:6])  # call and na.action won't match
  [1] TRUE
  > 
  > u.inst <- sort(unique(tdata$inst))
  > temp <- matrix(0, length(u.inst), 5)
  > for (i in 1:length(u.inst)) {
  +     temp[i,] <- with(subset(tdata, inst==u.inst[i]),        
  +                  allpair(time, status-1, -tpred, wt))
  + }
  > aeq(temp, cm4$count)
  [1] TRUE
  >     
  > eps <- 1e-6
  > keep <- (1:nrow(tdata))[-tfit$na.action]  # the obs that are not tossed
  > lmat <- matrix(0., length(keep), 5)
  > for (i in 1:length(keep)) {
  +     wt2 <- tdata$wt
  +     wt2[keep[i]] <- wt2[keep[i]] + eps
  + 
  +     test <- concordance(Surv(time, status) ~ predict(tfit) + strata(inst),
  +                    data=tdata, weight=wt2, group=group, reverse=TRUE,
  +                    keepstrata=TRUE)
  +     lmat[i,] <- colSums(test$count - cm4$count)/eps
  + }
  > aeq(lmat, cm4$influence, tolerance=eps)
  [1] TRUE
  > 
  > # Check that keepstrata gives the correct sum
  > cm4b <- concordance(tfit, keepstrata=FALSE)
  > aeq(cm4b$count, colSums(cm4$count))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     4.70    0.20    4.89 
Running the tests in 'tests/concordance3.R' failed.
Complete output:
  > library(survival)
  > options(na.action=na.exclude)
  > aeq <- function(x,y,...)  all.equal(as.vector(x), as.vector(y), ...)
  > 
  > #  Make sure strata is retained, and that the overall variance is correct
  > fit1 <- coxph(Surv(time, status) ~ age + offset(ph.ecog*0) +strata(sex), lung)
  > fit2 <- coxph(Surv(time, status) ~ age + ph.ecog +strata(sex), lung)
  > 
  > test <- concordance(fit1, fit2, influence=1)
  > 
  > ksex <- model.frame(fit1)[["strata(sex)"]]
  > test1 <- concordance(fit1$y ~ fit1$linear.predictors + strata(ksex), 
  +                      reverse=TRUE, influence=1)
  > test2 <- concordance(fit1$y ~ fit2$linear.predictors + strata(ksex), 
  +                      reverse=TRUE, influence=1)
  > aeq(test$concordance, c(test1$concordance, test2$concordance))
  [1] TRUE
  > aeq(diag(test$var), c(test1$var[1], test2$var[1]))
  [1] TRUE
  > aeq(test$dfbeta, cbind(test1$dfbeta, test2$dfbeta))
  [1] TRUE
  > 
  > cvec <- c(-1, 1)
  > aeq(cvec %*% test$var %*% cvec, sum((test1$dfbeta - test2$dfbeta)^2))
  [1] TRUE
  > 
  > # Time weights
  > # Start with a very small data set: aml has 23 subjects
  > #
  > atest1 <- concordance(Surv(time, status) ~ x, aml, ranks=TRUE)
  > atest2 <- concordance(Surv(time, status) ~ x, aml, ranks=TRUE, timewt='S')
  > atest3 <- concordance(Surv(time, status) ~ x, aml, ranks=TRUE, timewt='S/G')
  > atest4 <- concordance(Surv(time, status) ~ x, aml, ranks=TRUE, timewt='n/G2')
  > # The ranks data frame agrees for all but weights
  > all.equal(atest1$ranks[, -3], atest2$ranks[, -3])
  [1] TRUE
  > all.equal(atest1$ranks[, -3], atest3$ranks[, -3])
  [1] TRUE
  > all.equal(atest1$ranks[, -3], atest4$ranks[, -3])
  [1] TRUE
  > 
  > wt1 <- cbind(atest1$ranks[,"timewt"], atest2$ranks[,"timewt"],
  +              atest3$ranks[,"timewt"], atest4$ranks[,"timewt"])
  >             
  > # survfit0 adds time 0 to the curves
  > # to break ties between censor/death for G, we need to add an offset to
  > #  the censoring times.  Since time is integer, .1 works nicely
  > s1 <- survfit0(survfit(Surv(time, status) ~ 1, aml))
  > g1 <- survfit0(survfit(Surv(time + .1*(1-status), 1-status) ~1, aml))
  > 
  > # The ingredients of the weights
  > indx <- match(atest1$ranks[,"time"], s1$time)
  > nrisk  <- s1$n.risk[indx]
  > sminus <- s1$surv[indx-1]
  > gminus <- g1$surv[findInterval(atest1$ranks[,"time"], g1$time)]
  > n <- nrow(aml)
  > 
  > wt2 <- cbind(nrisk, n*sminus, n*sminus/gminus, nrisk/gminus^2)
  > aeq(wt1, wt2)
  [1] TRUE
  > 
  > # The sum of weighted ranks should equal (C-D) for a Cox model fit
  > tfun <- function(cfit, reverse=FALSE) {
  +     t1 <- sum(cfit$ranks$timewt * cfit$ranks$rank)
  +     t2 <- cfit$count[1] - cfit$count[2]
  +     all.equal(unname(t1), unname(t2))
  + }
  > tfun(atest1)
  [1] TRUE
  > tfun(atest2)
  [1] TRUE
  > tfun(atest3)
  [1] TRUE
  > tfun(atest4)
  [1] TRUE
  > 
  > # The nafld data set has strong and early censoring (one of the only ones
  > # in the package that does.) So it is a good check of time weights.
  > #
  > nfit <- coxph(Surv(futime, status) ~ male + pspline(age), nafld1)
  > cn1 <- concordance(nfit, timewt='n', ranks=TRUE)
  > cn2 <- concordance(nfit, timewt='S', ranks=TRUE)
  > cn3 <- concordance(nfit, timewt='S/G', ranks=TRUE)
  > cn4 <- concordance(nfit, timewt='n/G2', ranks=TRUE)
  > 
  > sfit <- survfit0(survfit(Surv(futime, status) ~ 1, nafld1))
  > gfit <- survfit0(survfit(Surv(futime + .1*(status==0), 1-status) ~0, nafld1))
  > 
  > # The ingredients of the weights
  > dtime <- cn1$ranks[, "time"]
  > indx <- match(dtime, sfit$time)
  > nrisk  <- sfit$n.risk[indx]
  > sminus <- sfit$surv[indx-1]
  > gminus <- gfit$surv[findInterval(dtime, gfit$time)]
  > n <- nrow(nafld1)
  > 
  > wt1 <- cbind(cn1$ranks[, "timewt"], cn2$ranks[,"timewt"],
  +              cn3$ranks[, "timewt"], cn4$ranks[,"timewt"])
  > wt2 <- cbind(nrisk, n*sminus, n*sminus/gminus, nrisk/gminus^2)
  > aeq(wt1, wt2)
  [1] TRUE
  > 
  > rd1 <- cn1$ranks
  > rd2 <- cn2$ranks
  > rd3 <- cn3$ranks
  > all.equal(rd1[c('time', 'rank', 'casewt')], rd2[c('time', 'rank', 'casewt')])
  [1] TRUE
  > all.equal(rd1[c('time', 'rank', 'casewt')], rd3[c('time', 'rank', 'casewt')])
  [1] TRUE
  > 
  > tfun(cn1)
  [1] TRUE
  > tfun(cn2)
  [1] TRUE
  > tfun(cn3)
  [1] TRUE
  > tfun(cn4)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.82    0.25    2.01 
Running the tests in 'tests/coxsurv.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # Test out subscripting in the case of a coxph survival curve
  > #
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > fit <- coxph(Surv(time, status) ~ age + sex + meal.cal + strata(ph.ecog),
  + 		data=lung)
  > surv1 <- survfit(fit)
  > temp <- surv1[2:3]
  > 
  > which <- cumsum(surv1$strata)
  > zed   <- (which[1]+1):(which[3])
  > aeq(surv1$surv[zed], temp$surv)
  [1] TRUE
  > aeq(surv1$time[zed], temp$time)
  [1] TRUE
  > 
  > # This call should not create a model frame in the code -- so same
  > #  answer but a different path through the underlying code
  > fit <- coxph(Surv(time, status) ~ age + sex + meal.cal + strata(ph.ecog),
  + 		x=T, data=lung)
  > surv2 <- survfit(fit)
  > all.equal(surv1, surv2)
  [1] TRUE
  > 
  > #
  > # Now a result with a matrix of survival curves
  > #
  > dummy <- data.frame(age=c(30,40,60), sex=c(1,2,2), meal.cal=c(500, 1000, 1500))
  > surv2 <- survfit(fit, newdata=dummy)
  > 
  > zed <- 1:which[1]
  > aeq(surv2$surv[zed,1], surv2[1,1]$surv)
  [1] TRUE
  > aeq(surv2$surv[zed,2], surv2[1,2]$surv)
  [1] TRUE
  > aeq(surv2$surv[zed,3], surv2[1,3]$surv)
  [1] TRUE
  > aeq(surv2$surv[zed, ], surv2[1,1:3]$surv)
  [1] TRUE
  > aeq(surv2$surv[zed], (surv2[1])$surv)
  [1] TRUE
  > aeq(surv2$surv[zed, ], surv2[1, ]$surv)
  [1] TRUE
  > 
  > # And the depreciated form - call with a named vector as 'newdata'
  > #  the resulting $call component  won't match so delete it before comparing
  > # newdata will have mismatched row names due to subscripting
  > surv3 <- survfit(fit, c(age=40, sex=2, meal.cal=1000))
  > keep <- which(!(names(surv3) %in% c("newdata", "call")))
  > all.equal(unclass(surv2[,2])[keep], unclass(surv3)[keep])
  [1] TRUE
  > 
  > # Test out offsets, which have recently become popular due to a Langholz paper
  > fit1 <- coxph(Surv(time, status) ~ age + ph.ecog, lung)
  > fit2 <- coxph(Surv(time, status) ~ age + offset(ph.ecog * fit1$coef[2]), lung)
  >  
  > surv1 <- survfit(fit1, newdata=data.frame(age=50, ph.ecog=1))
  > surv2 <- survfit(fit2, newdata=data.frame(age=50, ph.ecog=1))
  > all.equal(surv1$surv, surv2$surv)
  [1] TRUE
  > 
  > # And a model with only offsets.
  > eta <- cbind(lung$age, lung$ph.ecog) %*% coef(fit1)
  > fit3 <- coxph(Surv(time, status) ~ offset(eta), lung)
  > aeq(fit3$log, fit1$log[2])
  [1] TRUE
  > 
  > surv3 <- survfit(fit3, newdata=data.frame(eta= 50*fit1$coef[1] + fit1$coef[2]))
  > all.equal(surv3$surv, surv1$surv)
  [1] TRUE
  > 
  > #
  > # Check out the start.time option
  > #
  > surv3 <- survfit(fit1, newdata=data.frame(age=50, ph.ecog=1),
  +                  start.time=100)
  > index <- match(surv3$time, surv1$time)
  > rescale <- summary(surv1, time=100)$surv
  > all.equal(surv3$surv, surv1$surv[index]/rescale)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.32    0.12    1.39 
Running the tests in 'tests/coxsurv2.R' failed.
Complete output:
  > library(survival)
  > #
  > # Check that the survival curves from a Cox model with beta=0 
  > #  match ordinary survival
  > #
  > #  Aalen
  > surv1 <- survfit(Surv(time,status) ~ sex, data=lung, stype=2)
  > fit1 <- coxph(Surv(time, status) ~ age + strata(sex), data=lung, iter=0,
  +               ties='breslow')
  > fit1$var <- 0*fit1$var   #sneaky, causes the extra term in the Cox variance
  >                          # calculation to be zero
  > surv2 <- survfit(fit1, stype=2)
  > surv3 <- survfit(fit1)
  > 
  > arglist <- c('n', 'time', 'n.risk','n.event', 'n.censor', 'surv', 'strata',
  +              'std.err', 'upper', 'lower')
  > all.equal(unclass(surv1)[arglist], unclass(surv2)[arglist])
  [1] TRUE
  > all.equal(unclass(surv1)[arglist], unclass(surv3)[arglist])
  [1] TRUE
  > 
  > 
  > # Efron method
  > surv1 <- survfit(Surv(time,status) ~ sex, data=lung, stype=2, ctype=2)
  > surv2 <- survfit(fit1, ctype=2)
  > all.equal(unclass(surv1)[arglist], unclass(surv2)[arglist])
  [1] TRUE
  > 
  > # Kaplan-Meier
  > surv1 <- survfit(Surv(time,status) ~ sex, data=lung)
  > surv2 <- survfit(fit1, stype=1)
  > all.equal(unclass(surv1)[arglist], unclass(surv2)[arglist])
  [1] TRUE
  > 
  > 
  > # Now add some random weights
  > rwt <- runif(nrow(lung), .5, 3)
  > surv1 <- survfit(Surv(time,status) ~ sex, data=lung, stype=2, weight=rwt,
  +                  robust=FALSE)
  > fit1 <- coxph(Surv(time, status) ~ age + strata(sex), data=lung, iter=0,
  +               ties='breslow', weight=rwt, robust=FALSE)
  > fit1$var <- 0*fit1$var   #sneaky
  > surv2 <- survfit(fit1, stype=2, ctype=1)
  > surv3 <- survfit(fit1)
  > 
  > all.equal(unclass(surv1)[arglist], unclass(surv2)[arglist])
  [1] TRUE
  > all.equal(unclass(surv1)[arglist], unclass(surv3)[arglist])
  [1] TRUE
  > 
  > 
  > # Efron method
  > surv1 <- survfit(Surv(time,status) ~ sex, data=lung, stype=2, ctype=2,
  +                  weight=rwt, robust=FALSE)
  > surv2 <- survfit(fit1, ctype=2, robust=FALSE)
  > all.equal(unclass(surv1)[arglist], unclass(surv2)[arglist])
  [1] TRUE
  > 
  > # Kaplan-Meier
  > surv1 <- survfit(Surv(time,status) ~ sex, data=lung, weight=rwt, robust=FALSE)
  > surv2 <- survfit(fit1, stype=1, robust=FALSE)
  > all.equal(unclass(surv1)[arglist], unclass(surv2)[arglist])
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.40    0.18    1.56 
Running the tests in 'tests/coxsurv4.R' failed.
Complete output:
  > library(survival)
  > 
  > # Strata by covariate interactions, a case pointed out in early 2011
  > #  by Frank Harrell, which as it turns out had never been computed
  > #  correctly by any version of the package.  Which shows how often this
  > #  case arises in practice.
  > #
  > aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y))
  > fit1 <- coxph(Surv(time, status) ~ wt.loss + age*strata(sex) + strata(ph.ecog),
  +               data=lung)
  > tdata <- data.frame(wt.loss=c(10,5,0,10, 15,20,25),
  +                     age    =c(50,60,50,60,70,40,21),
  +                     sex    =c(1,1,2,2,1,1,1),
  +                     ph.ecog=c(0,0,1,1,2,2,2))
  > surv1 <- survfit(fit1, newdata=tdata)
  > 
  > fit2 <- coxph(Surv(time, status) ~ wt.loss + age + I(age*0), data=lung, 
  +               init=fit1$coef, iter=0, subset=(sex==1 & ph.ecog==0))
  > fit2$var <- fit1$var
  > 
  > surv2 <- survfit(fit2, newdata=list(wt.loss=c(10,5), age=c(50,60)))
  > s1 <- surv1[1:2]
  > aeq(s1$surv, surv2$surv) #first a vector, second a matrix
  [1] TRUE
  > aeq(s1$std.err, surv2$std.err)
  [1] TRUE
  > aeq(s1[1]$time, surv2$time)
  [1] TRUE
  > aeq(s1[1]$n.event, surv2$n.event)
  [1] TRUE
  > 
  > fit3 <- coxph(Surv(time, status) ~ wt.loss + age + I(age*1),
  +               data=lung, init=fit1$coef, iter=0,
  +               subset=(sex==2 & ph.ecog==1))
  > fit3$var <- fit1$var
  > surv3 <- survfit(fit3, newdata=list(wt.loss=c(0,10), age=c(50,60)))
  > aeq(surv1[3:4]$surv, surv3$surv)
  [1] TRUE
  > aeq(surv1[3:4]$std, surv3$std)
  [1] TRUE
  > 
  > fit4 <- coxph(Surv(time, status) ~ wt.loss + age + I(age*0),
  +               data=lung, init=fit1$coef, iter=0,
  +               subset=(sex==1 & ph.ecog==2))
  > fit4$var <- fit1$var
  > surv4 <- survfit(fit4, newdata=list(wt.loss=c(15,20,25), age=c(70,40,21)))
  > 
  > aeq(surv1[5:7]$surv, surv4$surv)
  [1] TRUE
  > aeq(surv1[5:7]$std.err, surv4$std.err)
  [1] TRUE
  > aeq(surv1[5]$n.risk, surv4$n.risk)
  [1] TRUE
  >              
  > 
  > proc.time()
     user  system elapsed 
     1.29    0.18    1.46 
Running the tests in 'tests/detail.R' failed.
Complete output:
  > # A short test on coxph.detail, to ensure that the computed hazard is
  > #  equal to the theoretical value
  > library(survival)
  > aeq <- function(a,b) all.equal(as.vector(a), as.vector(b))
  > 
  > # taken from book4.R
  > test2 <- data.frame(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
  +                     stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
  +                     event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
  +                     x    =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0) )
  > 
  > byhand <- function(beta, newx=0) {
  +     r <- exp(beta)
  +     loglik <- 4*beta - (log(r+1) + log(r+2) + 2*log(3*r+2) + 2*log(3*r+1) +
  +         log(2*r +2))
  +     u <- 1/(r+1) +  1/(3*r+1) + 2*(1/(3*r+2) + 1/(2*r+2)) -
  +                  ( r/(r+2) +3*r/(3*r+2) + 3*r/(3*r+1))
  +     imat <- r*(1/(r+1)^2 + 2/(r+2)^2 + 6/(3*r+2)^2 +
  +             6/(3*r+1)^2 + 6/(3*r+2)^2 + 4/(2*r +2)^2)
  + 
  +     hazard <-c( 1/(r+1), 1/(r+2), 1/(3*r+2), 1/(3*r+1), 1/(3*r+1),
  +                1/(3*r+2), 1/(2*r +2) )
  + 
  + 
  +     # The matrix of weights, one row per obs, one col per time
  +     #   deaths at 2,3,6,7,8,9
  +     wtmat <- matrix(c(1,0,0,0,1, 0, 0,0,0,0,
  +                       0,1,0,1,1, 0, 0,0,0,0,
  +                       0,0,1,1,1, 0, 1,1,0,0,
  +                       0,0,0,1,1, 0, 1,1,0,0,
  +                       0,0,0,0,1, 1, 1,1,0,0,
  +                       0,0,0,0,0, 1, 1,1,1,1,
  +                       0,0,0,0,0,.5,.5,1,1,1), ncol=7)
  +     wtmat <- diag(c(r,1,1,r,1,r,r,r,1,1)) %*% wtmat
  + 
  +     x      <- c(1,0,0,1,0,1,1,1,0,0)
  +     status <- c(1,1,1,1,1,1,1,0,0,0)
  +     xbar <- colSums(wtmat*x)/ colSums(wtmat)
  +     n <- length(x)
  + 
  +    # Table of sums for score and Schoenfeld resids
  +     hazmat <- wtmat %*% diag(hazard) #each subject's hazard over time
  +     dM <- -hazmat  #Expected part
  +     for (i in 1:5) dM[i,i] <- dM[i,i] +1  #observed
  +     dM[6:7,6:7] <- dM[6:7,6:7] +.5  # observed
  +     mart <- rowSums(dM)
  + 
  +     # Table of sums for score and Schoenfeld resids
  +     #  Looks like the last table of appendix E.2.1 of the book
  +     resid <- dM * outer(x, xbar, '-')
  +     score <- rowSums(resid)
  +     scho <- colSums(resid)
  + 
  +     # We need to add the ties back up (they are symmetric)
  +     scho[6:7] <- rep(mean(scho[6:7]), 2)
  + 
  +     list(loglik=loglik, u=u, imat=imat, xbar=xbar, haz=hazard* exp(beta*newx),
  + 	     mart=mart,  score=score, rmat=resid,
  + 		scho=scho)
  +     } 
  > 
  > # The actual coefficient of the fit is close to zero.  Using a larger
  > #  number pushes the test harder, but it should still work without
  > #  the init and iter arguments, i.e., for any coefficient.
  > fit1 <- coxph(Surv(start, stop, event) ~x, test2,init=-1, iter=0)
  > temp <- coxph.detail(fit1)
  > temp2 <- byhand(fit1$coef, fit1$means)
  > aeq(temp$haz, c(temp2$haz[1:5], sum(temp2$haz[6:7])))
  [1] TRUE
  >                   
  > 
  > proc.time()
     user  system elapsed 
     1.29    0.26    1.54 
Running the tests in 'tests/difftest.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # Test some more features of surv.diff
  > #
  > #  First, what happens when one group is a dummy
  > #
  > 
  > 
  > #
  > # The AML data, with a third group of early censorings "tacked on"
  > #
  > aml3 <- list(time=   c( 9, 13, 13, 18, 23, 28, 31, 34, 45, 48, 161,
  + 			   5,  5,  8,  8, 12, 16, 23, 27, 30, 33, 43, 45,
  + 			   1,  2,  2,  3,  3,  3,  4),
  + 	    status= c( 1,1,0,1,1,0,1,1,0,1,0, 1,1,1,1,1,0,1,1,1,1,1,1,
  + 				0,0,0,0,0,0,0),
  + 	    x     = as.factor(c(rep("Maintained", 11),
  + 				rep("Nonmaintained", 12), rep("Dummy",7) )))
  > 
  > aml3 <- data.frame(aml3)
  > 
  > # These should give the same result (chisq, df), but the second has an
  > #  extra group
  > survdiff(Surv(time, status) ~x, aml)
  Call:
  survdiff(formula = Surv(time, status) ~ x, data = aml)
  
                   N Observed Expected (O-E)^2/E (O-E)^2/V
  x=Maintained    11        7    10.69      1.27       3.4
  x=Nonmaintained 12       11     7.31      1.86       3.4
  
   Chisq= 3.4  on 1 degrees of freedom, p= 0.07 
  > survdiff(Surv(time, status) ~x, aml3)
  Call:
  survdiff(formula = Surv(time, status) ~ x, data = aml3)
  
                   N Observed Expected (O-E)^2/E (O-E)^2/V
  x=Dummy          7        0     0.00       NaN       NaN
  x=Maintained    11        7    10.69      1.27       3.4
  x=Nonmaintained 12       11     7.31      1.86       3.4
  
   Chisq= 3.4  on 1 degrees of freedom, p= 0.07 
  > 
  > 
  > #
  > # Now a test of the stratified log-rank
  > #   There are no tied times within institution, so the coxph program
  > #   can be used to give a complete test
  > #
  > fit <- survdiff(Surv(time, status) ~ pat.karno + strata(inst), lung)
  > 
  > cfit <- coxph(Surv(time, status) ~ factor(pat.karno) + strata(inst),
  + 		lung, iter=0)
  > 
  > tdata <- na.omit(lung[,c('time', 'status', 'pat.karno', 'inst')])
  > 
  > temp1 <- tapply(tdata$status-1, list(tdata$pat.karno, tdata$inst), sum)
  > temp1 <- ifelse(is.na(temp1), 0, temp1)
  > temp2 <- tapply(cfit$resid,  list(tdata$pat.karno, tdata$inst), sum)
  > temp2 <- ifelse(is.na(temp2), 0, temp2)
  > 
  > temp2 <- temp1 - temp2
  > 
  > #Now temp1=observed, temp2=expected
  > all.equal(c(temp1), c(fit$obs))
  [1] TRUE
  > all.equal(c(temp2), c(fit$exp))
  [1] TRUE
  > 
  > all.equal(fit$var[-1,-1], solve(cfit$var))
  [1] TRUE
  > 
  > rm(tdata, temp1, temp2)
  > 
  > proc.time()
     user  system elapsed 
     1.09    0.14    1.21 
Running the tests in 'tests/doaml.R' failed.
Complete output:
  > 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 
     1.15    0.17    1.31 
Running the tests in 'tests/expected.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #  Tests of expected survival
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > #
  > # This makes several scripts easier
  > # 
  > mdy.Date <- function(m, d, y) {
  +     y <- ifelse(y<100, y+1900, y)
  +     as.Date(paste(m,d,y, sep='/'), "%m/%d/%Y")
  +     }
  > 
  > # This function takes a single subject and walks down the rate table
  > # Input: the vector of starting points, futime, and a ratetable
  > # Output: the full history of walking through said table. Let n= #unique
  > #   rates that were used
  > #       cell = n by #dims of the table: index of the table cell
  > #       days = time spent in cell
  > #       hazard= accumulated hazard = days * rate
  > # This does not do date or factor conversions -- start has to be numeric
  > #
  > ratewalk <- function(start, futime, ratetable=survexp.us) {
  +     if (!is.ratetable(ratetable)) stop("Bad rate table")
  +     ratedim <- dim(ratetable)
  +     nvar <- length(ratedim)
  +     if (length(start) != nvar) stop("Wrong length for start")
  +     if (futime <=0) stop("Invalid futime")
  + 
  +     attR <- attributes(ratetable)
  +     discrete <- (attR$type ==1)  #discrete categories
  + 
  +     maxn <- sum(!discrete)*prod(ratedim[!discrete]) #most cells you can hit
  +     cell <- matrix(0, nrow=maxn, ncol=nvar)
  +     days <- hazard <- double(maxn)
  + 
  +     eps <- 1e-8  #Avoid round off error
  +     n <- 0
  +     while (futime >0) {
  +         n <- n+1
  +         #what cell am I in?
  +         # Note that at the edges of the rate table, we use the edge: if
  +         #   it only goes up the the year 2000, year 2000 is used for any
  +         #   dates beyond.  This effectively eliminates one boundary
  +         cell[n,discrete] <- start[discrete]
  +         edge <- futime  #time to nearest edge, or finish
  +         for (j in which(!discrete)) {
  +             indx <- sum(start[j] >= attR$cutpoints[[j]]-eps)
  +             cell[n, j] <- max(1, indx)
  +             if (indx < ratedim[j]) 
  +                 edge <- min(edge, (attR$cutpoints[[j]])[indx+1] - start[j])
  +             }
  +         days[n] <- edge  #this many days in the cell
  +         # using a matrix as a subscript is so handy sometimes
  +         hazard[n] <- edge * (as.matrix(ratetable))[cell[n,,drop=F]]
  +         futime <- futime - edge  #amount of time yet to account for
  +         start[!discrete] <- start[!discrete] + edge  #walk forward in time
  +         }
  +     list(cell=cell[1:n,], days=days[1:n], hazard=hazard[1:n])
  +     }
  > 
  > # Simple test of ratewalk: 20 years old, start on 7Sep 1960
  > #   116 days at the 1960, 20 year old male rate, through the end of the day
  > #     on 12/31/1960, then 84 days at the 1961 rate.  
  > #   The decennial q for 1960 males is .00169.
  > zz <- ratewalk(c(20.4*365.25, 1, as.Date("1960/09/07")), 200)
  > all.equal(zz$hazard[1], -(116/365.25)*log(1-.00169))
  [1] TRUE
  > all.equal(zz$days, c(116,84))
  [1] TRUE
  > 
  >         
  > #
  > # Simple case 1: a single male subject, born 1/1/36 and entered on study 1/2/55
  > #
  > #  Compute the 1, 5, 10 and 12 year expected survival
  > 
  > temp1 <- mdy.Date(1,1,36)
  > temp2 <- mdy.Date(1,2,55)
  > exp1 <- survexp(~1,  ratetable=survexp.usr,times=c(366, 1827, 3653, 4383),
  +                 rmap= list(year=temp2, age=(temp2-temp1), sex=1, race='white'))
  > 
  > t12 <- as.numeric(temp2-temp1)   # difftimes are a PITA
  > h1 <- ratewalk(c(t12, 1, 1, temp2), 366,  survexp.usr)
  > h2 <- ratewalk(c(t12, 1, 1, temp2), 1827, survexp.usr)
  > h3 <- ratewalk(c(t12, 1, 1, temp2), 3653, survexp.usr)
  > h4 <- ratewalk(c(t12, 1, 1, temp2), 4383, survexp.usr)
  > 
  > aeq(-log(exp1$surv), c(sum(h1$hazard), sum(h2$hazard), sum(h3$hazard),
  +                        sum(h4$hazard)))
  [1] TRUE
  > 
  > # pyears should give the same result
  > dummy <- data.frame(time = 4383,
  +                     year=temp2, sex = 1, age= temp2-temp1, race="white")
  > cuts <- tcut(0, c(0, 366, 1827, 3653, 4383))
  > exp1c <- pyears(time ~ cuts, data=dummy, ratetable=survexp.usr)
  > aeq(exp1$surv, exp(-cumsum(exp1c$expected)))
  [1] TRUE
  > 
  > 
  > # Just a little harder: 
  > #   Born 3/1/25 and entered the study on 6/10/55.  The code creates shifted
  > #   dates to align with US rate tables - entry is 59 days earlier (days from
  > #   1/1/1925 to 3/1/1925).
  > #
  > temp1 <- mdy.Date(3,1,25)
  > temp2 <- mdy.Date(6,10,55)
  > exp1 <- survexp(~1, ratetable=survexp.usr,times=c(366, 1827, 3653, 4383),
  +                 rmap= list(year=temp2, age=(temp2-temp1), sex=2, race='black'))
  > 
  > tyear <- temp2 - 59
  > t12 <-  as.numeric(temp2-temp1)
  > h1 <- ratewalk(c(t12, 2, 2, tyear), 366,  survexp.usr)
  > h2 <- ratewalk(c(t12, 2, 2, tyear), 1827, survexp.usr)
  > h3 <- ratewalk(c(t12, 2, 2, tyear), 3653, survexp.usr)
  > h4 <- ratewalk(c(t12, 2, 2, tyear), 4383, survexp.usr)
  > 
  > aeq(-log(exp1$surv), c(sum(h1$hazard), sum(h2$hazard), sum(h3$hazard),
  +                        sum(h4$hazard)))
  [1] TRUE
  > 
  > #
  > # Simple case 2: make sure that the averages are correct, for Ederer method
  > #
  > #  Compute the 1, 5, 10 and 12 year expected survival
  > 
  > temp1 <- mdy.Date(1:6,6:11,1890:1895)
  > temp2 <- mdy.Date(6:1,11:6,c(55:50))
  > temp3 <- c(1,2,1,2,1,2)
  > age <- temp2 - temp1
  > 
  > exp1 <- survexp(~1, rmap= list(year=temp2, age=(temp2-temp1), sex=temp3),
  + 		       times=c(366, 1827, 3653, 4383))
  > exp2 <- survexp(~ I(1:6), 
  +                 rmap= list(year=temp2, age=(temp2-temp1), sex=temp3),
  + 			times=c(366, 1827, 3653, 4383))
  > exp3 <- exp2$surv
  > for (i in 1:length(temp1)){
  +     exp3[,i] <- survexp(~ 1,
  +                         rmap = list(year=temp2, age=(temp2-temp1), sex=temp3),
  +                         times=c(366, 1827, 3653, 4383), subset=i)$surv
  +     }
  > 
  >              
  > print(aeq(exp2$surv, exp3))
  [1] TRUE
  > print(all.equal(exp1$surv, apply(exp2$surv, 1, mean)))
  [1] TRUE
  > 
  > # They agree, but are they right?
  > #
  > for (i in 1:length(temp1)) {
  +     offset <- as.numeric(temp1[i] - mdy.Date(1,1, 1889+i))
  +     tyear = temp2[i] - offset
  +     haz1 <- ratewalk(c(as.numeric(temp2-temp1)[i], temp3[i], tyear), 366)
  +     haz2 <- ratewalk(c(as.numeric(temp2-temp1)[i], temp3[i], tyear), 1827)
  +     haz3 <- ratewalk(c(as.numeric(temp2-temp1)[i], temp3[i], tyear), 3653)
  +     haz4 <- ratewalk(c(as.numeric(temp2-temp1)[i], temp3[i], tyear), 4383)
  +     print(aeq(-log(exp2$surv[,i]), c(sum(haz1$hazard), sum(haz2$hazard),
  +                                     sum(haz3$hazard), sum(haz4$hazard))))
  +     }
  [1] TRUE
  [1] TRUE
  [1] TRUE
  [1] TRUE
  [1] TRUE
  [1] TRUE
  > 
  > #
  > # Check that adding more time points doesn't change things
  > #
  > exp4 <- survexp(~ I(1:6),
  +           rmap= list(year=temp2, age=(temp2-temp1), sex=temp3),
  + 		times=sort(c(366, 1827, 3653, 4383, 30*(1:100))))
  > aeq(exp4$surv[match(exp2$time, exp4$time),], exp2$surv)
  [1] TRUE
  > 
  > exp4 <- survexp(~1,
  +                 rmap = list(year=temp2, age=(temp2-temp1), sex=temp3),
  + 		times=sort(c(366, 1827, 3653, 4383, 30*(1:100))))
  > aeq(exp1$surv, exp4$surv[match(exp1$time, exp4$time, nomatch=0)])
  [1] TRUE
  > 
  > 
  > #
  > # Now test Hakulinen's method, assuming an analysis date of 3/1/57
  > #
  > futime <- mdy.Date(3,1,57) - temp2
  > xtime  <- sort(c(futime, 30, 60, 185, 365))
  > 
  > exp1 <- survexp(futime ~ 1, rmap= list(year=temp2, age=(temp2-temp1), sex=1),
  + 		times=xtime, conditional=F)
  > exp2 <- survexp(~ I(1:6), times=futime, 
  +                 rmap= list(year=temp2, age=(temp2-temp1), sex=1))
  > 
  > wt <- rep(1,6)
  > con <- double(6)
  > for (i in 1:6) {
  +     con[i] <- sum(exp2$surv[i,i:6])/sum(wt[i:6])
  +     wt <- exp2$surv[i,]
  +     }
  > 
  > exp1$surv[match(futime, xtime)]
  [1] 0.9557362 0.9285840 0.9025661 0.8774220 0.8532489 0.8297416
  > aeq(exp1$surv[match(futime, xtime)], cumprod(con))
  [1] TRUE
  > 
  > 
  > #
  > # Now for the conditional method
  > #
  > exp1 <- survexp(futime ~ 1, rmap= list(year=temp2, age=(temp2-temp1), sex=1),
  + 		times=xtime, conditional=T)
  > 
  > cond <- exp2$surv
  > for (i in 6:2) cond[i,] <- (cond[i,]/cond[i-1,])  #conditional survival
  > for (i in 1:6) con[i] <- exp(mean(log(cond[i, i:6])))
  > 
  > all.equal(exp1$surv[match(futime, xtime)], cumprod(con))
  [1] TRUE
  > cumprod(con)
  [1] 0.9556656 0.9284398 0.9023612 0.8771798 0.8529944 0.8294940
  > 
  > #
  > # Test out expected survival, when the parent pop is another Cox model
  > #
  > test1 <- data.frame(time=  c(4, 3,1,1,2,2,3),
  +                     status=c(1,NA,1,0,1,1,0),
  +                     x=     c(0, 2,1,1,1,0,0))
  > 
  > fit <- coxph(Surv(time, status) ~x, test1, method='breslow')
  > 
  > dummy <- data.frame(time=c(.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5),
  + 		    status=c(1,0,1,0,1,0,1,1,1), x=(-4:4)/2)
  > 
  > efit <- survexp(time ~ 1, rmap= list(x=x), dummy, ratetable=fit, cohort=F)
  > 
  > #
  > # Now, compare to the true answer, which is known to us
  > #
  > ss <- exp(fit$coef)
  > haz <- c( 1/(3*ss+3), 2/(ss+3), 1) #truth at time 0,1,2,4+
  > chaz <- cumsum(c(0,haz))
  > chaz2 <- chaz[c(1,2,2,3,3,3,3,4,4)]
  > 
  > risk <- exp(fit$coef*dummy$x)
  > efit2 <- exp(-risk*chaz2)
  > 
  > all.equal(as.vector(efit), as.vector(efit2))  #ignore mismatched name attrib
  [1] TRUE
  > 
  > #
  > # Now test the direct-adjusted curve (Ederer)
  > #
  > efit <- survexp( ~ 1, dummy, ratetable=fit, se=F)
  > direct <- survfit(fit, newdata=dummy, censor=FALSE)$surv
  > 
  > chaz <- chaz[-1]                  #drop time 0
  > d2 <- exp(outer(-chaz, risk))
  > all.equal(as.vector(direct), as.vector(d2))   #this tests survfit
  [1] TRUE
  > 
  > all.equal(as.vector(efit$surv), as.vector(apply(direct,1,mean)))  #direct
  [1] TRUE
  > 
  > # Check out the "times" arg of survexp
  > efit2 <- survexp( ~1, dummy, ratetable=fit, se=F,
  +                   times=c(.5, 2, 3.5,6))
  > aeq(efit2$surv, c(1, efit$surv[c(2,2,3)]))
  [1] TRUE
  > 
  > #
  > # Now test out the Hakulinen method (Bonsel's method)
  > #  By construction, we have a large correlation between x and censoring
  > #
  > # In theory, hak1 and hak2 would be the same.  In practice, like a KM and
  > #   F-H, they differ when n is small.
  > #
  > efit <- survexp( time ~1, dummy, ratetable=fit, se=F)
  > 
  > surv  <- wt <- rep(1,9)
  > tt <- c(1,2,4)
  > hak1 <- hak2 <- NULL
  > for (i in 1:3) {
  +     wt[dummy$time < tt[i]]  <- 0
  +     hak1 <- c(hak1,  exp(-sum(haz[i]*risk*surv*wt)/sum(surv*wt)))
  +     hak2 <- c(hak2,  sum(exp(-haz[i]*risk)*surv*wt)/sum(surv*wt))
  +     surv <- surv * exp(-haz[i]*risk)
  +     }
  > 
  > all.equal(as.vector(efit$surv), as.vector(cumprod(hak1)))
  [1] TRUE
  > 
  > #
  > #  Now do the conditional estimate
  > #
  > efit <- survexp( time ~ 1, dummy, ratetable=fit, se=F,
  + 			conditional=T)
  > wt <- rep(1,9)
  > cond <- NULL
  > for (i in 1:3) {
  +     wt[dummy$time < tt[i]]  <- 0
  +     cond <- c(cond,  exp(-sum(haz[i]*risk*wt)/sum(wt)))
  +     }
  > 
  > all.equal(as.vector(efit$surv), as.vector(cumprod(cond)))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.87    0.15    2.03 
Running the tests in 'tests/factor2.R' failed.
Complete output:
  > library(survival)
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > options(na.action=na.exclude)
  > #
  > # More tests of factors in prediction, using a new data set
  > #
  > fit <- coxph(Surv(time, status) ~  factor(ph.ecog), lung)
  > 
  > tdata <- data.frame(ph.ecog = factor(0:3))
  > p1 <- predict(fit, newdata=tdata, type='lp')
  > p2 <- predict(fit, type='lp')
  > aeq(p1, p2[match(0:3, lung$ph.ecog)])
  [1] TRUE
  > 
  > fit2 <- coxph(Surv(time, status) ~ factor(ph.ecog) + factor(sex), lung)
  > tdata <- expand.grid(ph.ecog = factor(0:3), sex=factor(1:2))
  > p1 <- predict(fit2, newdata=tdata, type='risk')
  > 
  > xdata <- expand.grid(ph.ecog=factor(1:3), sex=factor(1:2))
  > p2 <- predict(fit2, newdata=xdata, type='risk')
  > all.equal(p2, p1[c(2:4, 6:8)], check.attributes=FALSE)
  [1] TRUE
  > 
  > 
  > fit3 <- survreg(Surv(time, status) ~ factor(ph.ecog) + age, lung)
  > tdata <- data.frame(ph.ecog=factor(0:3), age=50)
  > predict(fit, type='lp', newdata=tdata)
          1         2         3         4 
  0.0000000 0.3688401 0.9163870 2.2079803 
  > predict(fit3, type='lp', newdata=tdata)
         1        2        3        4 
  6.399571 6.142938 5.770523 4.916993 
  > 
  > proc.time()
     user  system elapsed 
     1.10    0.15    1.25 
Running the tests in 'tests/fr_cancer.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # Here is a test case with multiple smoothing terms
  > #
  > 
  > fit0 <- coxph(Surv(time, status) ~ ph.ecog + age, lung)
  > fit1 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,3), lung)
  > fit2 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,4), lung)
  > fit3 <- coxph(Surv(time, status) ~ ph.ecog + pspline(age,8), lung)
  > 
  > 
  > 
  > fit4 <- coxph(Surv(time, status) ~ ph.ecog + pspline(wt.loss,3), lung)
  > 
  > fit5 <-coxph(Surv(time, status) ~ ph.ecog + pspline(age,3) + 
  + 	     pspline(wt.loss,3), lung)
  > 
  > fit1
  Call:
  coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 3), 
      data = lung)
  
                              coef se(coef)      se2    Chisq   DF       p
  ph.ecog                  0.44802  0.11707  0.11678 14.64453 1.00 0.00013
  pspline(age, 3), linear  0.01126  0.00928  0.00928  1.47231 1.00 0.22498
  pspline(age, 3), nonlin                             2.07924 2.08 0.37143
  
  Iterations: 4 outer, 12 Newton-Raphson
       Theta= 0.861 
  Degrees of freedom for terms= 1.0 3.1 
  Likelihood ratio test=21.9  on 4.08 df, p=2e-04
  n= 227, number of events= 164 
     (1 observation deleted due to missingness)
  > fit2
  Call:
  coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 4), 
      data = lung)
  
                              coef se(coef)      se2    Chisq   DF       p
  ph.ecog                  0.45047  0.11766  0.11723 14.65751 1.00 0.00013
  pspline(age, 4), linear  0.01117  0.00927  0.00927  1.45195 1.00 0.22822
  pspline(age, 4), nonlin                             2.95816 3.08 0.41197
  
  Iterations: 4 outer, 11 Newton-Raphson
       Theta= 0.797 
  Degrees of freedom for terms= 1.0 4.1 
  Likelihood ratio test=22.7  on 5.07 df, p=4e-04
  n= 227, number of events= 164 
     (1 observation deleted due to missingness)
  > fit3
  Call:
  coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 8), 
      data = lung)
  
                              coef se(coef)      se2    Chisq   DF       p
  ph.ecog                  0.47640  0.12024  0.11925 15.69732 1.00 7.4e-05
  pspline(age, 8), linear  0.01172  0.00923  0.00923  1.61161 1.00    0.20
  pspline(age, 8), nonlin                             6.93188 6.99    0.43
  
  Iterations: 5 outer, 15 Newton-Raphson
       Theta= 0.691 
  Degrees of freedom for terms= 1 8 
  Likelihood ratio test=27.6  on 8.97 df, p=0.001
  n= 227, number of events= 164 
     (1 observation deleted due to missingness)
  > fit4
  Call:
  coxph(formula = Surv(time, status) ~ ph.ecog + pspline(wt.loss, 
      3), data = lung)
  
                                coef se(coef)      se2    Chisq   DF     p
  ph.ecog                    0.51545  0.12960  0.12737 15.81939 1.00 7e-05
  pspline(wt.loss, 3), line -0.00702  0.00655  0.00655  1.14638 1.00  0.28
  pspline(wt.loss, 3), nonl                             2.44612 2.09  0.31
  
  Iterations: 3 outer, 10 Newton-Raphson
       Theta= 0.776 
  Degrees of freedom for terms= 1.0 3.1 
  Likelihood ratio test=21.1  on 4.06 df, p=3e-04
  n= 213, number of events= 151 
     (15 observations deleted due to missingness)
  > fit5
  Call:
  coxph(formula = Surv(time, status) ~ ph.ecog + pspline(age, 3) + 
      pspline(wt.loss, 3), data = lung)
  
                                coef se(coef)      se2    Chisq   DF       p
  ph.ecog                    0.47422  0.13495  0.13206 12.34842 1.00 0.00044
  pspline(age, 3), linear    0.01368  0.00976  0.00974  1.96406 1.00 0.16108
  pspline(age, 3), nonlin                               1.90116 2.07 0.40284
  pspline(wt.loss, 3), line -0.00717  0.00661  0.00660  1.17529 1.00 0.27832
  pspline(wt.loss, 3), nonl                             2.07729 2.03 0.35929
  
  Iterations: 4 outer, 12 Newton-Raphson
       Theta= 0.85 
       Theta= 0.779 
  Degrees of freedom for terms= 1.0 3.1 3.0 
  Likelihood ratio test=25.2  on 7.06 df, p=7e-04
  n= 213, number of events= 151 
     (15 observations deleted due to missingness)
  > 
  > rm(fit1, fit2, fit3, fit4, fit5)
  > 
  > proc.time()
     user  system elapsed 
     1.17    0.17    1.32 
Running the tests in 'tests/fr_lung.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # A test with the lung data
  > #  This caused problems in one release
  > 
  > #
  > #   First, get rid of some missings
  > #
  > lung2 <- na.omit(lung[c('time', 'status', 'wt.loss')])
  > 
  > #
  > # Test the logliklihoods
  > #
  > fit <- coxph(Surv(time, status) ~ pspline(wt.loss,3), lung2, x=T)
  > fit0<- coxph(Surv(time, status) ~ 1, lung2)
  > fit1<- coxph(Surv(time, status) ~ fit$x, lung2, iter=0, init=fit$coef)
  > 
  > all.equal(fit$loglik[1], fit0$loglik)
  [1] TRUE
  > all.equal(fit$loglik[2], fit1$loglik[2])
  [1] TRUE
  > 
  > #
  > # Check variances
  > #
  > imat <- solve(fit1$var)
  > var2 <- fit$var %*% imat %*% fit$var
  > all.equal(fit$var2, var2)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.15    0.17    1.31 
Running the tests in 'tests/fr_ovarian.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # Test on the ovarian data
  > 
  > fit1 <- coxph(Surv(futime, fustat) ~ rx + age, ovarian)
  > fit2 <- coxph(Surv(futime, fustat) ~ rx + pspline(age, df=2), 
  + 		data=ovarian)
  > fit2$iter
  [1] 2 8
  > 
  > fit2$df
  [1] 0.9426611 1.9293051
  > 
  > fit2$history
  $`pspline(age, df = 2)`
  $`pspline(age, df = 2)`$theta
  [1] 0.4468868
  
  $`pspline(age, df = 2)`$done
  [1] TRUE
  
  $`pspline(age, df = 2)`$history
          thetas      dfs
  [1,] 1.0000000 1.000000
  [2,] 0.0000000 5.000000
  [3,] 0.6000000 1.734267
  [4,] 0.4845205 1.929305
  
  $`pspline(age, df = 2)`$half
  [1] 0
  
  
  > 
  > fit4 <- coxph(Surv(futime, fustat) ~ rx + pspline(age, df=4), 
  + 		data=ovarian)
  > fit4
  Call:
  coxph(formula = Surv(futime, fustat) ~ rx + pspline(age, df = 4), 
      data = ovarian)
  
                              coef se(coef)    se2  Chisq   DF      p
  rx                        -0.373    0.761  0.749  0.241 1.00 0.6238
  pspline(age, df = 4), lin  0.139    0.044  0.044  9.978 1.00 0.0016
  pspline(age, df = 4), non                         2.592 2.93 0.4457
  
  Iterations: 3 outer, 14 Newton-Raphson
       Theta= 0.242 
  Degrees of freedom for terms= 1.0 3.9 
  Likelihood ratio test=19.4  on 4.9 df, p=0.001
  n= 26, number of events= 12 
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.37    0.17    1.51 
Running the tests in 'tests/fr_simple.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # Test the logic of the penalized code by fitting some no-frailty models
  > #  (theta=0).  It should give exactly the same answers as 'ordinary' coxph.
  > #
  > test1 <- data.frame(time=  c(4, 3,1,1,2,2,3),
  +                     status=c(1,NA,1,0,1,1,0),
  +                     x=     c(0, 2,1,1,1,0,0))
  > 
  > test2 <- data.frame(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
  +                     stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
  +                     event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
  +                     x    =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0) )
  > 
  > zz <- zz2 <- rep(0, nrow(test1))
  > tfit1 <- coxph(Surv(time,status) ~x, test1, eps=1e-7)
  > tfit2 <- coxph(Surv(time,status) ~x + frailty(zz, theta=0, sparse=T), test1)
  > tfit3 <- coxph(Surv(zz,time,status) ~x + frailty(zz2, theta=0, sparse=T), test1)
  > 
  > temp <- c('coefficients', 'var', 'loglik', 'linear.predictors',
  + 	  'means', 'n', 'concordance')
  > 
  > all.equal(tfit1[temp], tfit2[temp])
  [1] TRUE
  > all.equal(tfit2[temp], tfit3[temp])
  [1] TRUE
  > 
  > zz <- rep(0, nrow(test2))
  > tfit1 <- coxph(Surv(start, stop, event) ~x, test2, eps=1e-7)
  > tfit2 <- coxph(Surv(start, stop, event) ~ x + frailty(zz, theta=0, sparse=T),
  + 	       test2)
  > all.equal(tfit1[temp], tfit2[temp])
  [1] TRUE
  > 
  > 
  > #
  > # Repeat the above tests, but with a strata added
  > #  Because the data set is simply doubled, the loglik will double,
  > #   beta is the same, variance is halved.
  > #
  > test3 <- rbind(test1, test1)
  > test3$x2 <- rep(1:2, rep(nrow(test1),2))
  > zz <- zz2 <- rep(0, nrow(test3))
  > tfit1 <- coxph(Surv(time,status) ~x + strata(x2), test3, eps=1e-7)
  > tfit2 <- coxph(Surv(time,status) ~x + frailty(zz, theta=0, sparse=T)
  + 	       + strata(x2), test3)
  > tfit3 <- coxph(Surv(zz,time,status) ~x + frailty(zz2, theta=0, sparse=T)
  + 	         + strata(x2), test3)
  > 
  > all.equal(tfit1[temp], tfit2[temp])
  [1] TRUE
  > all.equal(tfit2[temp], tfit3[temp])
  [1] TRUE
  > 
  > 
  > test4 <- rbind(test2, test2)
  > test4$x2 <- rep(1:2, rep(nrow(test2),2))
  > zz <- rep(0, nrow(test4))
  > tfit1 <- coxph(Surv(start, stop, event) ~x, test4, eps=1e-7)
  > tfit2 <- coxph(Surv(start, stop, event) ~ x + frailty(zz, theta=0, sparse=T),
  + 	       test4)
  > all.equal(tfit1[temp], tfit2[temp])
  [1] TRUE
  > 
  > rm(test3, test4, tfit1, tfit2, tfit3, temp, zz, zz2)
  > 
  > proc.time()
     user  system elapsed 
     1.34    0.12    1.46 
Running the tests in 'tests/mstate.R' failed.
Complete output:
  > #
  > # A tiny multi-state example
  > #
  > library(survival)
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > mtest <- data.frame(id= c(1, 1, 1,  2,  3,  4, 4, 4,  5, 5),
  +                     t1= c(0, 4, 9,  0,  2,  0, 2, 8,  1, 3),
  +                     t2= c(4, 9, 10, 5,  9,  2, 8, 9,  3, 11),
  +                     st= c(1, 2,  1, 2,  3,  1, 3, 0,  2,  0))
  > 
  > mtest$state <- factor(mtest$st, 0:3, c("censor", "a", "b", "c"))
  > 
  > if (FALSE) {
  +     # this graph is very useful when debugging
  +     temp <- survcheck(Surv(t1, t2, state) ~1, mtest, id=id)
  +     plot(c(0,11), c(1,5.1), type='n', xlab="Time", ylab= "Subject")
  +     with(mtest, segments(t1+.1, id, t2, id, col=as.numeric(temp$istate)))
  +     event <- subset(mtest, state!='censor')
  +     text(event$t2, event$id+.2, as.character(event$state))
  + }
  > 
  > mtest <- mtest[c(1,3,2,4,5,7,6,10, 9, 8),]  #not in time order
  > 
  > mfit <- survfit(Surv(t1, t2, state) ~ 1, mtest, id=id)
  > 
  > # True results
  > #
  > #time       state                    probabilities
  > #         entry  a   b  c         entry  a    b     c
  > #
  > #0        124                      1     0    0     0
  > #1+       1245
  > #2+       1235   4                3/4   1/4   0     0    4 -> a, add 3
  > #3+       123    4   5            9/16  1/4  3/16   0    5 -> b
  > #4+        23    14  5            6/16  7/16 3/16   0    1 -> a
  > #5+        3     14  5            3/16  7/16 6/16   0    2 -> b, exits
  > #8+        3     1   5  4         3/16  7/32 6/16  7/32  4 -> c
  > #9+                  15            0     0  19/32 13/32  1->b, 3->c & exit
  > # 10+            1   5                19/64 19/64 13/32  1->a
  > 
  > # In mfit, the "entry" state is first in the matrices, when this function was
  > #  first created it was the last.
  > swap <- c(4,1,2,3)  # at one time it was last
  > all.equal(mfit$n.risk,  matrix(c(4,4,3,2,1,1,0,0,
  +                                  0,1,1,2,2,1,0,0,
  +                                  0,0,1,1,1,1,2,1,
  +                                  0,0,0,0,0,1,0,0), ncol=4))
  [1] TRUE
  > all.equal(mfit$pstate,  matrix(c(24, 18, 12,  6, 6, 0, 0,  0,
  +                                   8,  8, 14, 14, 7, 0,  9.5, 9.5, 
  +                                   0,  6,  6, 12, 12,19,9.5, 9.5, 
  +                                   0,  0,  0,  0, 7, 13, 13, 13)/32, ncol=4))
  [1] TRUE
  > all.equal(mfit$n.event, matrix(c(0,0,0,0,0,0,0,0,
  +                                  1,0,1,0,0,0,1,0,
  +                                  0,1,0,1,0,1,0,0,
  +                                  0,0,0,0,1,1,0,0), ncol=4))
  [1] TRUE
  > all.equal(mfit$time, c(2, 3, 4, 5, 8, 9, 10, 11))
  [1] TRUE
  > 
  > 
  > # Somewhat more complex.
  > #  Scramble the input data
  > #  Not everyone starts at the same time or in the same state
  > #  Case weights
  > #
  > tdata <- data.frame(id= c(1, 1, 1,  2,  3,  4, 4, 4,  5,  5),
  +                     t1= c(0, 4, 9,  1,  2,  0, 2, 8,  1,  3),
  +                     t2= c(4, 9, 10, 5,  9,  2, 8, 9,  3, 11),
  +                     st= c(1, 2,  1, 2,  3,  1, 3, 0,  3,  0),
  +                     i0= c(4, 1,  2, 1,  4,  4, 1, 3,  2,  3),
  +                     wt= 1:10)
  > 
  > tdata$st <- factor(tdata$st, c(0:3),
  +                     labels=c("censor", "1", "2", "3"))
  > tdata$i0 <- factor(tdata$i0, 1:4,
  +                     labels=c("1", "2", "3", "entry"))
  > 
  > tfun <- function(data=tdata) {
  +     reorder <- c(10, 9, 1, 2, 5, 4, 3, 7, 8, 6)
  +     new <- data[reorder,]
  +     new
  + }
  > 
  > # These weight vectors are in the order of tdata
  > # w[9] is the weight for subject 5 at time 1.5, for instance
  > # p0 is defined as all those at risk just before the first event, which in
  > #  this data set is entry:a at time 2 for id=4; id 1,2,4,5 at risk
  > p0 <- function(w) c(w[4], w[9], 0, w[1]+ w[6])/ (w[1]+ w[4] + w[6] + w[9])
  > 
  > #  aj2 = Aalen-Johansen H matrix at time 2, etc.
  > aj2 <- function(w) {
  +     rbind(c(1, 0, 0, 0),    # state a (1) stays put
  +           c(0, 1, 0, 0),
  +           c(0, 0, 1, 0),
  +           c(w[6], 0, 0, w[1])/(w[1] + w[6]))  #subject 4 moves to 'a'
  + }
  > aj3 <- function(w) rbind(c(1, 0, 0, 0),   
  +                          c(0, 0, 1, 0),  # 5 moves from b to c
  +                          c(0, 0, 1, 0),
  +                          c(0, 0, 0, 1))
  > aj4 <- function(w) rbind(c(1, 0, 0, 0),
  +                          c(0, 1, 0, 0),  
  +                          c(0, 0, 1, 0),
  +                          c(w[1], 0, 0, w[5])/(w[1] + w[5])) #1 moves from 4 to a
  > aj5 <- function(w) rbind(c(w[2]+w[7], w[4], 0, 0)/(w[2]+ w[4] + w[7]), #2 to b
  +                          c(0, 1, 0, 0),  
  +                          c(0, 0, 1, 0),
  +                          c(0, 0, 0, 1))
  > aj8 <- function(w) rbind(c(w[2], 0, w[7], 0)/(w[2]+ w[7]), # 4  to c
  +                          c(0, 1, 0, 0),  
  +                          c(0, 0, 1, 0),
  +                          c(0, 0, 0, 1))
  > aj9 <- function(w) rbind(c(0, 1, 0, 0), # 1  to b
  +                          c(0, 1, 0, 0),  
  +                          c(0, 0, 1, 0),
  +                          c(0, 0, 1 ,0)) # 3 to c
  > aj10 <- function(w)rbind(c(1, 0, 0, 0),
  +                          c(1, 0, 0, 0),  #1 back to a
  +                          c(0, 0, 1, 0),
  +                          c(0, 0, 0, 1))
  > 
  > #time       state               
  > #         a   b  c  entry
  > #
  > #1        2   5     14       initial distribution
  > #2        24  5     1        4 -> a, add 3
  > #3        24     5  13       5 from b to c
  > #4       124     5   3       1 -> a
  > #5        14     5   3       2 -> b, exits
  > #8        1      45  3       4 -> c
  > #9            1  45          1->b, 3->c & exit
  > #10       1      45          1->a
  > 
  > # P is a product of matrices
  > dopstate <- function(w) {
  +     p1 <- p0(w)
  +     p2 <- p1 %*% aj2(w)
  +     p3 <- p2 %*% aj3(w)
  +     p4 <- p3 %*% aj4(w)
  +     p5 <- p4 %*% aj5(w)
  +     p8 <- p5 %*% aj8(w)
  +     p9 <- p8 %*% aj9(w)
  +     p10<- p9 %*% aj10(w)
  +     rbind(p2, p3, p4, p5, p8, p9, p10, p10)
  + }
  > 
  > # Check the pstate estimate
  > w1 <- rep(1,10)
  > mtest2 <- tfun(tdata)  # scrambled order
  > mfit2 <- survfit(Surv(t1, t2, st) ~ 1, tdata, id=id, istate=i0) # ordered
  > aeq(mfit2$pstate, dopstate(w1)[,swap])
  [1] TRUE
  > aeq(mfit2$p0, p0(w1)[swap])
  [1] TRUE
  > 
  > mfit2b <- survfit(Surv(t1, t2, st) ~ 1, mtest2, id=id, istate=i0)#scrambled
  > aeq(mfit2b$pstate, dopstate(w1)[,swap])
  [1] TRUE
  > aeq(mfit2b$p0, p0(w1)[swap])
  [1] TRUE
  > 
  > mfit2b$call <- mfit2$call <- NULL
  > all.equal(mfit2b, mfit2) 
  [1] TRUE
  > aeq(mfit2$transitions, c(2,0,1,0, 0,2,0,0, 1,1,1,0, 0,0,0,2))
  [1] TRUE
  > 
  > # Now the harder one, where subjects change weights
  > mfit3  <- survfit(Surv(t1, t2, st) ~ 1, tdata, id=id, istate=i0,
  +                   weights=wt, influence=TRUE)
  > aeq(mfit3$p0, p0(1:10)[swap])
  [1] TRUE
  > aeq(mfit3$pstate, dopstate(1:10)[,swap])
  [1] TRUE
  >     
  > 
  > # The derivative of a matrix product AB is (dA)B + A(dB) where dA is the
  > #  elementwise derivative of A and etc for B.
  > # dp0 creates the derivatives of p0 with respect to each subject, a 5 by 4
  > #  matrix
  > # All the functions below are hand coded for a weight vector that is in
  > #  exactly the same order as the rows of mtest.
  > # Since p0 = (w[4], w[9], 0, w[1]+ w[6])/ (w[1]+ w[4] + w[6] + w[9])
  > #      and subject id is 1,1,1, 2, 3, 4,4,4, 5,5
  > #   we get the derivative below
  > # 
  > dp0 <- function(w) {
  +   p <- p0(w)
  +   w0 <- w[c(1,4,6,9)]  # the 4 obs at the start, subjects 1, 2, 4, 5
  +   rbind(c(0, 0, 0, 1) - p,   # subject 1 affects p[4]
  +         c(1, 0, 0, 0) - p,   # subject 2 affects p[1]
  +         0,                   # subject 3 affects none
  +         c(0, 0, 0, 1) - p,   # subject 4 affect p[4]
  +         c(0, 1, 0, 0) - p)/   # subject 5 affects p[2]
  +       sum(w0)
  + }
  >   
  > 
  > dp2 <- function(w) {
  +     h2 <- aj2(w)   # H matrix at time 2
  +     part1 <- dp0(w) %*% h2
  + 
  +     # 1 and 4 in state 4, obs 4 moves from entry to a
  +     mult  <- p0(w)[4]/(w[1] + w[6])  #p(t-) / weights in state
  +     part2 <- rbind((c(0,0,0,1)- h2[4,]) * mult,
  +                    0,
  +                    0,
  +                    (c(1,0,0,0) - h2[4,]) * mult,
  +                    0)
  +     part1 + part2
  + }
  > 
  > dp3 <- function(w) {
  +     dp2(w) %*% aj3(w)
  + }
  > 
  > dp4 <- function(w) {
  +     h4 <- aj4(w)   # H matrix at time 4
  +     part1 <- dp3(w) %*% h4
  + 
  +     # subjects 1 and 3 in state 4, obs 1 and 5, 1 moves to a
  +     mult <- dopstate(w)[2,4]/ (w[1] + w[5])   # p_4(time 4-0) / wt
  +     part2 <- rbind((c(1,0,0,0)- h4[4,]) * mult,
  +                    0,
  +                    (c(0,0,0,1)- h4[4,]) * mult,
  +                    0,
  +                    0)
  +     part1 + part2
  + }
  > dp5 <- function(w) {
  +     h5 <- aj5(w)   # H matrix at time 5
  +     part1 <- dp4(w) %*% h5
  + 
  +     # subjects 124 in state 1, obs 2,4,7, 2 goes to 2
  +     mult <- dopstate(w)[3,1]/ (denom <- w[2] + w[4] + w[7]) 
  +     part2 <- rbind((c(1,0,0,0)- h5[1,]) * mult,
  +                    (c(0,1,0,0)- h5[1,]) * mult,
  +                    0,
  +                    (c(1,0,0,0)- h5[1,]) * mult,
  +                    0)
  +     part1 + part2
  + }
  > dp8 <- function(w) {
  +     h8 <- aj8(w)   # H matrix at time 8
  +     part1 <- dp5(w) %*% h8
  + 
  +     # subjects 14 in state 1, obs 2 &7, 4 goes to c
  +     mult <- dopstate(w)[4, 1]/ (w[2] + w[7]) 
  +     part2 <- rbind((c(1,0,0,0)- h8[1,]) * mult,
  +                    0,
  +                    0,
  +                    (c(0,0,1,0)- h8[1,]) * mult,
  +                    0)
  +     part1 + part2
  + }
  > dp9 <- function(w) dp8(w) %*% aj9(w)
  > dp10<- function(w) dp9(w) %*% aj10(w)
  > 
  > w1 <- 1:10
  > aeq(mfit3$influence[,1,], dp0(w1)[,swap])
  [1] TRUE
  > aeq(mfit3$influence[,2,], dp2(w1)[,swap])
  [1] TRUE
  > aeq(mfit3$influence[,3,], dp3(w1)[,swap])
  [1] TRUE
  > aeq(mfit3$influence[,4,], dp4(w1)[,swap])
  [1] TRUE
  > aeq(mfit3$influence[,5,], dp5(w1)[,swap])
  [1] TRUE
  > aeq(mfit3$influence[,6,], dp8(w1)[,swap])
  [1] TRUE
  > aeq(mfit3$influence[,7,], dp9(w1)[,swap])
  [1] TRUE
  > aeq(mfit3$influence[,8,], dp10(w1)[,swap])
  [1] TRUE
  > aeq(mfit3$influence[,9,], dp10(w1)[,swap]) # no changes at time 11
  [1] TRUE
  > 
  > 
  > # The cumulative hazard at each time point is remapped from a matrix
  > #  into a vector (in survfit)
  > # First check out the names
  > nstate <- length(mfit3$states)
  > temp <- matrix(0, nstate, nstate)
  > indx1 <- match(rownames(mfit3$transitions), mfit3$states)
  > indx2 <- match(colnames(mfit3$transitions), mfit3$states, nomatch=0)
  > temp[indx1, indx2] <- mfit3$transitions[, indx2>0]
  > # temp is an nstate by nstate version of the transitions matrix
  > from <- row(temp)[temp>0]
  > to   <- col(temp)[temp>0]
  >  
  > all.equal(colnames(mfit3$cumhaz), paste(from, to, sep='.'))
  [1] TRUE
  > 
  > 
  > hazard <- function(fit, i, indx=which(temp>0)) {
  +     nstate <- length(fit$states)
  +     cmat <- matrix(0, nstate, nstate)
  +     if (i==1) cmat[indx] <- fit$cumhaz[i,]
  +     else cmat[indx] <- fit$cumhaz[i,] - fit$cumhaz[i-1,]
  + 
  +     diag(cmat) <- 1- rowSums(cmat)
  +     cmat
  + }
  > 
  > aeq(hazard(mfit3, 1), aj2(w1)[swap, swap])
  [1] TRUE
  > aeq(hazard(mfit3, 2), aj3(w1)[swap, swap])
  [1] TRUE
  > aeq(hazard(mfit3, 3), aj4(w1)[swap, swap])
  [1] TRUE
  > aeq(hazard(mfit3, 4), aj5(w1)[swap, swap])
  [1] TRUE
  > aeq(hazard(mfit3, 5), aj8(w1)[swap, swap])
  [1] TRUE
  > aeq(hazard(mfit3, 6), aj9(w1)[swap, swap])
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.39    0.20    1.53 
Running the tests in 'tests/multi3.R' failed.
Complete output:
  > library(survival)
  > aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > # Check that a multi-state model, correctly set up, gives the same
  > # solution as a time-dependent covariate.
  > # This is a stronger test than mstrata: there the covariate which was mapped
  > #  into a state was constant, here it is time-dependent.
  > #
  > # First build the TD data set from pbcseq, with a categorical bilirubin
  > pbc1 <- pbcseq
  > pbc1$bili4 <- cut(pbc1$bili, c(0,1, 2,4, 100), 
  +                   c("normal", "1-2x", "2-4x", ">4"))
  > ptemp <- subset(pbc1, !duplicated(id))  # first row of each
  > 
  > pbc2 <- tmerge(ptemp[, c("id", "age", "sex")], ptemp, id,
  +                death= event(futime, status==2))
  > 
  > pbc2 <- tmerge(pbc2, pbc1, id=id, bili = tdc(day, bili),
  +                  bili4 = tdc(day, bili4), bstat = event(day, as.numeric(bili4)))
  > btemp <- with(pbc2, ifelse(death, 5, bstat))
  > 
  > # a row with the same starting and ending bili4 level is not an event
  > b2 <- ifelse(((as.numeric(pbc2$bili4)) == btemp), 0, btemp)
  > pbc2$bstat <- factor(b2, 0:5,
  +                      c("censor", "normal", "1-2x", "2-4x", ">4", "death"))
  > check1 <- survcheck(Surv(tstart, tstop, bstat) ~ 1, istate= bili4,
  +                     id = id, data=pbc2)
  > check1$transitions
          to
  from     normal 1-2x 2-4x  >4 death (censored)
    normal      0   81   10   3     9         77
    1-2x       61    0   68  15     9         36
    2-4x        2   33    0  94    12         24
    >4          1    3   28   0   110         35
    death       0    0    0   0     0          0
  > all.equal(as.character(pbc2$bili4), as.character(check1$istate))
  [1] TRUE
  > # the above verifies that I created the data set correctly
  > 
  > # Standard coxph fit with a time dependent bili4 variable.
  > fit1 <- coxph(Surv(tstart, tstop, death) ~ age + bili4, pbc2)
  > 
  > # An additive multi-state fit, where bili4 is a state
  > #  The three forms below should all give identical models
  > fit2 <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
  +                    c(1:4):5 ~ age / common + shared), id= id, istate=bili4,
  +               data=pbc2)
  > fit2b <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
  +                    1:5 + 2:5 + 3:5 + 4:5 ~ age / common + shared), 
  +               id= id, istate=bili4, data=pbc2)
  > fit2c <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
  +                    0:5 ~ age / common + shared), 
  +               id= id, istate=bili4, data=pbc2)
  > 
  > # Make sure the names are correct and the coefficients match
  > aeq(coef(fit1), coef(fit2))
  [1] TRUE
  > aeq(names(coef(fit2)), c("age", "ph(2:5/1:5)", "ph(3:5/1:5)", "ph(4:5/1:5)"))
  [1] TRUE
  > all.equal(coef(fit2), coef(fit2b))
  [1] TRUE
  > all.equal(coef(fit2), coef(fit2c))
  [1] TRUE
  > 
  > # Now a model with a separate age effect for each bilirubin group
  > fit3  <- coxph(Surv(tstart, tstop, death) ~ age*bili4, pbc2)
  > fit3b <- coxph(Surv(tstart, tstop, death) ~ bili4/age, pbc2)
  > fit4 <-  coxph(list(Surv(tstart, tstop, bstat) ~ 1,
  +                    c(1:4):5 ~ age / shared), id= id, istate=bili4,
  +               data=pbc2)
  > all.equal(fit3$loglik, fit3b$loglik)
  [1] TRUE
  > all.equal(fit3$loglik, fit4$loglik)
  [1] TRUE
  > 
  > # The coefficients are quite different due to different codings for dummy vars
  > # Unpack the interaction, first 4 coefs will be the age effect within each
  > #  bilirubin group
  > temp <- c(coef(fit3)[1] + c(0, coef(fit3)[5:7]), coef(fit3)[2:4])
  > names(temp)[1:4] <- c("age1", "age2", "age3", "age4")
  > aeq(temp, coef(fit3b)[c(4:7, 1:3)])
  [1] TRUE
  > aeq(temp, coef(fit4))
  [1] TRUE
  > 
  > # Third, a model with separate baseline hazards for each bili group
  > fit5 <- coxph(Surv(tstart, tstop, death) ~ strata(bili4)/age, pbc2,
  +               cluster=id)
  > fit6 <- coxph(list(Surv(tstart, tstop, bstat) ~ 1, 0:5 ~ age),
  +                    id=id, istate=bili4, pbc2)
  > aeq(coef(fit5), coef(fit6))
  [1] TRUE
  > aeq(fit5$var, fit6$var)
  [1] TRUE
  > aeq(fit5$naive.var, fit6$naive.var)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.81    0.18    1.98 
Running the tests in 'tests/multistate.R' failed.
Complete output:
  > #
  > # Tests for multi-state Cox models
  > #
  > library(survival)
  > 
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > #  There are a few subjects with progression and death on the same day. In the
  > # usual multi-state data set only one will count. 
  > data1 <- mgus2
  > data1$etime <- with(data1, ifelse(pstat==1, ptime, futime))
  > data1$event <- factor(ifelse(data1$pstat==1, 1, 2L*data1$death),
  +                       0:2, c("censor", "PCM", "death"))
  > 
  > # direct data set with 2 rows per subject, much like mstate package would do
  > data2 <- mgus2[rep(1:nrow(mgus2) ,2), c("id", "age", "sex", "mspike")]
  > data2$time <- rep(data1$etime, 2)
  > data2$status <- 1* c(data1$event=="PCM", data1$event=="death")
  > data2$type   <- rep(c(2:3), each=nrow(mgus2))
  > 
  > fit1 <- coxph(Surv(etime, event) ~ age + sex + mspike, data1, id=id, x=TRUE,
  +               robust=FALSE)
  > fit1a <- coxph(Surv(etime, event=="PCM") ~ age + sex + mspike, data1)
  > fit1b <- coxph(Surv(etime, event=='death') ~ age + sex + mspike, data1)
  > fit1c <- coxph(Surv(time, status) ~ strata(type)/(age + sex+ mspike), 
  +                data2, x=TRUE)
  > 
  > aeq(fit1$loglik, fit1a$loglik + fit1b$loglik)
  [1] TRUE
  > aeq(fit1$coef, c(fit1a$coef, fit1b$coef))
  [1] TRUE
  > aeq(fit1$var[1:3, 1:3], fit1a$var)
  [1] TRUE
  > aeq(fit1$var[4:6, 4:6], fit1b$var)
  [1] TRUE
  > aeq(fit1$coef[c(1,4,2,5,3,6)], fit1c$coef)
  [1] TRUE
  > 
  > # force a common age effect across all states
  > fit2 <- coxph(list(Surv(etime, event) ~ sex,
  +                                   1:0 ~ age / common), 
  +               data1, id=id)
  > 
  > data2 <- rbind(cbind(data1, status= (data1$event=="PCM"), etype=1),
  +                cbind(data1, status= (data1$event=='death'), etype=2))
  > fit2a <- coxph(Surv(etime, status) ~ age + strata(etype)/sex, data2)
  > 
  > aeq(coef(fit2), coef(fit2a)[c(2,1,3)])  # not in the same order
  [1] TRUE
  > aeq(fit2$loglik, fit2a$loglik)
  [1] TRUE
  > 
  > #same fit in more complex ways
  > data1$entry <- "Entry"
  > fit2b <-  coxph(list(Surv(etime, event) ~ sex,
  +                      "Entry":"PCM" + "Entry":"death" ~ age / common),
  +                 istate=entry, data1, id=id)
  > fit2c <-  coxph(list(Surv(etime, event) ~ sex,
  +                      "Entry":state(c("PCM", "death")) ~ age / common),
  +                 istate=entry, data1, id=id)
  > 
  > aeq(fit2b$loglik, fit2$loglik)
  [1] TRUE
  > aeq(fit2c$coef, fit2$coef)
  [1] TRUE
  > 
  > # mspike size as a covariate for PCM only
  > # first, 4 different ways to write the same
  > fit3 <- coxph(list(Surv(etime, event) ~ age + sex,
  +                    1:state("PCM") ~ mspike),
  +               data1, id=id)
  > fit3b <- coxph(list(Surv(etime, event) ~ age + sex,
  +                    1:"PCM" ~ mspike),
  +               data1, id=id)
  > fit3c <- coxph(list(Surv(etime, event) ~ age + sex,
  +                    1:c("PCM") ~ mspike),
  +               data1, id=id)
  > fit3d <- coxph(list(Surv(etime, event) ~ age + sex + mspike,
  +                     1:3 ~ -mspike), data1, id=id)
  > 
  > aeq(fit3b$coef, fit3$coef)
  [1] TRUE
  > aeq(fit3c$coef, fit3$coef)
  [1] TRUE
  > aeq(fit3d$coef, fit3$coef)
  [1] TRUE
  > 
  > data3 <- data2
  > data3$mspike[data3$etype==2] <- 0
  > fit3a <-  coxph(Surv(etime, status) ~ strata(etype)/(age + sex + mspike), data3)
  > aeq(fit3$loglik, fit3a$loglik)
  [1] TRUE
  > aeq(fit3$coef, fit3a$coef[c(1,3,5,2,4)])
  [1] TRUE
  > 
  > # models with strata
  > test1 <-  coxph(Surv(etime, event=="PCM") ~ age + mspike + strata(sex), data1)
  > test2 <-  coxph(Surv(etime, event=="death") ~ age + strata(sex), data1)
  > 
  > sfit1 <-  coxph(list(Surv(etime, event) ~ age + strata(sex), 
  +                    1:state("PCM") ~ mspike),
  +               data1, id=id)
  > aeq(coef(sfit1), c(coef(test1), coef(test2)))
  [1] TRUE
  > 
  > test3 <- coxph(Surv(etime, event=="death") ~ age +sex, data1)
  > sfit2 <- coxph(list(Surv(etime, event) ~ age + sex,
  +                     1:2 ~ mspike + strata(sex) - sex),  data1, id=id)
  > aeq(coef(sfit2), c(coef(test1), coef(test3)))
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.25    0.35    2.61 
Running the tests in 'tests/neardate.R' failed.
Complete output:
  > library(survival)
  > # the second data set is not sorted by id/date, on purpose
  > 
  > df1 <- data.frame(id= 1:10,
  +                   y1= as.Date(c("1992-01-01", "1996-01-01", "1997-03-20",
  +                               "2000-01-01", "2001-01-01", "2004-01-01",
  +                               "2014-03-27", "2014-01-30", "2000-08-01",
  +                               "1997-04-29")))
  > 
  > df2 <- data.frame(id= c(1, 1, 2, 3, 4, 4, 5, 6, 7, 7, 8, 9, 9, 9, 10, 
  +                         3, 3, 6, 6, 8),
  +                   y2= as.Date(c("1998-04-30", "2004-07-01", "1999-04-14", 
  +                         "2001-02-22", "2003-11-19", "2005-02-15", "2006-06-22",
  +                         "2007-09-20", "2013-08-02", "2015-01-09", "2014-01-15",
  +                         "2006-12-06", "1999-10-20", "2010-06-30", "1997-04-28",
  +                         "1995-04-20", "1997-03-20", "1998-04-30", "1995-04-20",
  +                         "2006-12-06")))
  > 
  > if (FALSE) {  # plot for visual check
  +     plot(y2 ~ id, df2, ylim=range(c(df1$y1, df2$y2)), type='n')
  +     text(df2$id, df2$y2, as.numeric(1:nrow(df2)))
  +     points(y1~id, df1, col=2, pch='+')
  + }
  > 
  > i1 <- neardate(df1$id, df2$id, df1$y1, df2$y2)
  > all.equal(i1, c(1, 3, 17, 5, 7, 8, 10, NA, 12, NA))
  [1] TRUE
  > 
  > i2 <- neardate(df1$id, df2$id, df1$y1, df2$y2, best="prior")
  > all.equal(i2, c(NA, NA, 17, NA, NA, 18, 9, 11, 13, 15))
  [1] TRUE
  > 
  > indx <- order(df2$id, df2$y2)
  > df3 <- df2[indx,]
  > i3 <- neardate(df1$id, df3$id, df1$y1, df3$y2)
  > all.equal(indx[i3], i1)
  [1] TRUE
  > 
  > i4 <- neardate(df1$id, df3$id, df1$y1, df3$y2, best="prior")
  > all.equal(indx[i4], i2)
  [1] TRUE
  > 
  > indx <- c(2,3,10,9, 4,5, 7,8,1,6)
  > df4 <- df1[indx,]
  > i5 <-  neardate(df4$id, df2$id, df4$y1, df2$y2)
  > all.equal(i1[indx], i5)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.04    0.12    1.15 
Running the tests in 'tests/nsk.R' failed.
Complete output:
  > library(survival)
  > library(splines)
  > 
  > # the nsk function should give the same solution as ns, but with a different
  > #  parameterization
  > #
  > xx <- runif(500, 1, 100)
  > yy <- 10*log(xx) + rnorm(500, 0, 2)
  > tdata <- data.frame(xx=xx, yy=yy)
  > fit1 <- lm(yy ~ ns(xx, df=4), tdata, model=TRUE)
  > fit2 <- lm(yy ~ nsk(xx, df=4, b=0), tdata)
  > all.equal(predict(fit1), predict(fit2))  # same solution
  [1] TRUE
  > 
  > xattr <- attributes(fit1$model[[2]])
  > allknots <- sort(c(xattr$knots, xattr$Boundary.knots)) # knots that were used
  > pred.knot <- predict(fit1, newdata=list(xx=allknots))
  > all.equal(pred.knot[-1] - pred.knot[1], coef(fit2)[-1],
  +           check.attributes = FALSE)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.15    0.14    1.28 
Running the tests in 'tests/overlap.R' failed.
Complete output:
  > #
  > # Make sure that useless intervals do not cause issues, i.e., any that do
  > #  not overlap at least one event time
  > #
  > library(survival)
  > test2 <- data.frame(time1 =c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8, 3),
  +                     time2 =c(2, 3, 6, 7, 8, 9, 9, 9,14,17, 5),
  +                     event =c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0),
  +                     x     =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0, 500) )
  > 
  > # The data set is the same as book3.R, except for the wild observation
  > #  with x=500 whose time interval of (4,5) overlaps no events.
  > 
  > fit1 <- coxph(Surv(time1, time2, event) ~ x, test2, subset=(x<100))
  > fit2 <- coxph(Surv(time1, time2, event) ~ x, test2)
  > 
  > ii <- match(c("coefficients", "var", "loglik", "score", "iter", 
  +               "wald.test", "concordance"), names(fit1))
  > all.equal(fit1[ii], fit2[ii])
  [1] TRUE
  > all.equal(c(fit1$residuals,0), fit2$residuals, check.attributes=FALSE)
  [1] TRUE
  > 
  > # The mean differs condiderably, and so to the linear predictors
  > 
  > # Now the same with a penalized model
  > fit3 <- coxph(Surv(time1, time2, event) ~ ridge(x, theta=.1), test2,
  +               subset= (x< 100))
  > fit4 <-  coxph(Surv(time1, time2, event) ~ ridge(x, theta=.1), test2)
  > fit5 <- coxph(Surv(time1,time2, event) ~ x, test2,
  +               iter=0, init=fit4$coef)
  > 
  > all.equal(fit3[ii], fit4[ii])
  [1] TRUE
  > all.equal(c(fit3$residuals,0),  fit4$residuals, check.attributes=FALSE)
  [1] TRUE
  > all.equal(fit4$residuals, fit5$residuals)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.18    0.09    1.26 
Running the tests in 'tests/r_sas.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # Reproduce example 1 in the SAS lifereg documentation
  > #
  > 
  > # this fit doesn't give the same log-lik that they claim
  > fit1 <- survreg(Surv(time, status) ~ I(1000/(273.2+temp)), imotor,
  + 		subset=(temp>150), dist='lognormal')
  > summary(fit1)
  
  Call:
  survreg(formula = Surv(time, status) ~ I(1000/(273.2 + temp)), 
      data = imotor, subset = (temp > 150), dist = "lognormal")
                           Value Std. Error     z       p
  (Intercept)            -10.471      2.772 -3.78 0.00016
  I(1000/(273.2 + temp))   8.322      1.284  6.48 9.1e-11
  Log(scale)              -0.504      0.183 -2.75 0.00596
  
  Scale= 0.604 
  
  Log Normal distribution
  Loglik(model)= -145.9   Loglik(intercept only)= -155
  	Chisq= 18.3 on 1 degrees of freedom, p= 1.9e-05 
  Number of Newton-Raphson Iterations: 6 
  n= 30 
  
  > 
  > # This one, with the loglik on the transformed scale (the inappropriate
  > #   scale, Ripley & Venables would argue) does agree.
  > # All coefs are of course identical.
  > fit2 <- survreg(Surv(log(time), status) ~ I(1000/(273.2+temp)), imotor,
  + 		subset=(temp>150), dist='gaussian')
  > 
  > 
  > # Give the quantile estimates, which is the lower half of "output 48.1.5"
  > #  in the SAS 9.2 manual
  > 
  > pp1 <- predict(fit1, newdata=list(temp=c(130,150)), p=c(.1, .5, .9),
  + 		      type='quantile', se=T)
  > pp2 <- predict(fit1, newdata=list(temp=c(130,150)), p=c(.1, .5, .9),
  + 		      type='uquantile', se=T)
  > pp1
  $fit
            [,1]      [,2]     [,3]
  [1,] 12033.185 26095.677 56592.20
  [2,]  4536.877  9838.864 21336.98
  
  $se.fit
           [,1]      [,2]      [,3]
  [1,] 5482.338 11359.450 26036.917
  [2,] 1443.072  2901.155  7172.343
  
  > 
  > temp130 <- matrix(0, nrow=3, ncol=6)
  > temp130[,1] <- pp1$fit[1,]
  > temp130[,2] <- pp1$se.fit[1,]
  > temp130[,3] <- pp2$fit[1,]
  > temp130[,4] <- pp2$se.fit[1,]
  > temp130[,5] <- exp(pp2$fit[1,] - 1.64*pp2$se.fit[1,])
  > temp130[,6] <- exp(pp2$fit[1,] + 1.64*pp2$se.fit[1,])
  > dimnames(temp130) <- list(c("p=.1", "p=.2", "p=.3"),
  +      c("Time", "se(time)", "log(time)", "se[log(time)]", 
  +        "lower 90", "upper 90"))
  > print(temp130)
           Time  se(time) log(time) se[log(time)]  lower 90  upper 90
  p=.1 12033.18  5482.338  9.395424     0.4556015  5700.089  25402.68
  p=.2 26095.68 11359.450 10.169525     0.4353001 12779.950  53285.37
  p=.3 56592.20 26036.917 10.943626     0.4600796 26611.422 120349.71
  > 
  > # A set of examples, copied from the manual pages of SAS procedure
  > #  "reliability", which is part of their QC product.
  > #
  > 
  > color <- c("black", "red", "green", "blue", "magenta", "red4",
  +                         "orange", "DarkGreen", "cyan2", "DarkViolet")
  > palette(color)
  > pdf(file='reliability.pdf')
  > 
  > #
  > # Insulating fluids example
  > #
  > 
  > # Adding a -1 to the fit just causes the each group to have it's own
  > # intercept, rather than a global intercept + constrasts.  The strata
  > # statement allows each to have a separate scale
  > ffit <- survreg(Surv(time) ~ factor(voltage) + strata(voltage) -1, ifluid)
  > 
  > # Get predicted quantiles at each of the voltages
  > # By default predict() would give a line of results for each observation,
  > #  I only want the unique set of x's, i.e., only 4 cases
  > uvolt <- sort(unique(ifluid$voltage))      #the unique levels
  > plist <- c(1, 2, 5, 1:9 *10, 95, 99)/100
  > pred  <- predict(ffit, type='quantile', p=plist,
  +                  newdata=data.frame(voltage=uvolt))
  > tfun <- function(x) log(-log(1-x))
  > 
  > matplot(t(pred), tfun(plist), type='l', log='x', lty=1,
  +         col=1:4, yaxt='n',
  +         xlab="Predicted", ylab="Quantile")
  > axis(2, tfun(plist), format(100*plist), adj=1)
  > 
  > kfit <- survfit(Surv(time) ~ voltage, ifluid, type='fleming') #KM fit
  > for (i in 1:4) {
  +     temp <- kfit[i]
  +     points(temp$time, tfun(1-temp$surv), col=i, pch=i)
  +     }
  > 
  > # Now a table
  > temp <- array(0, dim=c(4,4,4))  #4 groups by 4 parameters by 4 stats
  > temp[,1,1] <- ffit$coef         # "EV Location" in SAS manual
  > temp[,2,1] <- ffit$scale        # "EV scale"
  > temp[,3,1] <- exp(ffit$coef)    # "Weibull Scale"
  > temp[,4,1] <- 1/ffit$scale      # "Weibull Shape"
  >  
  > temp[,1,2] <- sqrt(diag(ffit$var))[1:4]   #standard error
  > temp[,2,2] <- sqrt(diag(ffit$var))[5:8] * ffit$scale
  > temp[,3,2] <- temp[,1,2] * temp[,3,1]
  > temp[,4,2] <- temp[,2,2] / (temp[,2,1])^2
  > 
  > temp[,1,3] <- temp[,1,1] - 1.96*temp[,1,2]  #lower conf limits 
  > temp[,1,4]  <- temp[,1,1] + 1.96*temp[,1,2] # upper
  > # log(scale) is the natural parameter, in which the routine did its fitting
  > #  and on which the std errors were computed
  > temp[,2, 3] <- exp(log(ffit$scale) - 1.96*sqrt(diag(ffit$var))[5:8])
  > temp[,2, 4] <- exp(log(ffit$scale) + 1.96*sqrt(diag(ffit$var))[5:8])
  > 
  > temp[,3, 3:4] <- exp(temp[,1,3:4])
  > temp[,4, 3:4] <- 1/temp[,2,4:3]
  > 
  > dimnames(temp) <- list(uvolt, c("EV Location", "EV Scale", "Weibull scale",
  +                                 "Weibull shape"),
  +                        c("Estimate", "SE", "lower 95% CI", "uppper 95% CI"))
  > print(aperm(temp, c(2,3,1)), digits=5)
  , , 26
  
                 Estimate         SE lower 95% CI uppper 95% CI
  EV Location     6.86249    1.10404      4.69857        9.0264
  EV Scale        1.83423    0.96114      0.65677        5.1227
  Weibull scale 955.74665 1055.18620    109.78973     8320.0103
  Weibull shape   0.54519    0.28568      0.19521        1.5226
  
  , , 30
  
                Estimate       SE lower 95% CI uppper 95% CI
  EV Location    4.35133  0.30151      3.76037        4.9423
  EV Scale       0.94446  0.22544      0.59156        1.5079
  Weibull scale 77.58159 23.39176     42.96420      140.0911
  Weibull shape  1.05881  0.25274      0.66318        1.6904
  
  , , 34
  
                Estimate      SE lower 95% CI uppper 95% CI
  EV Location    2.50326 0.31476      1.88632        3.1202
  EV Scale       1.29732 0.22895      0.91796        1.8334
  Weibull scale 12.22222 3.84707      6.59509       22.6506
  Weibull shape  0.77082 0.13603      0.54542        1.0894
  
  , , 38
  
                  Estimate      SE lower 95% CI uppper 95% CI
  EV Location   0.00092629 0.27318     -0.53450       0.53635
  EV Scale      0.73367610 0.20380      0.42565       1.26460
  Weibull scale 1.00092672 0.27343      0.58596       1.70976
  Weibull shape 1.36299929 0.37861      0.79077       2.34933
  
  > 
  > rm(temp, uvolt, plist, pred, ffit, kfit) 
  > 
  > #####################################################################
  > # Turbine cracks data
  > crack2 <- with(cracks, data.frame(day1=c(NA, days), day2=c(days, NA),
  +                                   n=c(fail, 167-sum(fail))))
  > cfit <- survreg(Surv(day1, day2, type='interval2') ~1, 
  +                 dist='weibull', data=crack2, weight=n)
  > 
  > summary(cfit)
  
  Call:
  survreg(formula = Surv(day1, day2, type = "interval2") ~ 1, data = crack2, 
      weights = n, dist = "weibull")
                Value Std. Error      z       p
  (Intercept)  7.6880     0.0744 103.30 < 2e-16
  Log(scale)  -0.3953     0.0987  -4.01 6.2e-05
  
  Scale= 0.674 
  
  Weibull distribution
  Loglik(model)= -309.6   Loglik(intercept only)= -309.6
  Number of Newton-Raphson Iterations: 5 
  n= 9 
  
  > #Their output also has Wiebull scale = exp(cfit$coef), shape = 1/(cfit$scale)
  > 
  > # Draw the SAS plot
  > #  The "type=fleming" argument reflects that they estimate hazards rather than
  > #  survival, and forces a Nelson-Aalen hazard estimate
  > #
  > plist <-  c(1, 2, 5, 1:8 *10)/100
  > plot(qsurvreg(plist, cfit$coef, cfit$scale), tfun(plist), log='x',
  +      yaxt='n', type='l',
  +      xlab="Weibull Plot for Time", ylab="Percent")
  > axis(2, tfun(plist), format(100*plist), adj=1)
  > 
  > kfit <- survfit(Surv(day1, day2, type='interval2') ~1, data=crack2,
  +                 weight=n, type='fleming')
  > # Only plot point where n.event > 0 
  > # Why?  I'm trying to match them.  Personally, all should be plotted.
  > who <- (kfit$n.event > 0)
  > points(kfit$time[who], tfun(1-kfit$surv[who]), pch='+')
  > points(kfit$time[who], tfun(1-kfit$upper[who]), pch='-')
  > points(kfit$time[who], tfun(1-kfit$lower[who]), pch='-')
  > 
  > text(rep(3,6), seq(.5, -1.0, length=6), 
  +          c("Scale", "Shape", "Right Censored", "Left Censored", 
  +            "Interval Censored", "Fit"), adj=0)
  > text(rep(9,6), seq(.5, -1.0, length=6), 
  +          c(format(round(exp(cfit$coef), 2)),
  +            format(round(1/cfit$scale, 2)),
  +            format(tapply(crack2$n, cfit$y[,3], sum)), "ML"), adj=1)
  > 
  > # Now a portion of his percentiles table
  > #  I don't get the same SE as SAS, I haven't checked out why.  The
  > #  estimates and se for the underlying Weibull model are the same.
  > temp <- predict(cfit, type='quantile', p=plist, se=T)
  > tempse <- sqrt(temp$se[1,])
  > mat <- cbind(temp$fit[1,], tempse, 
  +              temp$fit[1,] -1.96*tempse, temp$fit[1,] + 1.96*tempse)
  > dimnames(mat) <- list(plist*100, c("Estimate", "SE", "Lower .95", "Upper .95"))
  > print(mat)
       Estimate        SE  Lower .95 Upper .95
  1    98.47183  5.321676   88.04134  108.9023
  2   157.59360  6.186246  145.46856  169.7186
  5   295.17088  7.377034  280.71189  309.6299
  10  479.31737  8.227489  463.19149  495.4432
  20  794.55351  8.951640  777.00829  812.0987
  30 1089.70373  9.404613 1071.27068 1108.1368
  40 1387.95516  9.983696 1368.38712 1407.5232
  50 1704.71015 10.888487 1683.36872 1726.0516
  60 2057.23906 12.217230 2033.29329 2081.1848
  70 2472.58593 14.035515 2445.07632 2500.0955
  80 3006.43572 16.510569 2974.07501 3038.7964
  > 
  > #
  > # The cracks data has a particularly easy estimate, so use
  > # it to double check code
  > time <- c(crack2$day2[1], (crack2$day1 + crack2$day2)[2:8]/2, 
  +           crack2$day1[9])
  > cdf  <- cumsum(crack2$n)/sum(crack2$n)
  > all.equal(kfit$time, time)
  [1] TRUE
  > all.equal(kfit$surv, 1-cdf[c(1:8,8)]) 
  [1] TRUE
  > rm(time, cdf, kfit)
  > 
  > 
  > #######################################################
  > #
  > # Replacement of valve seats in diesel engines
  > #   The input data has id, time, and an indicator of whether there was an
  > #   event at that time: 0=no, 1=yes.  No one has an event at their last time.
  > #  The input data has two engines with dual failures: 328 loses 2 valves at 
  > #    time 653, and number 402 loses 2 at time 139.  For each, fudge the first
  > #    time to be .1 days earlier. 
  > #
  > ties <- which(diff(valveSeat$time)==0 & diff(valveSeat$id)==0)
  > temp <- valveSeat$time
  > temp[ties] <- temp[ties] - .1
  > n <- length(temp)
  > first <- !duplicated(valveSeat$id)
  > vtemp <- with(valveSeat, data.frame(id =id, 
  +                                     time1= ifelse(first, 0, c(0, temp[-n])),
  +                                     time2= temp, status=status))
  > 
  > kfit <- survfit(Surv(time1, time2, status) ~1, vtemp, id=id)
  > 
  > plot(kfit, fun='cumhaz', ylab="Sample Mean Cumulative Failures", xlab='Time')
  > title("Valve replacement data")
  > 
  > # The summary.survfit function doesn't have an option for printing out
  > #   cumulative hazards instead of survival --- need to add that
  > #   so I just reprise the central code of print.summary.survfit
  > xx <- summary(kfit)
  > temp <- cbind(xx$time, xx$n.risk, xx$n.event, xx$cumhaz, 
  +               xx$std.chaz, -log(xx$upper), -log(xx$lower))
  > dimnames(temp) <- list(rep("", nrow(temp)),
  +                        c("time", "n.risk", "n.event", "Cum haz", "std.err",
  +                          "lower 95%", "upper 95%"))
  > print(temp, digits=2)
   time n.risk n.event Cum haz std.err lower 95% upper 95%
     61     41       1   0.024   0.024    0.0000     0.073
     76     41       1   0.049   0.034    0.0000     0.117
     84     41       1   0.073   0.041    0.0000     0.156
     87     41       1   0.098   0.046    0.0057     0.192
     92     41       1   0.122   0.051    0.0208     0.226
     98     41       1   0.146   0.055    0.0373     0.259
    120     41       1   0.171   0.059    0.0548     0.291
    139     41       1   0.195   0.062    0.0732     0.322
    139     41       1   0.220   0.073    0.0750     0.369
    165     41       1   0.244   0.075    0.0954     0.398
    166     41       1   0.268   0.077    0.1163     0.427
    202     41       1   0.293   0.079    0.1376     0.455
    206     41       1   0.317   0.088    0.1452     0.497
    249     41       1   0.341   0.089    0.1675     0.524
    254     41       1   0.366   0.090    0.1903     0.551
    258     41       1   0.390   0.090    0.2133     0.577
    265     41       1   0.415   0.091    0.2368     0.603
    276     41       1   0.439   0.098    0.2479     0.641
    298     41       1   0.463   0.110    0.2490     0.689
    323     41       1   0.488   0.110    0.2734     0.714
    326     41       1   0.512   0.110    0.2981     0.739
    328     41       1   0.537   0.115    0.3124     0.774
    344     41       1   0.561   0.115    0.3376     0.798
    348     41       1   0.585   0.124    0.3430     0.842
    349     41       1   0.610   0.124    0.3686     0.866
    367     41       1   0.634   0.123    0.3945     0.890
    377     41       1   0.659   0.132    0.4018     0.932
    404     40       1   0.684   0.136    0.4189     0.965
    408     40       1   0.709   0.140    0.4365     0.998
    410     40       1   0.734   0.139    0.4631     1.022
    449     40       1   0.759   0.143    0.4813     1.055
    479     40       1   0.784   0.146    0.4998     1.087
    497     40       1   0.809   0.149    0.5187     1.118
    538     40       1   0.834   0.152    0.5380     1.150
    539     40       1   0.859   0.155    0.5576     1.181
    561     40       1   0.884   0.162    0.5696     1.220
    563     40       1   0.909   0.164    0.5898     1.250
    570     40       1   0.934   0.170    0.6029     1.287
    573     40       1   0.959   0.169    0.6310     1.310
    581     38       1   0.985   0.171    0.6532     1.341
    586     34       1   1.014   0.174    0.6777     1.376
    604     22       1   1.060   0.185    0.7009     1.446
    621     17       1   1.119   0.208    0.7144     1.554
    635     16       1   1.181   0.207    0.7795     1.618
    640     16       1   1.244   0.226    0.8033     1.723
    646     13       1   1.320   0.229    0.8775     1.809
    653      9       1   1.432   0.252    0.9394     1.983
    653      9       1   1.543   0.312    0.9201     2.238
  > 
  > # Note that I have the same estimates but different SE's than SAS.  We are using 
  > #  a different estimator. It's a statistical argument as to which is
  > #  better (one could defend both sides): SAS the more standard estimate found
  > #  in the reliability literature, mine the estimate from statistics literature
  > rm(temp, kfit, xx)
  >                     
  > ######################################################
  > # Turbine data, lognormal fit
  > turbine <- read.table('data.turbine', 
  +                       col.names=c("time1", "time2", "n"))
  > 
  > tfit <- survreg(Surv(time1, time2, type='interval2') ~1, turbine,
  +                 dist='lognormal', weights=n, subset=(n>0))
  > 
  > summary(tfit)
  
  Call:
  survreg(formula = Surv(time1, time2, type = "interval2") ~ 1, 
      data = turbine, weights = n, subset = (n > 0), dist = "lognormal")
                Value Std. Error     z      p
  (Intercept)  3.6999     0.0708 52.23 <2e-16
  Log(scale)  -0.3287     0.1232 -2.67 0.0076
  
  Scale= 0.72 
  
  Log Normal distribution
  Loglik(model)= -190.7   Loglik(intercept only)= -190.7
  Number of Newton-Raphson Iterations: 6 
  n= 21 
  
  > 
  > # Now, do his plot, but put bootstrap confidence bands on it!
  > #  First, make a simple data set without weights
  > tdata <- turbine[rep(1:nrow(turbine), turbine$n),]
  > 
  > qstat <- function(data) {
  +     temp <- survreg(Surv(time1, time2, type='interval2') ~1, data=data,
  +                     dist='lognormal')
  +     qsurvreg(plist, temp$coef, temp$scale, dist='lognormal')
  +     }
  > 
  > {if (exists('bootstrap')) {
  +     set.seed(1953)  # a good year :-)
  +     bfit <- bootstrap(tdata, qstat, B=1000)
  +     bci <- limits.bca(bfit, probs=c(.025, .975))
  +     }
  + else {
  +     values <- matrix(0, nrow=1000, ncol=length(plist))
  +     n <- nrow(tdata)
  +     for (i in 1:1000) {
  +         subset <- sample(1:n, n, replace=T)
  +         values[i,] <- qstat(tdata[subset,])
  +         }
  +     bci <- t(apply(values,2, quantile, c(.05, .95)))
  +     }
  +  }
  > xmat <- cbind(qsurvreg(plist, tfit$coef, tfit$scale, dist='lognormal'),
  +               bci)
  > 
  > 
  > matplot(xmat, qnorm(plist), 
  +         type='l', lty=c(1,2,2), col=c(1,1,1), 
  +         log='x', yaxt='n', ylab='Percent', 
  +         xlab='Time of Cracking (Hours x 100)')
  > axis(2, qnorm(plist), format(100*plist), adj=1)
  > title("Turbine Data")
  > kfit <- survfit(Surv(time1, time2, type='interval2') ~1, data=tdata)
  > points(kfit$time, qnorm(1-kfit$surv), pch='+')
  > 
  > dev.off()  #close the plot file
  null device 
            1 
  > 
  > 
  > proc.time()
     user  system elapsed 
     4.81    0.60    5.43 
Running the tests in 'tests/r_tdist.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # Test out the t-distribution
  > #
  > # First, a t-dist with 500 df should be nearly identical to the Gaussian
  > 
  > fitig <- survreg(Surv(time, status)~voltage, 
  +         dist = "gaussian", data = capacitor)
  > fit1 <- survreg(Surv(time, status) ~ voltage,
  + 		 dist='t', parms=500, capacitor)
  > fitig
  Call:
  survreg(formula = Surv(time, status) ~ voltage, data = capacitor, 
      dist = "gaussian")
  
  Coefficients:
  (Intercept)     voltage 
  1917.431558   -4.087366 
  
  Scale= 314.4748 
  
  Loglik(model)= -249.2   Loglik(intercept only)= -258.7
  	Chisq= 18.97 on 1 degrees of freedom, p= 1.32e-05 
  n= 64 
  > summary(fit1, corr=F)
  
  Call:
  survreg(formula = Surv(time, status) ~ voltage, data = capacitor, 
      dist = "t", parms = 500)
                 Value Std. Error     z       p
  (Intercept) 1918.329    219.536  8.74 < 2e-16
  voltage       -4.091      0.778 -5.26 1.5e-07
  Log(scale)     5.749      0.133 43.12 < 2e-16
  
  Scale= 314 
  
  Student-t distribution: parmameters= 500
  Loglik(model)= -249.2   Loglik(intercept only)= -258.7
  	Chisq= 19.01 on 1 degrees of freedom, p= 1.3e-05 
  Number of Newton-Raphson Iterations: 5 
  n= 64 
  
  > 
  > # A more realistic fit
  > fit2 <- survreg(Surv(time, status) ~ voltage,
  + 		 dist='t', parms=5, capacitor)
  > print(fit2)
  Call:
  survreg(formula = Surv(time, status) ~ voltage, data = capacitor, 
      dist = "t", parms = 5)
  
  Coefficients:
  (Intercept)     voltage 
  1975.155951   -4.315492 
  
  Scale= 276.799 
  
  Loglik(model)= -250.2   Loglik(intercept only)= -261
  	Chisq= 21.65 on 1 degrees of freedom, p= 3.27e-06 
  n= 64 
  > 
  > if (FALSE) {
  + resid(fit2, type='response')
  + resid(fit2, type='deviance')
  + resid(fit2, type='working') 
  + resid(fit2, type='dfbeta')
  + resid(fit2, type='dfbetas')
  + resid(fit2, type='ldresp')
  + resid(fit2, type='ldshape')
  + resid(fit2, type='ldcase')
  + resid(fit2, type='matrix')
  + 
  + predict(fit2, type='link')
  + predict(fit2, type='terms')
  + predict(fit2, type='quantile')
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.21    0.20    1.40 
Running the tests in 'tests/rttright.R' failed.
Complete output:
  > library(survival)
  > 
  > # start with the example used in chapter 2 of the book
  > 
  > bdata <- data.frame(time =   c(1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 
  +                                9, 10,11, 12,14, 15, 16, 16, 18, 20),
  +                     status = c(1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1,
  +                                0, 0, 1, 0, 0, 1, 0, 1, 0))
  > 
  > # First check: verify that the the RTTR reproduces the KM
  > kfit <- survfit(Surv(time, status) ~1, bdata)
  > bwt  <- rttright(Surv(time, status) ~1, bdata, renorm= FALSE)
  > 
  > cdf <- cumsum(bwt)/nrow(bdata)  # weighted CDF
  > cdf <- cdf[!duplicated(bdata$time, fromLast=TRUE)]  # remove duplicates
  > all.equal(kfit$surv, 1-cdf)
  [1] TRUE
  > 
  > 
  > # A covariate divides both survfit and rttr into disjoint groups, so repeat
  > #  the above check on subsets of the aml data
  > afit <- survfit(Surv(time, status) ~x, aml)
  > awt <-  rttright(Surv(time, status) ~x, aml, renorm=TRUE)
  > 
  > igroup <- as.numeric(aml$x)
  > for (i in 1:2) {
  +     atemp <- awt[igroup ==i]   # subset for this curve
  +     index <- order(aml$time[igroup ==i])
  +     acdf <- cumsum(atemp[index])
  +     acdf <- acdf[!duplicated(aml$time[igroup ==i], fromLast=TRUE)]
  +     print(all.equal(afit[i]$surv, 1-acdf))
  + }
  [1] TRUE
  [1] TRUE
  > 
  > 
  > ###########
  > # Alternate computation using inverse prob of censoring weights.
  > # First shift the censorings to avoid ties: if there is a death and a censor
  > #   at time 10, say, the death was not at risk of censoring. Censoring weights
  > #   happen "later".  This also results in a left-continuous curve.
  > delta <- min(diff(sort(unique(bdata$time)))) /3
  > offset <- ifelse(bdata$status==1, 0, delta)
  > cfit <- survfit(Surv(time+ offset, 1-status) ~ 1, bdata)
  > 
  > # interpolate
  > indx <- findInterval(bdata$time, cfit$time)
  > cwt <- ifelse(bdata$status==0, 0, 1/cfit$surv[indx])
  > all.equal(bwt, cwt) 
  [1] TRUE
  > 
  > # Multiple time points, this example is used in the vignette
  > tdata <- data.frame(time=   c(1,2,2,3,4,4,5,5,8,9),
  +                     status= c(1,1,0,1,0,0,1,0,1,1))
  > fit1 <- rttright(Surv(time, status) ~ 1, tdata, times=2:6, renorm=FALSE)
  > fit2 <- rttright(Surv(time, status) ~ 1, tdata, times=2:6, renorm=TRUE)
  > all.equal(fit1, 10*fit2)
  [1] TRUE
  > all.equal(fit1, cbind(7, c(7,7,0,8,8,8,8,8,8,8), 
  +                          c(7,7,0,8,8,8,8,8,8,8), 
  +                          c(7,7,0,8,0,0,12,12,12,12),
  +                          c(7,7,0,8,0,0,12, 0, 18,18))/7, check.attributes=FALSE)
  [1] TRUE
  > 
  > # Now test with (start, stop] data, should get the same results
  > b2 <- survSplit(Surv(time, status) ~ 1, bdata, cut= c(3,5, 7, 14),
  +                 id = "subject")
  > indx <- c(seq(1, 65, by=2), seq(64, 2, by= -2))
  > b2 <- b2[indx,]    # not in time within subject order (stronger test)
  > 
  > b2wt <- rttright(Surv(tstart, time, status) ~1, b2, id=subject)
  > indx2 <- order(b2$time)
  > cdf2 <- cumsum(b2wt[indx2])
  > cdf2 <- cdf2[!duplicated(b2$time[indx2], fromLast=TRUE)] # remove duplicates
  > utime2 <- sort(unique(b2$time))   # will have an extra time 7
  > utime1 <- sort(unique(bdata$time))
  > all.equal(cdf2[match(utime1, utime2)], cdf)
  [1] TRUE
  > 
  > 
  > # Competing risks
  > mdata <- mgus2
  > mdata$etime <- with(mgus2, ifelse(pstat==1, ptime, futime))
  > mdata$estat <- with(mgus2, ifelse(pstat==1, 1, 2*death))
  > mdata$estat <- factor(mdata$estat, 0:2, c('censor', 'pcm', 'death'))
  > mfit <- survfit(Surv(etime, estat) ~1, mdata, id=id)
  > mwt1 <- rttright(Surv(etime, estat) ~1, mdata, id=id)
  > 
  > morder <- order(mdata$etime)
  > mdata2 <- mdata[morder,]
  > mwt2   <- rttright(Surv(etime,estat) ~1, mdata2, id=id)
  > all.equal(mwt1[morder], mwt2)
  [1] TRUE
  > 
  > keep <- !duplicated(mdata2$etime, fromLast=TRUE)
  > csum1 <- cumsum(ifelse(mdata2$estat=="pcm", mwt2, 0))
  > csum2 <- cumsum(ifelse(mdata2$estat=="death", mwt2, 0))
  > 
  > all.equal(mfit$pstate[,2], csum1[keep])
  [1] TRUE
  > all.equal(mfit$pstate[,3], csum2[keep])
  [1] TRUE
  > 
  > # Case weights, at multiple times
  > bwt <- rep(1:2, length=nrow(bdata))
  > tm <- c(2, 6, 10, 15, 18)
  > fit1 <- rttright(Surv(time, status) ~1, bdata, weights=bwt, times= tm)
  > casefit <- survfit(Surv(time, status) ~ 1, bdata, weights= bwt)
  > csum1 <- summary(casefit, censor=FALSE, times= tm)
  > for (i in 1:length(tm)) {
  +     c1 <- sum(fit1[bdata$status==1 & bdata$time <= tm[i], i])
  +     print(all.equal(c1, 1-csum1$surv[i]))
  + }
  [1] TRUE
  [1] TRUE
  [1] TRUE
  [1] TRUE
  [1] TRUE
  > 
  > 
  > ###
  > # Delayed entry, tiny data set
  > #  Still to be done
  > delay <- data.frame(t0=c(0,0,0,0,3,0),
  +                     t1=1:6,
  +                     status=c(1,0,1,0,0,1),
  +                     id=1:6)
  > # dwt <- rttright(Surv(t0, t1, status) ~ 1, delay, id=id, times=0:5 + .9)
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.62    0.07    1.68 
Running the tests in 'tests/surv.R' failed.
Complete output:
  > #
  > library(survival)
  > 
  > # Some simple tests of the Surv function
  > #  The first two are motivated by a bug, pointed out by Kevin Buhr,
  > #    where a mixture of NAs and invalid values didn't work right
  > #  Even for the simplest things a test case is good.
  > #  All but the third should produce warning messages
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > temp <- Surv(c(1, 10, 20, 30), c(2, NA, 0, 40), c(1,1,1,1))
  Warning message:
  In Surv(c(1, 10, 20, 30), c(2, NA, 0, 40), c(1, 1, 1, 1)) :
    Stop time must be > start time, NA created
  > aeq(temp, c(1,10,NA,30,  2,NA,0,40, 1,1,1,1))
  [1] TRUE
  > 
  > temp <- Surv(c(1, 10, 20, 30), c(2, NA, 0, 40), type='interval2')
  Warning message:
  In Surv(c(1, 10, 20, 30), c(2, NA, 0, 40), type = "interval2") :
    Invalid interval: start > stop, NA created
  > aeq(temp, c(1,10,20,30,  2,1,1,40, 3,0,NA,3))
  [1] TRUE
  > 
  > #No error
  > temp <- Surv(1:5)
  > aeq(temp, c(1:5, 1,1,1,1,1))
  [1] TRUE
  > 
  > temp1 <- Surv(c(1,10,NA, 30, 30), c(1,NA,10,20, 40), type='interval2')
  Warning message:
  In Surv(c(1, 10, NA, 30, 30), c(1, NA, 10, 20, 40), type = "interval2") :
    Invalid interval: start > stop, NA created
  > temp2 <- Surv(c(1,10,10,30,30), c(9, NA, 5, 20,40), c(1, 0, 2,3,3),
  +               type='interval')
  Warning message:
  In Surv(c(1, 10, 10, 30, 30), c(9, NA, 5, 20, 40), c(1, 0, 2, 3,  :
    Invalid interval: start > stop, NA created
  > aeq(temp1, temp2)
  [1] TRUE
  > aeq(temp1, c(1,10,10,30,30, 1,1,1,1, 40, 1,0,2,NA,3))
  [1] TRUE
  > 
  > # Use of inf
  > temp1 <- Surv(c(1,10,NA, 30, 30), c(1,NA,10,30, 40), type='interval2')
  > temp2 <- Surv(c(1,10,-Inf, 30, 30), c(1,Inf,10,30, 40), type='interval2')
  > aeq(temp1, temp2)
  [1] TRUE
  > 
  > # Verify sorting and order routines
  > #  These fail in 3.4, succeed in 3.5 due to a system change in how
  > #  xtfrm.Surv is used.  
  > x1 <- Surv(c(4, 6, 3, 2, 1, NA, 2), c(1,0, NA, 0,1,1,1))
  > all.equal(order(x1), c(5,7, 4, 1, 2, 3, 6))
  [1] TRUE
  > all.equal(order(x1, decreasing=TRUE), c(2,1,4,7,5, 3, 6))
  [1] TRUE
  > all.equal(sort(x1), x1[c(5,7,4,1,2)])
  [1] TRUE
  > 
  > x2 <- Surv(c(4, 6, 3, 2, 1, NA, 2), c(1,0, NA, 0,1,1,1), type='left')
  > all.equal(order(x2), c(5,4, 7, 1, 2, 3, 6))
  [1] TRUE
  > 
  > x3 <- Surv(c(1,5,NA,7, 9), c(6, 6, 4, NA, 9), type="interval2")
  > all.equal(sort(x3), x3[c(1,3,2,4,5)])
  [1] TRUE
  > 
  > x4 <- Surv(c(1,5,6,5,2, 4), c(3, 7, 7, 6, 3, NA), factor(c(1, 2, 0, 1, 1, 0)))
  > all.equal(sort(x4), x4[c(1, 5, 4, 2,3)])
  [1] TRUE
  > all.equal(sort(x4, na.last=FALSE), x4[c(6,1,5,4,2,3)])
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.20    0.20    1.39 
Running the tests in 'tests/survSplit.R' failed.
Complete output:
  > library(survival)
  > # Make sure that the old-style and new-style calls both work
  > 
  > # new style
  > vet2 <- survSplit(Surv(time, status) ~ ., data= veteran, cut=c(90, 180), 
  +                   episode= "tgroup", id="id")
  > vet2[1:7, c("id", "tstart", "time", "status", "tgroup", "age", "karno")]
    id tstart time status tgroup age karno
  1  1      0   72      1      1  69    60
  2  2      0   90      0      1  64    70
  3  2     90  180      0      2  64    70
  4  2    180  411      1      3  64    70
  5  3      0   90      0      1  38    60
  6  3     90  180      0      2  38    60
  7  3    180  228      1      3  38    60
  > 
  > # old style
  > vet3 <- survSplit(veteran, end='time', event='status', cut=c(90,180),
  +                   episode="tgroup", id="id")
  > all.equal(vet2, vet3)
  [1] TRUE
  > 
  > all.equal(nrow(vet2), nrow(veteran) + sum(veteran$time >90) + 
  +                       sum(veteran$time > 180))
  [1] TRUE
  > 
  > 
  > # Do a parallel computation using survSplit, and pyears/tcut.  We should get
  > #  the same answer.
  > # Break subjects up by year of entry and current time on study.  Most of
  > #  the deaths are within 3 months
  > # pyears complains (justifiably) about the obs with 0 days of fu, so we add 1
  > #
  > data1 <- jasa
  > data1$ayear <- as.numeric(substring(as.character(jasa$accept.dt), 1,4))
  > temp <- round(c(0, .25, 1,2,5)*365.25)   # years
  > ftime <- tcut(rep(0, nrow(jasa)), temp, 
  +                   labels=paste(c(0, .25, 1:2), c(.25, 1,2,5), sep='-'))
  > pfit <- pyears(Surv(futime +1,fustat) ~ ayear + ftime, data1,
  +                scale=1)
  > 
  > data2 <- survSplit(Surv(futime+1, fustat) ~ ., cut=temp, data=data1,
  +                    episode = "tgroup")
  > 
  > tab1 <- with(data2, tapply(fustat, list(ayear, tgroup), sum))
  > tab1 <- ifelse(is.na(tab1), 0, tab1)
  > 
  > all.equal(as.vector(tab1), as.vector(pfit$event))  # ignore dimnames
  [1] TRUE
  > 
  > tab2 <- with(data2, tapply(tstop-tstart, list(ayear, tgroup), sum))
  > tab2 <- ifelse(is.na(tab2), 0, tab2)
  > 
  > all.equal(as.vector(tab2), as.vector(pfit$pyears))
  [1] TRUE
  > 
  > # double check that the "data" option gives the same values
  > pfit2 <- pyears(Surv(futime +1,fustat) ~ ayear + ftime, data1,
  +                scale=1, data.frame=TRUE)$data
  > all.equal(pfit2$pyears, pfit$pyears[pfit$pyears >0])
  [1] TRUE
  > all.equal(pfit2$event,  pfit$event[pfit$pyears >0])
  [1] TRUE
  > 
  > # and that the rows of data2 have the right labels
  > keep <- which(pfit$pyears >0)  # these are not in the data
  > rname <- rownames(pfit$pyears)[row(pfit$pyears)[keep]]
  > all.equal(rname, as.character(pfit2$ayear))
  [1] TRUE
  > 
  > cname <- colnames(pfit$pyears)[col(pfit$pyears)[keep]]
  > all.equal(cname, as.character(pfit2$ftime))
  [1] TRUE
  > 
  >                    
  > 
  > proc.time()
     user  system elapsed 
     1.09    0.12    1.21 
Running the tests in 'tests/survfit1.R' failed.
Complete output:
  > #
  > # Check out the survfit routine on the simple AML data set.
  > #  The leverage validation makes use of the fact that when all
  > #  weights are 1 and there is 1 obs per subject, the IJ variance is
  > #  equal to the Greenwood.
  > # There are 8 choices in the C code:  Nelson-Aalen or Fleming-Harrington
  > #  estimate of cumulative hazard,  KM or exp(cumhaz) estimate of survival,
  > #  regular or robust variance.  This tries to exercise them all.
  > 
  > library(survival)
  > aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > set.seed(1953)  # used only to reorder the data
  > adata <- aml
  > adata$id <- sample(LETTERS, nrow(aml)) # labels are not in time or data order
  > adata <- adata[sample(1:nrow(aml), nrow(aml)),] # data is unordered
  > adata$wt <- sample((2:30)/10, nrow(aml))  # non-integer weights
  > 
  > byhand <- function(time, status, weights, id) {
  +     # for a single curve
  +     utime <- sort(unique(time))
  +     ntime <- length(utime)
  +     n <- length(time)
  +     if (missing(weights)) weights <- rep(1.0, n)
  +     if (missing(id)) id <- seq(along=time)
  + 
  +     uid <- sort(unique(id))
  +     nid <- length(uid)
  +     id <- match(id, uid)  # change it to 1:nid
  + 
  +     n.risk <- n.event <- surv <- cumhaz <- double(ntime)
  +     KM <- 1; nelson <-0; 
  +     kvar <- 0; hvar<-0;
  +         
  +     U <- matrix(0, nid, 2)  # the two robust influence estimates
  +     V <- matrix(0, ntime, 4)  # variances
  +     usave <- array(0., dim=c(nid, 2, ntime))
  +     estimate <- matrix(0, ntime, 2)
  + 
  +     for (i in 1:ntime) {
  +         atrisk <- (time >= utime[i])
  +         n.risk[i] <- sum(weights[atrisk])
  +         deaths <- (time==utime[i] & status==1)                 
  +         n.event[i] <- sum(weights[deaths])
  +  
  +         haz <- n.event[i]/n.risk[i]
  +         dhaz <- (ifelse(deaths,1,0) - ifelse(atrisk, haz, 0))/n.risk[i]
  +         U[,1] <- U[,1]*(1-haz) - KM*tapply(dhaz*weights, id, sum)
  +         V[i,1] <- sum(U[,1]^2)
  +             
  +         U[,2] <- U[,2] + tapply(dhaz* weights, id, sum) #result in 'id' order
  +         V[i,2] <- sum(U[,2]^2)
  +         usave[,,i] <- U
  + 
  +         if (n.event[i] >0 ) {
  +             KM <- KM*(1-haz)
  +             nelson <- nelson + haz
  +             kvar <- kvar + n.event[i]/(n.risk[i] * (n.risk[i] - n.event[i]))
  +             hvar <- hvar + n.event[i]/(n.risk[i]^2)
  +         }
  +             
  +         V[i,3] <- kvar   # var of log(S)
  +         V[i,4] <- hvar
  +         estimate[i,] <- c(KM, nelson)
  +         }
  +     dimnames(usave) <- list(uid, c("KM", "chaz"), utime)
  +     list(time=utime, n.risk=n.risk, n.event=n.event, estimate=estimate,
  +          std = sqrt(V), influence=usave)
  + }
  > 
  > # the byhand function can only handle one group at a time
  > true1a <- with(subset(adata, x=="Maintained"), byhand(time, status, id=id))
  > true1b <- with(subset(adata, x!="Maintained"), byhand(time, status, id=id))
  > 
  > # The Greenwood and IJ estimates agree, except for a last point with
  > #  variance of zero.  These next few lines verify the byhand() function
  > aeq(true1a$std[,1], true1a$estimate[,1]*true1a$std[,3])
  [1] TRUE
  > aeq(true1b$std[1:9,1], true1b$estimate[1:9,1]*true1b$std[1:9,3])
  [1] TRUE
  > aeq(true1b$std[10,1], 0)   # variance of zero for jackknife
  [1] TRUE
  > !is.finite(true1b$std[10,3])   # Inf for Greenwood
  [1] TRUE
  > temp <- with(subset(adata, x=="Maintained"), byhand(time, status, id=id,
  +                                                     weights=rep(3,11)))
  > aeq(temp$std[,1:2], true1a$std[,1:2])  # IJ estimates should be invariant
  [1] TRUE
  > 
  > # fit1 uses the standard formulas: NA hazard, KM survival
  > fit1 <- survfit(Surv(time, status) ~ x, data=adata)
  > aeq(fit1$surv, c(true1a$estimate[,1], true1b$estimate[,1]))
  [1] TRUE
  > aeq(fit1$cumhaz, c(true1a$estimate[,2], true1b$estimate[,2]))
  [1] TRUE
  > aeq(fit1$std.err, c(true1a$std[,3], true1b$std[,3]))
  [1] TRUE
  > aeq(fit1$std.chaz, c(true1a$std[,4], true1b$std[,4]))
  [1] TRUE
  > aeq(fit1$n.risk, c(true1a$n.risk, true1b$n.risk))
  [1] TRUE
  > aeq(fit1$n.event, c(true1a$n.event, true1b$n.event))
  [1] TRUE
  > fit1$logse   # logse should be TRUE
  [1] TRUE
  > 
  > # fit2 will use the IJ method
  > fit2 <- survfit(Surv(time, status) ~ x, data=adata, id=id, robust=TRUE)
  > aeq(fit2$surv, c(true1a$estimate[,1], true1b$estimate[,1]))
  [1] TRUE
  > aeq(fit2$cumhaz, c(true1a$estimate[,2], true1b$estimate[,2]))
  [1] TRUE
  > aeq(fit2$std.err, c(true1a$std[,1], true1b$std[,1]))
  [1] TRUE
  > aeq(fit2$std.chaz, c(true1a$std[,2], true1b$std[,2]))
  [1] TRUE
  > aeq(fit2$n.risk, c(true1a$n.risk, true1b$n.risk))
  [1] TRUE
  > aeq(fit2$n.event, c(true1a$n.event, true1b$n.event))
  [1] TRUE
  > !fit2$logse  # logse should be FALSE
  [1] TRUE
  > 
  > # look at the leverage values
  > fit3 <- survfit(Surv(time, status) ~ x, data=adata, id=id, influence=3)
  > aeq(fit3$influence.surv[[1]], true1a$influence[,1,])
  [1] TRUE
  > aeq(fit3$influence.surv[[2]], true1b$influence[,1,])
  [1] TRUE
  > aeq(fit3$influence.chaz[[1]], true1a$influence[,2,])
  [1] TRUE
  > aeq(fit3$influence.chaz[[2]], true1b$influence[,2,])
  [1] TRUE
  > 
  > # compute the influence by brute force
  > tdata <- subset(adata, x != "Maintained")
  > tdata <- tdata[order(tdata$id),]   # easier to compare if it's in order
  > eps <- 1e-8
  > imat1 <- imat2 <-  matrix(0., 12, 10)
  > t1 <- survfit(Surv(time, status) ~x, data=tdata) 
  > for (i in 1:12) {
  +     wtemp <- rep(1.0, 12)
  +     wtemp[i] <- 1 + eps
  +     tfit <-survfit(Surv(time, status) ~x, data=tdata, weights=wtemp) 
  +     imat2[i,] <- (tfit$cumhaz - t1$cumhaz)/eps
  +     imat1[i,] <- (tfit$surv - t1$surv)/eps
  + }
  > aeq(imat1, true1b$influence[,1,], tol= sqrt(eps))
  [1] TRUE
  > aeq(imat2, true1b$influence[,2,], tol= sqrt(eps))
  [1] TRUE
  > 
  > # Repeat using the Nelson-Aalen hazard and exp(NA) for survival
  > fit1 <- survfit(Surv(time, status) ~ x, adata, stype=2)
  > aeq(fit1$surv, exp(-c(true1a$estimate[,2], true1b$estimate[,2])))
  [1] TRUE
  > aeq(fit1$cumhaz, c(true1a$estimate[,2], true1b$estimate[,2]))
  [1] TRUE
  > aeq(fit1$std.err, c(true1a$std[,4], true1b$std[,4]))
  [1] TRUE
  > aeq(fit1$std.chaz, c(true1a$std[,4], true1b$std[,4]))
  [1] TRUE
  > aeq(fit1$n.risk, c(true1a$n.risk, true1b$n.risk))
  [1] TRUE
  > 
  > # Nelson-Aalen + exp() surv, along with IJ variance
  > fit2 <- survfit(Surv(time, status) ~ x, data=adata, id=id, stype=2,
  +                 influence=3)
  > aeq(fit2$surv, exp(-c(true1a$estimate[,2], true1b$estimate[,2])))
  [1] TRUE
  > aeq(fit2$cumhaz, c(true1a$estimate[,2], true1b$estimate[,2]))
  [1] TRUE
  > aeq(fit2$std.err, c(true1a$std[,2], true1b$std[,2]))
  [1] TRUE
  > aeq(fit2$std.chaz, c(true1a$std[,2], true1b$std[,2]))
  [1] TRUE
  > aeq(fit2$n.risk, c(true1a$n.risk, true1b$n.risk))
  [1] TRUE
  > aeq(fit2$influence.chaz[[1]], true1a$influence[,2,])
  [1] TRUE
  > aeq(fit2$influence.chaz[[2]], true1b$influence[,2,])
  [1] TRUE
  > aeq(fit2$influence.surv[[2]], -true1b$influence[,2,]%*% diag(fit2[2]$surv))
  [1] TRUE
  > 
  > # Cumulative hazard is the same for fit1 and fit2
  > all.equal(fit2$influence.chaz, fit3$influence.chaz)
  [1] TRUE
  > 
  > # Weighted fits
  > true2a <- with(subset(adata, x=="Maintained"), byhand(time, status, id=id,
  +                                                       weights= wt))
  > true2b <- with(subset(adata, x!="Maintained"), byhand(time, status, id=id,
  +                                                       weights=wt))
  > fit3 <- survfit(Surv(time, status) ~ x, data=adata, id=id, weights=wt,
  +                  influence=TRUE)  
  > aeq(fit3$influence.surv[[1]], true2a$influence[,1,])
  [1] TRUE
  > aeq(fit3$influence.surv[[2]], true2b$influence[,1,])
  [1] TRUE
  > aeq(fit3$influence.chaz[[1]], true2a$influence[,2,])
  [1] TRUE
  > aeq(fit3$influence.chaz[[2]], true2b$influence[,2,])
  [1] TRUE
  > aeq(fit3$surv, c(true2a$estimate[,1], true2b$estimate[,1]))
  [1] TRUE
  > aeq(fit3$cumhaz, c(true2a$estimate[,2], true2b$estimate[,2]))
  [1] TRUE
  > aeq(fit3$std.err, c(true2a$std[,1], true2b$std[,1]))
  [1] TRUE
  > aeq(fit3$std.chaz, c(true2a$std[,2], true2b$std[,2]))
  [1] TRUE
  > aeq(fit3$n.risk, c(true2a$n.risk, true2b$n.risk))
  [1] TRUE
  > aeq(fit3$n.event, c(true2a$n.event, true2b$n.event))
  [1] TRUE
  > 
  > # Different survival, same hazard
  > fit3b <- survfit(Surv(time, status) ~ x, data=adata, id=id, weights=wt,
  +                  influence=2, stype=2) 
  > temp <- c("n", "time", "cumhaz", "std.chaz", "influence.chaz", "n.risk",
  +           "n.event")
  > aeq(unclass(fit3b)[temp], unclass(fit3)[temp])  # unclass avoids [.survfit
  [1] TRUE
  > aeq(fit3b$surv, exp(-c(true2a$estimate[,2], true2b$estimate[,2])))
  [1] TRUE
  > aeq(fit3b$std.err, fit3b$std.chaz)
  [1] TRUE
  > aeq(fit3b$logse, FALSE)
  [1] TRUE
  > aeq(fit3b$n.risk, c(true2a$n.risk, true2b$n.risk))
  [1] TRUE
  > aeq(fit3b$n.event, c(true2a$n.event, true2b$n.event))
  [1] TRUE
  > 
  > # The grouped jackknife
  > group <- rep("", nrow(adata))
  > temp <- table(adata$x)
  > group[adata$x == "Maintained"] <- rep(letters[1:4], length=temp[1])
  > group[adata$x != "Maintained"] <- rep(letters[4:7], length=temp[2])
  > adata$group <- group
  > fit4 <-  survfit(Surv(time, status) ~ x, data=adata, id=id, weights=wt,
  +                  influence=TRUE, cluster=group)
  > g1 <- adata$group[match(rownames(true2a$influence[,1,]), adata$id)]
  > g2 <- adata$group[match(rownames(true2b$influence[,1,]), adata$id)] 
  > aeq(fit4$influence.surv[[1]], rowsum(true2a$influence[,1,], g1))
  [1] TRUE
  > aeq(fit4$influence.surv[[2]], rowsum(true2b$influence[,1,], g2))
  [1] TRUE
  > aeq(fit4$influence.chaz[[1]], rowsum(true2a$influence[,2,], g1))
  [1] TRUE
  > aeq(fit4$influence.chaz[[2]], rowsum(true2b$influence[,2,], g2))
  [1] TRUE
  > 
  > aeq(c(colSums(fit4$influence.surv[[1]]^2), colSums(fit4$influence.surv[[2]]^2)),
  +     fit4$std.err^2)
  [1] TRUE
  > aeq(c(colSums(fit4$influence.chaz[[1]]^2), colSums(fit4$influence.chaz[[2]]^2)),
  +     fit4$std.chaz^2)
  [1] TRUE
  > 
  > # The Fleming-Harrington is a more complex formula.  Start with weights of
  > #   1.
  > fit5 <- survfit(Surv(time, status) ~x, adata, ctype=2)
  > nrisk <- c(11,10,8,7, 5,4,2, 12, 11, 10, 9, 8, 6:1)
  > chaz <- c(cumsum(1/nrisk[1:7])[c(1:4,4, 5,6,6,7,7)], 
  +           cumsum(1/nrisk[8:18])[c(2,4,5,5,6:11)])
  > aeq(fit5$cumhaz, chaz)
  [1] TRUE
  > aeq(fit5$std.chaz, sqrt(c(cumsum(1/nrisk[1:7]^2)[c(1:4,4, 5,6,6,7,7)], 
  +                           cumsum(1/nrisk[8:18]^2)[c(2,4,5,5,6:11)])))
  [1] TRUE
  > 
  > # We can compute the FH using a fake data set where each tie is spread out
  > #  over a set of fake times.
  > # 
  > fh <- function(time, status, weights, id) {
  +     counts <- table(time, status)
  +     utime <-  sort(unique(time))
  +     tied <- counts[,2] > 1
  + 
  +     if (missing(weights)) weights <- rep(1.0, length(time))
  +     if (missing(id))  id <- 1:length(time)
  + 
  +     # build the expanded data set
  +     delta <- min(diff(utime))/(2*max(counts[,2]))
  +     efun <- function(x) {
  +         who <- which(time==x & status==1)
  +         ntie <- length(who)
  +         data.frame(time = rep(x - (1:ntie -1)*delta, each=ntie),
  +                    id = rep(id[who], ntie),
  +                    status = rep(1, ntie^2),
  +                    weight = rep(weights[who]/ntie, ntie),
  +                    stringsAsFactors=FALSE
  +                    )
  +     }
  + 
  +     temp <- do.call(rbind, lapply(utime[tied], efun))
  +     notie <- (status==0 | !(time %in% utime[tied]))
  + 
  +     bfit <- byhand(time = c(time[notie], temp$time), 
  +                    status = c(status[notie], temp$status),
  +                    id = c(id[notie], temp$id),
  +                    weights = c(weights[notie], temp$weight)
  +                    )
  +     keep <- match(utime, bfit$time)  # the real time points
  + 
  +     list(time=bfit$time[keep], 
  +          n.risk=bfit$n.risk[keep - pmax(0, counts[,2]-1)],
  +          n.event = bfit$n.event[keep]* counts[,2],  
  +          estimate=bfit$estimate[keep,],
  +          std = bfit$std[keep,], influence=bfit$influence[,,keep])
  + }
  > 
  > # Case weights
  > true6a <- with(subset(adata, x=="Maintained"), fh(time, status, wt, id))
  > true6b <- with(subset(adata, x!="Maintained"), fh(time, status, wt, id))
  > 
  > fit6 <- survfit(Surv(time, status) ~ x, weight=wt, data=adata, stype=2, 
  +                 ctype=2, robust=FALSE)
  > aeq(fit6$cumhaz, c(true6a$estimate[,2], true6b$estimate[,2]))
  [1] TRUE
  > aeq(fit6$surv, exp(-c(true6a$estimate[,2], true6b$estimate[,2])))
  [1] TRUE
  > aeq(fit6$std.chaz, c(true6a$std[,4], true6b$std[,4]))
  [1] TRUE
  > aeq(fit6$n.risk, c(true6a$n.risk, true6b$n.risk))
  [1] TRUE
  > aeq(fit6$n.event, c(true6a$n.event, true6b$n.event))
  [1] TRUE
  > 
  > # Robust variance
  > fit7 <- survfit(Surv(time, status) ~ x, weight=wt, data=adata, stype=2,ctype=2, 
  +                 id=id, influence=2, robust=TRUE)
  > aeq(fit7$cumhaz, c(true6a$estimate[,2], true6b$estimate[,2]))
  [1] TRUE
  > aeq(fit7$surv, exp(-c(true6a$estimate[,2], true6b$estimate[,2])))
  [1] TRUE
  > aeq(fit7$std.chaz, c(true6a$std[,2], true6b$std[,2]))
  [1] TRUE
  > aeq(fit7$n.risk, c(true6a$n.risk, true6b$n.risk))
  [1] TRUE
  > aeq(fit7$n.event, c(true6a$n.event, true6b$n.event))
  [1] TRUE
  > aeq(fit7$influence.chaz[[1]], true6a$influence[,2,])
  [1] TRUE
  > aeq(fit7$influence.chaz[[2]], true6b$influence[,2,])
  [1] TRUE
  > 
  > # compute the influence by brute force
  > tdata <- subset(adata, x != "Maintained")
  > tdata <- tdata[order(tdata$id),]
  > eps <- 1e-8
  > imat <- matrix(0., 12, 10)
  > t1 <- survfit(Surv(time, status) ~x, data=tdata, ctype=2, weights=wt) 
  > for (i in 1:12) {
  +     wtemp <- tdata$wt
  +     wtemp[i] <- wtemp[i] + eps
  +     tfit <-survfit(Surv(time, status) ~x, data=tdata, ctype=2, 
  +               weights=wtemp)
  +     imat[i,] <- tdata$wt[i] * (tfit$cumhaz - t1$cumhaz)/eps
  + }
  > aeq(fit7$influence.chaz[[2]], imat, tol=sqrt(eps))
  [1] TRUE
  > 
  > #
  > # verify that the times and scale arguments work as expected.  They
  > #  are in the summary and print.survfit functions.
  > #
  > s1 <- summary(fit1, scale=1)
  > s2 <- summary(fit1, scale=2)
  > aeq(s1$time/2, s2$time)  #times change
  [1] TRUE
  > aeq(s1$surv, s2$surv)
  [1] TRUE
  > tscale <- rep(c(1,1,1,1, 2,2,2,2,2), each=2)  
  > aeq(s1$table, s2$table *tscale)
  [1] TRUE
  > 
  > s3 <- summary(fit1, scale=1, times=c(9, 18, 23, 33, 34))
  > s4 <- summary(fit1, scale=2, times=c(9, 18, 23, 33, 34))
  > aeq(s3$time, s4$time*2)
  [1] TRUE
  > aeq(s3$surv, s4$surv)
  [1] TRUE
  > 
  > print(fit1, rmean='common')
  Call: survfit(formula = Surv(time, status) ~ x, data = adata, stype = 2)
  
                   n events rmean* se(rmean) median 0.95LCL 0.95UCL
  x=Maintained    11      7   60.3     25.60     34      23      NA
  x=Nonmaintained 12     11   30.1      9.14     27       8      NA
      * restricted mean with upper limit =  161 
  > print(fit1, rmean='common', scale=2)
  Call: survfit(formula = Surv(time, status) ~ x, data = adata, stype = 2)
  
                   n events rmean* se(rmean) median 0.95LCL 0.95UCL
  x=Maintained    11      7   30.2     12.80   17.0    11.5      NA
  x=Nonmaintained 12     11   15.1      4.57   13.5     4.0      NA
      * restricted mean with upper limit =  80.5 
  > 
  > proc.time()
     user  system elapsed 
     1.25    0.26    1.50 
Running the tests in 'tests/survfit2.R' failed.
Complete output:
  > library(survival)
  > #
  > # Check out the Dory&Korn confidence interval option
  > #
  > tdata <- data.frame(time= 1:10,
  +                     status=c(1,0,1,0,1,0,0,0,1,0))
  > 
  > fit1 <- survfit(Surv(time, status) ~1, tdata, conf.lower='modified')
  > fit2 <- survfit(Surv(time, status) ~1, tdata)
  > 
  > stdlow <- fit2$std.err * sqrt(c(1, 10/9, 1, 8/7, 1, 6/5, 6/4, 6/3, 1, 2/1))
  > lower <- exp(log(fit2$surv) - qnorm(.975)*stdlow)
  > all.equal(fit1$lower, lower, check.attributes=FALSE)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.26    0.18    1.42 
Running the tests in 'tests/survtest.R' failed.
Complete output:
  > options(na.action=na.exclude) # preserve missings
  > options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
  > library(survival)
  > 
  > #
  > # Simple test of (start, stop] Kaplan-Meier curves, using the test2 data
  > #   set
  > #
  > test1 <- data.frame(time=  c(9, 3,1,1,6,6,8),
  +                     status=c(1,NA,1,0,1,1,0),
  +                     x=     c(0, 2,1,1,1,0,0))
  > test2 <- data.frame(start=c(1, 2, 5, 2, 1, 7, 3, 4, 8, 8),
  +                     stop =c(2, 3, 6, 7, 8, 9, 9, 9,14,17),
  +                     event=c(1, 1, 1, 1, 1, 1, 1, 0, 0, 0),
  +                     x    =c(1, 0, 0, 1, 0, 1, 1, 1, 0, 0),
  +                     wt   = 1:10)
  > aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > test2 <- test2[c(1,6,2,7,3,8,4,9,5,10),]  # unsorted data is a harder test
  > 
  > fit1 <- survfit(Surv(start, stop, event) ~1, test2, type='fh2',
  +                 error='tsiatis')
  > fit2 <- survfit(Surv(start, stop, event) ~x, test2, start.time=3,
  + 		type='fh2')
  > 
  > cfit1<- survfit(coxph(Surv(start, stop, event)~1, test2))
  > cfit2<- survfit(coxph(Surv(start, stop, event) ~ strata(x), test2, subset=-1))
  > 
  > deaths <- (fit1$n.event + fit1$n.censor)>0
  > aeq(fit1$time[deaths], cfit1$time)
  [1] TRUE
  > aeq(fit1$n.risk[deaths], cfit1$n.risk)
  [1] TRUE
  > aeq(fit1$n.event[deaths], cfit1$n.event)
  [1] TRUE
  > aeq(fit1$surv[deaths], cfit1$surv)
  [1] TRUE
  > aeq(fit1$std.err[deaths], cfit1$std.err)
  [1] TRUE
  > 
  > deaths <- (fit2$n.event + fit2$n.censor)>0
  > aeq(fit2$time[deaths], cfit2$time)
  [1] TRUE
  > aeq(fit2$n.risk[deaths], cfit2$n.risk)
  [1] TRUE
  > aeq(fit2$n.event[deaths], cfit2$n.event)
  [1] TRUE
  > aeq(fit2$surv[deaths], cfit2$surv)
  [1] TRUE
  > 
  > fit3 <- survfit(Surv(start, stop, event) ~1, test2) #Kaplan-Meier
  > aeq(fit3$n, 10)
  [1] TRUE
  > aeq(fit3$time, sort(unique(test2$stop)))
  [1] TRUE
  > aeq(fit3$n.risk, c(2,3,5,4,4,5,2,1))
  [1] TRUE
  > aeq(fit3$n.event,c(1,1,1,1,1,2,0,0))
  [1] TRUE
  > aeq(fit3$surv[fit3$n.event>0], c(.5, 1/3, 4/15, 1/5, 3/20, 9/100))
  [1] TRUE
  > temp <- with(fit3, n.event/(n.risk * (n.risk - n.event)))
  > aeq(fit3$std.err, sqrt(cumsum(temp)))
  [1] TRUE
  > 
  > #
  > #  Verify that both surv AND n.risk are right between time points.
  > #
  > fit <- survfit(Surv(time, status) ~1, test1)
  > temp <- summary(fit, time=c(.5,1, 1.5, 6, 7.5, 8, 8.9, 9, 10), extend=TRUE)
  > 
  > aeq(temp$n.risk, c(6,6,4,4,2,2,1,1,0))
  [1] TRUE
  > aeq(temp$surv, c(1, fit$surv[c(1,1,2,2,3,3,4,4)]))
  [1] TRUE
  > aeq(temp$n.event, c(0,1,0,2,0,0,0,1,0))
  [1] TRUE
  > aeq(temp$std.err, c(0, (fit$surv*fit$std.err)[c(1,1,2,2,3,3,4,4)]))
  [1] TRUE
  > 
  > 
  > fit <- survfit(Surv(start, stop, event) ~1, test2)
  > temp <- summary(fit, times=c(.5, 1.5, 2.5, 3, 6.5, 14.5, 16.5))
  > aeq(temp$surv, c(1, 1, fit$surv[c(1,2,3,6,6)]))
  [1] TRUE
  > 
  > # This next fails.  With start-stop data the number at risk at intermediate
  > #  endpoints is not known precisely, since the underlying routine does not report
  > #  time points at which only an addition occured. 
  > if (FALSE) aeq(temp$n.risk, c(0, 2, 3, 3, 4, 1,1))
  > 
  > # compute conditional survival
  > fit1 <- survfit(Surv(start, stop, event)~1, test2, weights=wt)
  > fit2 <- survfit(Surv(start, stop, event)~1, test2, weights=wt, start.time=5)
  > 
  > aeq(fit1$surv[2], summary(fit1, time=5)$surv)  # verify my subscript
  [1] TRUE
  > aeq(fit2$surv, fit1$surv[3:8]/fit1$surv[2])
  [1] TRUE
  > aeq(fit2$std.err^2,   fit1$std.err[3:8]^2 - fit1$std.err[2]^2)
  [1] TRUE
  > aeq(fit2$cumhaz, fit1$cumhaz[3:8] - fit1$cumhaz[2])
  [1] TRUE
  > aeq(fit2$std.chaz^2,   fit1$std.chaz[3:8]^2 - fit1$std.chaz[2]^2)
  [1] TRUE
  > 
  > # Now with a Cox model
  > cfit <- coxph(Surv(start, stop, event)~1, test2, weights=wt)
  > fit1 <- survfit(cfit)
  > fit2 <- survfit(cfit, start.time=5)
  > aeq(fit2$surv, fit1$surv[3:8]/fit1$surv[2])
  [1] TRUE
  > aeq(fit2$std.err^2,   fit1$std.err[3:8]^2 - fit1$std.err[2]^2)
  [1] TRUE
  > aeq(fit2$cumhaz, fit1$cumhaz[3:8] - fit1$cumhaz[2])
  [1] TRUE
  > aeq(fit2$std.chaz^2,   fit1$std.chaz[3:8]^2 - fit1$std.chaz[2]^2)
  [1] TRUE
  > 
  > 
  > # bigger data set, with covariates and some tied event times
  > mfit <- coxph(Surv(age, age+futime/12, death) ~ sex + mspike, mgus2)
  > dummy <- data.frame(sex='F', mspike=1.3)
  > 
  > msurv1 <- survfit(mfit, newdata=dummy)
  > msurv2 <- survfit(mfit, newdata=dummy, start.time=80)
  > j <- max(which(msurv1$time < 80))
  > k <- seq(j+1, length(msurv1$time))
  > aeq(msurv2$surv, msurv1$surv[k] / msurv1$surv[j])
  [1] TRUE
  > aeq(msurv2$cumhaz, msurv1$cumhaz[k] - msurv1$cumhaz[j])
  [1] TRUE
  > # standard errors now have a term due to vmat(mfit), so don't factor
  > #  into a simple sum
  > 
  > proc.time()
     user  system elapsed 
     1.15    0.14    1.25 
Running the tests in 'tests/tiedtime.R' failed.
Complete output:
  > library(survival)
  > 
  > #
  > # The survival code was failing for certain data sets when called as
  > #    survfit(Surv(time2-time1, status) ~ ......
  > # The issue was how tied floating point numbers are handled, and the
  > #    fact that unique(x), factor(x) and tapply(x) are not guarranteed to
  > #    all be the same.  
  > # This test fails in survival 2.36-5, fixed in 2.36-6.  Data sets that
  > #    can cause it are few and far between.
  > #
  > 
  > load('ties.rda')
  > x <- time2 -time1
  > 
  > #  Here is the heart of the old problem
  > #       length(unique(x))== length(table(x))
  > #  And the prior fix which worked ALMOST always
  > #       x <- round(x, 15)
  > #       length(unique(round(x,15)))== length(table(round(x,15)))
  > 
  > fit1 <- survfit(Surv(x) ~1)
  > length(fit1$time) == length(fit1$surv)
  [1] TRUE
  > 
  > 
  > # a second test, once "rounding.R"
  > 
  > tdata <- data.frame(time=c(1,2, sqrt(2)^2, 2, sqrt(2)^2),
  +                     status=rep(1,5), 
  +                     group=c(1,1,1,2,2))
  > fit <- survfit(Surv(time, status) ~ group, data=tdata)
  > 
  > all.equal(sum(fit$strata), length(fit$time))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.01    0.06    1.06 
Running the tests in 'tests/tmerge2.R' failed.
Complete output:
  > library(survival)
  > 
  > # This test is based on a user report that a 0/1 variable would not reset
  > #  to zero.  It turned out to be a bug when data2 was not sorted
  > 
  > baseline <- data.frame(idd=1:5,  futime=c(20, 30, 40, 30, 20),
  +                        status= c(0, 1, 0, 1, 0))
  > tests <- data.frame(idd = c(2,3,3,3,4,4,5),
  +                     date = c(25, -1, 15, 23, 17, 19, 14),
  +                     onoff= c( 1, 1, 0, 1, 1, 0, 1))
  > tests <- tests[c(7,2,6,3,4,1,5),]  #scramble data2
  > 
  > mydata <- tmerge(baseline, baseline, id=idd, death=event(futime, status))
  > mydata <- tmerge(mydata, tests, id=idd, ondrug=tdc(date, onoff))
  > 
  > all.equal(mydata$ondrug, c(NA, NA,1, 1,0,1, NA, 1,0, NA, 1))
  [1] TRUE
  > 
  > 
  > # Check out addition of a factor, character, and logical
  > tests$ff <- factor(tests$onoff, 0:1, letters[4:5])
  > tests$fchar <- as.character(tests$ff)
  > tests$logic <- as.logical(tests$onoff)
  > tests$num <- rep(1:3, length=nrow(tests))
  > 
  > mydata <- tmerge(mydata, tests, id=idd, fgrp= tdc(date, ff),
  +                  chgrp = tdc(date, fchar), 
  +                  options=list(tdcstart="new"))
  > all.equal(mydata$fgrp, 
  +           factor(c(3,3,2,2,1,2,3,2,1,3,2), labels=c("d", "e", "new")))
  [1] TRUE
  > all.equal(mydata$chgrp, 
  +           c("d", "e", "new")[c(3,3,2,2,1,2,3,2,1,3,2)])
  [1] TRUE
  > 
  > mydat2  <-  tmerge(mydata, tests, id=idd, 
  +                  logic1 = tdc(date, logic), logic2= event(date, logic))
  > all.equal(mydat2$logic1, c(FALSE, TRUE, NA)[as.numeric(mydat2$fgrp)])
  [1] TRUE
  > all.equal(mydat2$logic2, c(FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, TRUE,
  +                            FALSE, FALSE, TRUE, FALSE))
  [1] TRUE
  > 
  > mydat3 <- tmerge(mydata, tests, id=idd,
  +                  xx = tdc(date, num), options=list(tdcstart=5))
  > all.equal(mydat3$xx, c(5,5,3,2,1,2,5,1,3,5,1))
  [1] TRUE
  > temp <- tmerge(mydata, tests, id=idd, xx=tdc(date, num, 5)) # alternate default
  > all.equal(mydat3$xx, temp$xx)
  [1] TRUE
  > 
  > # Multiple chained calls. 
  > temp <- outer(cgd0$id, 100*0:6, "+")
  > colnames(temp) <- paste0("x", 1:7)   # add a time dependent covariate too 
  > test <- cbind(cgd0, temp)
  > 
  > newcgd <- tmerge(data1=cgd0[, 1:13], data2=cgd0, id=id, tstop=futime)
  > newcgd <- tmerge(newcgd, test, id=id, infect = event(etime1), xx=cumtdc(etime1, x1, 0)) 
  > newcgd <- tmerge(newcgd, test, id=id, infect = event(etime2), xx=cumtdc(etime2, x2)) 
  > newcgd <- tmerge(newcgd, test, id=id, infect = event(etime3), xx=cumtdc(etime3, x3)) 
  > newcgd <- tmerge(newcgd, test, id=id, infect = event(etime4), xx=cumtdc(etime4, x4)) 
  > newcgd <- tmerge(newcgd, test, id=id, infect = event(etime5), xx=cumtdc(etime5, x5)) 
  > newcgd <- tmerge(newcgd, test, id=id, infect = event(etime6), xx=cumtdc(etime6, x6)) 
  > newcgd <- tmerge(newcgd, test, id=id, infect = event(etime7), xx=cumtdc(etime7, x7)) 
  > newcgd <- tmerge(newcgd, newcgd, id, enum=cumtdc(tstart))
  > all.equal(dim(newcgd), c(203,18))
  [1] TRUE
  > all.equal(as.vector(table(newcgd$infect)), c(127, 76))
  [1] TRUE
  > temp <- with(newcgd, ifelse(enum==1, 0, id + (enum-2)*100))
  > temp2 <- tapply(temp, newcgd$id, cumsum)
  > all.equal(newcgd$xx, unlist(temp2), check.attributes=FALSE)
  [1] TRUE
  > tcount <- attr(newcgd, "tcount")
  > all(tcount[,1:3] ==0)  # no early, late, or gap
  [1] TRUE
  > 
  > # table with number of subjects who have etime1 < futime (row 1)
  > #  and  etime1==futime (row 2)
  > # the table command ignores the missings
  > temp <- subset(cgd0, select=etime1:etime7)
  > counts <- sapply(temp, function(x) 
  +     as.vector(table(factor(x>= cgd0$futime, c(FALSE, TRUE)))))
  > 
  > all(tcount[c(1,3,5,7,9,11,13), c("within", "trailing")] == t(counts))
  [1] TRUE
  > 
  > 
  > #
  > # Merging with a date as the time variable.  In this case tstart/tstop are required
  > #  A default start of 0 has no meaning
  > #
  > base2 <- baseline
  > base2$date1 <- as.Date("1953-03-10")   # everyone enrolled that day
  > base2$date2 <- as.Date("1953-03-10") + base2$futime
  > base2$futime <- NULL
  > test2 <- tests
  > test2$date <- as.Date("1953-03-10") + test2$date
  > 
  > mydata2 <- tmerge(base2, base2, id=idd, death=event(date2, status), 
  +                  tstart = date1, tstop= date2,
  +                  options=list(tstartname="date1", tstopname="date2"))
  > mydata2 <- tmerge(mydata2, test2, id=idd, ondrug=tdc(date, onoff))
  > all.equal(mydata$ondrug, c(NA, NA,1, 1,0,1, NA, 1,0, NA, 1))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.06    0.07    1.12 
Running the tests in 'tests/tmerge3.R' failed.
Complete output:
  > # Understanding edge cases
  > library(survival)
  > 
  > #
  > # this is from a user report of a problem with cumevents.  When there is
  > #  a row merged in that is a censor, don't mark it as a cumevent.
  > #
  > base <- data.frame(
  +   id = 1:2, tstart = c(0, 0), tstop  = c(10, 10), got_flu = c(0, 0), 
  +   has_flu = factor(c("no", "no"), levels = c("no", "yes")))
  > base <- tmerge(base, base, id = id, got_flu = event(tstop, got_flu))
  > 
  > # add time-varying covariates
  > vars <- data.frame(id = c(1, rep(2, 5)), time = c(0, (0:4) * 2), x = rnorm(6))
  > base <- tmerge(base, vars, id = id, x = tdc(time, x))
  > 
  > # add cumevents, using a covariate
  > events <- data.frame(
  +   id = c(2, 2, 2), 
  +   # notice the zero -- the second row should not add an event
  +   got_flu = c(1,0,2), 
  +   has_flu = c("yes", "no", "yes"),
  +   time = c(3, 5, 8))
  > b2 <- tmerge(base, events, id = id, got_flu = cumevent(time, got_flu),
  +                has_flu = tdc(time, has_flu))
  Warning message:
  In tmerge(base, events, id = id, got_flu = cumevent(time, got_flu),  :
    replacement of variable 'has_flu'
  > 
  > all.equal(b2$got_flu, c(0,0,1,0,0,0,3,0))
  [1] TRUE
  > 
  > 
  > # Tied times in the merger data set
  > # for all of them missings are essentially ignored
  > # last obs wins for tdc and event
  > tiedat <- data.frame(id=c(1, 1, 1, 2,2,2), time=c(3,4, 4, 3, 5, 5), 
  +                      x=c(1, NA,0, 2,3,4))
  > b3 <- tmerge(base, tiedat, id=id, x1= tdc(time, x), x2=cumtdc(time, x),
  +              x3= event(time, x), x4 = cumevent(time, x))
  > all.equal(b3$x1, c(NA, 1, 0, NA, NA, 2,2, 4,4,4))
  [1] TRUE
  > all.equal(b3$x2, c(NA, 1, 1, NA, NA, 2,2, 9,9,9))
  [1] TRUE
  > all.equal(b3$x3, c(1,0,0,0,2,0,4,0,0,0))
  [1] TRUE
  > all.equal(b3$x4, c(1,0,0,0,2,0,9,0,0,0))
  [1] TRUE
  > 
  > # Multiple overlapping time windows in the first step.
  > #  Should generate an error message
  > test <- tryCatch(
  +             {tmerge(pbcseq[, c("id", "trt", "age", "sex")], pbcseq, id,
  +                death = event(futime, status==2))},
  +             error= function(cond) {
  +                 if (grepl("duplicate identifiers", cond)) 
  +                     cat("successful tmerge error test\n")
  +             }
  + )
  successful tmerge error test
  > 
  > # Using a tdc that depends on more than one variable.  If they are not
  > #  exactly the same class, tmerge should fail.
  > # Happens with wide data sets
  > 
  > tdata <- data.frame(id= 1:3, age=c(40,44,38), dtime=c(700, 600, 500),  
  +                     t1 = c(111, 211, 311), x1= as.integer(c(4, 5, 6)),
  +                     t2 = c(120, 240, 400.3), x2=c( 9, 8, 7),
  +                     t3 = c(400, 500, 450), x3=c(12,2, 0))
  > # This works
  > wide1 <- tmerge(tdata[,1:2], tdata, id=id, death= event(dtime),
  +                 x = tdc(t1, x1),  x= tdc(t2, x2), x= tdc(t3, x3))
  > 
  > r1 <- data.frame(id=rep(1:3, each=4), 
  +                  age= tdata$age[rep(1:3, each=4)],
  +                  tstart=c(0,111, 120, 400, 0, 211, 240, 500, 0, 311,400.3, 450),
  +                  tstop =c(111, 120, 400, 700, 211, 240, 500, 600, 
  +                           311, 400.3, 450, 500),
  +                  death= rep(c(0,0,0,1), 3),
  +                  x= c(NA,4, 9,12, NA, 5, 8, 2, NA, 6,7, 0))
  > all.equal(r1, wide1, check.attributes=FALSE)
  [1] TRUE
  > 
  > tdata$x2[2] <- 'c'  # different data type
  > test <- tryCatch(
  +             {tmerge(tdata[,1:2], tdata, id=id, death= event(dtime),
  +                 x = tdc(t1, x1),  x= tdc(t2, x2), x= tdc(t3, x3))},
  +              error= function(cond) {
  +                 if (grepl("tdc update does not match prior variable type: x", cond)) 
  +                     cat("successful tmerge error test\n")
  +             }
  + )
  successful tmerge error test
  > 
  > proc.time()
     user  system elapsed 
     0.84    0.06    0.89 
Running the tests in 'tests/tt2.R' failed.
Complete output:
  > # A reprise of tt.R, using (time1, time2) data.
  > library(survival)
  > library(splines)
  > aeq <- function(x, y) all.equal(as.vector(x), as.vector(y))
  > 
  > # A contrived example for the tt function
  > #
  > mkdata <- function(n, beta) {
  +     age <- round(runif(n, 20, 60))
  +     x <- rbinom(n, 1, .5)
  + 
  +     futime <- rep(40, n)   # everyone has 40 years of follow-up
  +     entry  <- pmax(0, seq(-10, 30, length=n))  # 1/4 enter at 0
  +     entry  <- round(entry)
  +     status <- rep(0, n)
  +     dtime <-  runif(n/2, 1, 40)  # 1/2 of them die
  +     dtime <- sort(dtime)
  + 
  +     # The risk is set to beta[1]*x + beta[2]* f(current_age)
  +     #   where f= 0 up to age 40, rises linear to age 70, flat after that
  +     for (i in 1:length(dtime)) {
  +         atrisk <- (futime >= dtime[i] & entry < dtime[i])
  +         c.age <- age + dtime
  +         age2 <- pmin(30, pmax(0, c.age-40))
  +         xbeta <- beta[1]*x + beta[2]*age2
  +         
  +         # Select a death according to risk
  +         risk <- ifelse(atrisk, exp(xbeta), 0)
  +         dead <- sample(1:n, 1, prob=risk/sum(risk))
  +         
  +         futime[dead] <- dtime[i]
  +         status[dead] <- 1
  +     }
  +     out <- data.frame(time1= entry, time2=round(futime,1), status=status, 
  +                       age=age, x=x, risk=risk,
  +                casewt = sample(1:5, n, replace=TRUE),
  +                grp = sample(1:15, n, replace=TRUE), id= 1:n)
  +     subset(out, time1 < time2)
  + }
  > 
  > set.seed(1953)  # a good year
  > # Make n larger for the (time1, time2) case; more stress.
  > tdata <- mkdata(250, c(log(1.5), 2/30))   # data set has many ties
  > #tdata <- mkdata(100, c(log(1.5), 2/30))   # data set has many ties
  > tdata$strat <- floor(tdata$grp/10)
  > 
  > dtime <- sort(unique(tdata$time2[tdata$status==1]))
  > data2 <- survSplit(Surv(time1, time2, status) ~., tdata, cut=dtime)
  > data2$c.age <- data2$age + data2$time2  # current age
  > 
  > # fit1 uses data at the event times, fit2$c.age might have a 
  > #  wider range due to censorings.  To make the two fits agree
  > #  fix the knots.  I know a priori that 20 to 101 will cover it.
  > ns2 <- function(x) ns(x, Boundary.knots=c(20, 101), knots=c(45, 60, 75))
  > 
  > fit1 <- coxph(Surv(time1, time2, status)~ x + tt(age), tdata,
  +               tt= function(x, t, ...) ns2(x+t))
  > 
  > fit2 <- coxph(Surv(time1, time2, status) ~ x + ns2(c.age), data2)
  > 
  > aeq(coef(fit1), coef(fit2))
  [1] TRUE
  > aeq(vcov(fit1), vcov(fit2))
  [1] TRUE
  > 
  > #
  > # Check that cluster, weight, and offset were correctly expanded
  > #
  > fit3a <- coxph(Surv(time1, time2, status)~ x + tt(age), tdata, weights=casewt, 
  +               tt= function(x, t, ...) ns2(x+t), x=TRUE)
  > fit3b <-  coxph(Surv(time1, time2, status) ~ x + ns2(c.age), data2,
  +                 weights=casewt)
  > aeq(coef(fit3a), coef(fit3b))
  [1] TRUE
  > aeq(vcov(fit3a), vcov(fit3b))
  [1] TRUE
  > 
  > fit4a <- coxph(Surv(time1, time2, status)~ x + tt(age), tdata,
  +               tt= function(x, t, ...) ns2(x+t), cluster=grp)
  > fit4b <-  coxph(Surv(time1, time2, status) ~ x + ns2(c.age), data2,
  +                 cluster=grp)
  > fit4c <- coxph(Surv(time1, time2, status) ~ x + ns2(c.age) + cluster(grp),
  +                data2)
  > aeq(coef(fit4a), coef(fit4b))
  [1] TRUE
  > aeq(vcov(fit4a), vcov(fit4b))
  [1] TRUE
  > aeq(coef(fit4a), coef(fit4c))
  [1] TRUE
  > aeq(vcov(fit4a), vcov(fit4c))
  [1] TRUE
  > 
  > fit5a <- coxph(Surv(time1, time2, status)~ x + tt(age) + offset(grp/10), tdata,
  +               tt= function(x, t, ...) ns2(x+t),)
  > fit5b <-  coxph(Surv(time1, time2, status) ~ x + ns2(c.age)+ offset(grp/10)
  +                 , data=data2)
  > aeq(coef(fit5a), coef(fit5b))
  [1] TRUE
  > aeq(vcov(fit5a), vcov(fit5b))
  [1] TRUE
  > 
  > # Check that strata is correct
  > fit6a <- coxph(Surv(time1, time2, status) ~ x + tt(age) + strata(strat), tdata,
  +                tt = function(x, t, ...) (x+t)^2, x=TRUE)
  > fit6b <- coxph(Surv(time1, time2, status) ~ x + I(c.age^2) +strata(strat), data2)
  > aeq(coef(fit6a), coef(fit6b))
  [1] TRUE
  > aeq(vcov(fit6a), vcov(fit6b))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.43    0.14    1.56 
Running the tests in 'tests/update.R' failed.
Complete output:
  > library(survival)
  > 
  > # the way a +cluster() term is handled in coxph has implications for update.
  > 
  > fit1 <- coxph(Surv(time, status) ~ age, cluster= inst, lung)
  > fit2 <- coxph(Surv(time, status) ~ age + cluster(inst), lung)
  > all.equal(fit1, fit2)
  [1] TRUE
  > 
  > fit3 <- coxph(Surv(time, status) ~ age + sex + cluster(inst), lung)
  > 
  > test1 <- update(fit1, .~ .+ sex)
  > all.equal(test1, fit3)
  [1] TRUE
  > 
  > # Gives a spurious warning message
  > test2 <- update(fit1, . ~ age + sex + cluster(inst), lung)
  Warning message:
  In coxph(formula = Surv(time, status) ~ age + sex + cluster(inst),  :
    cluster appears both in a formula and as an argument, formula term ignored
  > all.equal(test2, fit3)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.85    0.09    0.89 
Running the tests in 'tests/yates0.R' failed.
Complete output:
  > #
  > # This code verified some oddities about model frames
  > #
  > options(contrasts=c("contr.treatment", "contr.poly")) # clean slate
  > 
  > tdata <- data.frame(y=1:10, x1=letters[c(1,2,3,1,2,3,1,2,3,1)],
  +                             x2=LETTERS[c(1,2,3,4,1,2,3,4,1,2)],
  +                     stringsAsFactors=TRUE)
  > tdata$x3 <- as.character(tdata$x1)
  > 
  > fit1 <- lm(y ~ x1 + x2, tdata, x=TRUE)
  > m1 <- fit1$model
  > t1 <- terms(fit1)
  > 
  > # Lesson 1: xlev is ignored when the variable already has levels
  > #  as an attribute
  > 
  > temp <- list(x1=c("a", "c", "b"), x2=LETTERS[4:1])
  > x2 <- model.matrix(t1, m1, xlev=temp)
  > x3 <- model.matrix(t1, tdata, xlev=temp)
  > 
  > all.equal(x2, fit1$x)
  [1] TRUE
  > !is.logical(all.equal(x2, x3))  # x2 and x3 do not agree
  [1] TRUE
  > attributes(m1$x1)
  $levels
  [1] "a" "b" "c"
  
  $class
  [1] "factor"
  
  > 
  > 
  > # Lesson 2: character variables do not have their levels
  > # remembered as attributes in the model frame, but these are
  > #   found in fit$xlevels and fit$contrasts.
  > # However, the xlev argument is still ignored for a model frame!
  > fit2 <- lm(y ~ x3 + x2, tdata, x=TRUE)
  > m2 <- fit2$model
  > x3 <- model.matrix(terms(fit2), m2, xlev=fit2$xlevels)
  > x4 <- model.matrix(terms(fit2), m2, xlev=list(x3=letters[3:1]))
  > x5 <- model.matrix(terms(fit2), tdata, xlev=list(x3=letters[3:1]))
  > all.equal(fit2$x, x3)
  [1] TRUE
  > all.equal(x3, x4)
  [1] TRUE
  > all.equal(x3, x5)  # FALSE
  [1] "Attributes: < Component \"dimnames\": Component 2: 1 string mismatch >"
  [2] "Mean relative difference: 2.333333"                                    
  > 
  > 
  > # Lesson 3: contrasts.arg is relevant, even when the model frame
  > #  has a saved contrast
  > ctemp <- list(x1="contr.SAS", x2= contr.helmert(LETTERS[1:4]))
  > x4 <- model.matrix(t1, m1, contrasts.arg=ctemp)  # no saved contrast
  > 
  > fit3 <- lm(y ~ x1 + C(x2, contr.SAS), tdata)   
  > m3 <- fit3$model
  > attr(m3[[3]], 'contr')   # the contrast is saved
    A B C
  A 1 0 0
  B 0 1 0
  C 0 0 1
  D 0 0 0
  > c2 <- ctemp
  > names(c2) <- names(fit3$contrasts)
  > x5 <- model.matrix(terms(fit3), m3, contrasts.arg=c2)
  > all.equal(x4, x5, check.attributes=FALSE)
  [1] TRUE
  > 
  > # Lesson 4: and this holds for a character variable as well
  > fit4 <- lm(y ~ x3 + C(x2, contr.SAS), tdata, x=TRUE)   
  > c4 <- ctemp
  > names(c4) <- names(fit4$contrasts)
  > x6 <- model.matrix(terms(fit4), fit4$model, xlev=fit4$xlevels, contrasts.arg=c4)
  > x7 <- model.matrix(terms(fit4), fit4$model, contrasts.arg=c4)
  > all.equal(x7, x6)
  [1] TRUE
  > all.equal(x7, x5, check.attributes=FALSE)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     0.14    0.03    0.15 
Running the tests in 'tests/yates1.R' failed.
Complete output:
  > library(survival)
  > aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
  > 
  > fit1 <- lm(skips ~ Opening + Solder + Mask + PadType + Panel,
  +               data=solder)
  > y1 <- yates(fit1, "Opening")
  > 
  > temp <- levels(solder$Opening)
  > tpred <- matrix(0., nrow(solder), 3)
  > for (i in 1:3) {
  +     tdata <- solder
  +     tdata$Opening <- temp[i]
  +     tpred[,i] <- predict(fit1, newdata=tdata)
  +  }
  > all.equal(y1$estimate[,"pmm"], colMeans(tpred))
  [1] TRUE
  > 
  > # This fit is deficient: there are no Opening=L and Mask=A6 obs
  > # The MPV for Mask=A6 and Opening L will therefore be NA, as well
  > #   as for all levels of Solder, but we can compute the others.
  > # Solder will be NA for all levels
  > fit2 <- lm(skips ~ Opening*Mask + Solder,
  +               data=solder)
  > y2a <- yates(fit2, "Mask", population="factorial")
  > y2b <- yates(fit2, "Opening",  population="factorial")
  > y2c <- yates(fit2, "Solder",  population="factorial")
  > 
  > # The predict.lm function gives correct predictions for estimable
  > #  functions (all but L,A6) and nonsense for others.  It knows that
  > #  some are not estimable due to the NA coefficients, but not which ones,
  > #  so always prints a warning.  Hence the suppressWarnings call.
  > tdata <- do.call(expand.grid, fit2$xlevels[1:3])
  > temp <- levels(solder$Mask)
  > tpreda <- matrix(0., nrow(tdata), length(temp),
  +                  dimnames=list(NULL, temp))
  > for (i in seq(along=temp)) {
  +     tdata$Mask <- temp[i]
  +     suppressWarnings(tpreda[,i] <- predict(fit2, newdata=tdata))
  +  }
  > tpreda[,"A6"] <- NA  # the A6 estimate is deficient
  > aeq(y2a$estimate[,"pmm"], colMeans(tpreda))
  [1] TRUE
  > 
  > tdata <- do.call(expand.grid, fit2$xlevels[1:3])
  > temp <- levels(solder$Opening)
  > tpredb <- matrix(0., nrow(tdata), length(temp),
  +                  dimnames=list(NULL, temp))
  > for (i in seq(along=temp)) {
  +     tdata$Opening <- temp[i]
  +     suppressWarnings(tpredb[,i] <- predict(fit2, newdata=tdata))
  +  }
  > tpredb[,"L"] <- NA  
  > aeq(y2b$estimate[,"pmm"], colMeans(tpredb))
  [1] TRUE
  > 
  > # Solder should be all NA
  > all(is.na(y2c$estimate[,"pmm"]))
  [1] TRUE
  > 
  > # Tests for Solder are defined for a non-factorial population, however.
  > # the [] below retains the factor structure of the variable, where the
  > #  runs above did not.  R gets prediction correct  both ways.
  > y2d <- yates(fit2, ~Solder)
  > temp <- levels(solder$Solder)
  > tdata <- solder
  > tpredd <- matrix(0, nrow(tdata), length(temp),
  +                  dimnames=list(NULL, temp))
  > for (i in seq(along=temp)) {
  +     tdata$Solder[] <- temp[i]
  +     suppressWarnings(tpredd[,i] <- predict(fit2, newdata=tdata))
  + }
  > aeq(y2d$estimate$pmm, colMeans(tpredd))
  [1] TRUE
  > 
  > #
  > # Verify that the result is unchanged by how dummies are coded
  > #   The coefs move all over the map, but predictions are unchanged
  > fit3 <- lm(skips ~ C(Opening, contr.helmert)*Mask + C(Solder, contr.SAS),
  +            data=solder)
  > y3a <- yates(fit3, ~Mask, population='yates')
  > equal <- c("estimate", "test", "mvar")
  > all.equal(y3a[equal], y2a[equal])
  [1] TRUE
  > 
  > tdata <- do.call(expand.grid, fit2$xlevels[1:3]) # use orignal variable names
  > temp <- levels(solder$Mask)
  > cpred <- matrix(0., nrow(tdata), length(temp),
  +                  dimnames=list(NULL, temp))
  > for (i in seq(along=temp)) {
  +     tdata$Mask <- temp[i]
  +     suppressWarnings(cpred[,i] <- predict(fit3, newdata=tdata))
  +  }
  > aeq(cpred[, temp!="A6"], tpreda[, temp!= "A6"])  # same predictions
  [1] TRUE
  > all.equal(y3a$estimate, y2a$estimate)
  [1] TRUE
  > 
  > y3b <- yates(fit3, ~Opening, population='yates')
  > # column names will differ
  > all.equal(y3b$estimate, y2b$estimate, check.attributes=FALSE)
  [1] TRUE
  > 
  > y3d <- yates(fit3, ~Solder)
  > for (i in 1:3) {
  +     print(all.equal(y3d[[i]], y2d[[i]], check.attributes=FALSE))
  + }
  [1] TRUE
  [1] TRUE
  [1] TRUE
  > 
  > # Reprise this with a character variable in the model
  > sdata <- solder
  > sdata$Mask <- as.character(sdata$Mask)
  > fit4 <-  lm(skips ~ Opening*Mask + Solder, data=sdata)
  > y4a <- yates(fit4, ~ Mask, population= "yates")
  > y4b <- yates(fit4, ~ Opening, population= "yates")
  > y4d <- yates(fit4, ~ Solder)
  > equal <- c("estimate", "tests", "mvar", "cmat")
  > all.equal(y2a[equal], y4a[equal])  # the "call" component differs
  [1] TRUE
  > all.equal(y2b[equal], y4b[equal])
  [1] TRUE
  > all.equal(y2d[equal], y4d[equal])
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     0.93    0.14    1.04 
Running the tests in 'tests/yates2.R' failed.
Complete output:
  > library(survival)
  > # Tests for glm
  > 
  > gfit1 <- glm(skips ~ Mask* Opening + Solder, data=solder, poisson)
  > yg1 <- yates(gfit1, ~Mask)
  > yg2 <- yates(gfit1, ~ Mask, predict='response')
  > 
  > proc.time()
     user  system elapsed 
     0.85    0.15    1.01 
Running the tests in 'tests/zph.R' failed.
Complete output:
  > library(survival)
  > aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
  > 
  > test1 <- data.frame(time=  c(9, 3,1,1,6,6,8),
  +                     status=c(1,NA,1,0,1,1,0),
  +                     x=     c(0, 2,1,1,1,0,0))
  > 
  > # Verify that cox.zph computes a score test
  > # First for the Breslow estimate
  > r <- (3 + sqrt(33))/2   # actual MLE for log(beta)
  > U <- c(1/(r+1), 3/(r+3), -r/(r+3), 0)   # score statistic
  > imat <- c(r/(r+1)^2, 3*r/(r+3)^2, 3*r/(r+3)^2, 0)  # information matrix
  > g = c(1, 6, 6, 9)  # death times
  > 
  > u2 <- c(sum(U), sum(g*U))  # first derivative
  > i2 <- matrix(c(sum(imat), sum(g*imat), sum(g*imat), sum(g^2*imat)),
  +                2,2)  # second derivative
  > sctest <- solve(i2, u2) %*% u2
  > 
  > # Verify that centering makes no difference for the test (though i2 changes)
  > g2 <- g - mean(g)
  > u2b <- c(sum(U), sum(g2*U))
  > i2b <- matrix(c(sum(imat), sum(g2*imat), sum(g2*imat), sum(g2^2*imat)),
  +                2,2) 
  > sctest2 <- solve(i2b, u2b) %*% u2b
  > all.equal(sctest, sctest2)
  [1] TRUE
  > 
  > # Now check the program
  > fit1 <- coxph(Surv(time, status) ~ x, test1, ties='breslow')
  > aeq(fit1$coef, log(r))
  [1] TRUE
  > zp1 <- cox.zph(fit1, transform='identity', global=FALSE)
  > aeq(zp1$table[,1], sctest)
  [1] TRUE
  > aeq(zp1$y, resid(fit1, type="scaledsch"))
  [1] TRUE
  > 
  > dummy <- rep(0, nrow(test1))
  > fit1b <- coxph(Surv(dummy, time, status) ~ x, test1, ties='breslow')
  > aeq(fit1b$coef, log(r))
  [1] TRUE
  > zp1b <- cox.zph(fit1b, transform='identity', global=FALSE)
  > aeq(zp1b$table[,1], sctest)
  [1] TRUE
  > # the pair of tied times gets reversed in the zph result
  > #  but since the 'y' values are only used to plot it doesn't matter
  > aeq(zp1b$y[c(1,3,2,4)], resid(fit1b, type="scaledsch"))  
  [1] TRUE
  > 
  > # log time check
  > g3 <- log(g) - mean(log(g))
  > u3 <- c(sum(U), sum(g3*U))  # first derivative
  > i3 <- matrix(c(sum(imat), sum(g3*imat), sum(g3*imat), sum(g3^2*imat)),
  +                2,2)  # second derivative
  > sctest3 <- solve(i3, u3) %*% u3
  > zp3 <- cox.zph(fit1, transform='log', global=FALSE)
  > aeq(zp3$table[,1], sctest3)
  [1] TRUE
  > 
  > # Efron approximation
  > phi <- acos((45/23)*sqrt(3/23))
  > r <- 2*sqrt(23/3)* cos(phi/3)   # actual MLE for log(beta)
  > U <- c(1/(r+1), 3/(r+3), -r/(r+5), 0)   # score statistic
  > imat <- c(r/(r+1)^2, 3*r/(r+3)^2, 5*r/(r+5)^2, 0)  # information matrix
  > 
  > u4 <- c(sum(U), sum(g3*U))  # first derivative
  > i4 <- matrix(c(sum(imat), sum(g3*imat), sum(g3*imat), sum(g3^2*imat)),
  +                2,2)  # second derivative
  > sctest4 <- solve(i4, u4) %*% u4
  > 
  > fit4 <- coxph(Surv(time, status) ~ x, test1, ties='efron')
  > aeq(fit4$coef, log(r))
  [1] TRUE
  > zp4 <- cox.zph(fit4, transform='log', global=FALSE)
  > aeq(zp4$table[,1], sctest4)
  [1] TRUE
  > aeq(zp4$y, resid(fit4, type="scaledsch"))
  [1] TRUE
  > 
  > fit5 <-  coxph(Surv(dummy, time, status) ~ x, test1, ties="efron")
  > aeq(fit5$coef, log(r))
  [1] TRUE
  > zp5 <- cox.zph(fit5, transform="log", global=FALSE)
  > aeq(zp5$table[,1], sctest4)
  [1] TRUE
  > 
  > # Artificial stratification
  > test2 <- rbind(test1, test1)
  > test2$group <- rep(letters[1:2], each=nrow(test1))
  > # U, imat, and sctest will all double
  > dummy <- c(dummy, dummy)
  > fit6 <- coxph(Surv(dummy, time, status) ~ x + strata(group), test2)
  > aeq(fit6$coef, log(r))
  [1] TRUE
  > zp6 <- cox.zph(fit6, transform="log", globa=FALSE)
  > aeq(zp6$table[,1], 2*sctest4)
  [1] TRUE
  > 
  > # A multi-state check, 2 covariates
  > #  Verify that the multi-state result = the single state Cox models
  > etime <- with(mgus2, ifelse(pstat==0, futime, ptime))
  > event <- with(mgus2, ifelse(pstat==0, 2*death, 1))
  > event <- factor(event, 0:2, labels=c("censor", "pcm", "death"))
  > table(event)
  event
  censor    pcm  death 
     409    115    860 
  > 
  > ct1 <- coxph(Surv(etime, event) ~ sex + age, mgus2, id=id)
  > ct2 <- coxph(Surv(etime, event=='pcm') ~ sex + age, mgus2)
  > ct3 <- coxph(Surv(etime, event=='death') ~ sex + age, mgus2)
  > 
  > zp1 <- cox.zph(ct1, transform='identity')
  > zp2 <- cox.zph(ct2, transform='identity')
  > zp3 <- cox.zph(ct3, transform='identity')
  > aeq(zp1$table[1:2,], zp2$table[1:2,])
  [1] TRUE
  > aeq(zp1$table[3:4,], zp3$table[1:2,])
  [1] TRUE
  > 
  > # Now add a starting time of zero
  > dummy <- rep(0, nrow(mgus2))
  > ct4 <- coxph(Surv(dummy, etime, event) ~ sex + age, mgus2, id=id)
  > ct5 <- coxph(Surv(dummy, etime, event=='pcm') ~ sex + age, mgus2)
  > ct6 <- coxph(Surv(dummy, etime, event=='death') ~ sex + age, mgus2)
  > zp4 <- cox.zph(ct4, transform='identity')
  > zp5 <- cox.zph(ct5, transform='identity')
  > zp6 <- cox.zph(ct6, transform='identity')
  > aeq(zp4$table[1:2,], zp5$table[1:2,])
  [1] TRUE
  > aeq(zp4$table[3:4,], zp6$table[1:2,])
  [1] TRUE
  > 
  > 
  > # Direct check of a multivariate model with start, stop data
  > p1 <- pbcseq[!duplicated(pbcseq$id), 1:6]
  > pdata <- tmerge(p1[, c("id", "trt", "age", "sex")], p1, id=id,
  +                 death = event(futime, status==2))
  > pdata <- tmerge(pdata, pbcseq, id=id, bili=tdc(day, bili),
  +                 edema = tdc(day, edema), albumin=tdc(day, albumin),
  +                 protime = tdc(day, protime))
  > pfit <- coxph(Surv(tstart, tstop, death) ~ log(bili) + albumin + edema +
  +                   age + log(protime), data = pdata, ties='efron')
  > zp7  <- cox.zph(pfit, transform='log', global=FALSE)
  > 
  > direct <- function(fit) {
  +     nvar <- length(fit$coef)
  +     dt <- coxph.detail(fit)
  +     gtime <- log(dt$time) - mean(log(dt$time))
  +     # key idea: at any event time I have a first deriviative vector
  +     #     c(dt$score[i,],  gtime[i]* dt$score[i,]) 
  +     # and second derivative matrix
  +     #     dt$imat[,,i]              gtime[i]   * dt$imat[,,i]
  +     #     gtime[i]*dt$imat[,,i]     gtime[i]^2 * dt$imat[,,i]
  +     # for the expanded model, where imat[,,i] is symmetric, 
  +     # and colSums(dt$score) =0  (since the model converged)
  +     # 
  +     # Create score tests for adding one time-dependent variable 
  +     #  gtime * x[,j] at a time: first derivative of this test is
  +     #     c(dt$score[i,],  gtime[i]* dt$score[i,j]) 
  +     #  and etc.
  +     unew  <-   colSums(gtime * dt$score)
  +     temp1 <- apply(dt$imat, 1:2, sum)
  +     temp2 <- apply(dt$imat, 1:2, function(x) sum(x*gtime))
  +     temp3 <- apply(dt$imat, 1:2, function(x) sum(x * gtime^2))
  + 
  +     score <- double(nvar)
  +     smat  <- matrix(0., nvar+1, nvar+1)  # second deriv matrix for the test
  +     smat[1:nvar, 1:nvar] <- temp1
  +     for (i in 1:nvar) {
  +         smat[nvar+1,] <- c(temp2[i,], temp3[i,i])
  +         smat[,nvar+1] <- c(temp2[,i], temp3[i,i])
  +         utemp <- c(rep(0,nvar), unew[i])
  +         score[i] <- solve(smat, utemp) %*% utemp
  +     }  
  +     list(sctest = score, u= c(colSums(dt$score), unew), 
  +          imat=cbind(rbind(temp1, temp2), rbind(temp2, temp3)))
  + }   
  > 
  > aeq(zp7$table[,1], direct(pfit)$sctest)
  [1] TRUE
  > 
  > # Last, make sure that NA coefficients are ignored
  > d1 <- survSplit(Surv(time, status) ~ ., veteran, cut=150, episode="epoch")
  > fit <- coxph(Surv(tstart, time, status) ~ celltype:strata(epoch) + age, d1)
  > zz <- cox.zph(fit)
  > 
  > fit2 <- coxph(Surv(tstart, time, status) ~ celltype:strata(epoch) + age, d1,
  +               x=TRUE)
  > zz2 <- cox.zph(fit2)
  > 
  > x2 <- fit2$x[, !is.na(fit$coefficients)][,-1]
  > fit3 <- coxph(Surv(tstart, time, status) ~ age + x2, d1)
  > all.equal(fit3$loglik, fit2$loglik)
  [1] TRUE
  > zz3 <- cox.zph(fit3)
  > 
  > all.equal(unclass(zz)[1:7], unclass(zz2)[1:7]) #ignore the call component
  [1] TRUE
  > all.equal(as.vector(zz$table), as.vector(zz3$table)) # variable names change
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.54    0.15    1.68