### R code from vignette source 'multi.Rnw'

###################################################
### code chunk number 1: multi.Rnw:32-75
###################################################
require(survival)
#require(Rcolorbrewer)
#brewer.pal(5, "Dark2")
palette(c("#000000", "#D95F02", "#1B9E77", "#7570B3", "#E7298A", "#66A61E"))
options(continue=' ')

# These functions are used in the document, but not discussed until the end
crisk <- function(what, horizontal = TRUE, ...) {
    nstate <- length(what)
    connect <- matrix(0, nstate, nstate,
                      dimnames=list(what, what))
    connect[1,-1] <- 1  # an arrow from state 1 to each of the others
    if (horizontal) statefig(c(1, nstate-1),  connect, ...)
    else statefig(matrix(c(1, nstate-1), ncol=1), connect, ...)
}

state3 <- function(what, horizontal=TRUE, ...) {
    if (length(what) != 3) stop("Should be 3 states")
    connect <- matrix(c(0,0,0, 1,0,0, 1,1,0), 3,3,
                      dimnames=list(what, what))
    if (horizontal) statefig(1:2, connect, ...)
    else statefig(matrix(1:2, ncol=1), connect, ...)
}

state4 <- function() {
    sname <- c("Entry", "CR", "Transplant", "Transplant")
    layout <- cbind(c(1/2, 3/4, 1/4, 3/4),
                    c(5/6, 1/2, 1/2, 1/6))
    connect <- matrix(0,4,4, dimnames=list(sname, sname))
    connect[1, 2:3] <- 1
    connect[2,4] <- 1
    statefig(layout, connect)
}

state5 <- function(what, ...) {
    sname <- c("Entry", "CR", "Tx", "Rel", "Death")
    connect <- matrix(0, 5, 5, dimnames=list(sname, sname))
    connect[1, -1] <- c(1,1,1, 1.4)
    connect[2, 3:5] <- c(1, 1.4, 1)
    connect[3, c(2,4,5)] <- 1
    connect[4, c(3,5)]  <- 1
    statefig(matrix(c(1,3,1)), connect, cex=.8,...)
}


###################################################
### code chunk number 2: multi.Rnw:81-82 (eval = FALSE)
###################################################
## curves <- survfit(Surv(time, status) ~ group, data=mydata)


###################################################
### code chunk number 3: simple1
###################################################
set.seed(1952)
crdata <- data.frame(time=1:11, 
                     endpoint=factor(c(1,1,2,0,1,1,3,0,2,3,0),
                                     labels=c("censor", "a", "b", "c")))
tfit  <- survfit(Surv(time, endpoint) ~ 1, data=crdata)
dim(tfit)
summary(tfit)


###################################################
### code chunk number 4: multi.Rnw:112-113
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(tfit, col=1:3, lwd=2, ylab="Probability in state")


###################################################
### code chunk number 5: overall
###################################################
myeloid[1:5,]


###################################################
### code chunk number 6: sfit0
###################################################
getOption("SweaveHooks")[["fig"]]()
sfit0 <- survfit(Surv(futime, death) ~ trt, myeloid)
plot(sfit0, xscale=365.25, xaxs='r', col=1:2, lwd=2,
     xlab="Years post enrollment", ylab="Survival")
legend(20, .4, c("Arm A", "Arm B"),
       col=1:2, lwd=2, bty='n')


###################################################
### code chunk number 7: data1
###################################################
data1 <- myeloid
data1$crstat <- factor(with(data1, ifelse(is.na(crtime), death, 2)),
                        labels=c("censor", "death", "CR"))
data1$crtime <- with(data1, ifelse(crstat=="CR", crtime, futime))

data1$txstat <- factor(with(data1, ifelse(is.na(txtime), death, 2)),
                        labels=c("censor", "death", "transplant"))
data1$txtime <- with(data1, ifelse(txstat=="transplant", txtime, futime))
for (i in c("futime", "crtime", "txtime", "rltime"))
    data1[[i]] <- data1[[i]] * 12/365.25  #rescale to months


###################################################
### code chunk number 8: curve1
###################################################
getOption("SweaveHooks")[["fig"]]()
sfit1 <- survfit(Surv(futime, death) ~ trt, data1) #survival
sfit2 <- survfit(Surv(crtime, crstat) ~ trt, data1) # CR
sfit3 <- survfit(Surv(txtime, txstat) ~ trt, data1)

layout(matrix(c(1,1,1,2,3,4), 3,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1.1, .1))
plot(sfit2[,2], mark.time=FALSE, fun='event', xmax=48,
         lty=3, lwd=2, col=1:2, xaxt='n',
     xlab="Months post enrollment", ylab="Events")
