The R Project SVN
R-packages
(root)/trunk/Old/survival/Examples/test/Rtest – Rev 7500
Rev 1816 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
## Set up for the test##dyn.load("../loadmod.o")#attach("../.Data")library(survival)options(na.action="na.exclude")## Compute some answers "by hand"# For test1 and test2 databyhand1 <- function(coef, newx) {s <- exp(coef)loglik <- 2*coef - (log(3*s+3) + 2*log(s+3))u <- (6 + 3*s - s^2) / ((s+1)*(s+3))imat <- s/(s+1)^2 + 6*s/(s+3)^2x <- c(1,1,1,0,0,0)status <- c(1,0,1,1,0,1)xbar <- c(s/(s+1), s/(s+3), 0, 0)haz <- c(1/(3*s+3), 2/(s+3), 0, 1 )ties <- c(1,1,2,2,3,4)wt <- c(s,s,s,1,1,1)mart <- c(1,0,1,1,0,1) - wt* (cumsum(haz))[ties]score <- rep(6,0)for (i in 1:6) {j <- ties[i]score[i] <- -wt[i]*(cumsum((x[i]-xbar) * haz))[j]score[i] <- score[i] + wt[i]* ((x[i]-xbar)*status[i])[j]}scho <- c(1/(s+1), (3-s)/(3+s), 0)surv <- exp(-cumsum(haz)* exp(coef*newx))varhaz.g <- cumsum(c(1/(3*s+3)^2, 2/(s+3)^2, 0, 1 ))varhaz.d <- cumsum((newx-xbar) * haz)varhaz <- (varhaz.g + varhaz.d0^2/ imat) * exp(2*coef*newx)names(xbar) <- names(haz) <- 1:4names(surv) <- names(varhaz) <- 1:4list(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)}byhand2 <- function(coef, newx) {s <- exp(coef)loglik <- 4*coef - log(s+1) - log(s+2) - 3*log(3*s+2) - 2*log(3*s+1)u <- 1/(s+1) + 1/(3*s+1) + 4/(3*s+2) -( s/(s+2) +3*s/(3*s+2) + 3*s/(3*s+1))imat <- s/(s+1)^2 + 2*s/(s+2)^2 + 6*s/(3*s+2)^2 +3*s/(3*s+1)^2 + 3*s/(3*s+1)^2 + 12*s/(3*s+2)^2hazard <-c( 1/(s+1), 1/(s+2), 1/(3*s+2), 1/(3*s+1), 1/(3*s+1), 2/(3*s+2) )xbar <- c(s/(s+1), s/(s+2), 3*s/(3*s+2), 3*s/(3*s+1), 3*s/(3*s+1),3*s/(3*s+2))var.g <- cumsum(hazard*hazard /c(1,1,1,1,1,2))var.d <- cumsum( (xbar-newx)*hazard)surv <- exp(-cumsum(hazard) * exp(coef*newx))varhaz <- (var.g + var.d^2/imat)* exp(2*coef*newx)list(loglik=loglik, u=u, imat=imat, hazard=hazard,xbar=xbar, surv=surv, varhaz=varhaz, var.g=var.g, var.d=var.d)}byhand3 <- function(coef) {#Hard coded -- what is found in the Agreg.3 comments files <- as.vector(exp(coef)) #kill the names attrimat <- s/(s+1)^2 + 2*s/(s+2)^2 + 6*s/(3*s+2)^2 +3*s/(3*s+1)^2 + 3*s/(3*s+1)^2 + 12*s/(3*s+2)^2hazard <-c( 1/(s+1), 1/(s+2), 1/(3*s+2), 1/(3*s+1), 1/(3*s+1), 2/(3*s+2) )xbar <- c(s/(s+1), s/(s+2), 3*s/(3*s+2), 3*s/(3*s+1), 3*s/(3*s+1),3*s/(3*s+2))newx <- c(0,0,1,1,1,0,0, 2,2,2,2)wt <- exp(coef*newx)indx <- c(1,2,4,5,6,1,2,3,4,5,6)var.g <- hazard*hazard /c(1,1,1,1,1,2)surv <- exp(-cumsum(hazard[indx]*wt))var1 <- cumsum(var.g[indx]*wt*wt)d <- cumsum( (xbar[indx] - newx)* hazard[indx] * wt)var2 <- d^2/imatnames(surv) <- names(var1) <-names(var2) <- NULLlist(time= c(2,3,7,8,9,12,13,16,17,18,19),surv=surv, std= sqrt(var1 + var2), var.g=var1, var.d=var2, d=d,hazard=hazard, wt=wt)}# Create the simplest test data set#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))## Do the test on the simple data set#fit <-coxph(Surv(time, status) ~x, test1, method='breslow')fitfit0 <-coxph(Surv(time, status) ~x, test1, iter=0)fit0$coefcoxph(Surv(time, status) ~x, test1, iter=1, method='breslow')$coefcoxph(Surv(time, status) ~x, test1, iter=2, method='breslow')$coefcoxph(Surv(time, status) ~x, test1, iter=3, method='breslow')$coefcoxph(Surv(time, status) ~ x, test1, method='efron')coxph(Surv(time, status) ~ x, test1, method='exact')resid(fit0)resid(coxph(Surv(time, status) ~1, test1))resid(fit0, 'scor')resid(fit0, 'scho')resid(fit)resid(fit, 'scor')resid(fit, 'scho')predict(fit, type='lp')predict(fit, type='risk')predict(fit, type='expected')predict(fit, type='terms')predict(fit, type='lp', se.fit=T)predict(fit, type='risk', se.fit=T)predict(fit, type='expected', se.fit=T)predict(fit, type='terms', se.fit=T)summary(survfit(fit))summary(survfit(fit, list(x=2)))## Create the simple data for agreg#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) )## Do the test on the simple agreg data set#fit <-coxph(Surv(start, stop, event)~ x, test2, method='breslow')fitfit0 <-coxph(Surv(start, stop, event)~ x, test2, iter=0)fit0$coefcoxph(Surv(start, stop, event)~ x, test2, iter=1, method='breslow')$coefcoxph(Surv(start, stop, event)~ x, test2, iter=2, method='breslow')$coefcoxph(Surv(start, stop, event)~ x, test2, iter=3, method='breslow')$coefcoxph(Surv(start, stop, event) ~ x, test2, method='efron')coxph(Surv(start, stop, event) ~ x, test2, method='exact')resid(fit0)resid(fit0, 'scor')resid(fit0, 'scho')resid(fit)resid(fit, 'scor')resid(fit, 'scho')resid(coxph(Surv(start, stop, event)~ x, test2, iter=0, init=log(2)), 'score')sfit <-survfit(fit)sfitsummary(sfit)sfit.km <- survfit(Surv(start, stop, event)~ x, data=test2)sfit.kmsummary(sfit.km)# make a doubled data settemp <- rbind(test2, test2)temp <- data.frame(temp, x2=c(test2$x, test2$x^2),ss=c(rep(0, nrow(test2)), rep(1, nrow(test2))))fitx <- coxph(Surv(start, stop, event) ~ x2 * strata(ss), data=temp,method='breslow')sfit <- survfit(fitx, c(fitx$means[1], 0) )sfitsummary(sfit)## Even though everyone in strata 1 has x2==0, I won't get the same survival# curve above if survfit is called without forcing predicted x2 to be# zero-- otherwise I am asking for a different baseline than the# simple model did. In this particular case coef[2] is nearly zero, so# the curves are the same, but the variances differ.## This mimics 'byhand3' and the documentationfit <- coxph(Surv(start, stop, event) ~x, test2, method='breslow')tdata <- data.frame( start=c(0,20, 6,0,5),stop =c(5,23,10,5,15),event=rep(0,5),x=c(0,0,1,0,2) )temp <- survfit(fit, tdata, individual=T)temp2 <- byhand3(fit$coef)all.equal(temp$surv, temp2$surv)all.equal(temp2$std, temp$std.err)temp2## Test out the survival curve and variance, in gory detail#xdata <- data.frame(rbind(test1,test1), ss = rep(1:2, rep(nrow(test1),2)))fit <- coxph(Surv(time, status)~x + strata(ss), xdata, method='breslow')sfit <- survfit(fit, list(x=0)) #type='aalen' is default# From the hand worked notesbb <- as.vector(exp(fit$coef))realhaz <- cumsum(c(1/(3*bb+3), 2/(bb+3), 1))all.equal(sfit$surv, exp(-rep(realhaz,2)))realvar <- cumsum(c( (1/(3*bb+3))^2, 2/(bb+3)^2, 1))dd <- cumsum(c( (bb/(bb+1))*(1/(3*bb+3)),(bb/(bb+3))*(2/(bb+3)), 0*1))realvar <- realvar + dd^2 * fit$varall.equal(sfit$std, sqrt(rep(realvar,2)))summary(sfit)# Get the Kalbfleisch-Prentice survival with Greenwood variance#sfit <- survfit(fit, list(x=0), type='kalb')# second term is the solution to an equation, in the neighborhood of .6tfun <- function(alpha) (bb/(1-alpha^bb) + 1/(1-alpha) - (bb+3))^2## not in R## temp <- nlminb(0.6, tfun, lower=.1, upper=.9)temp <- optimize(tfun, lower=.1, upper=.9,tol=sqrt(.Machine$double.eps))realkm <- cumprod(c((1- bb/(3*bb+3))^(1/bb), temp$minimum, 0))all.equal(sfit$surv, rep(realkm,2))realvar <- cumsum(c( 1/((3*bb+3)*(2*bb+3)), 2/((bb+3)*2), 1))dd <- cumsum(c( (bb/(bb+1))*(1/(3*bb+3)),(bb/(bb+3))*(2/(bb+3)), 0*1))realvar <- realvar + dd^2 * fit$varall.equal(sfit$std, sqrt(rep(realvar,2)))summary(sfit)## Repeat with the Efron approximation for Cox, Efron estimates#fit <- coxph(Surv(time, status)~x + strata(ss), xdata)sfit <- survfit(fit, list(x=0))bb <- as.vector(exp(fit$coef))realhaz <- cumsum(c(1/(3*bb+3), 1/(bb+3) + 2/(bb+5), 1))all.equal(sfit$surv, exp(-rep(realhaz,2)))realvar <- cumsum(c( (1/(3*bb+3))^2, 1/(bb+3)^2 + 4/(bb+5)^2, 1))dd <- cumsum(c( (bb/(bb+1))*(1/(3*bb+3)),(bb/(bb+3))*(1/(bb+3))+ (bb/(bb+5))*(2/(bb+5)), 0*1))realvar <- realvar + dd^2 * fit$varall.equal(sfit$std, sqrt(rep(realvar,2)))summary(sfit)rm(bb, sfit, realhaz, realvar, fit, realkm, xdata)## The AML data, from Miller, "Survival Analysis", page 49.#aml <- list(time= c( 9, 13, 13, 18, 23, 28, 31, 34, 45, 48, 161,5, 5, 8, 8, 12, 16, 23, 27, 30, 33, 43, 45),status= c( 1,1,0,1,1,0,1,1,0,1,0, 1,1,1,1,1,0,1,1,1,1,1,1),x = as.factor(c(rep("Maintained", 11),rep("Nonmaintained", 12) )))aml <- data.frame(aml)## These results can be found in Miller#fit <- coxph(Surv(aml$time, aml$status) ~ aml$x, method='breslow')fitresid(fit, type='mart')resid(fit, type='score')resid(fit, type='scho')fit <- survfit(Surv(aml$time, aml$status) ~ aml$x)fitsummary(fit)survdiff(Surv(aml$time, aml$status)~ aml$x)## Test out the weighted K-M## First, equal case weights- shouldn't change the survival, but will# halve the variancetemp2 <-survfit(Surv(aml$time, aml$status), type='kaplan', weight=rep(2,23))temp <-survfit(Surv(time, status)~1, aml)temp$surv/temp2$surv(temp$std.err/temp2$std.err)^2# Risk weights-- use a null Cox modeltfit <- coxph(Surv(aml$time, aml$status) ~ offset(log(1:23)))sfit <- survfit(tfit, type='aalen')# Now compute it by hand# Ties are a challengeatime <- sort(aml$time)denom <- rev(cumsum(rev((1:23)[order(aml$time)])))denom <- denom[match(unique(atime), atime)]deaths <- tapply(aml$status, aml$time, sum)chaz <- cumsum(deaths/denom)all.equal(sfit$surv, as.vector(exp(-chaz[deaths>0])))cvar <- cumsum(deaths/denom^2)all.equal(sfit$std^2, as.vector(cvar[deaths>0]))summary(sfit)# And the Efron resultsummary(survfit(tfit))# Lots of ties, so its a good test casex1 <- coxph(Surv(time, status)~x, aml, method='efron')x1x2 <- coxph(Surv(rep(0,23),time, status) ~x, aml, method='efron')x1$coef - x2$coefrm(x1, x2, atime, denom, deaths, chaz,cvar, tfit, sfit, temp, temp2, fit)## Trivial test of stratified residuals# Make a second strata = replicate of the first, and I should get the# exact same answerstemp <- 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.frametstrat <- 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])all.equal(resid(fit1, type='score') , (resid(fit2, type='score'))[1:n])all.equal(resid(fit1, type='schoe') , (resid(fit2, type='schoe'))[1:ndead])#AG modeltemp <- as.matrix(test2)n <- nrow(temp)ndead<- sum(test2$event[!is.na(test2$event)])temp <- data.frame(rbind(temp, temp))tstrat <- rep(1:2, c(n,n))fit1 <- coxph(Surv(start, stop, event) ~x, test2)fit2 <- coxph(Surv(start, stop, event) ~x + strata(tstrat), temp)all.equal(resid(fit1) , (resid(fit2))[1:n])all.equal(resid(fit1, type='score') , (resid(fit2, type='score'))[1:n])all.equal(resid(fit1, type='schoe') , (resid(fit2, type='schoe'))[1:ndead])## A test to exercise the "infinity" check on 2 variables#test3 <- data.frame(futime=1:12, fustat=c(1,0,1,0,1,0,0,0,0,0,0,0),x1=rep(0:1,6), x2=c(rep(0,6), rep(1,6)))coxph(Surv(futime, fustat) ~ x1 + x2, test3)# Create a "counting process" version of the simplest test data set#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))## 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$coef*x), test1)fitd <- coxph(Surv(start, stop, status) ~ offset(fit$coef*x), test1b)resid(fit0)resid(fit0b, collapse=test1b$id)resid(fit)resid(fitb, collapse=test1b$id)resid(fitc)resid(fitd, collapse=test1b$id)all.equal(resid(fitc), resid(fit))resid(fit0, type='score')resid(fit0b, type='score', collapse=test1b$id)resid(fit, type='score')resid(fitb, type='score', collapse=test1b$id)resid(fit0, type='scho')resid(fit0b, type='scho', collapse=test1b$id)resid(fit, type='scho')resid(fitb, type='scho', collapse=test1b$id)summary(survfit(fit0, list(x=0)))summary(survfit(fit, list(x=0)))summary(survfit(fitb,list(x=0)))summary(survfit(fit))# 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.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))fit0 <- coxph(Surv(time, status) ~x, testw1, weights=wt,method='breslow', iter=0)fit <- coxph(Surv(time, status) ~x, testw1, weights=wt, method='breslow')fit0summary(fit)resid(fit0, type='mart')resid(fit0, type='score')resid(fit0, type='scho')fit0 <- coxph(Surv(time, status) ~x,testw1, weights=wt, iter=0)resid(fit0, 'mart')resid(coxph(Surv(time, status) ~1, testw1, weights=wt)) #Null modelresid(fit, type='mart')resid(fit, type='score')resid(fit, type='scho')fit <- coxph(Surv(time, status) ~x, testw1, weights=wt, method='efron')fitresid(fit, type='mart')resid(fit, type='score')resid(fit, type='scho')# Tests of the weighted Cox model, AG form of the data# Same solution as doweight1.s#testw2 <- data.frame(id = c( 1, 1, 2, 3, 3, 3, 4, 5, 5, 6, 7, 8, 8, 9),begin= c( 0, 5, 0, 0,10,15, 0, 0,14, 0, 0, 0,23, 0),time= c( 5,10,10,10,15,20,20,14,20,20,30,23,40,50),status= c( 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0),x= c( 2, 2, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0),wt = c( 1, 1, 2, 3, 3, 3, 4, 3, 3, 2, 1, 2, 2, 1))fit0 <- coxph(Surv(begin,time, status) ~x, testw2, weights=wt,method='breslow', iter=0)fit <- coxph(Surv(begin,time, status) ~x, testw2, weights=wt, method='breslow')fit0summary(fit)resid(fit0, type='mart', collapse=testw2$id)resid(fit0, type='score', collapse=testw2$id)resid(fit0, type='scho')resid(fit, type='mart', collapse=testw2$id)resid(fit, type='score', collapse=testw2$id)resid(fit, type='scho')fit0 <- coxph(Surv(begin, time, status) ~x,testw2, weights=wt, iter=0)resid(fit0, 'mart', collapse=testw2$id)resid(coxph(Surv(begin, time, status) ~1, testw2, weights=wt), collapse=testw2$id) #Null modelfit <- coxph(Surv(begin,time, status) ~x, testw2, weights=wt, method='efron')fitresid(fit, type='mart', collapse=testw2$id)resid(fit, type='score', collapse=testw2$id)resid(fit, type='scho')## The test data set from Turnbull, JASA 1974, 169-73.## status 0=right censored# 1=exact# 2=left censored#turnbull <- data.frame( time =c( 1,1,1, 2,2,2, 3,3,3, 4,4,4),status=c( 1,0,2, 1,0,2, 1,0,2, 1,0,2),n =c(12,3,2, 6,2,4, 2,0,2, 3,3,5))## Compute the K-M for the Turnbull data# via a slow EM calculation#emsurv <- function(time, status, wt, verbose=T) {left.cen <- (status==2)if (!any(left.cen)) stop("No left censored data!")if (!any(status==1))stop("Must have some exact death times")tempy <- Surv(time[!left.cen], status[!left.cen])ww <- wt[!left.cen]tempx <- factor(rep(1, sum(!left.cen)))tfit <- survival:::survfit.km(tempx, tempy, casewt=ww)if (verbose)cat("Iteration 0, survival=", format(round(tfit$surv[tfit$n.event>0],3)),"\n")stimes <- tfit$time[tfit$n.event>0]ltime <- time[left.cen]lwt <- wt[left.cen]tempx <- factor(rep(1, length(stimes) + sum(!left.cen)))tempy <- Surv(c(time[!left.cen], stimes),c(status[!left.cen], rep(1, length(stimes))))for (iter in 1:4) {wt2 <- stimes*0ssurv <- tfit$surv[tfit$n.event>0]sjump <- diff(c(1, ssurv))for (j in 1:(length(ltime))) {k <- sum(ltime[j]>=stimes) #index of the death timeif (k==0)stop("Left censored observation before the first death")wt2[1:k] <- wt2[1:k] + lwt[j]*sjump[1:k] /(ssurv[k]-1)}tfit <- survival:::survfit.km(tempx, tempy, casewt=c(ww, wt2))if (verbose) {cat("Iteration", iter, "survival=",format(round(tfit$surv[tfit$n.event>0],3)), "\n")cat(" weights=", format(round(wt2,3)), "\n")}}survfit(tempy ~ tempx, weights=c(ww, wt2))}temp <-emsurv(turnbull$time, turnbull$status, turnbull$n)print(summary(temp))## Read in the ovarian data#ovarian <- read.table("data.ovarian", row.names=NULL,col.names= c("futime", "fustat", "age", "resid.ds", "rx", "ecog.ps"))## Test the coxph program on the Ovarian data#xx <- order(ovarian$futime) #put data in same order as SAS green booktemp <- ovarian[xx,]attach(temp)# List the datatempsummary(survfit(Surv(futime, fustat)), censor=T)# Various modelscoxph(Surv(futime, fustat)~ age)coxph(Surv(futime, fustat)~ resid.ds)coxph(Surv(futime, fustat)~ rx)coxph(Surv(futime, fustat)~ ecog.ps)coxph(Surv(futime, fustat)~ resid.ds + rx + ecog.ps)coxph(Surv(futime, fustat)~ age + rx + ecog.ps)coxph(Surv(futime, fustat)~ age + resid.ds + ecog.ps)coxph(Surv(futime, fustat)~ age + resid.ds + rx)# Residualsfit <- coxph(Surv(futime, fustat)~ age + resid.ds + rx + ecog.ps )resid(fit)resid(fit, 'dev')resid(fit, 'scor')resid(fit, 'scho')fit <- coxph(Surv(futime, fustat) ~ age + ecog.ps + strata(rx))summary(fit)summary(survfit(fit))sfit <- survfit(fit, list(age=c(30,70), ecog.ps=c(2,3))) #two columnssfitsummary(sfit)detach(2)# Test the robust=T option of coxphfit <- coxph(Surv(futime, fustat) ~ age + ecog.ps + rx, ovarian, robust=T)rr <- resid(fit, type='dfbeta')all.equal(as.vector(t(rr) %*% rr), as.vector(fit$var))temp <- scan("data.jasa", what=list(id=0, b.mo=0, b.d=0, b.y=0,a.mo=0, a.d=0, a.y=0,t.mo=0, t.d=0, t.y=0,f.mo=0, f.d=0, f.y=0,fu.stat=0, surg=0, mismatch=0,hla.a2=0, mscore=0, reject=0))temp3 <- mdy.date(temp$b.mo, temp$b.d, temp$b.y,nineteen=TRUE)temp4 <- mdy.date(temp$a.mo, temp$a.d, temp$a.y,nineteen=TRUE)temp5 <- mdy.date(temp$t.mo, temp$t.d, temp$t.y,nineteen=TRUE)temp6 <- mdy.date(temp$f.mo, temp$f.d, temp$f.y,nineteen=TRUE)# Make sure that a particular idiocy is turned off: turning logicals into# factors!as.data.frame.logical <- as.data.frame.vectorjasa <- data.frame( birth.dt=temp3, accept.dt=temp4, tx.date=temp5,fu.date=temp6, fustat=temp$fu.stat,surgery = temp$surg, age= temp4-temp3,futime = 1+temp6 - temp4, wait.time= 1+temp5 - temp4,transplant = c(!is.na(temp5)),mismatch=temp$mismatch, hla.a2=temp$hla.a2,mscore = temp$mscore, reject=temp$reject)row.names(jasa) <- temp$id# The "1+" above causes us to match the analysis in Kalbfleisch and Prentice.# Someone accepted and transplanted on the same day is assumed to have one# day under the non-transplant risk.rm(temp, temp3, temp4, temp5, temp6)data(ratetables)expect <- survexp(futime ~ ratetable(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))# Do a Stanford heart transplant data the way that K&P do## Input - a data frame containing the raw data# Output- a data frame with multiple obs for the transplanted subjects## There are more efficient ways to do this than the "for" loop below, but# this script is much more readable than any of them.#stan1 <- function(jasa) {id <- row.names(jasa)tx <- 1*(jasa$transplant)covar <- cbind(jasa$age/365.25 -48,(jasa$accept.dt - mdy.date(10,1,1967))/365.25,jasa$surgery)n <- length(tx) #number in studyntx <- sum(tx) #number transplanted# Per paragraph 2, p138, the patient who died on the day of transplant is# treated differently, for all others deaths are assumed to occur "earlier# in the day" than transplants.special <- id ==38wait <- jasa$wait.timewait[special] <- wait[special] - .5age <- year <- surgery <- transplant <- id2 <- double(n+ntx)start <- stop <- event <- double(n+ntx)ii <- 1for (i in 1:n) {age[ii] <- covar[i,1]year[ii] <- covar[i,2]surgery[ii] <- covar[i,3]transplant[ii] <- 0id2[ii] <- id[i]if (tx[i]) { #transplanted - 2 lines if datastart[ii] <- 0stop[ii] <- wait[i]event[ii] <- 0ii <- ii+1start[ii] <- wait[i]stop[ii] <- jasa$futime[i]event[ii] <- jasa$fustat[i]age[ii] <- covar[i,1]year[ii] <- covar[i,2]surgery[ii] <- covar[i,3]transplant[ii] <- 1id2[ii] <- id[i]}else { # one line of datastart[ii] <-0stop[ii] <- jasa$futime[i]event[ii]<- jasa$fustat[i]}ii <- ii+1}data.frame(start, stop, event, age, year, surgery,transplant=factor(transplant), id=id2)}jasa1 <- stan1(jasa)attach(jasa1)# Now fit the 6 models found in Kalbfleisch and Prentice, p139#options(contrasts=c("contr.treatment", "contr.treatment"))sfit.1 <- coxph(Surv(start, stop, event)~ (age + surgery)*transplant,method='breslow')sfit.2 <- coxph(Surv(start, stop, event)~ year*transplant,method='breslow')sfit.3 <- coxph(Surv(start, stop, event)~ (age + year)*transplant,method='breslow')sfit.4 <- coxph(Surv(start, stop, event)~ (year +surgery) *transplant,method='breslow')sfit.5 <- coxph(Surv(start, stop, event)~ (age + surgery)*transplant + year ,method='breslow')sfit.6 <- coxph(Surv(start, stop, event)~ age*transplant + surgery + year,method='breslow')summary(sfit.1)sfit.2summary(sfit.3)sfit.4sfit.5sfit.6# Survival curve for the "average" subjectsummary(survfit(sfit.1))# Survival curve for a subject of age 50, with prior surgery, tx at 6 monthsdata <- data.frame(start=c(0,183), stop=c(183,3*365), event=c(1,1),age=c(50,50), surgery=c(1,1), transplant=c(0,1))summary(survfit(sfit.1, data, individual=T))# These should all give the same answerj.age <- jasa$age/365.25 -48fit1 <- coxph(Surv(futime, fustat) ~ j.age, data=jasa)fit2 <- coxph(Surv(futime, fustat) ~ j.age, jasa, init=fit1$coef, iter=0)fit3 <- coxph(Surv(start, stop, event) ~ age)fit4 <- coxph(Surv(start, stop, event) ~ offset((age-fit3$means)*fit1$coef))s1 <- survfit(fit1, fit3$means)s2 <- survfit(fit2, fit3$means)s3 <- survfit(fit3)s4 <- survfit(fit4)all.equal(s1$surv, s2$surv)all.equal(s1$surv, s3$surv)all.equal(s1$surv, s4$surv)# Still the same answer, fit multiple strata at once# Strata 1 has independent coefs of strata 2, so putting in# the other data should not affect itll <- length(start)ss <- rep(0:1, c(ll,ll))fit <- coxph(Surv(rep(start,2), rep(stop,2), rep(event,2)) ~rep(age,2)*strata(ss) + I(rep(age,2)^2*ss) )fit$coef[1] - fit3$coefs4 <- survfit(fit, c(fit3$means, 0,0))all.equal(s4$surv[1:(s4$strata[1])], s3$surv)detach("jasa1")# A subset of some local data on lung cancer patients. A useful internal test# because it has multiple strata and lots of missing values.## inst = enrolling institution# sex 1=male 2=female# ph.ecog physician's estimate of the ECOG performace score. 0=fully active,# 4=bedridden# ph.karno physician's estimate of the Karnofsky score, a competitor to the# ECOG ps.# pat.karno patient's assesment of his/her Karnofsky score# meal.cal # calories consumed at meals (exclude beverages and snacks)# wt.loss weight loss in the last 6 months# 12/98: this is now part of the Splus distribution, but they call it# "lung"#cancer <- lung## In R it's data(cancer) or data(lung)data(cancer)# Test out all of the routines on a more complex data set#temp <- survfit(Surv(time, status) ~ ph.ecog, cancer)summary(temp, times=c(30*1:11, 365*1:3))print(temp[2:3])temp <- survfit(Surv(time, status), cancer, type='fleming',conf.int=.9, conf.type='log-log', error='tsiatis')summary(temp, times=30 *1:5)temp <- survdiff(Surv(time, status) ~ inst, cancer, rho=.5)print(temp, digits=6)temp <- coxph(Surv(time, status) ~ ph.ecog + ph.karno + pat.karno + wt.loss+ sex + age + meal.cal + strata(inst), cancer)summary(temp)cox.zph(temp)cox.zph(temp, transform='identity')coxph(Surv(rep(0,length(time)), time, status) ~ ph.ecog + ph.karno + pat.karno+ wt.loss + sex + age + meal.cal + strata(inst), cancer)bladder <- read.table('data.bladder4',col.names=c('id', 'rx', 'number', 'size', 'start', 'stop','event', 'enum'))bladder2 <- bladder[bladder$start< bladder$stop, ]bladder$start <- NULLbladder <- bladder[bladder$enum<5,]## Fit the models found in Wei et. al.##options(contrasts='contr.treatment')wfit <- coxph(Surv(stop, event) ~ (rx + size + number)* strata(enum) +cluster(id), bladder, method='breslow')wfit# Check the rx coefs versus Wei, et al, JASA 1989rx <- c(1,4,5,6) # the treatment coefs abovecmat <- diag(4); cmat[1,] <- 1; #contrast matrixwfit$coef[rx] %*% cmat # the coefs in their paper (table 5)t(cmat) %*% wfit$var[rx,rx] %*% cmat # var matrix (eqn 3.2)# Anderson-Gill fitfita <- coxph(Surv(start, stop, event) ~ rx + size + number + cluster(id),bladder2, method='breslow')fita# Prentice fits. Their model 1 a and b are the samefit1p <- coxph(Surv(stop, event) ~ rx + size + number, bladder2,subset=(enum==1), method='breslow')fit2pa <- coxph(Surv(stop, event) ~ rx + size + number, bladder2,subset=(enum==2), method='breslow')fit2pb <- coxph(Surv(stop-start, event) ~ rx + size + number, bladder2,subset=(enum==2), method='breslow')fit3pa <- coxph(Surv(stop, event) ~ rx + size + number, bladder2,subset=(enum==3), method='breslow')#and etc.fit1pfit2pafit2pbfit3pa## A test of NULL models#fit1 <- coxph(Surv(stop, event) ~ rx + strata(number), bladder, iter=0)fit2 <- coxph(Surv(stop, event) ~ strata(number), bladder)all.equal(fit1$loglik[2], fit2$loglik)all.equal(fit1$resid, fit2$resid)fit1 <- coxph(Surv(start, stop, event) ~ rx + strata(number), bladder2, iter=0)fit2 <- coxph(Surv(start, stop, event) ~ strata(number), bladder2)all.equal(fit1$loglik[2], fit2$loglik[2])all.equal(fit1$resid, fit2$resid)# Tests of expected survival## Simple case 1: a single male subject, born 6/6/36 and entered on study 6/6/55.## Compute the 1, 5, 10 and 12 year expected survivaltemp1 <- mdy.date(6,6,36)temp2 <- mdy.date(6,6,55)##exp1 <- survexp(~ratetable(year=temp2, age=(temp2-temp1), sex=1),## ratetable=survexp.uswhite,times=c(366, 1827, 3653, 4383))exp1 <- survexp(~ratetable(year=temp2, age=(temp2-temp1), sex=1,race="white"),ratetable=survexp.usr,times=c(366, 1827, 3653, 4383))# Well, almost easy. The uswhite table assumes that someone age 19 will have# seen 5 leap years-- but this lad has only seen 4. Thus the routine sees# him as 1 day shy of his birthday.# (The age 19 category starts at 6940, but he is 6939 days old at entry).# Thus his passage through the table is a bit more complicated# First epoch: 1 day at age 18, 1954 6/6/55# 365 days at age 19, 1955 6/7/55 - 6/5/56# Second 365 days at age 20, 1956 6/6/56 - 6/5/57# 365 days at age 21, 1957 6/6/57 - 6/5/58# 366 days at age 22, 1958 6/6/58 - 6/6/59# 365 days at age 23, 1959 6/7/59 - 6/5/60# Third 365 days at age 24, 1960 6/6/60 - 6/5/61# 365 days at age 25, 1961 6/6/61 - 6/5/62# 366 days at age 26, 1962 6/6/62 - 6/6/63# 365 days at age 27, 1963 6/7/63 - 6/5/64# 365 days at age 28, 1964 6/6/64 - 6/5/64# Fourth 365 days at age 29, 1965 6/6/65 - 6/5/66# 365 days at age 30, 1966 6/6/66 - 6/5/67# remember, a first subscript of "1" is for age 0##xx <- survexp.uswhite[,1,]## Edited to account for change to survexp.usr, which has more categories## before age 1.xx <- survexp.usr[,1,"white",]check <- c( (.6*xx["18","1950"] + .4*xx["18","1960"]) +365* (.5*xx["19","1950"] + .5*xx["19","1960"]) ,365* (.4*xx["20","1950"] + .6*xx["20","1960"]) +365* (.3*xx["21","1950"] + .7*xx["21","1960"]) +366* (.2*xx["22","1950"] + .8*xx["23","1960"]) +365* (.1*xx["23","1950"] + .9*xx["24","1960"]) ,365* ( xx["24","1960"] ) +365* (.9*xx["25","1960"] + .1*xx["25","1970"]) +366* (.8*xx["26","1960"] + .2*xx["26","1970"]) +365* (.7*xx["27","1960"] + .3*xx["27","1970"]) +365* (.6*xx["28","1960"] + .4*xx["28","1970"]) ,365* (.5*xx["29","1960"] + .5*xx["29","1970"]) +365* (.4*xx["30","1960"] + .6*xx["30","1970"]) )print(exp1$surv)print(exp(-cumsum(check))) ###?FIXME# This does not pass the "all.equal" test. Because of leap year (again),# the internal S code does not do exactly what is stated above. US# rate tables are special: the entry for age 20 in 1950 is the probability# that someone who becomes 20 years old in 1950 will reach his 21st birthday.# In order to fit this into the general "cutpoints" calculations of# person-years, dates are adjusted so that everyone appears to have a# birthday on Jan 1. But because of leap years, a birthday then moves to# Dec 31 sometimes -- this happens in each of the 366 day intervals above,# e.g., 1 day at age 22 with 1957 rates, 365 at age 22 with 1958 rates.# Upshot: they only agree to 6 decimals. Close enough for me!rm(temp1, temp2, exp1, check, xx)# Tests of expected survival## Simple case 2: make sure that the averages are correct, for Ederer method## Compute the 1, 5, 10 and 12 year expected survivaltemp1 <- mdy.date(1:6,6:11,36:41)temp2 <- mdy.date(6:1,11:6,55:50)age <- temp2 - temp1exp1 <- survexp(~ratetable(year=temp2, age=(temp2-temp1), sex=1),times=c(366, 1827, 3653, 4383))exp2 <- survexp(~ratetable(year=temp2, age=(temp2-temp1), sex=1) + I(1:6),times=c(366, 1827, 3653, 4383))print(all.equal(exp1$surv, apply(exp2$surv, 1, mean)))## Check that adding more time points doesn't change things#exp3 <- survexp(~ratetable(year=temp2, age=(temp2-temp1), sex=1),times=sort(c(366, 1827, 3653, 4383, 30*(1:100))))exp1$survexp3$surv[match(exp1$time, exp3$time, nomatch=0)]## Now test Hakulinen's method, assuming an analysis date of 3/1/57#futime <- mdy.date(3,1,57) - temp2xtime <- sort(c(futime, 30, 60, 185, 365))exp1 <- survexp(futime ~ ratetable(year=temp2, age=(temp2-temp1), sex=1),times=xtime, conditional=F)exp2 <- survexp(~ratetable(year=temp2, age=(temp2-temp1), sex=1) + I(1:6),times=futime)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)]cumprod(con) #should be equal## Now for the conditional method#exp1 <- survexp(futime ~ ratetable(year=temp2, age=(temp2-temp1), sex=1),times=xtime, conditional=T)cond <- exp2$survfor (i in 6:2) cond[i,] <- (cond[i,]/cond[i-1,]) #conditional survivalfor (i in 1:6) con[i] <- exp(mean(log(cond[i, i:6])))all.equal(exp1$surv[match(futime, xtime)], cumprod(con))cumprod(con)rm(con, cond, exp1, exp2, wt, temp1, temp2, age)rm(exp3, futime, xtime)## Test out expected survival, when the parent pop is another Cox model#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 ~ ratetable(x=x), dummy, ratetable=fit, cohort=F)## Now, compare to the true answer, which is known to us#ss <- exp(fit$coef)haz <- c( 1/(3*ss+3), 2/(ss+3), 1) #truth at time 0,1,2,4+chaz <- cumsum(c(0,haz))chaz2 <- chaz[c(1,2,2,3,3,3,3,4,4)]risk <- exp(fit$coef*dummy$x)efit2 <- exp(-risk*chaz2)all.equal(as.vector(efit), as.vector(efit2)) #ignore mismatched name attrib## Now test the direct-adjusted curve (Ederer)#efit <- survexp( ~ ratetable(x=x), dummy, ratetable=fit, se=F)direct <- survfit(fit, newdata=dummy)$survchaz <- chaz[-1] #drop time 0d2 <- exp(outer(-chaz, risk))all.equal(as.vector(direct), as.vector(d2)) #this tests survfitall.equal(as.vector(efit$surv), as.vector(apply(direct,1,mean))) #direct## Now test out the Hakulinen method (Bonsel's method)# By construction, we have a large correlation between x and censoring## In theory, hak1 and hak2 would be the same. In practice, like a KM and# F-H, they differ when n is small.#efit <- survexp( time ~ ratetable(x=x), dummy, ratetable=fit, se=F)surv <- wt <- rep(1,9)tt <- c(1,2,4)hak1 <- hak2 <- NULLfor (i in 1:3) {wt[dummy$time < tt[i]] <- 0hak1 <- c(hak1, exp(-sum(haz[i]*risk*surv*wt)/sum(surv*wt)))hak2 <- c(hak2, sum(exp(-haz[i]*risk)*surv*wt)/sum(surv*wt))surv <- surv * exp(-haz[i]*risk)}all.equal(as.vector(efit$surv), as.vector(cumprod(hak2)))## Now do the conditional estimate#efit <- survexp( time ~ ratetable(x=x), dummy, ratetable=fit, se=F,conditional=T)wt <- rep(1,9)cond <- NULLfor (i in 1:3) {wt[dummy$time < tt[i]] <- 0cond <- c(cond, exp(-sum(haz[i]*risk*wt)/sum(wt)))}all.equal(as.vector(efit$surv), as.vector(cumprod(cond)))rm(wt, cond, efit, tt, surv, hak1, hak2)rm(fit, dummy, ss, efit2, chaz, chaz2, risk)rm(d2, direct)## Run a test that can be verified using SAS's LIFEREG#fit1w <- survreg(Surv(time, status) ~x, test1, dist='weibull')fit1wsummary(fit1w)fit1e <- survreg(Surv(time, status) ~x, test1, dist='exp')fit1esummary(fit1e)fit1l <- survreg(Surv(time, status) ~x, test1, dist='loglogistic')fit1lsummary(fit1l)fit1g <- survreg(Surv(time, status) ~x, test1, dist='lognormal')summary(fit1g)## Do a test with the ovarian data#fitfw <- survreg.old(Surv(futime, fustat) ~ age + ecog.ps, ovarian,link='log', dist='extreme')fitfwfitfl <- survreg.old(Surv(futime, fustat) ~ age + ecog.ps, ovarian,link='log', dist='logistic')fitflflem2 <- scan("fleming.data2", what=list(ltime=0, rtime=0))flsurv<- Surv(flem2$ltime, flem2$rtime, type='interval2')fitfw2 <- survreg(flsurv ~ age + ecog.ps, ovarian, dist='weibull')summary(fitfw2)fitfl2 <- survreg(flsurv ~ age + ecog.ps, ovarian, dist='loglogistic')summary(fitfl2)fitfg2 <- survreg(flsurv ~ age + ecog.ps, ovarian, dist='lognormal')summary(fitfg2)## A simple test of the singular=ok option#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))temp <- rep(0:3, rep(7,4))stest <- data.frame(start = 10*temp,stop = 10*temp + test1$time,status = rep(test1$status,4),x = c(test1$x+ 1:7, rep(test1$x,3)),epoch = rep(1:4, rep(7,4)))fit1 <- coxph(Surv(start, stop, status) ~ x * factor(epoch), stest)fit1$coef # elements 2:4 should be NA## Test some more features of surv.diff## First, what happens when one group is a dummy### The AML data, with a third group of early censorings "tacked on"#aml3 <- list(time= c( 9, 13, 13, 18, 23, 28, 31, 34, 45, 48, 161,5, 5, 8, 8, 12, 16, 23, 27, 30, 33, 43, 45,1, 2, 2, 3, 3, 3, 4),status= c( 1,1,0,1,1,0,1,1,0,1,0, 1,1,1,1,1,0,1,1,1,1,1,1,0,0,0,0,0,0,0),x = as.factor(c(rep("Maintained", 11),rep("Nonmaintained", 12), rep("Dummy",7) )))aml3 <- data.frame(aml3)# These should give the same result (chisq, df), but the second has an# extra groupsurvdiff(Surv(time, status) ~x, aml)survdiff(Surv(time, status) ~x, aml3)## Now a test of the stratified log-rank# There are no tied times within institution, so the coxph program# can be used to give a complete test#fit <- survdiff(Surv(time, status) ~ pat.karno + strata(inst), cancer)#options(contrasts='contr.treatment')cfit <- coxph(Surv(time, status) ~ factor(pat.karno) + strata(inst),cancer, iter=0)tdata <- na.omit(cancer[,c('time', 'status', 'pat.karno', 'inst')])temp1 <- tapply(tdata$status-1, list(tdata$pat.karno, tdata$inst), sum)temp1 <- ifelse(is.na(temp1), 0, temp1)temp2 <- tapply(cfit$resid, list(tdata$pat.karno, tdata$inst), sum)temp2 <- ifelse(is.na(temp2), 0, temp2)temp2 <- temp1 - temp2#Now temp1=observed, temp2=expectedall.equal(c(temp1), c(fit$obs))all.equal(c(temp2), c(fit$exp))all.equal(fit$var[-1,-1], solve(cfit$var))rm(tdata, temp1, temp2)## Simple case: a single male subject, born 6/6/36 and entered on study 6/6/55.#temp1 <- mdy.date(6,6,36)temp2 <- mdy.date(6,6,55)# Now compare the results from person-years#temp.age <- tcut(temp2-temp1, floor(c(-1, (18:31 * 365.24))),labels=c('0-18', paste(18:30, 19:31, sep='-')))temp.yr <- tcut(temp2, mdy.date(1,1,1954:1965), labels=1954:1964)temp.time <- 3700 #total days of fupy1 <- pyears(temp.time ~ temp.age + temp.yr, scale=1) #output in days# The subject should appear in 20 cells# 6/6/55 - 12/31/55, 209 days, age 19-20, 1955# 1/1/56 - 6/ 4/56, 156 days, age 19-20, 1956# 6/5/56 - 12/31/56, 210 days, age 20-21, 1956 (a leap year, and his# birthday computes one day earlier)# 1/1/57 - 6/ 5/57, 156 days, age 20-21, 1957# 6/6/57 - 12/31/57, 209 days, age 21-22, 1957# and etc# with 203 days "off table", ie, beyond the last cell of the table## It is a nuisance, but tcut follows 'cut' in that we give the ENDS of# the intervals, whereas the survival tables use the starts of intervals.# Thus this breakdown does not match that in doexpect.s#xx <- matrix(0, nrow=14, ncol=11)xx[cbind(3:11, 3:11)] <- 156xx[cbind(3:12, 2:11)] <- c(209, 210, rep(c(209, 209, 209, 210),2))dimnames(xx) <- list(c('0-18', paste(18:30, 19:31, sep='-')), 1954:1964)all.equal(xx, py1$pyears)all.equal(203, py1$offtable)all.equal(1*(xx>0), py1$n)## Now with expecteds#py2 <- pyears(temp.time ~ temp.age + temp.yr+ ratetable(age=temp2-temp1, year=temp2, sex=1,race="white"),scale=1, ratetable=survexp.usr ) #output in daysall.equal(xx, py2$pyears)all.equal(203, py2$offtable)all.equal(1*(xx>0), py2$n)py3 <- pyears(temp.time ~ temp.age + temp.yr+ ratetable(age=temp2-temp1, year=temp2, sex=1,race="white"),scale=1, ratetable=survexp.usr , expect='pyears')all.equal(py2$n, py3$n)all.equal(py2$pyear, py3$pyear)all.equal(py3$n, 1*(py3$expect>0))# Now, compute the py3 result "by hand". Since there is only one person# it can be derived from py2.#xx1 <- py2$expect[py2$n>0] # the hazard over each intervalcumhaz <- cumsum(c(0, xx1[-length(xx1)])) # the cumulative hazardxx2 <- py3$expect[py3$n>0] # the expected number of person daysxx3 <- py3$pyears[py3$n>0] # the potential number of person days# This is the integral of the curve "exp(-haz *t)" over the intervalintegral <- xx3 * exp(-cumhaz)* (1- exp(-xx1))/ xx1# They might not be exactly equal, since the C code tracks changes in the# rate tables that occur -within- an interval. So try for 7 digitsall.equal(round(integral,4)/10000, round(xx2,4)/10000)# Cut off the bottom of the table, instead of the sidetemp.age <- tcut(temp2-temp1, floor(c(-1, (18:27 * 365.24))),labels=c('0-18', paste(18:26, 19:27, sep='-')))py4 <- eval(py3$call)all.equal(py4$pyear, py3$pyear[1:10,])all.equal(py4$expect, py3$expect[1:10,])rm(xx1, xx2, xx3, cumhaz, temp1, temp2)rm(temp.age, temp.yr, temp.time, xx)## Test out subscripting in the case of a coxph survival curve#fit <- coxph(Surv(time, status) ~ age + sex + meal.cal + strata(ph.ecog),data=cancer)surv1 <- survfit(fit)temp <- surv1[2:3]which <- cumsum(surv1$strata)zed <- (which[1]+1):(which[3])all.equal(surv1$surv[zed], temp$surv)all.equal(surv1$time[zed], temp$time)## Now a result with a matrix of survival curves#dummy <- data.frame(age=c(30,40,60), sex=c(1,2,2), meal.cal=c(500, 1000, 1500))surv2 <- survfit(fit, newdata=dummy)zed <- 1:which[1]all.equal(surv2$surv[zed,1], surv2[1,1]$surv)all.equal(surv2$surv[zed,2], surv2[1,2]$surv)all.equal(surv2$surv[zed,3], surv2[1,3]$surv)all.equal(surv2$surv[zed, ], surv2[1,1:3]$surv)all.equal(surv2$surv[zed], (surv2[1]$surv)[,1])all.equal(surv2$surv[zed, ], surv2[1, ]$surv)