File: C:/Program Files/R/R-4.4.0/library/survival/tests/survfit1.Rout.save
R Under development (unstable) (2021-06-11 r80486) -- "Unsuffered Consequences"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> #
> # Check out the survfit routine on the simple AML data set.
> # The leverage validation makes use of the fact that when all
> # weights are 1 and there is 1 obs per subject, the IJ variance is
> # equal to the Greenwood.
> # There are 8 choices in the C code: Nelson-Aalen or Fleming-Harrington
> # estimate of cumulative hazard, KM or exp(cumhaz) estimate of survival,
> # regular or robust variance. This tries to exercise them all.
>
> library(survival)
> aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)
>
> set.seed(1953) # used only to reorder the data
> adata <- aml
> adata$id <- sample(LETTERS, nrow(aml)) # labels are not in time or data order
> adata <- adata[sample(1:nrow(aml), nrow(aml)),] # data is unordered
> adata$wt <- sample((2:30)/10, nrow(aml)) # non-integer weights
>
> byhand <- function(time, status, weights, id) {
+ # for a single curve
+ utime <- sort(unique(time))
+ ntime <- length(utime)
+ n <- length(time)
+ if (missing(weights)) weights <- rep(1.0, n)
+ if (missing(id)) id <- seq(along=time)
+
+ uid <- sort(unique(id))
+ nid <- length(uid)
+ id <- match(id, uid) # change it to 1:nid
+
+ n.risk <- n.event <- surv <- cumhaz <- double(ntime)
+ KM <- 1; nelson <-0;
+ kvar <- 0; hvar<-0;
+
+ U <- matrix(0, nid, 2) # the two robust influence estimates
+ V <- matrix(0, ntime, 4) # variances
+ usave <- array(0., dim=c(nid, 2, ntime))
+ estimate <- matrix(0, ntime, 2)
+
+ for (i in 1:ntime) {
+ atrisk <- (time >= utime[i])
+ n.risk[i] <- sum(weights[atrisk])
+ deaths <- (time==utime[i] & status==1)
+ n.event[i] <- sum(weights[deaths])
+
+ haz <- n.event[i]/n.risk[i]
+ dhaz <- (ifelse(deaths,1,0) - ifelse(atrisk, haz, 0))/n.risk[i]
+ U[,1] <- U[,1]*(1-haz) - KM*tapply(dhaz*weights, id, sum)
+ V[i,1] <- sum(U[,1]^2)
+
+ U[,2] <- U[,2] + tapply(dhaz* weights, id, sum) #result in 'id' order
+ V[i,2] <- sum(U[,2]^2)
+ usave[,,i] <- U
+
+ if (n.event[i] >0 ) {
+ KM <- KM*(1-haz)
+ nelson <- nelson + haz
+ kvar <- kvar + n.event[i]/(n.risk[i] * (n.risk[i] - n.event[i]))
+ hvar <- hvar + n.event[i]/(n.risk[i]^2)
+ }
+
+ V[i,3] <- kvar # var of log(S)
+ V[i,4] <- hvar
+ estimate[i,] <- c(KM, nelson)
+ }
+ dimnames(usave) <- list(uid, c("KM", "chaz"), utime)
+ list(time=utime, n.risk=n.risk, n.event=n.event, estimate=estimate,
+ std = sqrt(V), influence=usave)
+ }
>
> # the byhand function can only handle one group at a time
> true1a <- with(subset(adata, x=="Maintained"), byhand(time, status, id=id))
> true1b <- with(subset(adata, x!="Maintained"), byhand(time, status, id=id))
>
> # The Greenwood and IJ estimates agree, except for a last point with
> # variance of zero. These next few lines verify the byhand() function
> aeq(true1a$std[,1], true1a$estimate[,1]*true1a$std[,3])
[1] TRUE
> aeq(true1b$std[1:9,1], true1b$estimate[1:9,1]*true1b$std[1:9,3])
[1] TRUE
> aeq(true1b$std[10,1], 0) # variance of zero for jackknife
[1] TRUE
> !is.finite(true1b$std[10,3]) # Inf for Greenwood
[1] TRUE
> temp <- with(subset(adata, x=="Maintained"), byhand(time, status, id=id,
+ weights=rep(3,11)))
> aeq(temp$std[,1:2], true1a$std[,1:2]) # IJ estimates should be invariant
[1] TRUE
>
> # fit1 uses the standard formulas: NA hazard, KM survival
> fit1 <- survfit(Surv(time, status) ~ x, data=adata)
> aeq(fit1$surv, c(true1a$estimate[,1], true1b$estimate[,1]))
[1] TRUE
> aeq(fit1$cumhaz, c(true1a$estimate[,2], true1b$estimate[,2]))
[1] TRUE
> aeq(fit1$std.err, c(true1a$std[,3], true1b$std[,3]))
[1] TRUE
> aeq(fit1$std.chaz, c(true1a$std[,4], true1b$std[,4]))
[1] TRUE
> aeq(fit1$n.risk, c(true1a$n.risk, true1b$n.risk))
[1] TRUE
> aeq(fit1$n.event, c(true1a$n.event, true1b$n.event))
[1] TRUE
> fit1$logse # logse should be TRUE
[1] TRUE
>
> # fit2 will use the IJ method
> fit2 <- survfit(Surv(time, status) ~ x, data=adata, id=id, robust=TRUE)
> aeq(fit2$surv, c(true1a$estimate[,1], true1b$estimate[,1]))
[1] TRUE
> aeq(fit2$cumhaz, c(true1a$estimate[,2], true1b$estimate[,2]))
[1] TRUE
> aeq(fit2$std.err, c(true1a$std[,1], true1b$std[,1]))
[1] TRUE
> aeq(fit2$std.chaz, c(true1a$std[,2], true1b$std[,2]))
[1] TRUE
> aeq(fit2$n.risk, c(true1a$n.risk, true1b$n.risk))
[1] TRUE
> aeq(fit2$n.event, c(true1a$n.event, true1b$n.event))
[1] TRUE
> !fit2$logse # logse should be FALSE
[1] TRUE
>
> # look at the leverage values
> fit3 <- survfit(Surv(time, status) ~ x, data=adata, id=id, influence=3)
> aeq(fit3$influence.surv[[1]], true1a$influence[,1,])
[1] TRUE
> aeq(fit3$influence.surv[[2]], true1b$influence[,1,])
[1] TRUE
> aeq(fit3$influence.chaz[[1]], true1a$influence[,2,])
[1] TRUE
> aeq(fit3$influence.chaz[[2]], true1b$influence[,2,])
[1] TRUE
>
> # compute the influence by brute force
> tdata <- subset(adata, x != "Maintained")
> tdata <- tdata[order(tdata$id),] # easier to compare if it's in order
> eps <- 1e-8
> imat1 <- imat2 <- matrix(0., 12, 10)
> t1 <- survfit(Surv(time, status) ~x, data=tdata)
> for (i in 1:12) {
+ wtemp <- rep(1.0, 12)
+ wtemp[i] <- 1 + eps
+ tfit <-survfit(Surv(time, status) ~x, data=tdata, weights=wtemp)
+ imat2[i,] <- (tfit$cumhaz - t1$cumhaz)/eps
+ imat1[i,] <- (tfit$surv - t1$surv)/eps
+ }
> aeq(imat1, true1b$influence[,1,], tol= sqrt(eps))
[1] TRUE
> aeq(imat2, true1b$influence[,2,], tol= sqrt(eps))
[1] TRUE
>
> # Repeat using the Nelson-Aalen hazard and exp(NA) for survival
> fit1 <- survfit(Surv(time, status) ~ x, adata, stype=2)
> aeq(fit1$surv, exp(-c(true1a$estimate[,2], true1b$estimate[,2])))
[1] TRUE
> aeq(fit1$cumhaz, c(true1a$estimate[,2], true1b$estimate[,2]))
[1] TRUE
> aeq(fit1$std.err, c(true1a$std[,4], true1b$std[,4]))
[1] TRUE
> aeq(fit1$std.chaz, c(true1a$std[,4], true1b$std[,4]))
[1] TRUE
> aeq(fit1$n.risk, c(true1a$n.risk, true1b$n.risk))
[1] TRUE
>
> # Nelson-Aalen + exp() surv, along with IJ variance
> fit2 <- survfit(Surv(time, status) ~ x, data=adata, id=id, stype=2,
+ influence=3)
> aeq(fit2$surv, exp(-c(true1a$estimate[,2], true1b$estimate[,2])))
[1] TRUE
> aeq(fit2$cumhaz, c(true1a$estimate[,2], true1b$estimate[,2]))
[1] TRUE
> aeq(fit2$std.err, c(true1a$std[,2], true1b$std[,2]))
[1] TRUE
> aeq(fit2$std.chaz, c(true1a$std[,2], true1b$std[,2]))
[1] TRUE
> aeq(fit2$n.risk, c(true1a$n.risk, true1b$n.risk))
[1] TRUE
> aeq(fit2$influence.chaz[[1]], true1a$influence[,2,])
[1] TRUE
> aeq(fit2$influence.chaz[[2]], true1b$influence[,2,])
[1] TRUE
> aeq(fit2$influence.surv[[2]], -true1b$influence[,2,]%*% diag(fit2[2]$surv))
[1] TRUE
>
> # Cumulative hazard is the same for fit1 and fit2
> all.equal(fit2$influence.chaz, fit3$influence.chaz)
[1] TRUE
>
> # Weighted fits
> true2a <- with(subset(adata, x=="Maintained"), byhand(time, status, id=id,
+ weights= wt))
> true2b <- with(subset(adata, x!="Maintained"), byhand(time, status, id=id,
+ weights=wt))
> fit3 <- survfit(Surv(time, status) ~ x, data=adata, id=id, weights=wt,
+ influence=TRUE)
> aeq(fit3$influence.surv[[1]], true2a$influence[,1,])
[1] TRUE
> aeq(fit3$influence.surv[[2]], true2b$influence[,1,])
[1] TRUE
> aeq(fit3$influence.chaz[[1]], true2a$influence[,2,])
[1] TRUE
> aeq(fit3$influence.chaz[[2]], true2b$influence[,2,])
[1] TRUE
> aeq(fit3$surv, c(true2a$estimate[,1], true2b$estimate[,1]))
[1] TRUE
> aeq(fit3$cumhaz, c(true2a$estimate[,2], true2b$estimate[,2]))
[1] TRUE
> aeq(fit3$std.err, c(true2a$std[,1], true2b$std[,1]))
[1] TRUE
> aeq(fit3$std.chaz, c(true2a$std[,2], true2b$std[,2]))
[1] TRUE
> aeq(fit3$n.risk, c(true2a$n.risk, true2b$n.risk))
[1] TRUE
> aeq(fit3$n.event, c(true2a$n.event, true2b$n.event))
[1] TRUE
>
> # Different survival, same hazard
> fit3b <- survfit(Surv(time, status) ~ x, data=adata, id=id, weights=wt,
+ influence=2, stype=2)
> temp <- c("n", "time", "cumhaz", "std.chaz", "influence.chaz", "n.risk",
+ "n.event")
> aeq(unclass(fit3b)[temp], unclass(fit3)[temp]) # unclass avoids [.survfit
[1] TRUE
> aeq(fit3b$surv, exp(-c(true2a$estimate[,2], true2b$estimate[,2])))
[1] TRUE
> aeq(fit3b$std.err, fit3b$std.chaz)
[1] TRUE
> aeq(fit3b$logse, FALSE)
[1] TRUE
> aeq(fit3b$n.risk, c(true2a$n.risk, true2b$n.risk))
[1] TRUE
> aeq(fit3b$n.event, c(true2a$n.event, true2b$n.event))
[1] TRUE
>
> # The grouped jackknife
> group <- rep("", nrow(adata))
> temp <- table(adata$x)
> group[adata$x == "Maintained"] <- rep(letters[1:4], length=temp[1])
> group[adata$x != "Maintained"] <- rep(letters[4:7], length=temp[2])
> adata$group <- group
> fit4 <- survfit(Surv(time, status) ~ x, data=adata, id=id, weights=wt,
+ influence=TRUE, cluster=group)
> g1 <- adata$group[match(rownames(true2a$influence[,1,]), adata$id)]
> g2 <- adata$group[match(rownames(true2b$influence[,1,]), adata$id)]
> aeq(fit4$influence.surv[[1]], rowsum(true2a$influence[,1,], g1))
[1] TRUE
> aeq(fit4$influence.surv[[2]], rowsum(true2b$influence[,1,], g2))
[1] TRUE
> aeq(fit4$influence.chaz[[1]], rowsum(true2a$influence[,2,], g1))
[1] TRUE
> aeq(fit4$influence.chaz[[2]], rowsum(true2b$influence[,2,], g2))
[1] TRUE
>
> aeq(c(colSums(fit4$influence.surv[[1]]^2), colSums(fit4$influence.surv[[2]]^2)),
+ fit4$std.err^2)
[1] TRUE
> aeq(c(colSums(fit4$influence.chaz[[1]]^2), colSums(fit4$influence.chaz[[2]]^2)),
+ fit4$std.chaz^2)
[1] TRUE
>
> # The Fleming-Harrington is a more complex formula. Start with weights of
> # 1.
> fit5 <- survfit(Surv(time, status) ~x, adata, ctype=2)
> nrisk <- c(11,10,8,7, 5,4,2, 12, 11, 10, 9, 8, 6:1)
> chaz <- c(cumsum(1/nrisk[1:7])[c(1:4,4, 5,6,6,7,7)],
+ cumsum(1/nrisk[8:18])[c(2,4,5,5,6:11)])
> aeq(fit5$cumhaz, chaz)
[1] TRUE
> aeq(fit5$std.chaz, sqrt(c(cumsum(1/nrisk[1:7]^2)[c(1:4,4, 5,6,6,7,7)],
+ cumsum(1/nrisk[8:18]^2)[c(2,4,5,5,6:11)])))
[1] TRUE
>
> # We can compute the FH using a fake data set where each tie is spread out
> # over a set of fake times.
> #
> fh <- function(time, status, weights, id) {
+ counts <- table(time, status)
+ utime <- sort(unique(time))
+ tied <- counts[,2] > 1
+
+ if (missing(weights)) weights <- rep(1.0, length(time))
+ if (missing(id)) id <- 1:length(time)
+
+ # build the expanded data set
+ delta <- min(diff(utime))/(2*max(counts[,2]))
+ efun <- function(x) {
+ who <- which(time==x & status==1)
+ ntie <- length(who)
+ data.frame(time = rep(x - (1:ntie -1)*delta, each=ntie),
+ id = rep(id[who], ntie),
+ status = rep(1, ntie^2),
+ weight = rep(weights[who]/ntie, ntie),
+ stringsAsFactors=FALSE
+ )
+ }
+
+ temp <- do.call(rbind, lapply(utime[tied], efun))
+ notie <- (status==0 | !(time %in% utime[tied]))
+
+ bfit <- byhand(time = c(time[notie], temp$time),
+ status = c(status[notie], temp$status),
+ id = c(id[notie], temp$id),
+ weights = c(weights[notie], temp$weight)
+ )
+ keep <- match(utime, bfit$time) # the real time points
+
+ list(time=bfit$time[keep],
+ n.risk=bfit$n.risk[keep - pmax(0, counts[,2]-1)],
+ n.event = bfit$n.event[keep]* counts[,2],
+ estimate=bfit$estimate[keep,],
+ std = bfit$std[keep,], influence=bfit$influence[,,keep])
+ }
>
> # Case weights
> true6a <- with(subset(adata, x=="Maintained"), fh(time, status, wt, id))
> true6b <- with(subset(adata, x!="Maintained"), fh(time, status, wt, id))
>
> fit6 <- survfit(Surv(time, status) ~ x, weight=wt, data=adata, stype=2,
+ ctype=2, robust=FALSE)
> aeq(fit6$cumhaz, c(true6a$estimate[,2], true6b$estimate[,2]))
[1] TRUE
> aeq(fit6$surv, exp(-c(true6a$estimate[,2], true6b$estimate[,2])))
[1] TRUE
> aeq(fit6$std.chaz, c(true6a$std[,4], true6b$std[,4]))
[1] TRUE
> aeq(fit6$n.risk, c(true6a$n.risk, true6b$n.risk))
[1] TRUE
> aeq(fit6$n.event, c(true6a$n.event, true6b$n.event))
[1] TRUE
>
> # Robust variance
> fit7 <- survfit(Surv(time, status) ~ x, weight=wt, data=adata, stype=2,ctype=2,
+ id=id, influence=2, robust=TRUE)
> aeq(fit7$cumhaz, c(true6a$estimate[,2], true6b$estimate[,2]))
[1] TRUE
> aeq(fit7$surv, exp(-c(true6a$estimate[,2], true6b$estimate[,2])))
[1] TRUE
> aeq(fit7$std.chaz, c(true6a$std[,2], true6b$std[,2]))
[1] TRUE
> aeq(fit7$n.risk, c(true6a$n.risk, true6b$n.risk))
[1] TRUE
> aeq(fit7$n.event, c(true6a$n.event, true6b$n.event))
[1] TRUE
> aeq(fit7$influence.chaz[[1]], true6a$influence[,2,])
[1] TRUE
> aeq(fit7$influence.chaz[[2]], true6b$influence[,2,])
[1] TRUE
>
> # compute the influence by brute force
> tdata <- subset(adata, x != "Maintained")
> tdata <- tdata[order(tdata$id),]
> eps <- 1e-8
> imat <- matrix(0., 12, 10)
> t1 <- survfit(Surv(time, status) ~x, data=tdata, ctype=2, weights=wt)
> for (i in 1:12) {
+ wtemp <- tdata$wt
+ wtemp[i] <- wtemp[i] + eps
+ tfit <-survfit(Surv(time, status) ~x, data=tdata, ctype=2,
+ weights=wtemp)
+ imat[i,] <- tdata$wt[i] * (tfit$cumhaz - t1$cumhaz)/eps
+ }
> aeq(fit7$influence.chaz[[2]], imat, tol=sqrt(eps))
[1] TRUE
>
> #
> # verify that the times and scale arguments work as expected. They
> # are in the summary and print.survfit functions.
> #
> s1 <- summary(fit1, scale=1)
> s2 <- summary(fit1, scale=2)
> aeq(s1$time/2, s2$time) #times change
[1] TRUE
> aeq(s1$surv, s2$surv)
[1] TRUE
> tscale <- rep(c(1,1,1,1, 2,2,2,2,2), each=2)
> aeq(s1$table, s2$table *tscale)
[1] TRUE
>
> s3 <- summary(fit1, scale=1, times=c(9, 18, 23, 33, 34))
> s4 <- summary(fit1, scale=2, times=c(9, 18, 23, 33, 34))
> aeq(s3$time, s4$time*2)
[1] TRUE
> aeq(s3$surv, s4$surv)
[1] TRUE
>
> print(fit1, rmean='common')
Call: survfit(formula = Surv(time, status) ~ x, data = adata, stype = 2)
n events rmean* se(rmean) median 0.95LCL 0.95UCL
x=Maintained 11 7 60.3 25.60 34 23 NA
x=Nonmaintained 12 11 30.1 9.14 27 8 NA
* restricted mean with upper limit = 161
> print(fit1, rmean='common', scale=2)
Call: survfit(formula = Surv(time, status) ~ x, data = adata, stype = 2)
n events rmean* se(rmean) median 0.95LCL 0.95UCL
x=Maintained 11 7 30.2 12.80 17.0 11.5 NA
x=Nonmaintained 12 11 15.1 4.57 13.5 4.0 NA
* restricted mean with upper limit = 80.5
>
> proc.time()
user system elapsed
0.944 0.056 0.993