lines(sfit1, mark.time=FALSE, xmax=48, fun='event', col=1:2, lwd=2)
lines(sfit3[,2], mark.time=FALSE, xmax=48, fun='event', col=1:2, 
          lty=2, lwd=2)

xtime <- c(0, 6, 12, 24, 36, 48)
axis(1, xtime, xtime) #marks every year rather than 10 months
temp <- outer(c("A", "B"), c("death", "transplant", "CR"),  paste)
temp[7] <- ""
legend(25, .3, temp[c(1,2,7,3,4,7,5,6,7)], lty=c(1,1,1, 2,2,2 ,3,3,3),
       col=c(1,2,0), bty='n', lwd=2)
abline(v=2, lty=2, col=3)

# add the state space diagrams
par(mar=c(4,.1,1,1))
crisk(c("Entry","Death", "CR"), alty=3)
crisk(c("Entry","Death", "Tx"), alty=2)
crisk(c("Entry","Death"))
par(oldpar)


###################################################
### code chunk number 9: badfit
###################################################
getOption("SweaveHooks")[["fig"]]()
badfit <- survfit(Surv(txtime, txstat=="transplant") ~ trt, data1)

layout(matrix(c(1,1,1,2,3,4), 3,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1.1, .1))
plot(badfit, fun="event", xmax=48, xaxt='n', col=1:2, lty=2, lwd=2,
     xlab="Months from enrollment", ylab="P(state)")
axis(1, xtime, xtime)
lines(sfit3[,2], fun='event',  xmax=48, col=1:2, lwd=2)
legend(24, .3, c("Arm A", "Arm B"), lty=1, lwd=2,
       col=1:2, bty='n', cex=1.2)

par(mar=c(4,.1,1,1))
crisk(c("Entry", "transplant"), alty=2, cex=1.2)
crisk(c("Entry","transplant", "Death"), cex=1.2)
par(oldpar)


###################################################
### code chunk number 10: <data2
###################################################
temp <- myeloid
id <- which(temp$crtime == temp$txtime) # the one special person
temp$crtime[id] <- temp$crtime[id] -1   # move their CR back by 1 day
data2 <- tmerge(myeloid[, c('id', 'trt')], temp,
                 id=id, death=event(futime, death),
                        transplant = event(txtime),
                        response   = event(crtime),
                        relapse    = event(rltime),
                        priortx    = tdc(txtime),
                        priorcr    = tdc(crtime))
attr(data2, "tcount")

data2$event <- with(data2, factor(death + 2*response + 3*transplant + 
                                4*relapse, 0:4,
                                labels=c("censor", "death", "CR", 
                                         "transplant", "relapse")))
data2[1:10,c(1:4, 11, 9, 10)]


###################################################
### code chunk number 11: data2b
###################################################
for (i in c("tstart", "tstop"))
    data2[[i]] <- data2[[i]] *12/365.25  #scale to months

ctab <- table(table(data2$id))
ctab
with(data2, table( table(id, event)))

etab <- table(data2$event, useNA="ifany")
etab


###################################################
### code chunk number 12: cr2
###################################################
getOption("SweaveHooks")[["fig"]]()
crstat <- data2$event
crstat[crstat=="transplant"] <- "censor" # ignore transplants
crsurv <- survfit(Surv(tstart, tstop, crstat) ~ trt,
                  data= data2, id=id, influence=TRUE)

layout(matrix(c(1,1,2,3), 2,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1.1, .1))
plot(sfit2[,2], lty=3, lwd=2, col=1:2, xmax=12, 
     xlab="Months", ylab="CR")
lines(crsurv[,2], lty=1, lwd=2, col=1:2, xmax=12)
par(mar=c(4, .1, 1, 1))
crisk( c("Entry","CR", "Death"), alty=3)
state3(c("Entry", "CR", "Death/Relapse"))

par(oldpar)


###################################################
### code chunk number 13: cr2b
###################################################
print(crsurv, rmean=48, digits=2)


###################################################
### code chunk number 14: cr2c
###################################################
temp <- summary(crsurv, rmean=48)$table
delta <- round(temp[4,3] - temp[3,3], 2)


###################################################
### code chunk number 15: txsurv
###################################################
getOption("SweaveHooks")[["fig"]]()
event2 <- with(data2, ifelse(event=="transplant" & priorcr==1, 6,
               as.numeric(event)))
