Last updated on 2026-10-24 12:51:49 CEST.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 3.8-12 | 59.70 | 585.94 | 645.64 | OK | |
| r-devel-linux-x86_64-debian-gcc | 3.8-12 | 44.90 | 340.38 | 385.28 | ERROR | |
| r-devel-linux-x86_64-fedora-clang | 3.8-12 | 43.00 | 519.70 | 562.70 | OK | |
| r-devel-linux-x86_64-fedora-gcc | 3.8-12 | 43.00 | 414.06 | 457.06 | OK | |
| r-devel-windows-x86_64 | 3.8-12 | 89.00 | 858.00 | 947.00 | NOTE | |
| r-patched-linux-x86_64 | 3.8-12 | 63.75 | 558.41 | 622.16 | OK | |
| r-release-linux-x86_64 | 3.8-12 | 56.35 | 560.39 | 616.74 | OK | |
| r-release-macos-arm64 | 3.8-12 | 15.00 | 154.00 | 169.00 | OK | |
| r-release-macos-x86_64 | 3.8-12 | 48.00 | 674.00 | 722.00 | OK | |
| r-release-windows-x86_64 | 3.8-12 | 93.00 | 746.00 | 839.00 | OK | |
| r-oldrel-macos-arm64 | 3.8-12 | 14.00 | 155.00 | 169.00 | OK | |
| r-oldrel-macos-x86_64 | 3.8-12 | 46.00 | 784.00 | 830.00 | OK | |
| r-oldrel-windows-x86_64 | 3.8-12 | 112.00 | 790.00 | 902.00 | OK |
Version: 3.8-12
Check: examples
Result: ERROR
Running examples in ‘survival-Ex.R’ failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: aggregate.survfit
> ### Title: Average survival curves
> ### Aliases: aggregate.survfit
> ### Keywords: survival
>
> ### ** Examples
>
> cfit <- coxph(Surv(futime, death) ~ sex + age*hgb, data=mgus2)
> # marginal effect of sex, after adjusting for the others
> dummy <- rbind(mgus2, mgus2)
> dummy$sex <- rep(c("F", "M"), each=nrow(mgus2)) # population data set
> dummy <- na.omit(dummy) # don't count missing hgb in our "population
> csurv <- survfit(cfit, newdata=dummy)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 3.8-12
Check: tests
Result: ERROR
Running ‘aareg.R’ [2s/2s]
Running ‘anova.R’ [2s/2s]
Running ‘badname.R’ [1s/2s]
Running ‘bladder.R’ [1s/2s]
Comparing ‘bladder.Rout’ to ‘bladder.Rout.save’ ... OK
Running ‘book1.R’ [1s/2s]
Running ‘book2.R’ [1s/2s]
Running ‘book3.R’ [1s/2s]
Running ‘book4.R’ [1s/2s]
Running ‘book5.R’ [1s/2s]
Running ‘book6.R’ [1s/2s]
Running ‘book7.R’ [2s/2s]
Comparing ‘book7.Rout’ to ‘book7.Rout.save’ ... OK
Running ‘brier.R’ [1s/2s]
Running ‘cancer.R’ [2s/2s]
Running ‘checkSurv2.R’ [2s/2s]
Comparing ‘checkSurv2.Rout’ to ‘checkSurv2.Rout.save’ ... OK
Running ‘clogit.R’ [1s/2s]
Comparing ‘clogit.Rout’ to ‘clogit.Rout.save’ ... OK
Running ‘concordance.R’ [2s/2s]
Comparing ‘concordance.Rout’ to ‘concordance.Rout.save’ ... OK
Running ‘concordance2.R’ [2s/2s]
Running ‘concordance3.R’ [2s/3s]
Comparing ‘concordance3.Rout’ to ‘concordance3.Rout.save’ ... OK
Running ‘counting.R’ [1s/2s]
Running ‘coxsurv.R’ [1s/2s]
Running ‘coxsurv2.R’ [1s/2s]
Running ‘coxsurv3.R’ [1s/2s]
Running ‘coxsurv4.R’ [1s/2s]
Running ‘coxsurv5.R’ [2s/2s]
Running ‘coxsurv6.R’ [1s/2s]
Running ‘detail.R’ [1s/1s]
Running ‘difftest.R’ [1s/1s]
Comparing ‘difftest.Rout’ to ‘difftest.Rout.save’ ... OK
Running ‘doaml.R’ [1s/2s]
Running ‘doublecolon.R’ [1s/2s]
Comparing ‘doublecolon.Rout’ to ‘doublecolon.Rout.save’ ... OK
Running ‘doweight.R’ [1s/2s]
Running ‘dropspecial.R’ [1s/2s]
Comparing ‘dropspecial.Rout’ to ‘dropspecial.Rout.save’ ... OK
Running ‘ekm.R’ [1s/2s]
Comparing ‘ekm.Rout’ to ‘ekm.Rout.save’ ... OK
Running ‘expected.R’ [2s/2s]
Running ‘expected2.R’ [1s/2s]
Running ‘factor.R’ [1s/2s]
Running ‘factor2.R’ [1s/2s]
Running ‘finegray.R’ [2s/2s]
Comparing ‘finegray.Rout’ to ‘finegray.Rout.save’ ... OK
Running ‘fr_cancer.R’ [2s/2s]
Comparing ‘fr_cancer.Rout’ to ‘fr_cancer.Rout.save’ ... OK
Running ‘fr_kidney.R’ [2s/2s]
Comparing ‘fr_kidney.Rout’ to ‘fr_kidney.Rout.save’ ... OK
Running ‘fr_lung.R’ [1s/2s]
Comparing ‘fr_lung.Rout’ to ‘fr_lung.Rout.save’ ... OK
Running ‘fr_ovarian.R’ [1s/2s]
Comparing ‘fr_ovarian.Rout’ to ‘fr_ovarian.Rout.save’ ... OK
Running ‘fr_rat1.R’ [1s/2s]
Comparing ‘fr_rat1.Rout’ to ‘fr_rat1.Rout.save’ ... OK
Running ‘fr_resid.R’ [1s/1s]
Running ‘fr_simple.R’ [1s/2s]
Comparing ‘fr_simple.Rout’ to ‘fr_simple.Rout.save’ ... OK
Running ‘frailty.R’ [1s/1s]
Running ‘frank.R’ [1s/2s]
Running ‘infcox.R’ [1s/2s]
Comparing ‘infcox.Rout’ to ‘infcox.Rout.save’ ... OK
Running ‘jasa.R’ [1s/2s]
Running ‘model.matrix.R’ [1s/1s]
Running ‘mstate.R’ [2s/2s]
Comparing ‘mstate.Rout’ to ‘mstate.Rout.save’ ... OK
Running ‘mstate2.R’ [2s/2s]
Comparing ‘mstate2.Rout’ to ‘mstate2.Rout.save’ ... OK
Running ‘mstrata.R’ [1s/1s]
Running ‘multi2.R’ [2s/2s]
Running ‘multi3.R’ [2s/2s]
Comparing ‘multi3.Rout’ to ‘multi3.Rout.save’ ... OK
Running ‘multistate.R’ [2s/2s]
Comparing ‘multistate.Rout’ to ‘multistate.Rout.save’ ... OK
Running ‘multistrata.R’ [2s/2s]
Running ‘neardate.R’ [1s/2s]
Comparing ‘neardate.Rout’ to ‘neardate.Rout.save’ ... OK
Running ‘nested.R’ [1s/2s]
Running ‘ovarian.R’ [2s/2s]
Running ‘overlap.R’ [1s/2s]
Comparing ‘overlap.Rout’ to ‘overlap.Rout.save’ ... OK
Running ‘prednew.R’ [1s/2s]
Running ‘predsurv.R’ [1s/2s]
Running ‘pseudo.R’ [2s/2s]
Comparing ‘pseudo.Rout’ to ‘pseudo.Rout.save’ ... OK
Running ‘pspline.R’ [1s/1s]
Running ‘pyear.R’ [1s/2s]
Comparing ‘pyear.Rout’ to ‘pyear.Rout.save’ ... OK
Running ‘quantile.R’ [1s/2s]
Running ‘r_lung.R’ [1s/2s]
Comparing ‘r_lung.Rout’ to ‘r_lung.Rout.save’ ... OK
Running ‘r_resid.R’ [1s/2s]
Comparing ‘r_resid.Rout’ to ‘r_resid.Rout.save’ ... OK
Running ‘r_sas.R’ [5s/5s]
Comparing ‘r_sas.Rout’ to ‘r_sas.Rout.save’ ... OK
Running ‘r_scale.R’ [1s/2s]
Comparing ‘r_scale.Rout’ to ‘r_scale.Rout.save’ ... OK
Running ‘r_stanford.R’ [2s/2s]
Comparing ‘r_stanford.Rout’ to ‘r_stanford.Rout.save’ ... OK
Running ‘r_strata.R’ [1s/2s]
Comparing ‘r_strata.Rout’ to ‘r_strata.Rout.save’ ... OK
Running ‘r_tdist.R’ [1s/1s]
Comparing ‘r_tdist.Rout’ to ‘r_tdist.Rout.save’ ... OK
Running ‘r_user.R’ [1s/2s]
Comparing ‘r_user.Rout’ to ‘r_user.Rout.save’ ... OK
Running ‘ratetable.R’ [1s/1s]
Comparing ‘ratetable.Rout’ to ‘ratetable.Rout.save’ ... OK
Running ‘residms.R’ [1s/2s]
Running ‘residsf.R’ [1s/2s]
Comparing ‘residsf.Rout’ to ‘residsf.Rout.save’ ... OK
Running ‘royston.R’ [1s/2s]
Running ‘rttright.R’ [2s/2s]
Comparing ‘rttright.Rout’ to ‘rttright.Rout.save’ ... OK
Running ‘singtest.R’ [1s/1s]
Comparing ‘singtest.Rout’ to ‘singtest.Rout.save’ ... OK
Running ‘strata2.R’ [1s/2s]
Running ‘stratatest.R’ [1s/2s]
Running ‘summary_survfit.R’ [1s/2s]
Comparing ‘summary_survfit.Rout’ to ‘summary_survfit.Rout.save’ ... OK
Running ‘summarydf.R’ [1s/2s]
Running ‘surv.R’ [1s/2s]
Comparing ‘surv.Rout’ to ‘surv.Rout.save’ ... OK
Running ‘survSplit.R’ [1s/2s]
Comparing ‘survSplit.Rout’ to ‘survSplit.Rout.save’ ... OK
Running ‘survcheck.R’ [2s/2s]
Comparing ‘survcheck.Rout’ to ‘survcheck.Rout.save’ ... OK
Running ‘survexpm.R’ [1s/2s]
Comparing ‘survexpm.Rout’ to ‘survexpm.Rout.save’ ... OK
Running ‘survfit1.R’ [2s/2s]
Comparing ‘survfit1.Rout’ to ‘survfit1.Rout.save’ ... OK
Running ‘survfit2.R’ [1s/1s]
Comparing ‘survfit2.Rout’ to ‘survfit2.Rout.save’ ... OK
Running ‘survfit3.R’ [1s/2s]
Comparing ‘survfit3.Rout’ to ‘survfit3.Rout.save’ ... OK
Running ‘survreg1.R’ [1s/2s]
Comparing ‘survreg1.Rout’ to ‘survreg1.Rout.save’ ... OK
Running ‘survreg2.R’ [2s/3s]
Comparing ‘survreg2.Rout’ to ‘survreg2.Rout.save’ ... OK
Running ‘survtest.R’ [1s/2s]
Running ‘testci.R’ [2s/2s]
Comparing ‘testci.Rout’ to ‘testci.Rout.save’ ... OK
Running ‘testci2.R’ [1s/2s]
Comparing ‘testci2.Rout’ to ‘testci2.Rout.save’ ... OK
Running ‘testnull.R’ [1s/2s]
Comparing ‘testnull.Rout’ to ‘testnull.Rout.save’ ... OK
Running ‘testreg.R’ [1s/2s]
Comparing ‘testreg.Rout’ to ‘testreg.Rout.save’ ... OK
Running ‘tiedtime.R’ [2s/2s]
Comparing ‘tiedtime.Rout’ to ‘tiedtime.Rout.save’ ... OK
Running ‘timeline.R’ [2s/2s]
Running ‘tmerge.R’ [1s/2s]
Comparing ‘tmerge.Rout’ to ‘tmerge.Rout.save’ ... OK
Running ‘tmerge2.R’ [1s/2s]
Comparing ‘tmerge2.Rout’ to ‘tmerge2.Rout.save’ ... OK
Running ‘tmerge3.R’ [1s/2s]
Comparing ‘tmerge3.Rout’ to ‘tmerge3.Rout.save’ ... OK
Running ‘tmerge4.R’ [1s/1s]
Comparing ‘tmerge4.Rout’ to ‘tmerge4.Rout.save’ ... OK
Running ‘tt.R’ [2s/2s]
Comparing ‘tt.Rout’ to ‘tt.Rout.save’ ... OK
Running ‘tt2.R’ [2s/2s]
Comparing ‘tt2.Rout’ to ‘tt2.Rout.save’ ... OK
Running ‘turnbull.R’ [2s/2s]
Comparing ‘turnbull.Rout’ to ‘turnbull.Rout.save’ ... OK
Running ‘update.R’ [1s/2s]
Comparing ‘update.Rout’ to ‘update.Rout.save’ ... OK
Running ‘yates.R’ [0s/0s]
Running ‘yates0.R’ [0s/0s]
Comparing ‘yates0.Rout’ to ‘yates0.Rout.save’ ... OK
Running ‘yates1.R’ [2s/2s]
Comparing ‘yates1.Rout’ to ‘yates1.Rout.save’ ... OK
Running ‘yates2.R’ [1s/2s]
Running ‘zph.R’ [1s/1s]
Running the tests in ‘tests/aareg.R’ failed.
Complete output:
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
> library(survival)
>
> #
> # Test aareg, for some simple data where the answers can be computed
> # in closed form
> #
> aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> 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),
+ wt= c(1, 1:6))
>
> tfit <- aareg(Surv(time, status) ~ x, test1)
> aeq(tfit$times, c(1,2,2))
[1] TRUE
> aeq(tfit$nrisk, c(6,4,4))
[1] TRUE
> aeq(tfit$coefficient, matrix(c(0,0,1/3, 1/3, 1, -1/3), ncol=2))
[1] TRUE
> aeq(tfit$tweight, matrix(c(3,3,3, 3/2, 3/4, 3/4), ncol=2))
[1] TRUE
> aeq(tfit$test.statistic, c(1,1))
[1] TRUE
> aeq(tfit$test.var, c(1, -1/4, -1/4, 1/4 + 9/16 + 1/16))
[1] TRUE
>
> tfit <- aareg(Surv(time, status) ~ x, test1, test='nrisk')
> aeq(tfit$tweight, matrix(c(3,3,3, 3/2, 3/4, 3/4), ncol=2)) #should be as before
[1] TRUE
> aeq(tfit$test.statistic, c(4/3, 6/3+ 4 - 4/3))
[1] TRUE
> aeq(tfit$test.var, c(16/9, -16/9, -16/9, 36/9 + 16 + 16/9))
[1] TRUE
>
> # In the 1-variable case, this is the same as the default Aalen weight
> tfit <- aareg(Surv(time, status) ~ x, test1, test='variance')
> aeq(tfit$test.statistic, c(1,1))
[1] TRUE
> aeq(tfit$test.var, c(1, -1/4, -1/4, 1/4 + 9/16 + 1/16))
[1] TRUE
>
> #
> # Repeat the above, with case weights
> #
> tfit <- aareg(Surv(time, status) ~x, test1, weights=wt)
> aeq(tfit$times, c(1,2,2))
[1] TRUE
> aeq(tfit$nrisk, c(21,16,16))
[1] TRUE
> aeq(tfit$coefficient, matrix(c(0,0,5/12, 2/9, 1, -5/12), ncol=2))
[1] TRUE
> aeq(tfit$tweight, matrix(c(12,12,12, 36/7, 3,3), ncol=2))
[1] TRUE
> aeq(tfit$test.statistic, c(5, 72/63 + 3 - 15/12))
[1] TRUE
> aeq(tfit$test.var, c(25, -25/4, -25/4, (72/63)^2 + 9 + (5/4)^2))
[1] TRUE
>
> tfit <- aareg(Surv(time, status) ~x, test1, weights=wt, test='nrisk')
> aeq(tfit$test.statistic, c(20/3, 42/9 + 16 - 16*5/12))
[1] TRUE
> aeq(tfit$test.var, c(400/9, -400/9, -400/9,
+ (42/9)^2 + 16^2 + (16*5/12)^2))
[1] TRUE
>
> #
> # Make a test data set with no NAs, in sorted order, no ties,
> # 15 observations
> tdata <- lung[15:29, c('time', 'status', 'age', 'sex', 'ph.ecog')]
> tdata$status <- tdata$status -1
> tdata <- tdata[order(tdata$time, tdata$status),]
> row.names(tdata) <- 1:15
> tdata$status[8] <- 0 #for some variety
>
> afit <- aareg(Surv(time, status) ~ age + sex + ph.ecog, tdata, nmin=6)
> #
> # Now, do it "by hand"
> cfit <- coxph(Surv(time, status) ~ age + sex + ph.ecog, tdata, iter=0,
+ method='breslow')
> dt1 <- coxph.detail(cfit)
Error in a$specials[[i]] : subscript out of bounds
Calls: coxph.detail ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: anova ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/badname.R’ failed.
Complete output:
> # check for acceptance of non-syntatic variable names
> # survfit was the one with issues
> library(survival)
> aeq <- function(x,y,...) all.equal(as.vector(x), as.vector(y))
>
> lung2 <- lung
> lung2["bad sex"] <- lung2$sex
>
> fit1a <- coxph(Surv(time, status) ~ ph.ecog + sex, lung2)
> fit1b <- coxph(Surv(time, status) ~ ph.ecog + `bad sex`, lung2)
> aeq(coef(fit1a), coef(fit1b))
[1] TRUE
> aeq(fit1a$loglik, fit1b$loglik)
[1] TRUE
>
> dummy <- expand.grid(ph.ecog=0:2, sex=1:2)
> dummy["bad sex"] <- dummy$sex
> surv1a <- survfit(fit1a, newdata= dummy)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
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$residuals[c(2:6,1)])
[1] TRUE
> aeq(truth0$scho, resid(fit0, 'schoen'))
Error in a$specials[[i]] : subscript out of bounds
Calls: aeq ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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$residuals[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'))
Error in a$specials[[i]] : subscript out of bounds
Calls: aeq ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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$residuals)
[1] TRUE
> aeq(truth0$scho, resid(fit0, 'schoen'))
Error in a$specials[[i]] : subscript out of bounds
Calls: aeq ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/book4.R’ failed.
Complete output:
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
> library(survival)
>
> #
> # Tests from the appendix of Therneau and Grambsch
> # d. Data set 2 and Efron 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) + 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,
+ mart=mart, score=score, rmat=resid,
+ scho=scho)
+ }
>
>
> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
>
> fit0 <-coxph(Surv(start, stop, event) ~x, test2, 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$residuals)
[1] TRUE
> aeq(truth0$scho, resid(fit0, 'schoen'))
Error in a$specials[[i]] : subscript out of bounds
Calls: aeq ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/book5.R’ failed.
Complete output:
> library(survival)
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
>
> # Tests of the weighted Cox model
> # This is section 1.3 of my appendix -- not yet found in the book
> # though, it awaits the next edition
> #
> # Similar data set to test1, but add weights,
> # a double-death/censor tied time
> # a censored last subject
> # The latter two are cases covered only feebly elsewhere.
> #
> # The data set testw2 has the same data, but done via replication
> #
> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
>
> testw1 <- data.frame(time= c(1,1,2,2,2,2,3,4,5),
+ status= c(1,0,1,1,1,0,0,1,0),
+ x= c(2,0,1,1,0,1,0,1,0),
+ wt = c(1,2,3,4,3,2,1,2,1),
+ id = 1:9)
> # Expanded data set
> testw2 <- testw1[rep(1:9, testw1$wt), -4]
> row.names(testw2) <- NULL
> indx <- match(1:9, testw2$id)
>
> # Breslow estimate
> byhand <- function(beta, newx=0) {
+ r <- exp(beta)
+ loglik <- 11*beta - (log(r^2 + 11*r +7) + 10*log(11*r +5) +2*log(2*r+1))
+ hazard <- c(1/(r^2 + 11*r +7), 10/(11*r +5), 2/(2*r+1))
+ xbar <- c((2*r^2 + 11*r)*hazard[1], 11*r/(11*r +5), r*hazard[3])
+ U <- 11- (xbar[1] + 10*xbar[2] + 2*xbar[3])
+ imat <- (4*r^2 + 11*r)*hazard[1] - xbar[1]^2 +
+ 10*(xbar[2] - xbar[2]^2) + 2*(xbar[3] - xbar[3]^2)
+
+ temp <- cumsum(hazard)
+ risk <- c(r^2, 1,r,r,1,r,1,r,1)
+ expected <- risk* temp[c(1,1,2,2,2,2,2,3,3)]
+
+ # The matrix of weights, one row per obs, one col per death
+ # deaths at 1,2,2,2, and 4
+ riskmat <- matrix(c(1,1,1,1,1,1,1,1,1,
+ 0,0,1,1,1,1,1,1,1,
+ 0,0,1,1,1,1,1,1,1,
+ 0,0,1,1,1,1,1,1,1,
+ 0,0,0,0,0,0,0,1,1), ncol=5)
+ wtmat <- diag(c(r^2, 2, 3*r, 4*r, 3, 2*r, 1, 2*r, 1)) %*% riskmat
+
+ x <- c(2,0,1,1,0,1,0,1,0)
+ status <- c(1,0,1,1,1,0,0,1,0)
+ wt <- c(1,2,3,4,3,2,1,2,1)
+ # Table of sums for score and Schoenfeld resids
+ hazmat <- riskmat %*% diag(c(1,3,4,3,2)/colSums(wtmat))
+ dM <- -risk*hazmat #Expected part
+ dM[1,1] <- dM[1,1] +1 # deaths at time 1
+ for (i in 2:4) dM[i+1, i] <- dM[i+1,i] +1
+ dM[8,5] <- dM[8,5] +1
+ mart <- rowSums(dM)
+ resid <-dM * outer(x, xbar[c(1,2,2,2,3)] ,'-')
+
+ # Increments to the variance of the hazard
+ var.g <- cumsum(hazard^2/ c(1,10,2))
+ var.d <- cumsum((xbar-newx)*hazard)
+
+ list(loglik=loglik, U=U, imat=imat, hazard=hazard, xbar=xbar,
+ mart=c(1,0,1,1,1,0,0,1,0)-expected, expected=expected,
+ score=rowSums(resid), schoen=c(2,1,1,0,1) - xbar[c(1,2,2,2,3)],
+ varhaz=(var.g + var.d^2/imat)* exp(2*beta*newx))
+ }
>
> aeq(byhand(0)$expected, c(1/19, 1/19, rep(103/152, 5), rep(613/456,2))) #verify
[1] TRUE
>
> fit0 <- coxph(Surv(time, status) ~x, testw1, weights=wt,
+ method='breslow', iter=0)
> fit0b <- coxph(Surv(time, status) ~x, testw2, method='breslow', iter=0)
> fit <- coxph(Surv(time, status) ~x, testw1, weights=wt, method='breslow')
> fitb <- coxph(Surv(time, status) ~x, testw2, method='breslow')
>
> aeq(resid(fit0, type='mart'), (resid(fit0b, type='mart'))[indx])
[1] TRUE
> aeq(resid(fit0, type='scor'), (resid(fit0b, type='scor'))[indx])
Error in a$specials[[i]] : subscript out of bounds
Calls: aeq ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/book6.R’ failed.
Complete output:
> library(survival)
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
>
> # Tests of the weighted Cox model
> # This is section 1.3 of my appendix -- no yet found in any of the
> # printings though, it awaits the next edition
> #
> # Efron approximation
> #
> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
>
> testw1 <- data.frame(time= c(1,1,2,2,2,2,3,4,5),
+ status= c(1,0,1,1,1,0,0,1,0),
+ x= c(2,0,1,1,0,1,0,1,0),
+ wt = c(1,2,3,4,3,2,1,2,1))
> xx <- testw1$wt
>
> # Efron estimate
> byhand <- function(beta, newx=0) {
+ r <- exp(beta)
+ a <- 7*r +3; b<- 4*r+2
+ loglik <- 11*beta - (log(r^2 + 11*r +7) + 10*log(11*r +5)/3 +
+ 10*log(a*2/3 +b)/3 + 10*log(a/3 +b)/3 +2*log(2*r+1))
+
+ hazard <- c(1/(r^2 + 11*r +7),
+ 10/(3*c(11*r +5, a*2/3 +b, a/3+b)), 2/(2*r+1))
+ temp <- c(hazard[1], hazard[1]+hazard[2] + hazard[3]*2/3 + hazard[4]/3,
+ cumsum(hazard)[4:5])
+ risk <- c(r^2, 1,r,r,1,r,1,r,1)
+ expected <- risk* temp[c(1,1,2,2,2,3,3,4,4)]
+
+ # The matrix of weights, one row per obs, one col per death
+ # deaths at 1,2,2,2, and 4
+ riskmat <- matrix(c(1,1,1,1,1,1,1,1,1,
+ 0,0,1,1,1,1,1,1,1,
+ 0,0,2/3,2/3,2/3,1,1,1,1,
+ 0,0,1/3,1/3,1/3,1,1,1,1,
+ 0,0,0,0,0,0,0,1,1), ncol=5)
+ wtmat <- diag(c(r^2, 2, 3*r, 4*r, 3, 2*r, 1, 2*r, 1)) %*% riskmat
+
+ x <- c(2,0,1,1,0,1,0,1,0)
+ xbar <- colSums(x*wtmat)/ colSums(wtmat)
+ imat <- (4*r^2 + 11*r)*hazard[1] - xbar[1]^2 +
+ 10* mean(xbar[2:4] - xbar[2:4]^2) + 2*(xbar[5] - xbar[5]^2)
+
+ status <- c(1,0,1,1,1,0,0,1,0)
+ wt <- c(1,2,3,4,3,2,1,2,1)
+ # Table of sums for score resids
+ hazmat <- riskmat %*% diag(c(1,10/3,10/3, 10/3,2)/colSums(wtmat))
+ dM <- -risk*hazmat #Expected part
+ dM[1,1] <- dM[1,1] +1 # deaths at time 1
+ for (i in 2:4) dM[3:5, i] <- dM[3:5,i] + 1/3
+ dM[8,5] <- dM[8,5] +1
+ mart <- rowSums(dM)
+ resid <-dM * outer(x, xbar ,'-')
+
+ # Increments to the variance of the hazard
+ var.g <- cumsum(hazard^2* c(1,3/10, 3/10, 3/10, 1/2))
+ var.d <- cumsum((xbar-newx)*hazard)
+
+ sxbar <- c(xbar[1], mean(xbar[2:4]), xbar[5]) #xbar for Schoen
+ list(loglik=loglik, imat=imat, hazard=hazard, xbar=xbar,
+ mart=status-expected, expected=expected,
+ score=rowSums(resid), schoen=c(2,1,1,0,1) - sxbar[c(1,2,2,2,3)],
+ varhaz=((var.g + var.d^2/imat)* exp(2*beta*newx))[c(1,4,5)])
+ }
>
> # Verify
> temp <- byhand(0,0)
> aeq(temp$xbar, c(13/19, 11/16, 26/38, 19/28, 2/3))
[1] TRUE
> aeq(temp$hazard, c(1/19, 5/24, 5/19, 5/14, 2/3))
[1] TRUE
>
> fit0 <- coxph(Surv(time, status) ~x, testw1, weights=wt, iter=0)
> fit <- coxph(Surv(time, status) ~x, testw1, weights=wt)
>
> truth0 <- byhand(0,pi)
> aeq(fit0$loglik[1], truth0$loglik)
[1] TRUE
> aeq(1/truth0$imat, fit0$var)
[1] TRUE
> aeq(truth0$mart, fit0$residuals)
[1] TRUE
> aeq(truth0$schoen, resid(fit0, 'schoen'))
Error in a$specials[[i]] : subscript out of bounds
Calls: aeq ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: brier ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/cancer.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), ...)
>
> #
> # Test out all of the routines on a more complex data set
> #
> temp <- survfit(Surv(time, status) ~ ph.ecog, lung)
> summary(temp, times=c(30*1:11, 365*1:3))
Call: survfit(formula = Surv(time, status) ~ ph.ecog, data = lung)
1 observation deleted due to missingness
ph.ecog=0
time n.risk n.event survival std.err lower 95% CI upper 95% CI
30 60 3 0.952 0.0268 0.9012 1.000
60 58 2 0.921 0.0341 0.8562 0.990
90 56 2 0.889 0.0396 0.8146 0.970
120 56 0 0.889 0.0396 0.8146 0.970
150 55 1 0.873 0.0419 0.7946 0.959
180 52 2 0.841 0.0461 0.7553 0.936
210 48 2 0.808 0.0498 0.7164 0.912
240 45 0 0.808 0.0498 0.7164 0.912
270 38 2 0.770 0.0543 0.6709 0.884
300 33 2 0.727 0.0591 0.6203 0.853
330 29 2 0.681 0.0637 0.5670 0.818
365 22 6 0.535 0.0728 0.4100 0.699
730 5 11 0.193 0.0707 0.0943 0.396
ph.ecog=1
time n.risk n.event survival std.err lower 95% CI upper 95% CI
30 111 2 0.982 0.0124 0.9583 1.000
60 110 3 0.956 0.0193 0.9186 0.994
90 104 4 0.920 0.0255 0.8718 0.972
120 99 5 0.876 0.0310 0.8174 0.939
150 93 6 0.823 0.0359 0.7556 0.896
180 82 8 0.751 0.0407 0.6756 0.836
210 68 9 0.666 0.0450 0.5831 0.760
240 57 6 0.604 0.0474 0.5176 0.704
270 53 4 0.561 0.0487 0.4729 0.665
300 46 3 0.527 0.0495 0.4384 0.633
330 40 4 0.480 0.0504 0.3903 0.589
365 34 4 0.431 0.0509 0.3417 0.543
730 7 21 0.114 0.0388 0.0582 0.222
ph.ecog=2
time n.risk n.event survival std.err lower 95% CI upper 95% CI
30 46 5 0.9000 0.0424 0.82057 0.987
60 43 2 0.8600 0.0491 0.76900 0.962
90 40 3 0.8000 0.0566 0.69647 0.919
120 34 4 0.7174 0.0641 0.60216 0.855
150 31 3 0.6541 0.0680 0.53342 0.802
180 26 6 0.5275 0.0719 0.40385 0.689
210 21 4 0.4431 0.0717 0.32266 0.608
240 17 3 0.3766 0.0705 0.26100 0.543
270 17 0 0.3766 0.0705 0.26100 0.543
300 13 3 0.3102 0.0677 0.20223 0.476
330 11 2 0.2624 0.0651 0.16135 0.427
365 9 2 0.2147 0.0614 0.12258 0.376
730 1 6 0.0371 0.0345 0.00601 0.229
ph.ecog=3
time n.risk n.event survival std.err lower 95% CI upper 95% CI
30 1 0 1 0 1 1
60 1 0 1 0 1 1
90 1 0 1 0 1 1
> print(temp[2:3])
Call: survfit(formula = Surv(time, status) ~ ph.ecog, data = lung)
n events median 0.95LCL 0.95UCL
ph.ecog=1 113 82 306 268 429
ph.ecog=2 50 44 199 156 288
>
> temp <- survfit(Surv(time, status)~1, lung, type='fleming',
+ conf.int=.9, conf.type='log-log', error='tsiatis')
> summary(temp, times=30 *1:5)
Call: survfit(formula = Surv(time, status) ~ 1, data = lung, error = "tsiatis",
type = "fleming", conf.int = 0.9, conf.type = "log-log")
time n.risk n.event survival std.err lower 90% CI upper 90% CI
30 219 10 0.956 0.0135 0.928 0.974
60 213 7 0.926 0.0173 0.891 0.950
90 201 10 0.882 0.0213 0.842 0.913
120 189 10 0.838 0.0244 0.793 0.874
150 179 10 0.794 0.0268 0.745 0.834
>
> temp <- survdiff(Surv(time, status) ~ inst, lung, rho=.5)
> print(temp, digits=6)
Call:
survdiff(formula = Surv(time, status) ~ inst, data = lung, rho = 0.5)
n=227, 1 observation deleted due to missingness.
N Observed Expected (O-E)^2/E (O-E)^2/V
inst=1 36 21.190058 17.455181 0.799149708 1.171232977
inst=2 5 3.173330 1.964395 0.744007932 0.860140808
inst=3 19 10.663476 11.958755 0.140294489 0.200472362
inst=4 4 2.245347 3.559344 0.485085848 0.677874608
inst=5 9 5.010883 4.500982 0.057765161 0.077128402
inst=6 14 8.862602 7.078516 0.449665221 0.582743947
inst=7 8 4.445647 4.416133 0.000197254 0.000253632
inst=10 4 2.901923 2.223283 0.207150016 0.249077097
inst=11 18 7.807867 9.525163 0.309611863 0.422142221
inst=12 23 14.009656 12.216768 0.263117640 0.365712493
inst=13 20 9.140983 11.863298 0.624699853 0.874238212
inst=15 6 3.170744 3.558447 0.042241456 0.057938955
inst=16 16 8.870360 9.992612 0.126038005 0.175170113
inst=21 13 9.263733 4.460746 5.171484268 6.149354145
inst=22 17 8.278566 11.971473 1.139171459 1.645863937
inst=26 6 1.627074 3.542694 1.035821659 1.286365543
inst=32 7 1.792468 2.679904 0.293869782 0.343966668
inst=33 2 0.929177 0.416202 0.632249272 0.676682390
Chisq= 15.1 on 17 degrees of freedom, p= 0.5904
>
> # verify that the zph routine does the actual score test
> dtime <- lung$time[lung$status==2]
> lung2 <- survSplit(Surv(time, status) ~ ., lung, cut=dtime)
>
> cfit1 <-coxph(Surv(time, status) ~ ph.ecog + ph.karno + pat.karno + wt.loss
+ + sex + age + strata(inst), lung)
> cfit2 <-coxph(Surv(tstart, time, status) ~ ph.ecog + ph.karno + pat.karno +
+ wt.loss + sex + age + strata(inst), lung2)
> all.equal(cfit1$loglik, cfit2$loglik)
[1] TRUE
> all.equal(coef(cfit1), coef(cfit2))
[1] TRUE
>
> # the above verifies that the data set is correct
> zp1 <- cox.zph(cfit1, transform="log")
Error in a$specials[[i]] : subscript out of bounds
Calls: cox.zph ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: concordance ... eval -> predict -> predict.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/counting.R’ failed.
Complete output:
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
> library(survival)
>
> # Create a "counting process" version of the simplest test 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))
>
> test1b<- list(start= c(0, 3, 0, 0, 5, 0, 6,14, 0, 0, 10,20,30, 0),
+ stop = c(3,10, 10, 5,20, 6,14,20, 30, 10,20,30,40, 10),
+ status=c(0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0),
+ x= c(1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, NA),
+ id = c(3, 3, 4, 5, 5, 6, 6, 6, 7, 1, 1, 1, 1, 2))
>
> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
> #
> # Check out the various residuals under an Efron approximation
> #
> fit0 <- coxph(Surv(time, status)~ x, test1, iter=0)
> fit <- coxph(Surv(time, status) ~x, test1)
> fit0b <- coxph(Surv(start, stop, status) ~ x, test1b, iter=0)
> fitb <- coxph(Surv(start, stop, status) ~x, test1b)
> fitc <- coxph(Surv(time, status) ~ offset(fit$coefficients*x), test1)
> fitd <- coxph(Surv(start, stop, status) ~ offset(fit$coefficients*x), test1b)
>
> aeq(fit0b$coefficients, fit0$coefficients)
[1] TRUE
>
> aeq(resid(fit0), resid(fit0b, collapse=test1b$id))
[1] TRUE
> aeq(resid(fit), resid(fitb, collapse=test1b$id))
[1] TRUE
> aeq(resid(fitc), resid(fitd, collapse=test1b$id))
[1] TRUE
> aeq(resid(fitc), resid(fit))
[1] TRUE
>
> aeq(resid(fit0, type='score'), resid(fit0b, type='score', collapse=test1b$id))
Error in a$specials[[i]] : subscript out of bounds
Calls: aeq ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/coxsurv3.R’ failed.
Complete output:
> library(survival)
> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
> # One more test on coxph survival curves, to test out the individual
> # option. First fit a model with a time dependent covariate.
> # This is test data 2 of section 4 of the validation vignette (appendix E2 of
> # Therneau and Grambsch), i.e. all the results are known in closed form.
> #
> 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) )
>
> # True hazard function, and components of the variance
> lambda <- function(beta, x=0, method="efron") {
+ time <- c(2,3,6,7,8,9)
+ r <- exp(beta)
+ lambda <- 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))
+ xbar <- c(r/(r+1), r/(r+2), 3*r/(3*r +2), 3*r/(3*r+1),
+ 3*r/(3*r+1), (1.5*r)/(3*r +2) + r/(2*r+2))
+
+ if (method == "breslow") {
+ lambda[6] <- 2/(3*r +2)
+ xbar[6] <- 3*r/(3*r+2)
+ }
+
+ list(time=time, lambda=lambda, xbar=xbar)
+ }
>
> fit <- coxph(Surv(start, stop, event) ~x, test2)
> # A curve for someone who never changes
> surv1 <-survfit(fit, newdata=list(x=0), censor=FALSE)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/coxsurv5.R’ failed.
Complete output:
> library(survival)
> aeq <- function(x, y) all.equal(as.vector(x), as.vector(y))
> #
> # Compute the hazard functions for a multi-state Cox model
> #
> coxhaz <- function(y, id, risk, wt, expm=TRUE) {
+ # y should be a multi-state survival
+ if (!inherits(y, "Surv") || attr(y, "type") != "mcounting")
+ stop("invalid response")
+
+ n <- nrow(y)
+ if (missing(id) || length(id) !=n) stop("invalid id")
+
+ if (missing(wt)) wt <- rep(1.0, n)
+ else if (length(wt) !=n || any(wt <=0)) stop("invalid wt")
+
+ # get the current state, and the list of transtions
+ # transitions to censor don't count
+ mcheck <- survcheck(y~1, id= id)
+ states <- mcheck$states
+ nstate <- length(states)
+ istate <- mcheck$istate
+
+ event <- y[,3] > 0
+ temp <- attr(y, 'states')[y[event,3]]
+ tmat <- table(y[event,2], from=istate[event], to=temp)
+ tmat2 <- tapply(wt[event], list(y[event,2], from=mcheck$istate[event],
+ to=temp), sum)
+ tmat2 <- ifelse(is.na(tmat2), 0, tmat2)
+
+
+ # Hazards can be done one at a time. For each of them the risk
+ # weight vector for the subjects can be different.
+ # First organize the material as a 2 dim matrix
+ temp <- apply(tmat, 2:3, sum) # count of transtions
+ keep <- which(temp>0)
+ from <- states[row(temp)[keep]]
+ hlab <- outer(match(rownames(temp), states),
+ match(colnames(temp), states), paste, sep=':')[keep]
+ nhaz <- length(keep)
+ nevent <- matrix(tmat2, nrow(tmat2))[,keep]
+ dtime <- sort(unique(y[event,2]))
+ ntime <- length(dtime)
+ dimnames(nevent) <- list(NULL, hlab)
+ aindex <- cbind(as.numeric(substring(hlab,1,1)),
+ as.numeric(substring(hlab,3,3)))
+
+ if (missing(risk)) risk <- matrix(1, nrow=n, ncol=nhaz)
+ if (!is.matrix(risk) || nrow(risk) != n || ncol(risk) != nhaz)
+ stop("invalid risk matrix")
+ risk <- risk * wt
+
+ # get the weighted at risk at each time
+ wtrisk <- matrix(, length(dtime), nhaz)
+ statematch <- outer(istate, from, function(x, y) x==y)
+ risk <- ifelse(statematch, risk, 0)
+ for (i in 1:ntime) {
+ atrisk <- (y[,1]< dtime[i] & y[,2] >= dtime[i])
+ wtrisk[i,] <- colSums(risk[atrisk,, drop=FALSE])
+ }
+
+ haz <- nevent/ifelse(wtrisk==0, 1, wtrisk) # avoid 0/0
+ chaz<- apply(haz, 2, cumsum)
+
+ # compute the probability in state, with p(0)= 1,0, ..
+ pstate <- matrix(0, ntime+1, nstate)
+ pstate[1,1] <- 1
+ for (i in 1:ntime) {
+ tmat <- matrix(0, nstate, nstate)
+ tmat[aindex] <- haz[i,]
+ if (expm) {
+ diag(tmat) <- -rowSums(tmat)
+ pstate[i+1,] <- pstate[i,] %*% as.matrix(Matrix::expm(tmat))
+ } else {
+ diag(tmat) <- 1-rowSums(tmat)
+ pstate[i+1,] <- pstate[i,] %*% tmat
+ }
+ }
+
+ list(time=dtime, nrisk=wtrisk, nevent=nevent,
+ haz=haz, cumhaz=chaz, states=states, pstate= pstate[-1,])
+ }
>
> 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),
+ state= c(1, 2, 1, 2, 3, 1, 3, 0, 2, 0),
+ x = c(0, 0, 0, 1, 1, 0, 0, 0, 2, 2))
> mtest$state <- factor(mtest$state, 0:3, c("censor", "a", "b", "c"))
>
> # True results
> #
> #time at risk events
> # entry a b c
> #
> #2 1245 4 -> a
> #3 1235 4 5 -> b
> #4 123 4 5 1 -> a
> #5 23 14 5 2 -> b, exits
> #8 3 14 5 4 -> c
> #9 3 1 5 4 1->b, 3->c & exit, 4 censored
> #10 15 1->a, exit
> #11 5 censor
>
> # with all coefficients =0
> check1 <- with(mtest, coxhaz(Surv(t1, t2, state), id))
> fit1 <- survfit(Surv(t1, t2, state) ~1, mtest, id=id)
> aeq(check1$cumhaz, fit1$cumhaz[match(check1$time, fit1$time),])
[1] TRUE
>
> dummy <- data.frame(x=1:2)
> cox0 <- coxph(Surv(t1, t2, state) ~x, iter=0, mtest, id=id)
> cfit0 <- survfit(cox0, newdata=dummy)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/coxsurv6.R’ failed.
Complete output:
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> # Test the survival curve for a fit with shared hazards.
> # Use the pbcseq data set, and turn bilirubin into a time-dependent state with
> # 4 levels, and a shared baseline hazard for the 4 transitions to death.
> # The subtlety is that coefficients for a shared (proportional) baseline hazard
> # are attached to a state, not to an observation.
> # (A bilirubin value of <1 is normal.)
> pbc1 <- pbcseq
> pbc1$bili4 <- cut(pbc1$bili, c(0,1, 2,4, 100),
+ c("normal", "1-2", "2-4", ">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-2", "2-4", ">4", "death"))
> check1 <- survcheck(Surv(tstart, tstop, bstat) ~ 1, istate= bili4,
+ id = id, data=pbc2)
> check1$transitions
to
from normal 1-2 2-4 >4 death (censor)
normal 0 81 10 3 9 77
1-2 61 0 68 15 9 36
2-4 2 33 0 94 12 24
>4 1 3 28 0 110 35
death 0 0 0 0 0 0
>
> fit2 <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
+ c(1:4):5 ~ age / common + shared), id= id, istate=bili4,
+ data=pbc2)
>
> # Before we tackle fit2, start small with just 9 subjects, coefs fixed to
> # simple values to make hand computation easier. There are no transitions
> # from state 3 to death in this subset. The design of coxph is such that
> # the set of (observed transitions) is assumed to be the set of (possible
> # transitions): it does not know that "death" is something special. However,
> # the above call explicitly included 3:5 via the shared argument. This leads
> # to some 'interesting' special cases, not all of which have yet been
> # worked through. More discussion can be found in the "shared coefficients
> # and shared baseline" vignette (survivalVignettes package).
> #
> # The parsecovar3 has an option to only count transitions with at least one
> # actual case --- if we comment out selected lines.
> #
> pbc3 <- subset(pbc2, id < 10)
> pbc3$age <- round(pbc3$age) # easier to do "by hand" sums
> survcheck(Surv(tstart, tstop, bstat) ~1, pbc3, id=id, istate=bili4)$transitions
to
from normal 1-2 2-4 >4 death (censor)
normal 0 4 0 0 1 0
1-2 2 0 3 1 1 1
2-4 0 1 0 5 0 0
>4 0 0 1 0 4 2
death 0 0 0 0 0 0
>
> init3 <- c(.05, .7, 1.5, 3.4)
> fit3 <- coxph(list(Surv(tstart, tstop, bstat) ~ 1,
+ c(1:4):5 ~ age / common + shared), x=TRUE,
+ id= id, istate=bili4, data=pbc3, init= init3, iter=0)
> # a mixed p0 gives a stronger test than our usual (1, 0,0,0,0)
> surv3 <- survfit(fit3, newdata=list(age=50), p0=c(.4, .3, .2, .1, 0))
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: coxph.detail ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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')
Error in a$specials[[i]] : subscript out of bounds
Calls: resid ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/doweight.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 the weighted Cox model
> #
> # Similar data set to test1, but add weights,
> # a double-death/censor tied time
> # a censored last subject
> # The latter two are cases covered only feebly elsewhere.
> #
> # The data set testw2 has the same data, but done via replication
> #
> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
>
> testw1 <- data.frame(time= c(1,1,2,2,2,2,3,4,5),
+ status= c(1,0,1,1,1,0,0,1,0),
+ x= c(2,0,1,1,0,1,0,1,0),
+ wt = c(1,2,3,4,3,2,1,2,1))
> xx <- c(1,2,3,4,3,2,1,2,1)
> testw2 <- data.frame(time= rep(c(1,1,2,2,2,2,3,4,5), xx),
+ status= rep(c(1,0,1,1,1,0,0,1,0), xx),
+ x= rep(c(2,0,1,1,0,1,0,1,0), xx),
+ id= rep(1:9, xx))
> indx <- match(1:9, testw2$id)
> testw2 <- data.frame(time= rep(c(1,1,2,2,2,2,3,4,5), xx),
+ status= rep(c(1,0,1,1,1,0,0,1,0), xx),
+ x= rep(c(2,0,1,1,0,1,0,1,0), xx),
+ id= rep(1:9, xx))
> indx <- match(1:9, testw2$id)
>
> fit0 <- coxph(Surv(time, status) ~x, testw1, weights=wt,
+ method='breslow', iter=0)
> fit0b <- coxph(Surv(time, status) ~x, testw2, ties='breslow', iter=0)
> fit <- coxph(Surv(time, status) ~x, testw1, weights=wt, ties='breslow')
> fitb <- coxph(Surv(time, status) ~x, testw2, ties='breslow')
>
> texp <- function(beta) { # expected, Breslow estimate
+ r <- exp(beta)
+ temp <- cumsum(c(1/(r^2 + 11*r +7), 10/(11*r +5), 2/(2*r+1)))
+ c(r^2, 1,r,r,1,r,1,r,1)* temp[c(1,1,2,2,2,2,2,3,3)]
+ }
> aeq(texp(0), c(1/19, 1/19, rep(103/152, 5), rep(613/456,2))) #verify texp()
[1] TRUE
>
> xbar <- function(beta) { # xbar, Breslow estimate
+ r <- exp(beta)
+ temp <- r* rep(c(2*r + 11, 11/10, 1), c(2, 5, 2))
+ temp * texp(beta)
+ }
>
> fit0
Call:
coxph(formula = Surv(time, status) ~ x, data = testw1, weights = wt,
method = "breslow", iter = 0)
coef exp(coef) se(coef) z p
x 0.0000 1.0000 0.5858 0 1
Likelihood ratio test=0 on 1 df, p=1
n= 9, number of events= 5
> summary(fit)
Call:
coxph(formula = Surv(time, status) ~ x, data = testw1, weights = wt,
ties = "breslow")
n= 9, number of events= 5
coef exp(coef) se(coef) z Pr(>|z|)
x 0.8596 2.3621 0.7131 1.205 0.228
exp(coef) exp(-coef) lower .95 upper .95
x 2.362 0.4233 0.5839 9.556
Concordance= 0.637 (se = 0.161 )
Likelihood ratio test= 1.69 on 1 df, p=0.2
Wald test = 1.45 on 1 df, p=0.2
Score (logrank) test = 1.52 on 1 df, p=0.2
> aeq(resid(fit0), testw1$status - texp(0))
[1] TRUE
> resid(fit0, type='score')
Error in a$specials[[i]] : subscript out of bounds
Calls: resid ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: survexp -> all.vars -> delete.response -> off
Execution halted
Running the tests in ‘tests/expected2.R’ failed.
Complete output:
> library(survival)
> #
> # A Cox model with a factor, followed by survexp.
> #
> pfit2 <- coxph(Surv(time, status > 0) ~ trt + log(bili) +
+ log(protime) + age + platelet + sex, data = pbc)
> esurv <- survexp(~ trt, ratetable = pfit2, data = pbc)
Error in a$specials[[i]] : subscript out of bounds
Calls: survexp -> all.vars -> delete.response -> off
Execution halted
Running the tests in ‘tests/factor.R’ failed.
Complete output:
> #
> # Ensure that factors work in prediction
> #
> 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), ...)
>
> tfit <- coxph(Surv(time, status) ~ age + factor(ph.ecog), lung)
> p1 <- predict(tfit, type='risk')
Error in a$specials[[i]] : subscript out of bounds
Calls: predict -> predict.coxph -> delete.response -> off
Execution halted
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')
Error in a$specials[[i]] : subscript out of bounds
Calls: predict -> predict.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/fr_resid.R’ failed.
Complete output:
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
> library(survival)
>
> #
> # The residual methods treat a sparse frailty as a fixed offset with
> # no variance
> #
> aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> kfit1 <- coxph(Surv(time, status) ~ age + sex +
+ frailty(id, dist='gauss'), kidney)
> tempf <- predict(kfit1, type='terms')[,3]
Error in a$specials[[i]] : subscript out of bounds
Calls: predict ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/frailty.R’ failed.
Complete output:
> library(survival)
> #
> # The constuction of a survival curve with sparse frailties
> #
> # In this case the coefficient vector is kept in two parts, the
> # fixed coefs and the (often very large) random effects coefficients
> # The survfit function treats the second set of coefficients as fixed
> # values, to avoid an unmanagable variance matrix, and behaves like
> # the second fit below.
>
> fit1 <- coxph(Surv(time, status) ~ age + frailty(inst), lung)
> sfit1 <- survfit(fit1)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/frank.R’ failed.
Complete output:
> library(survival)
> #
> # Check out intercept/interaction for Frank H
> #
> age2 <- lung$age - 50
> fit1 <- coxph(Surv(time, status) ~ age * strata(sex), lung)
> fit2 <- coxph(Surv(time, status) ~ age2*strata(sex), lung)
>
> tdata <- data.frame(age=50:60, age2=0:10, sex=c(1,2,1,2,1,2,1,2,1,2,1))
>
> surv1 <- survfit(fit1, tdata)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/jasa.R’ failed.
Complete output:
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
> library(survival)
>
> expect <- survexp(futime ~ 1,
+ rmap = list(age=(accept.dt - birth.dt), sex=1,
+ year=accept.dt, race='white'), jasa, cohort=F,
+ ratetable=survexp.usr)
>
> survdiff(Surv(jasa$futime, jasa$fustat) ~ offset(expect))
Call:
survdiff(formula = Surv(jasa$futime, jasa$fustat) ~ offset(expect))
Observed Expected Z p
75.000 0.587 -97.119 0.000
> # Now fit the 6 models found in Kalbfleisch and Prentice, p139
> sfit.1 <- coxph(Surv(start, stop, event)~ (age + surgery)*transplant,
+ jasa1, method='breslow')
> sfit.2 <- coxph(Surv(start, stop, event)~ year*transplant,
+ jasa1, method='breslow')
> sfit.3 <- coxph(Surv(start, stop, event)~ (age + year)*transplant,
+ jasa1, method='breslow')
> sfit.4 <- coxph(Surv(start, stop, event)~ (year +surgery) *transplant,
+ jasa1, method='breslow')
> sfit.5 <- coxph(Surv(start, stop, event)~ (age + surgery)*transplant + year ,
+ jasa1, method='breslow')
> sfit.6 <- coxph(Surv(start, stop, event)~ age*transplant + surgery + year,
+ jasa1, method='breslow')
>
> summary(sfit.1)
Call:
coxph(formula = Surv(start, stop, event) ~ (age + surgery) *
transplant, data = jasa1, method = "breslow")
n= 170, number of events= 75
coef exp(coef) se(coef) z Pr(>|z|)
age 0.01386 1.01395 0.01813 0.765 0.445
surgery -0.54652 0.57896 0.61091 -0.895 0.371
transplant 0.11572 1.12268 0.32729 0.354 0.724
age:transplant 0.03473 1.03534 0.02725 1.274 0.202
surgery:transplant -0.29037 0.74799 0.75819 -0.383 0.702
exp(coef) exp(-coef) lower .95 upper .95
age 1.014 0.9862 0.9786 1.051
surgery 0.579 1.7272 0.1748 1.917
transplant 1.123 0.8907 0.5911 2.132
age:transplant 1.035 0.9659 0.9815 1.092
surgery:transplant 0.748 1.3369 0.1692 3.306
Concordance= 0.595 (se = 0.037 )
Likelihood ratio test= 12.45 on 5 df, p=0.03
Wald test = 11.62 on 5 df, p=0.04
Score (logrank) test = 12.02 on 5 df, p=0.03
> sfit.2
Call:
coxph(formula = Surv(start, stop, event) ~ year * transplant,
data = jasa1, method = "breslow")
coef exp(coef) se(coef) z p
year -0.2652 0.7671 0.1051 -2.522 0.0117
transplant -0.2871 0.7504 0.5137 -0.559 0.5762
year:transplant 0.1371 1.1469 0.1409 0.973 0.3306
Likelihood ratio test=8.61 on 3 df, p=0.03488
n= 170, number of events= 75
> summary(sfit.3)
Call:
coxph(formula = Surv(start, stop, event) ~ (age + year) * transplant,
data = jasa1, method = "breslow")
n= 170, number of events= 75
coef exp(coef) se(coef) z Pr(>|z|)
age 0.01558 1.01571 0.01734 0.899 0.36887
year -0.27413 0.76023 0.10588 -2.589 0.00962 **
transplant -0.59388 0.55218 0.54222 -1.095 0.27339
age:transplant 0.03380 1.03438 0.02795 1.209 0.22653
year:transplant 0.20228 1.22419 0.14247 1.420 0.15566
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0157 0.9845 0.9818 1.0508
year 0.7602 1.3154 0.6178 0.9356
transplant 0.5522 1.8110 0.1908 1.5981
age:transplant 1.0344 0.9668 0.9792 1.0926
year:transplant 1.2242 0.8169 0.9259 1.6185
Concordance= 0.63 (se = 0.036 )
Likelihood ratio test= 14.85 on 5 df, p=0.01
Wald test = 13.77 on 5 df, p=0.02
Score (logrank) test = 14.06 on 5 df, p=0.02
> sfit.4
Call:
coxph(formula = Surv(start, stop, event) ~ (year + surgery) *
transplant, data = jasa1, method = "breslow")
coef exp(coef) se(coef) z p
year -0.2541 0.7756 0.1077 -2.360 0.0183
surgery -0.2365 0.7894 0.6282 -0.377 0.7065
transplant -0.2969 0.7431 0.5054 -0.588 0.5568
year:transplant 0.1653 1.1797 0.1416 1.167 0.2431
surgery:transplant -0.5499 0.5770 0.7759 -0.709 0.4785
Likelihood ratio test=12.36 on 5 df, p=0.03021
n= 170, number of events= 75
> sfit.5
Call:
coxph(formula = Surv(start, stop, event) ~ (age + surgery) *
transplant + year, data = jasa1, method = "breslow")
coef exp(coef) se(coef) z p
age 0.01504 1.01516 0.01760 0.855 0.3928
surgery -0.42025 0.65689 0.61564 -0.683 0.4948
transplant 0.07412 1.07694 0.33113 0.224 0.8229
year -0.13631 0.87258 0.07097 -1.921 0.0548
age:transplant 0.02691 1.02728 0.02712 0.992 0.3210
surgery:transplant -0.29656 0.74337 0.75800 -0.391 0.6956
Likelihood ratio test=16.2 on 6 df, p=0.0127
n= 170, number of events= 75
> sfit.6
Call:
coxph(formula = Surv(start, stop, event) ~ age * transplant +
surgery + year, data = jasa1, method = "breslow")
coef exp(coef) se(coef) z p
age 0.01527 1.01539 0.01750 0.873 0.3828
transplant 0.04457 1.04558 0.32171 0.139 0.8898
surgery -0.62113 0.53734 0.36786 -1.688 0.0913
year -0.13608 0.87278 0.07090 -1.919 0.0550
age:transplant 0.02703 1.02740 0.02714 0.996 0.3193
Likelihood ratio test=16.06 on 5 df, p=0.006687
n= 170, number of events= 75
>
> # Survival curve for an "average" subject,
> # done once as overall, once via individual method
> surv1 <- survfit(sfit.1, newdata=list(age=-2, surgery=0, transplant=0))
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/model.matrix.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), ...)
>
> #
> # Test out the revised model.matrix code
> #
> 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),
+ z= factor(c('a', 'a', 'b', 'b', 'c', 'c', 'a')),
+ stringsAsFactors=FALSE)
>
> fit1 <- coxph(Surv(time, status) ~ z, test1, iter=1)
> fit2 <- coxph(Surv(time, status) ~z, test1, x=T, iter=1)
> all.equal(model.matrix(fit1), fit2$x)
Error in a$specials[[i]] : subscript out of bounds
Calls: all.equal ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/mstrata.R’ failed.
Complete output:
> #
> # Verify that using multiple states + proportional baselines
> # will mimic a factor covariate
> #
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> tdata <- subset(lung, ph.ecog < 3) # there is only one row with ph.ecog=3
> tdata$state <- factor(tdata$status, 1:2, c("censor", "death"))
> tdata$cstate<- factor(tdata$ph.ecog, 0:2, c("ph0", "ph1", "ph2"))
> tdata$id <- 1:nrow(tdata)
> survcheck(Surv(time, state) ~ 1, id=id, istate=cstate, tdata)
Call:
survcheck(formula = Surv(time, state) ~ 1, data = tdata, id = id,
istate = cstate)
Unique identifiers Observations Transitions
226 226 163
Transitions table:
to
from death (censor)
ph0 37 26
ph1 82 31
ph2 44 6
death 0 0
Number of subjects with 0, 1, ... transitions to each state:
0 1
death 63 163
>
> # standard coxph fit, stratified by the ph0/1/2 groups
> fit1 <- coxph(Surv(time, status) ~ age + sex + factor(ph.ecog), tdata,
+ ties="breslow")
> # multi-state fit, where ph0/1/2 are states with a shared hazard
> fit2 <- coxph(list(Surv(time, state) ~1,
+ 1:4 + 2:4 + 3:4~ age + sex/ common + shared),
+ id=id, istate=cstate, data= tdata, ties="breslow")
>
> aeq(coef(fit1), coef(fit2)) # the names are quite different, values the same
[1] TRUE
> all.equal(fit1$loglik, fit2$loglik)
[1] TRUE
>
> # Three curves in the usual way: ph0, 1, or 2 for all time, common baseline
> csurv1 <- survfit(fit1, newdata=expand.grid(age=65, sex=1, ph.ecog=0:2))
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/multi2.R’ failed.
Complete output:
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> # Check that estimates from a multi-state model agree with single state models
> # Use a simplified version of the myeloid data set
> tdata <- tmerge(myeloid[,1:3], myeloid, id=id, death=event(futime,death),
+ priortx = tdc(txtime), sct= event(txtime))
> tdata$event <- factor(with(tdata, sct + 2*death), 0:2,
+ c("censor", "sct", "death"))
> fit <- coxph(Surv(tstart, tstop, event) ~ trt + sex, tdata, id=id,
+ iter=4, x=TRUE, robust=FALSE)
>
> # Multi-state now defaults to breslow rather than efron
> fit12 <- coxph(Surv(tstart, tstop, event=='sct') ~ trt + sex, tdata,
+ subset=(priortx==0), iter=4, x=TRUE, method='breslow')
> fit13 <- coxph(Surv(tstart, tstop, event=='death') ~ trt + sex, tdata,
+ subset=(priortx==0), iter=4, x=TRUE, method= 'breslow')
> fit23 <- coxph(Surv(tstart, tstop, event=='death') ~ trt + sex, tdata,
+ subset=(priortx==1), iter=4, x=TRUE, method="breslow")
> aeq(coef(fit), c(coef(fit12), coef(fit13), coef(fit23)))
[1] TRUE
> aeq(fit$loglik, fit12$loglik + fit13$loglik + fit23$loglik)
[1] TRUE
> temp <- matrix(0, 6,6)
> temp[1:2, 1:2] <- fit12$var
> temp[3:4, 3:4] <- fit13$var
> temp[5:6, 5:6] <- fit23$var
> aeq(fit$var, temp)
[1] TRUE
>
> # check out model.frame
> fita <- coxph(Surv(tstart, tstop, event) ~ trt, tdata, id=id)
> fitb <- coxph(Surv(tstart, tstop, event) ~ trt, tdata, id=id, model=TRUE)
> all.equal(model.frame(fita), fitb$model)
[1] "Component \"trt\": 'current' is not a factor"
> # model.frame fails due to an internal rule in R, factors vs characters
> # result when the xlev arg is in the call. So model.frame(fita) has trt
> # as a factor, not character.
>
> #check residuals. As of 3.8-4 the martingale resids are a matrix with one
> # row per row of input data, one col per transition, and structural zeros
> # for obs that are not at risk
> scheck <- survcheck(Surv(tstart, tstop, event) ~1, tdata, id=id)
> istate <- as.numeric(scheck$istate) # istate==1: obs at risk for 1:2 or 1:3
> aeq(residuals(fit)[istate==1, 1], residuals(fit12))
[1] TRUE
> aeq(residuals(fit)[istate==1, 2], residuals(fit13))
[1] TRUE
> aeq(residuals(fit)[istate==2, 3], residuals(fit23))
[1] TRUE
>
> # score residuals are an array with a row per obs, one col per coef
> temp <- residuals(fit, type='score')
Error in a$specials[[i]] : subscript out of bounds
Calls: residuals ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/nested.R’ failed.
Complete output:
> library(survival)
> #
> # A test of nesting. It makes sure the model.frame is built correctly
> #
> tfun <- function(fit, mydata) {
+ survfit(fit, newdata=mydata)
+ }
>
> myfit <- coxph(Surv(time, status) ~ age + factor(sex), lung)
>
> temp1 <- tfun(myfit, lung[1:5,])
Error in a$specials[[i]] : subscript out of bounds
Calls: tfun ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/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 the coxph program on the Ovarian data
> summary(survfit(Surv(futime, fustat)~1, data=ovarian), censor=TRUE)
Call: survfit(formula = Surv(futime, fustat) ~ 1, data = ovarian)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
59 26 1 0.962 0.0377 0.890 1.000
115 25 1 0.923 0.0523 0.826 1.000
156 24 1 0.885 0.0627 0.770 1.000
268 23 1 0.846 0.0708 0.718 0.997
329 22 1 0.808 0.0773 0.670 0.974
353 21 1 0.769 0.0826 0.623 0.949
365 20 1 0.731 0.0870 0.579 0.923
377 19 0 0.731 0.0870 0.579 0.923
421 18 0 0.731 0.0870 0.579 0.923
431 17 1 0.688 0.0919 0.529 0.894
448 16 0 0.688 0.0919 0.529 0.894
464 15 1 0.642 0.0965 0.478 0.862
475 14 1 0.596 0.0999 0.429 0.828
477 13 0 0.596 0.0999 0.429 0.828
563 12 1 0.546 0.1032 0.377 0.791
638 11 1 0.497 0.1051 0.328 0.752
744 10 0 0.497 0.1051 0.328 0.752
769 9 0 0.497 0.1051 0.328 0.752
770 8 0 0.497 0.1051 0.328 0.752
803 7 0 0.497 0.1051 0.328 0.752
855 6 0 0.497 0.1051 0.328 0.752
1040 5 0 0.497 0.1051 0.328 0.752
1106 4 0 0.497 0.1051 0.328 0.752
1129 3 0 0.497 0.1051 0.328 0.752
1206 2 0 0.497 0.1051 0.328 0.752
1227 1 0 0.497 0.1051 0.328 0.752
>
> # Various models
> coxph(Surv(futime, fustat)~ age, data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ age, data = ovarian)
coef exp(coef) se(coef) z p
age 0.16162 1.17541 0.04974 3.249 0.00116
Likelihood ratio test=14.29 on 1 df, p=0.0001564
n= 26, number of events= 12
> coxph(Surv(futime, fustat)~ resid.ds, data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ resid.ds, data = ovarian)
coef exp(coef) se(coef) z p
resid.ds 1.2092 3.3507 0.6724 1.798 0.0721
Likelihood ratio test=3.76 on 1 df, p=0.05251
n= 26, number of events= 12
> coxph(Surv(futime, fustat)~ rx, data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ rx, data = ovarian)
coef exp(coef) se(coef) z p
rx -0.5964 0.5508 0.5870 -1.016 0.31
Likelihood ratio test=1.05 on 1 df, p=0.3052
n= 26, number of events= 12
> coxph(Surv(futime, fustat)~ ecog.ps, data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ ecog.ps, data = ovarian)
coef exp(coef) se(coef) z p
ecog.ps 0.3984 1.4894 0.5864 0.679 0.497
Likelihood ratio test=0.47 on 1 df, p=0.4935
n= 26, number of events= 12
>
> coxph(Surv(futime, fustat)~ resid.ds + rx + ecog.ps, data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ resid.ds + rx + ecog.ps,
data = ovarian)
coef exp(coef) se(coef) z p
resid.ds 1.3466 3.8444 0.6803 1.980 0.0478
rx -0.7493 0.4727 0.5949 -1.260 0.2078
ecog.ps 0.4528 1.5727 0.5903 0.767 0.4431
Likelihood ratio test=6.03 on 3 df, p=0.1101
n= 26, number of events= 12
> coxph(Surv(futime, fustat)~ age + rx + ecog.ps, data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ age + rx + ecog.ps, data = ovarian)
coef exp(coef) se(coef) z p
age 0.1470 1.1583 0.0463 3.175 0.0015
rx -0.8146 0.4428 0.6342 -1.285 0.1990
ecog.ps 0.1032 1.1087 0.6064 0.170 0.8649
Likelihood ratio test=15.92 on 3 df, p=0.00118
n= 26, number of events= 12
> coxph(Surv(futime, fustat)~ age + resid.ds + ecog.ps, data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ age + resid.ds + ecog.ps,
data = ovarian)
coef exp(coef) se(coef) z p
age 0.14240 1.15303 0.05196 2.740 0.00614
resid.ds 0.66261 1.93984 0.75045 0.883 0.37727
ecog.ps 0.16634 1.18097 0.61486 0.271 0.78675
Likelihood ratio test=15.1 on 3 df, p=0.001732
n= 26, number of events= 12
> coxph(Surv(futime, fustat)~ age + resid.ds + rx, data=ovarian)
Call:
coxph(formula = Surv(futime, fustat) ~ age + resid.ds + rx, data = ovarian)
coef exp(coef) se(coef) z p
age 0.1285 1.1372 0.0473 2.718 0.00657
resid.ds 0.6964 2.0065 0.7585 0.918 0.35858
rx -0.8489 0.4279 0.6392 -1.328 0.18416
Likelihood ratio test=16.77 on 3 df, p=0.0007889
n= 26, number of events= 12
>
> # Residuals
> fit <- coxph(Surv(futime, fustat)~ age + resid.ds + rx + ecog.ps, data=ovarian)
> resid(fit)
1 2 3 4 5 6
0.84103277 0.54424388 0.59670824 -0.11281376 0.75111588 -0.32609026
7 8 9 10 11 12
0.59998927 0.29570718 -2.15325805 0.76243469 0.06474272 -0.11680752
13 14 15 16 17 18
-1.22562781 -0.63474839 -0.07535824 -0.17058905 -0.22986038 -0.14654862
19 20 21 22 23 24
-0.18762920 -0.12771548 -0.53373114 -0.65480022 0.95866131 0.82111675
25 26
0.55136554 -0.09154014
> resid(fit, 'dev')
1 2 3 4 5 6
1.41281595 0.69505907 0.78916003 -0.47500266 1.13106322 -0.80757694
7 8 9 10 11 12
0.79532966 0.33122166 -2.07521471 1.16179002 0.06619519 -0.48333740
13 14 15 16 17 18
-1.56564862 -1.12671948 -0.38822221 -0.58410453 -0.67802711 -0.54138455
19 20 21 22 23 24
-0.61258338 -0.50540178 -1.03318066 -0.54976346 2.11059000 1.34157009
25 26
0.70736314 -0.42787881
> resid(fit, 'scor')
Error in a$specials[[i]] : subscript out of bounds
Calls: resid ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/prednew.R’ failed.
Complete output:
> #
> # Make sure that the newdata argument works for various
> # predictions
> # We purposely use a subset of the lung data that has only some
> # of the levels of ph.ecog
> library(survival)
> options(na.action=na.exclude, contrasts=c('contr.treatment', 'contr.poly'))
> aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
>
> myfit <- coxph(Surv(time, status) ~ age + factor(ph.ecog) + strata(sex), lung)
>
> keep <- which(lung$inst<13 & (lung$ph.ecog==1 | lung$ph.ecog==2))
> p1 <- predict(myfit, type='lp')
Error in a$specials[[i]] : subscript out of bounds
Calls: predict -> predict.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/predsurv.R’ failed.
Complete output:
> #
> # Verify that predict(coxfit) agrees with survfit for type= expected and survival
> # This also acts as a check on summary.survfit
> #
> library(survival)
> aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> fit1 <- coxph(Surv(time, status) ~ age + ph.ecog, lung)
>
> d1 <- data.frame(age= 40 * 1:8*5, ph.ecog= rep(0:2, length=8))
> curves <- survfit(fit1, newdata= d1)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/pspline.R’ failed.
Complete output:
> library(survival)
> #
> # Tests with the pspline function, to verify the prediction aspects
> #
> options(na.action=na.exclude)
> aeq <- function(x,y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> spfit <- coxph(Surv(time, status) ~ pspline(age) + ph.ecog, lung)
>
> spfit2 <- coxph(Surv(time, status) ~ pspline(age) + ph.ecog, lung, x=TRUE)
> x2 <- model.matrix(spfit)
Error in a$specials[[i]] : subscript out of bounds
Calls: model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/quantile.R’ failed.
Complete output:
> #
> # Formal test of the quantile routine for survfit
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> # There are 8 cases: strata Y/N, ncol(surv) >1, conf.int = T/F
> # Subcase: the quantile exactly agrees with a horizontal segment of
> # the curve or not.
> # First do the 4 cases where fit$surv is a vector
> #
> test1 <- data.frame(time= c(9, 3,1,1,6,6,8, 10),
+ status=c(1,NA,1,0,1,1,0, 0),
+ x= c(0, 2,1,1,1,0,0, 0))
>
> # True survival = (6/7) * (3/5) * (1/2) for overall
> # The q's are chosen to include a point < first jump, mid, after last jump,
> # and exact intersections with the "flats" of the curve.
> #
> qq <- c(13/14, 6/7, 2/3, .5, 9/35, .1)
>
> # Nothing on the right hand side, simple survival (no strata)
> fit1 <- survfit(Surv(time, status) ~ 1, test1, conf.type='none')
> aeq(quantile(fit1, 1-qq), c(1, 3.5, 6, 9, 9.5, NA)) #without conf.int
[1] TRUE
>
> fit2 <- survfit(Surv(time, status) ~ 1, test1) #with conf.int
> aeq(quantile(fit2, 1-qq),
+ list(quantile = c(1, 3.5, 6, 9, 9.5, NA),
+ lower = c(1,1,1,6,6,9),
+ upper = rep(as.numeric(NA), 6)), check.attributes=FALSE)
[1] TRUE
> aeq(quantile(fit2, 1-qq, FALSE), c(1, 3.5, 6, 9, 9.5, NA))
[1] TRUE
>
>
> # Now a variable on the right (strata in the result)
> # curve 0: (t=6, S=3/4), (t=9, S=3/8)
> # curve 1: (t=1, S=2/3), (t=6, S= 0)
> fit1 <- survfit(Surv(time, status) ~ x, test1, conf.type='none')
> aeq(quantile(fit1, 1-qq),
+ matrix(c(6,6,9,9,NA,NA, 1,1,3.5, 6,6,6), nrow=2, byrow=T))
[1] TRUE
>
> fit2 <- survfit(Surv(time, status) ~ x, test1)
> aeq(quantile(fit2, 1-qq, FALSE),
+ matrix(c(6,6,9,9,NA,NA, 1,1,3.5, 6,6,6), nrow=2, byrow=T))
[1] TRUE
>
> temp <- quantile(fit2, 1-qq)
> aeq(temp$quantile, matrix(c(6,6,9,9,NA,NA, 1,1,3.5, 6,6,6), nrow=2, byrow=T))
[1] TRUE
> aeq(temp$lower, matrix(c(6,6,6,6,9,9, 1,1,1,1, NA,NA), nrow=2, byrow=T))
[1] TRUE
> aeq(temp$upper, rep(as.numeric(NA), 12))
[1] TRUE
>
> # Second major case set -- a survfit object where fit$surv is a matrix
> # This arises from coxph models
> # There is only 1 subject with ph.ecog=3 which is a nice edge case
> cfit <- coxph(Surv(time, status) ~ age + strata(ph.ecog), lung)
> sfit <- survfit(cfit, newdata=data.frame(age=c(50, 70)))
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/residms.R’ failed.
Complete output:
> #
> # Test residuals for a coxphms model.
> # This is done by explicitly fitting the three individual transitions and
> # computing the residuals on those.
> # Use the same myeloid data set as multi2.R
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
> table2 <- function(...) table(..., useNA= "ifany")
> tdata <- tmerge(myeloid[,1:4], myeloid, id=id, death=event(futime,death),
+ priortx = tdc(txtime), sct= event(txtime))
> tdata$event <- factor(with(tdata, sct + 2*death), 0:2,
+ c("censor", "sct", "death"))
>
> # Make the processing a little harder by making sex missing for 3 people,
> # use it as a covariate for transitions to death, but not entry to sct
> # flt3 missing for 3, a covariate for entry:sct and sct:death, but not 1:3.
> #
> # Subject 275 for instance has two rows (427,428): entry to sct, then sct to
> # censor, missing sex. They will not be omitted from the data frame since
> # they are at risk for entry:sct, but not sct:death. Subject 271 is missing
> # flt for both entry:sct and sct:death, but will be counted in the risk set
> # for an entry:death transtion. Only one subject, 273, is deleted completely.
> # The 646 subject/1009 obs data set is first thinned down to 645/1006; this is
> # what shows as 'deleted due to missing' in the printout. Rows with a
> # prior stem cell transplant (priortx) and missing sex or flt have
> # nowhere to go, nor id 273 who is missing both.
> tdata$sex[tdata$id %in% 273:275] <- NA # obs 425 to 428
> tdata$flt3[tdata$id %in% 271:273] <- NA # obs 422 to 425
> subset(tdata, id %in% 270:275)
id trt sex flt3 tstart tstop death priortx sct event
420 270 A m C 0 121 0 0 1 sct
421 270 A m C 121 191 1 1 0 death
422 271 A f <NA> 0 223 0 0 1 sct
423 271 A f <NA> 223 707 1 1 0 death
424 272 B m <NA> 0 103 1 0 0 death
425 273 A <NA> <NA> 0 283 1 0 0 death
426 274 B <NA> C 0 293 1 0 0 death
427 275 B <NA> A 0 102 0 0 1 sct
428 275 B <NA> A 102 1442 0 1 0 censor
>
> # Now scramble the data set, which tests the logic in multimiss
> set.seed(1960)
> index <- sample(1:nrow(tdata), nrow(tdata))
> tdata <- tdata[index,]
>
> check <- survcheck(Surv(tstart, tstop, event) ~ 1, tdata, id=id)
> fit <- coxph(list(Surv(tstart, tstop, event) ~ trt,
+ 1:3 + 2:3 ~ sex,
+ 1:2 + 2:3 ~ flt3), tdata, id=id)
>
> aeq(check$transitions, fit$transitions + c(1,0,0, 2,1,0, 0,0,0))
[1] TRUE
> # The above is due to the difference between coxph, which has the more
> # sophisticated list form of a formula, and survfit/survcheck which do not.
> # survcheck with ~1 on the right uses all 1009 obs, 3 more obs than fit, but
> # with ~sex on the right it will use fewer rows than fit
>
> # Multi-state now defaults to breslow rather than efron
> # The id option allows for collapse=TRUE later
> fit12 <- coxph(Surv(tstart, tstop, event=='sct') ~ trt + flt3, tdata, id=id,
+ subset=(priortx==0), iter=4, x=TRUE, method='breslow')
> fit13 <- coxph(Surv(tstart, tstop, event=='death') ~ trt + sex, tdata, id =id,
+ subset=(priortx==0), iter=4, x=TRUE, method= 'breslow')
> fit23 <- coxph(Surv(tstart, tstop, event=='death') ~ trt + sex + flt3,
+ tdata, id=id,
+ subset=(priortx==1), iter=4, x=TRUE, method="breslow")
>
> # martingale residuals
> # one row per retained obs, one col per transition
> rr1 <- resid(fit)
> aeq(dim(rr1), c(nrow(tdata)-3, 3))
[1] TRUE
>
> #
> # Obs 1-2 start in state (s0) so contribute to the 1:2 and 1:3 transitions,
> # columns 1 and 2 of rr1, while rr1[1:2, 3] =0.
> all(rr1[1:2,3] ==0)
[1] TRUE
>
> # Comparing the overall fit to the per-transaction ones is a bit of nuisance
> # rr1 = resid(fit) has 1006 rows and 3 columns.
> # Rows where priortx=0 are not at risk for a 2:3 transition, the 2:3 colum will
> # be deterministicly 0 in rr1. Where priortx==1 the 1:2 and 1:3 columns are 0.
> mf2 <- tdata[-fit$na.action, ] # the 1006 subset
> dim(mf2)
[1] 1006 10
> all(rr1[mf2$priortx==0, 3] ==0)
[1] TRUE
> all(rr1[mf2$priortx==1, 1:2] ==0)
[1] TRUE
>
> # fit12 has 643 obs while mf2 has 645 with priortx==0, two of those have missing
> # flt. To match the fit residuals with fit12 we need to pick off the proper
> # 643 rows, indx1 below. The 1:3 transition is similar, but with sex.
> with(mf2, table2(priortx))
priortx
0 1
645 361
> with(mf2, table2(priortx, flt3))
flt3
priortx A B C <NA>
0 149 317 177 2
1 71 182 108 0
> indx1 <- (mf2$priortx==0 & !is.na(mf2$flt3))
> aeq(rr1[indx1, 1], resid(fit12))
[1] TRUE
> indx2 <- (mf2$priortx==0 & !is.na(mf2$sex))
> aeq(rr1[indx2, 2], resid(fit13))
[1] TRUE
> indx3 <- (mf2$priortx==1 & !is.na(mf2$sex) & !is.na(mf2$flt3))
> aeq(rr1[indx3, 3], resid(fit23))
[1] TRUE
>
> # dfbeta residuals have dim of (data row, coefficient)
> rr3 <- resid(fit, type='dfbeta')
Error in a$specials[[i]] : subscript out of bounds
Calls: resid ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/royston.R’ failed.
Complete output:
> # Verify the values found in the Royston paper
> library(survival)
>
> pbc2 <- na.omit(pbc[,-1]) # no id variable, no missings
>
> pfit1 <- coxph(Surv(time, status==2) ~ . + log(bili) - bili, pbc2,
+ ties="breslow")
> # backwards elimination was used to eliminate all but 8
> pfit2 <- coxph(Surv(time, status==2) ~ age + log(bili) + edema + albumin +
+ stage + copper, data=pbc2, ties="breslow")
>
> temp <- rbind(royston(pfit1), royston(pfit1, adjust=TRUE),
+ royston(pfit2), royston(pfit2, adjust=TRUE))
Error in a$specials[[i]] : subscript out of bounds
Calls: rbind ... royston -> predict -> predict.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/strata2.R’ failed.
Complete output:
> #
> # New tests 4/2010 to validate strata by covariate interactions
> #
> 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))
>
> tdata <- lung
> tdata$sex <- lung$sex +3
>
> # Both of these should produce warning messages about singular X, since there
> # are ph.ecog=3 subjects in only 1 of the strata.
> # Does not affect the test
> fit1 <- coxph(Surv(time, status) ~ age + sex:strata(ph.ecog), lung)
> fit2 <- coxph(Surv(time, status) ~ age + sex:strata(ph.ecog), tdata)
>
> aeq(fit1$coef, fit2$coef)
[1] TRUE
> aeq(fit1$var, fit2$var)
[1] TRUE
> aeq(predict(fit1), predict(fit2))
Error in a$specials[[i]] : subscript out of bounds
Calls: aeq ... as.vector -> predict -> predict.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/stratatest.R’ failed.
Complete output:
> options(na.action=na.exclude) # preserve missings
> options(contrasts=c('contr.treatment', 'contr.poly')) #ensure constrast type
> library(survival)
>
> #
> # Trivial test of stratified residuals
> # Make a second strata = replicate of the first, and I should get the
> # exact same answers
> 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) )
>
> temp <- as.matrix(test1)
> n <- nrow(temp)
> ndead<- sum(test1$status[!is.na(test1$status)])
> temp <- data.frame(rbind(temp, temp)) #later releases of S have rbind.data.frame
> tstrat <- rep(1:2, c(n,n))
>
> fit1 <- coxph(Surv(time, status) ~x, test1)
> fit2 <- coxph(Surv(time, status) ~x + strata(tstrat), temp)
>
> all.equal(resid(fit1) , (resid(fit2))[1:n])
[1] TRUE
> all.equal(resid(fit1, type='score') , (resid(fit2, type='score'))[1:n])
Error in a$specials[[i]] : subscript out of bounds
Calls: all.equal ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/summarydf.R’ failed.
Complete output:
> #
> # Check out the data.frame option of summary.survfit
> #
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> surv1 <- survfit(Surv(time, status) ~ ph.ecog, lung)
> summ1a <- summary(surv1, times= round(30.5*1:32))
> summ1b <- summary(surv1, times= round(30.5*1:32), data.frame=TRUE)
> test <- with(summ1a, data.frame(time, n.risk, n.event, n.censor,
+ surv, cumhaz, std.err, std.chaz, lower, upper, strata))
> all.equal(test, summ1b)
[1] TRUE
>
> # more complex
> cox1 <- coxph(Surv(time, status) ~ ph.ecog + sex + strata(inst), lung)
> dummy <- expand.grid(ph.ecog=0:2, sex=1:2)
> surv2 <- survfit(cox1, newdata=dummy)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
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))
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... survfit.coxph -> model.matrix.coxph -> delete.response -> off
Execution halted
Running the tests in ‘tests/timeline.R’ failed.
Complete output:
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> #
> # Follow on to the checkSurv2.R tests
> # Multi-state data set with repeated states and covariates
> #
> data1 <- subset(pbcseq,, c(id, day, age, bili, albumin, edema, ast))
> data1$bstat <- as.numeric( cut(data1$bili, c(0, 1, 4, 100)))
> temp <- subset(pbcseq, !duplicated(id), c(id, futime, status))
> names(temp) <- c("id", "day", "bstat")
> temp$bstat <- 4*(temp$bstat==2)
> pdata <- merge(data1, temp, all=TRUE)
> pdata$bstat <- factor(pdata$bstat, 0:4,
+ c("censor", "normal", "1-4", "4+", "death"))
>
> # Throw in some missings: id 1 is missing edema, id 2 missing even visit
> # albumins, id 3 missing the first ast, bstat2 is missing on row 2.
> # Models with edema should throw out subject 1, those with albumin will have
> # albumin filled in, and those with ast will have delayed entry of subject 3.
> # A model using bstat2 should fail survcheck since id 1 will have a "hole"
> pdata$edema[pdata$id==1] <- NA
> indx <- which(pdata$id ==2)
> indx <- indx[indx%%2== 1]
> pdata$albumin[indx] <- NA
> pdata$ast[which(pdata$id==3)[1]] <- NA
> pdata$bstat2 <- pdata$bstat
> pdata$bstat2[2] <- NA
>
> subset(pdata, id<4) # the printout looks good
id day bstat age bili albumin edema ast bstat2
1 1 0 4+ 58.76523 14.5 2.60 NA 138.0 4+
2 1 192 4+ 58.76523 21.3 2.94 NA 6.2 <NA>
3 1 400 death NA NA NA NA NA death
4 2 0 1-4 56.44627 1.1 4.14 0.0 113.5 1-4
5 2 182 normal 56.44627 0.8 NA 0.0 139.5 normal
6 2 365 normal 56.44627 1.0 3.55 0.0 144.2 normal
7 2 768 1-4 56.44627 1.9 NA 0.0 144.2 1-4
8 2 1790 1-4 56.44627 2.6 3.32 0.5 131.8 1-4
9 2 2151 1-4 56.44627 3.6 NA 1.0 131.8 1-4
10 2 2515 4+ 56.44627 4.2 2.73 1.0 145.7 4+
11 2 2882 1-4 56.44627 3.6 NA 1.0 119.0 1-4
12 2 3226 4+ 56.44627 4.6 2.67 1.0 88.0 4+
13 2 5169 censor NA NA NA NA NA censor
14 3 0 1-4 70.07255 1.4 3.48 0.5 NA 1-4
15 3 176 1-4 70.07255 1.1 3.29 0.0 69.8 1-4
16 3 364 1-4 70.07255 1.5 3.57 0.5 57.4 1-4
17 3 743 1-4 70.07255 1.8 3.25 0.5 88.4 1-4
18 3 1012 death NA NA NA NA NA death
> cdata <- fromtimeline(Surv(day, bstat) ~ ., pdata, id=id)
> subset(cdata, id<4)
id age bili albumin edema ast bstat2 istate day1 day2 bstat
1 1 58.76523 14.5 2.60 NA 138.0 4+ 4+ 0 192 censor
2 1 58.76523 21.3 2.94 NA 6.2 4+ 4+ 192 400 death
4 2 56.44627 1.1 4.14 0.0 113.5 1-4 1-4 0 182 normal
5 2 56.44627 0.8 4.14 0.0 139.5 normal normal 182 365 censor
6 2 56.44627 1.0 3.55 0.0 144.2 normal normal 365 768 1-4
7 2 56.44627 1.9 3.55 0.0 144.2 1-4 1-4 768 1790 censor
8 2 56.44627 2.6 3.32 0.5 131.8 1-4 1-4 1790 2151 censor
9 2 56.44627 3.6 3.32 1.0 131.8 1-4 1-4 2151 2515 4+
10 2 56.44627 4.2 2.73 1.0 145.7 4+ 4+ 2515 2882 1-4
11 2 56.44627 3.6 2.73 1.0 119.0 1-4 1-4 2882 3226 4+
12 2 56.44627 4.6 2.67 1.0 88.0 4+ 4+ 3226 5169 censor
14 3 70.07255 1.4 3.48 0.5 NA 1-4 1-4 0 176 censor
15 3 70.07255 1.1 3.29 0.0 69.8 1-4 1-4 176 364 censor
16 3 70.07255 1.5 3.57 0.5 57.4 1-4 1-4 364 743 censor
17 3 70.07255 1.8 3.25 0.5 88.4 1-4 1-4 743 1012 death
> # in the above, bstat2 also gets the lvcf treatment -- it's just a covariate
>
> fit1a <- coxph(list(Surv2(day, bstat) ~ 1,
+ (1:3):4 ~ edema + albumin + ast), pdata, id=id)
>
> fit1b <- coxph(list(Surv(day1, day2, bstat) ~ 1,
+ (1:3):4 ~ edema + albumin + ast),
+ cdata, id=id, istate=istate)
> all.equal(fit1a[c("coefficients", "var", "loglik", "score")],
+ fit1b[c("coefficients", "var", "loglik", "score")])
[1] TRUE
>
> # reconstruct the model frame
> fit1c <- coxph(list(Surv2(day, bstat) ~ 1,
+ (1:3):4 ~ edema + albumin + ast), pdata, id=id,
+ model=TRUE)
> m1 <- model.frame(fit1a)
> m2 <- model.frame(fit1b)
> m3 <- model.frame(fit1c)
> all.equal(m1,m3)
[1] TRUE
> # m2 will have a different name for the first column, and that first col
> # has a different "inputAttributes" attribute
> all.equal(m2[,-1], m3[,-1])
[1] TRUE
> aeq(m2[,1], m3[,1])
[1] TRUE
>
> dummy <- expand.grid(edema=0:1, albumin=3, ast=c(75, 150))
> s1 <- survfit(fit1a, newdata=dummy)
Error in a$specials[[i]] : subscript out of bounds
Calls: survfit ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
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')
>
> # Fit a model with an NA coef
> # there is only one subject with ph.ecog=3, so the interaction of
> # ph.ecog=3 with sex is not estimable.
> fit2 <- coxph(Surv(time, status) ~ factor(ph.ecog) *sex + age, lung)
> fit2
Call:
coxph(formula = Surv(time, status) ~ factor(ph.ecog) * sex +
age, data = lung)
coef exp(coef) se(coef) z p
factor(ph.ecog)1 0.335017 1.397964 0.609291 0.550 0.5824
factor(ph.ecog)2 0.618092 1.855384 0.686960 0.900 0.3683
factor(ph.ecog)3 1.941570 6.969687 1.033409 1.879 0.0603
sex -0.641345 0.526583 0.385890 -1.662 0.0965
age 0.009667 1.009714 0.009601 1.007 0.3140
factor(ph.ecog)1:sex 0.063350 1.065400 0.451918 0.140 0.8885
factor(ph.ecog)2:sex 0.221865 1.248403 0.507016 0.438 0.6617
factor(ph.ecog)3:sex NA NA 0.000000 NA NA
Likelihood ratio test=31.09 on 7 df, p=5.976e-05
n= 227, number of events= 164
(1 observation deleted due to missingness)
>
> # The marginal estimate for that group, averaged over sex, is not available
> # so becomes NA, as does the test for "all 4 ph.ecog population averages are
> # the same
> yf1 <- yates(fit2, ~ph.ecog)
Error in a$specials[[i]] : subscript out of bounds
Calls: yates -> delete.response -> off
Execution halted
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)
Error in a$specials[[i]] : subscript out of bounds
Calls: cox.zph ... model.matrix -> model.matrix.coxph -> delete.response -> off
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 3.8-12
Check: re-building of vignette outputs
Result: ERROR
Error(s) in re-building vignettes:
...
--- re-building ‘adjcurve.Rnw’ using Sweave
Error: processing vignette 'adjcurve.Rnw' failed with diagnostics:
chunk 20 (label = flc8)
Error in a$specials[[i]] : subscript out of bounds
--- failed re-building ‘adjcurve.Rnw’
--- re-building ‘approximate.Rnw’ using Sweave
--- finished re-building ‘approximate.Rnw’
--- re-building ‘compete.Rnw’ using Sweave
Error: processing vignette 'compete.Rnw' failed with diagnostics:
chunk 12 (label = PCMcurve)
Error in a$specials[[i]] : subscript out of bounds
--- failed re-building ‘compete.Rnw’
--- re-building ‘concordance.Rnw’ using Sweave
Error: processing vignette 'concordance.Rnw' failed with diagnostics:
chunk 8
Error in a$specials[[i]] : subscript out of bounds
--- failed re-building ‘concordance.Rnw’
--- re-building ‘matrix.Rnw’ using Sweave
--- finished re-building ‘matrix.Rnw’
--- re-building ‘methods.Rnw’ using Sweave
Error: processing vignette 'methods.Rnw' failed with diagnostics:
chunk 8 (label = tdcoef)
Error in a$specials[[i]] : subscript out of bounds
--- failed re-building ‘methods.Rnw’
--- re-building ‘modelframe.Rnw’ using Sweave
--- finished re-building ‘modelframe.Rnw’
--- re-building ‘multi.Rnw’ using Sweave
--- finished re-building ‘multi.Rnw’
--- re-building ‘other.Rnw’ using Sweave
--- finished re-building ‘other.Rnw’
--- re-building ‘population.Rnw’ using Sweave
Error: processing vignette 'population.Rnw' failed with diagnostics:
chunk 17 (label = coxfit2)
Error in a$specials[[i]] : subscript out of bounds
--- failed re-building ‘population.Rnw’
--- re-building ‘redistribute.Rnw’ using Sweave
Error: processing vignette 'redistribute.Rnw' failed with diagnostics:
chunk 4 (label = rotterdam)
Error in a$specials[[i]] : subscript out of bounds
--- failed re-building ‘redistribute.Rnw’
--- re-building ‘splines.Rnw’ using Sweave
Error: processing vignette 'splines.Rnw' failed with diagnostics:
chunk 2 (label = mplot)
Error in a$specials[[i]] : subscript out of bounds
--- failed re-building ‘splines.Rnw’
--- re-building ‘survival.Rnw’ using Sweave
Error: processing vignette 'survival.Rnw' failed with diagnostics:
chunk 37 (label = lung1)
Error in a$specials[[i]] : subscript out of bounds
--- failed re-building ‘survival.Rnw’
--- re-building ‘tiedtimes.Rnw’ using Sweave
--- finished re-building ‘tiedtimes.Rnw’
--- re-building ‘timeline.Rnw’ using Sweave
--- finished re-building ‘timeline.Rnw’
--- re-building ‘validate.Rnw’ using Sweave
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Ran out of iterations and did not converge
Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
Ran out of iterations and did not converge
--- finished re-building ‘validate.Rnw’
SUMMARY: processing the following files failed:
‘adjcurve.Rnw’ ‘compete.Rnw’ ‘concordance.Rnw’ ‘methods.Rnw’
‘population.Rnw’ ‘redistribute.Rnw’ ‘splines.Rnw’ ‘survival.Rnw’
Error: Vignette re-building failed.
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 3.8-12
Check: compiled code
Result: NOTE
Note: information on .o files for x64 is not available
File 'd:/Rcompile/CRANpkg/lib/4.7/survival/libs/x64/survival.dll':
Found '_exit', possibly from '_exit' (C)
Found 'abort', possibly from 'abort' (C), 'runtime' (Fortran)
Compiled code should not call entry points which might terminate R nor
write to stdout/stderr instead of to the console, nor use Fortran I/O
nor system RNGs nor [v]sprintf. The detected symbols are linked into
the code but might come from libraries and not actually be called.
See 'Writing portable packages' in the 'Writing R Extensions' manual.
Flavor: r-devel-windows-x86_64