event2 <- factor(event2, 1:6, c(levels(data2$event), "tx after CR"))
txsurv <- survfit(Surv(tstart, tstop, event2) ~ trt, data2, id=id,
                  subset=(priortx ==0))

layout(matrix(c(1,1,1,2,2,0),3,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1,.1))
plot(txsurv[,c(3,5)], col=1:2, lty=c(1,1,2,2), lwd=2, xmax=48,
     xaxt='n', xlab="Months", ylab="Transplanted")
axis(1, xtime, xtime)
legend(15, .13, c("A, transplant without CR", "B, transplant without CR",
                 "A, transplant after CR", "B, transplant after CR"),
       col=1:2, lty=c(1,1,2,2), lwd=2, bty='n')
state4()  # add the state figure
par(oldpar)


###################################################
### code chunk number 16: sfit4
###################################################
getOption("SweaveHooks")[["fig"]]()
sfit4 <- survfit(Surv(tstart, tstop, event) ~ trt, data2, id=id)
sfit4$transitions
layout(matrix(1:2,1,2), widths=2:1)
oldpar <- par(mar=c(5.1, 4.1, 1,.1))
plot(sfit4, col=rep(1:4,each=2), lwd=2, lty=1:2, xmax=48, xaxt='n',
     xlab="Months", ylab="Current state")
axis(1, xtime, xtime)
text(c(40, 40, 40, 40), c(.51, .13, .32, .01),
     c("Death", "CR", "Transplant", "Recurrence"), col=1:4)

par(mar=c(5.1, .1, 1, .1))
state5()
par(oldpar)


###################################################
### code chunk number 17: reprise
###################################################
crsurv <- survfit(Surv(tstart, tstop, crstat) ~ trt,
                  data= data2, id=id, influence=TRUE)
curveA <- crsurv[1,]  # select treatment A
dim(curveA$pstate)    # P matrix for treatement A
dim(curveA$influence) # influence matrix for treatment A
table(data1$trt)
curveA$p0             # state distribution at time 0


###################################################
### code chunk number 18: meantime
###################################################
t48 <- pmin(48, curveA$time)   
delta <- diff(c(0, t48, 48))  # width of intervals
rfun <- function(pmat, delta) colSums(pmat * delta)  #area under the curve
rmean <- rfun(rbind(curveA$p0, curveA$pstate), delta)
round(rmean, 2)

# Apply the same calculation to each subject's influence slice
inf <- apply(curveA$influence, 1, rfun, delta=delta)
# inf is now a 5 state by 310 subject matrix, containing the IJ estimates
#  on the AUC or mean time.  The sum of squares is a variance.
se.rmean <- sqrt(rowSums(inf^2))
round(se.rmean, 2)


###################################################
### code chunk number 19: crisk
###################################################
crisk <- function(what, horizontal = TRUE, ...) {
    nstate <- length(what)
    connect <- matrix(0, nstate, nstate,
                      dimnames=list(what, what))
    connect[1,-1] <- 1  # an arrow from state 1 to each of the others
    if (horizontal) statefig(c(1, nstate-1),  connect, ...)
    else statefig(matrix(c(1, nstate-1), ncol=1), connect, ...)
}


###################################################
### code chunk number 20: state3
###################################################
state3 <- function(what, horizontal=TRUE, ...) {
    if (length(what) != 3) stop("Should be 3 states")
    connect <- matrix(c(0,0,0, 1,0,0, 1,1,0), 3,3,
                      dimnames=list(what, what))
    if (horizontal) statefig(1:2, connect, ...)
    else statefig(matrix(1:2, ncol=1), connect, ...)
}


###################################################
### code chunk number 21: state5
###################################################
state5 <- function(what, ...) {
    sname <- c("Entry", "CR", "Tx", "Rel", "Death")
    connect <- matrix(0, 5, 5, dimnames=list(sname, sname))
    connect[1, -1] <- c(1,1,1, 1.4)
    connect[2, 3:5] <- c(1, 1.4, 1)
    connect[3, c(2,4,5)] <- 1
    connect[4, c(3,5)]  <- 1
    statefig(matrix(c(1,3,1)), connect, cex=.8, ...)
}


###################################################
### code chunk number 22: state4
###################################################
state4 <- function() {
    sname <- c("Entry", "CR", "Transplant", "Transplant")
    layout <- cbind(x =c(1/2, 3/4, 1/4, 3/4),
                    y =c(5/6, 1/2, 1/2, 1/6))
    connect <- matrix(0,4,4, dimnames=list(sname, sname))
    connect[1, 2:3] <- 1
    connect[2,4] <- 1
    statefig(layout, connect)
}