Nothing
### R code from vignette source 'survival.Rnw'
###################################################
### code chunk number 1: survival.Rnw:39-83
###################################################
options(continue=" ", width=70)
options(SweaveHooks=list(fig=function() par(mar=c(4.1, 4.1, .3, 1.1))))
pdf.options(pointsize=10) #text in graph about the same as regular text
options(contrasts=c("contr.treatment", "contr.poly")) #ensure default
library("survival")
palette(c("#000000", "#D95F02", "#1B9E77", "#7570B3", "#E7298A", "#66A61E"))
# 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: states
###################################################
getOption("SweaveHooks")[["fig"]]()
oldpar <- par(mar=c(.1, .1, .1, .1), mfrow=c(2,2))
sname1 <- c("Alive", "Dead")
cmat1 <- matrix(c(0,0,1,0), nrow=2,
dimnames=list(sname1, sname1))
statefig(c(1,1), cmat1)
sname2 <- c("0", "1", "2", "...")
cmat2 <- matrix(0, 4,4, dimnames= list(sname2, sname2))
cmat2[1,2] <- cmat2[2,3] <- cmat2[3,4] <- 1
statefig(c(1,1,1,1), cmat2, bcol=c(1,1,1,0))
sname3 <- c("Entry", "Transplant", "Withdrawal", "Death")
cmat3 <- matrix(0, 4,4, dimnames=list(sname3, sname3))
cmat3[1, -1] <- 1
statefig(c(1,3), cmat3)
sname4 <- c("Health", "Illness", "Death")
cmat4 <- matrix(0, 3, 3, dimnames = list(sname4, sname4))
cmat4[1,2] <- cmat4[2,1] <- cmat4[-3, 3] <- 1
statefig(c(1,2), cmat4, offset=.03)
par(oldpar)
###################################################
### code chunk number 3: survfit1
###################################################
fit1 <- survfit(Surv(futime, fustat) ~ resid.ds, data=ovarian)
print(fit1, rmean= 730)
summary(fit1, times= (0:4)*182.5, scale=365)
###################################################
### code chunk number 4: survfit2
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(fit1, col=1:2, xscale=365.25, lwd=2, mark.time=TRUE,
xlab="Years since study entry", ylab="Survival")
legend(750, .9, c("No residual disease", "Residual disease"),
col=1:2, lwd=2, bty='n')
###################################################
### code chunk number 5: survfit3
###################################################
getOption("SweaveHooks")[["fig"]]()
fit2 <- survfit(Surv(time, status) ~ sex + ph.ecog, data=lung)
fit2
plot(fit2[1:3], lty=1:3, lwd=2, xscale=365.25, fun='event',
xlab="Years after enrollment", ylab="Survival")
legend(550, .6, paste("Performance Score", 0:2, sep=' ='),
lty=1:3, lwd=2, bty='n')
text(400, .95, "Males", cex=2)
###################################################
### code chunk number 6: survival.Rnw:576-578
###################################################
data.frame(id=rep(392,3), time1=c(0, 258, 328), time2=c(258, 328, 377),
status=c(1,1,0))
###################################################
### code chunk number 7: survival4
###################################################
vdata <- with(valveSeat, data.frame(id=id, time2=time, status=status))
first <- !duplicated(vdata$id)
vdata$time1 <- ifelse(first, 0, c(0, vdata$time[-nrow(vdata)]))
double <- which(vdata$time1 == vdata$time2)
vdata$time1[double] <- vdata$time1[double] -.01
vdata$time2[double-1] <- vdata$time1[double]
vdata[1:7, c("id", "time1", "time2", "status")]
survcheck(Surv(time1, time2, status) ~ 1, id=id, data=vdata)
###################################################
### code chunk number 8: survival5
###################################################
getOption("SweaveHooks")[["fig"]]()
vfit <- survfit(Surv(time1, time2, status) ~1, data=vdata, id=id)
plot(vfit, cumhaz=TRUE, xlab="Days", ylab="Cumulative hazard")
###################################################
### code chunk number 9: cgd1d
###################################################
getOption("SweaveHooks")[["fig"]]()
cgdsurv <- survfit(Surv(tstart, tstop, status) ~ treat, cgd, id=id)
plot(cgdsurv, cumhaz=TRUE, col=1:2, conf.times=c(100, 200, 300, 400),
xlab="Days since randomization", ylab="Cumulative hazard")
###################################################
### code chunk number 10: simple1
###################################################
crdata <- data.frame(time= c(1:8, 6:8),
endpoint=factor(c(1,1,2,0,1,1,3,0,2,3,0),
labels=c("censor", "a", "b", "c")),
istate=rep("entry", 11),
id= LETTERS[1:11])
tfit <- survfit(Surv(time, endpoint) ~ 1, data=crdata, id=id, istate=istate)
dim(tfit)
summary(tfit)
###################################################
### code chunk number 11: survival.Rnw:696-697
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(tfit, col=1:4, lty=1:4, lwd=2, ylab="Probability in state")
###################################################
### code chunk number 12: survival.Rnw:710-712
###################################################
dim(tfit)
tfit$states
###################################################
### code chunk number 13: mgus1
###################################################
mgus2[55:59, -(4:7)]
###################################################
### code chunk number 14: mgus2
###################################################
getOption("SweaveHooks")[["fig"]]()
event <- with(mgus2, ifelse(pstat==1, 1, 2*death))
event <- factor(event, 0:2, c("censored", "progression", "death"))
etime <- with(mgus2, ifelse(pstat==1, ptime, futime))
crfit <- survfit(Surv(etime, event) ~ sex, mgus2)
crfit
plot(crfit, col=1:2, noplot="",
lty=c(3,3,2,2,1,1), lwd=2, xscale=12,
xlab="Years post diagnosis", ylab="P(state)")
legend(240, .65, c("Female, death", "Male, death", "malignancy", "(s0)"),
lty=c(1,1,2,3), col=c(1,2,1,1), bty='n', lwd=2)
###################################################
### code chunk number 15: mgus3
###################################################
getOption("SweaveHooks")[["fig"]]()
pcmbad <- survfit(Surv(etime, pstat) ~ sex, data=mgus2)
plot(pcmbad[2], mark.time=FALSE, lwd=2, fun="event", conf.int=FALSE, xscale=12,
xlab="Years post diagnosis", ylab="Fraction with PCM")
lines(crfit[2,2], lty=2, lwd=2, mark.time=FALSE, conf.int=FALSE)
legend(0, .25, c("Males, PCM, incorrect curve", "Males, PCM, competing risk"),
col=1, lwd=2, lty=c(1,2), bty='n')
###################################################
### code chunk number 16: survival.Rnw:844-847
###################################################
dim(crfit)
crfit$strata
crfit$states
###################################################
### code chunk number 17: overall
###################################################
myeloid[1:5,]
###################################################
### code chunk number 18: 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 19: sfit0a
###################################################
mdata <- tmerge(myeloid[,1:2], myeloid, id=id, death= event(futime, death),
sct = event(txtime), cr = event(crtime),
relapse = event(rltime))
temp <- with(mdata, cr + 2*sct + 4*relapse + 8*death)
table(temp)
###################################################
### code chunk number 20: sfit0b
###################################################
tdata <- myeloid # temporary working copy
tied <- with(tdata, (!is.na(crtime) & !is.na(txtime) & crtime==txtime))
tdata$crtime[tied] <- tdata$crtime[tied] -1
mdata <- tmerge(tdata[,1:2], tdata, id=id, death= event(futime, death),
sct = event(txtime), cr = event(crtime),
relapse = event(rltime),
priorcr = tdc(crtime), priortx = tdc(txtime))
temp <- with(mdata, cr + 2*sct + 4*relapse + 8*death)
table(temp)
mdata$event <- factor(temp, c(0,1,2,4,8),
c("none", "CR", "SCT", "relapse", "death"))
mdata[1:7, c("id", "trt", "tstart", "tstop", "event", "priorcr", "priortx")]
###################################################
### code chunk number 21: survival.Rnw:982-983
###################################################
survcheck(Surv(tstart, tstop, event) ~1, mdata, id=id)
###################################################
### code chunk number 22: newevent
###################################################
levels(mdata$event)
temp1 <- with(mdata, ifelse(priorcr, 0, c(0,1,0,0,2)[event]))
mdata$crstat <- factor(temp1, 0:2, c("none", "CR", "death"))
temp2 <- with(mdata, ifelse(priortx, 0, c(0,0,1,0,2)[event]))
mdata$txstat <- factor(temp2, 0:2, c("censor", "SCT", "death"))
temp3 <- with(mdata, c(0,0,1,0,2)[event] + priortx)
mdata$tx2 <- factor(temp3, 0:3,
c("censor", "SCT", "death w/o SCT", "death after SCT"))
###################################################
### code chunk number 23: curve1
###################################################
getOption("SweaveHooks")[["fig"]]()
# I want to have the plots in months, it is simpler to fix time
# once rather than repeat xscale many times
tdata$futime <- tdata$futime * 12 /365.25
mdata$tstart <- mdata$tstart * 12 /365.25
mdata$tstop <- mdata$tstop * 12 /365.25
sfit1 <- survfit(Surv(futime, death) ~ trt, tdata) # survival
sfit2 <- survfit(Surv(tstart, tstop, crstat) ~ trt,
data= mdata, id = id) # CR
sfit3 <- survfit(Surv(tstart, tstop, txstat) ~ trt,
data= mdata, id =id) # SCT
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))
mlim <- c(0, 48) # and only show the first 4 years
plot(sfit2[,"CR"], xlim=mlim,
lty=3, lwd=2, col=1:2, xaxt='n',
xlab="Months post enrollment", ylab="Fraction with the endpoint")
lines(sfit1, mark.time=FALSE, xlim=mlim,
fun='event', col=1:2, lwd=2)
lines(sfit3[,"SCT"], xlim=mlim, col=1:2,
lty=2, lwd=2)
xtime <- c(0, 6, 12, 24, 36, 48)
axis(1, xtime, xtime) #axis marks every year rather than 10 months
temp <- outer(c("A", "B"), c("CR", "transplant", "death"), paste)
temp[7] <- ""
legend(25, .3, temp[c(1,2,7,3,4,7,5,6,7)], lty=c(3,3,3, 2,2,2 ,1,1,1),
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", "CR", "Death"), alty=3)
crisk(c("Entry", "Tx", "Death"), alty=2)
crisk(c("Entry","Death"))
par(oldpar)
layout(1)
###################################################
### code chunk number 24: badfit
###################################################
getOption("SweaveHooks")[["fig"]]()
badfit <- survfit(Surv(tstart, tstop, event=="SCT") ~ trt,
id=id, mdata, subset=(priortx==0))
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], 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)
layout(1)
###################################################
### code chunk number 25: cr2
###################################################
getOption("SweaveHooks")[["fig"]]()
cr2 <- mdata$event
cr2[cr2=="SCT"] <- "none" # ignore transplants
crsurv <- survfit(Surv(tstart, tstop, cr2) ~ trt,
data= mdata, 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)
par(mar=c(4, .1, 1, 1))
crisk( c("Entry","CR", "Death"), alty=3)
state3(c("Entry", "CR", "Death/Relapse"))
par(oldpar)
layout(1)
###################################################
### code chunk number 26: cr2b
###################################################
print(crsurv, rmean=48, digits=2)
###################################################
### code chunk number 27: cr2c
###################################################
temp <- summary(crsurv, rmean=48)$table
delta <- round(temp[4,3] - temp[3,3], 2)
###################################################
### code chunk number 28: txsurv
###################################################
getOption("SweaveHooks")[["fig"]]()
event2 <- with(mdata, ifelse(event=="SCT" & priorcr==1, 6,
as.numeric(event)))
event2 <- factor(event2, 1:6, c(levels(mdata$event), "SCT after CR"))
txsurv <- survfit(Surv(tstart, tstop, event2) ~ trt, mdata, id=id,
subset=(priortx ==0))
dim(txsurv) # number of strata by number of states
txsurv$states # Names of states
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,6)], 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 29: sfit4
###################################################
getOption("SweaveHooks")[["fig"]]()
sfit4 <- survfit(Surv(tstart, tstop, event) ~ trt, mdata, 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=c(4,1,2,3))
par(mar=c(5.1, .1, 1, .1))
state5()
par(oldpar)
###################################################
### code chunk number 30: reprise
###################################################
crsurv <- survfit(Surv(tstart, tstop, cr2) ~ trt,
data= mdata, id=id, influence=TRUE)
curveA <- crsurv[1,] # select treatment A
dim(curveA) # P matrix for treatement A
curveA$states
dim(curveA$pstate) # 426 time points, 5 states
dim(curveA$influence) # influence matrix for treatment A
table(myeloid$trt)
###################################################
### code chunk number 31: meantime
###################################################
t48 <- pmin(48, curveA$time)
delta <- diff(c(t48, 48)) # width of intervals
rfun <- function(pmat, delta) colSums(pmat * delta) #area under the curve
rmean <- rfun(curveA$pstate, delta)
# 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(rbind(rmean, se.rmean), 2)
print(curveA, rmean=48, digits=2)
###################################################
### code chunk number 32: survdiff
###################################################
survdiff(Surv(time, status) ~ x, aml)
###################################################
### code chunk number 33: 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 34: 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 35: 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 36: 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)
}
###################################################
### code chunk number 37: lung1
###################################################
options(show.signif.stars=FALSE) # display statistical intelligence
cfit1 <- coxph(Surv(time, status) ~ age + sex + wt.loss, data=lung)
print(cfit1, digits=3)
summary(cfit1, digits=3)
anova(cfit1)
###################################################
### code chunk number 38: na.action
###################################################
cfit1a <- coxph(Surv(time, status) ~ age + sex + wt.loss, data=lung,
na.action = na.omit)
cfit1b <- coxph(Surv(time, status) ~ age + sex + wt.loss, data=lung,
na.action = na.exclude)
r1 <- residuals(cfit1a)
r2 <- residuals(cfit1b)
length(r1)
length(r2)
###################################################
### code chunk number 39: cox12
###################################################
cfit2 <- coxph(Surv(time, status) ~ age + sex + wt.loss + strata(inst),
data=lung)
round(cbind(simple= coef(cfit1), stratified=coef(cfit2)), 4)
###################################################
### code chunk number 40: cox13
###################################################
getOption("SweaveHooks")[["fig"]]()
dummy <- expand.grid(age=c(50, 60), sex=1, wt.loss=5)
dummy
csurv1 <- survfit(cfit1, newdata=dummy)
csurv2 <- survfit(cfit2, newdata=dummy)
dim(csurv1)
dim(csurv2)
plot(csurv1, col=1:2, xscale=365.25, xlab="Years", ylab="Survival")
dummy2 <- data.frame(age=c(50, 60), sex=1:2, wt.loss=5, inst=c(6,11))
csurv3 <- survfit(cfit2, newdata=dummy2)
dim(csurv3)
###################################################
### code chunk number 41: lung2
###################################################
getOption("SweaveHooks")[["fig"]]()
zp1 <- cox.zph(cfit1)
zp1
plot(zp1[2], resid=FALSE)
abline(coef(cfit1)[2] ,0, lty=3)
###################################################
### code chunk number 42: lung3
###################################################
getOption("SweaveHooks")[["fig"]]()
cfit3 <- coxph(Surv(time, status) ~ pspline(age) + sex + wt.loss, lung)
print(cfit3, digits=2)
termplot(cfit3, term=1, se=TRUE)
cfit4 <- update(cfit1, . ~ . + age*sex)
anova(cfit1, cfit4)
###################################################
### code chunk number 43: cgd1
###################################################
cfit1 <- coxph(Surv(tstart, tstop, status) ~ treat + inherit + steroids +
age + strata(hos.cat), data=cgd)
print(cfit1, digits=2)
###################################################
### code chunk number 44: cgd1b
###################################################
cfit2 <- coxph(Surv(tstart, tstop, status) ~ treat + inherit+
age + strata(hos.cat), data=cgd)
print(cfit2, digits=2)
###################################################
### code chunk number 45: cgd3
###################################################
getOption("SweaveHooks")[["fig"]]()
dummy <- expand.grid(age=c(6,12), inherit='X-linked',
treat=levels(cgd$treat))
dummy
csurv <- survfit(cfit2, newdata=dummy)
dim(csurv)
plot(csurv[1,], fun="event", col=1:2, lty=c(1,1,2,2),
xlab="Days on study", ylab="Pr( any infection )")
###################################################
### code chunk number 46: cfit4
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(csurv[1,], cumhaz=TRUE, col=1:2, lty=c(1,1,2,2), lwd=2,
xlab="Days on study", ylab="E( number of infections )")
legend(20, 1.5, c("Age 6, control", "Age 12, control",
"Age 6, gamma interferon", "Age 12, gamma interferon"),
lty=c(2,2,1,1), col=c(1,2,1,2), lwd=2, bty='n')
###################################################
### code chunk number 47: survfit-mgus1
###################################################
getOption("SweaveHooks")[["fig"]]()
mgus2[56:59,]
sname <- c("MGUS", "Malignancy", "Death")
smat <- matrix(c(0,0,0, 1,0,0, 1,1,0), 3, 3,
dimnames = list(sname, sname))
statefig(c(1,2), smat)
###################################################
### code chunk number 48: survfit-mgus2
###################################################
crdata <- mgus2
crdata$etime <- pmin(crdata$ptime, crdata$futime)
crdata$event <- ifelse(crdata$pstat==1, 1, 2*crdata$death)
crdata$event <- factor(crdata$event, 0:2, c("censor", "PCM", "death"))
quantile(crdata$age, na.rm=TRUE)
table(crdata$sex)
quantile(crdata$mspike, na.rm=TRUE)
cfit <- coxph(Surv(etime, event) ~ I(age/10) + sex + mspike,
id = id, crdata)
print(cfit, digits=1) # narrow the printout a bit
###################################################
### code chunk number 49: PCMcurve
###################################################
getOption("SweaveHooks")[["fig"]]()
dummy <- expand.grid(sex=c("F", "M"), age=c(60, 80), mspike=1.2)
csurv <- survfit(cfit, newdata=dummy)
plot(csurv[,2], xmax=20*12, xscale=12,
xlab="Years after MGUS diagnosis", ylab="Pr(has entered PCM state)",
col=1:2, lty=c(1,1,2,2), lwd=2)
legend(100, .04, outer(c("female,", "male, "),
c("diagnosis at age 60", "diagnosis at age 80"),
paste),
col=1:2, lty=c(1,1,2,2), bty='n', lwd=2)
###################################################
### code chunk number 50: mrate
###################################################
mpfit <- glm(pstat ~ sex -1 + offset(log(ptime)), data=mgus2, poisson)
exp(coef(mpfit)) * 12 # rate per year
###################################################
### code chunk number 51: msingle
###################################################
getOption("SweaveHooks")[["fig"]]()
sfit <- coxph(Surv(etime, event=="PCM") ~ I(age/10) + sex + mspike, crdata)
rbind(single = coef(sfit),
multi = coef(cfit)[1:3])
#par(mfrow=c(1,2))
ssurv <- survfit(sfit, newdata=dummy)
plot(ssurv[3:4], col=1:2, lty=2, xscale=12, xmax=12*20, lwd=2, fun="event",
xlab="Years from diagnosis", ylab= "Pr(has entered PCM state)")
lines(csurv[3:4, 2], col=1:2, lty=1, lwd=2)
legend(20, .22, outer(c("80 year old female,", "80 year old male,"),
c("incorrect", "correct"), paste),
col=1:2, lty=c(2,2,1,1), lwd=2, bty='n')
###################################################
### code chunk number 52: state5
###################################################
getOption("SweaveHooks")[["fig"]]()
state5 <- c("0MC", "1MC", "2MC", "3MC", "death")
tmat <- matrix(0L, 5, 5, dimnames=list(state5, state5))
tmat[1,2] <- tmat[2,3] <- tmat[3,4] <- 1
tmat[-5,5] <- 1
statefig(rbind(4,1), tmat)
###################################################
### code chunk number 53: nafld1
###################################################
ndata <- tmerge(nafld1[,1:8], nafld1, id=id, death= event(futime, status))
ndata <- tmerge(ndata, subset(nafld3, event=="nafld"), id,
nafld= tdc(days))
ndata <- tmerge(ndata, subset(nafld3, event=="diabetes"), id = id,
diabetes = tdc(days), e1= cumevent(days))
ndata <- tmerge(ndata, subset(nafld3, event=="htn"), id = id,
htn = tdc(days), e2 = cumevent(days))
ndata <- tmerge(ndata, subset(nafld3, event=="dyslipidemia"), id=id,
lipid = tdc(days), e3= cumevent(days))
ndata <- tmerge(ndata, subset(nafld3, event %in% c("diabetes", "htn",
"dyslipidemia")),
id=id, comorbid= cumevent(days))
summary(ndata)
###################################################
### code chunk number 54: survival.Rnw:2034-2037
###################################################
tc <- attr(ndata, 'tcount') # shorter name for use in Sexpr below
icount <- table(table(nafld3$id)) #number with 1, 2, ... intervals
ncount <- sum(nafld3$event=="nafld")
###################################################
### code chunk number 55: nafld2
###################################################
with(ndata, if (any(e1>1 | e2>1 | e3>1)) stop("multiple events"))
ndata$cstate <- with(ndata, factor(diabetes + htn + lipid, 0:3,
c("0mc", "1mc", "2mc", "3mc")))
temp <- with(ndata, ifelse(death, 4, comorbid))
ndata$event <- factor(temp, 0:4,
c("censored", "1mc", "2mc", "3mc", "death"))
ndata$age1 <- ndata$age + ndata$tstart/365.25 # analysis on age scale
ndata$age2 <- ndata$age + ndata$tstop/365.25
check1 <- survcheck(Surv(age1, age2, event) ~ nafld + male, data=ndata,
id=id, istate=cstate)
check1
###################################################
### code chunk number 56: nafld3
###################################################
getOption("SweaveHooks")[["fig"]]()
states <- c("No comorbidity", "1 comorbidity", "2 comorbidities",
"3 comorbitities", "Death")
cmat <- matrix(0, 5,5)
cmat[,5] <- 1
cmat[1,2] <- cmat[2,3] <- cmat[3,4] <- 1
cmat[1,3] <- cmat[2,4] <- 1.6
cmat[1,4] <- 1.6
dimnames(cmat) <- list(states, states)
statefig(cbind(4,1), cmat, alty=c(1,2,1,2,2,1,1,1,1,1,1))
###################################################
### code chunk number 57: nafld4
###################################################
nfit1 <- coxph(list(Surv(age1, age2, event) ~ nafld + male,
"0mc":state("1mc", "2mc", "3mc") ~ nafld+ male / common,
2:3 + 2:4 ~ nafld + male / common,
0:"death" ~ male / common),
data=ndata, id=id, istate=cstate)
nfit1$states
nfit1$cmap
###################################################
### code chunk number 58: nafld5b
###################################################
print(coef(nfit1), digits=3)
print(coef(nfit1, matrix=TRUE), digits=3) # alternate form
print(nfit1)
###################################################
### code chunk number 59: survival.Rnw:2232-2234
###################################################
options(show.signif.stars = FALSE) # display statistical maturity
summary(nfit1, digits=3)
###################################################
### code chunk number 60: nafld5c
###################################################
nfit2 <- coxph(list(Surv(age1, age2, event) ~ nafld + male,
"0mc":state("1mc", "2mc", "3mc") ~ nafld+ male / common,
2:3 + 2:4 ~ nafld + male / common,
1:5 + 2:5 +3:5 ~ male / common + shared),
data=ndata, id=id, istate=cstate)
nfit2$cmap
###################################################
### code chunk number 61: timeline1
###################################################
ctime <- with(mgus2, ifelse(pstat==1, ptime, futime))
cstat <- with(mgus2, ifelse(pstat==1, 1, 2*death))
cstat <- factor(cstat, 0:2, c("censor", "PCM", "death"))
tdata <- data.frame(id=mgus2$id, days=ctime, cstat=cstat)
# counting process
mdata1 <- tmerge(mgus2[,1:7], tdata, id, state=event(days, cstat))
mfit1 <- coxph(Surv(tstart, tstop, state) ~ age + sex, id=id, mdata1)
# timeline
mdata2 <- data.frame(mgus2[,1:7], days=0)
mdata2 <- merge(mdata2, tdata, all=TRUE)
mfit2 <- coxph(Surv2(days, cstat) ~ age + sex, id=id, mdata2)
all.equal(coef(mfit1), coef(mfit2))
###################################################
### code chunk number 62: timeline1b
###################################################
mdata1[1:3,]
print(mdata2[1:6,], na.print='.')
###################################################
### code chunk number 63: timeline2
###################################################
tldata <- data.frame(nafld1[,1:7],
days= 0, death=0, iage=nafld1$age, nafld=0)
tldata <- merge(tldata, with(nafld1, data.frame(id=id, days=futime, death=status)),
all=TRUE)
# Add in the comorbidities of interest. None of these 4 happen to have
# duplicates (MI does, for instance).
# Start by removing the the 13 rows with a "confirmed NAFLD" (actual NAFLD + 1 year)
# that is after the actual last follow-up date.
# Treat diabetes before day 0 as diabetes on day 0.
badrow <- which(nafld3$days > nafld1$futime[match(nafld3$id, nafld1$id)])
fixnf3 <- nafld3[-badrow,]
tldata <- merge(tldata, with(subset(fixnf3, event=="diabetes"),
data.frame(id=id, days=pmax(0,days), diabetes=1)),
all=TRUE, by=c("id", "days"))
tldata <- merge(tldata, with(subset(fixnf3, event=="htn"),
data.frame(id=id, days=pmax(0,days), htn=1)),
all=TRUE, by=c("id", "days"))
tldata <- merge(tldata, with(subset(fixnf3, event=="dyslipidemia"),
data.frame(id=id, days= pmax(0, days), dyslipid=1)),
all=TRUE, by=c("id", "days"))
tldata <- merge(tldata, with(subset(fixnf3, event=="nafld"),
data.frame(id=id, days= pmax(0,days), nafld=1)),
by=c("id", "days"), all=TRUE)
tldata$nafld <- with(tldata, ifelse(is.na(nafld.y), nafld.x, nafld.y))
###################################################
### code chunk number 64: timeline3
###################################################
#
# For cumulative events within subject we use a helper function
cumevent <- function(id, time, status, istate) {
# do all the work on ordered data
ord <- order(id, time)
id2 <- id[ord]
time2 <- time[ord]
stat2 <- ifelse(is.na(status[ord]), 0, status[ord])
firstid <- !duplicated(id)
csum <- cumsum(stat2)
indx <- match(id2, id2)
cstat<- csum + stat2[indx] - csum[indx]
cstat[stat2==0] <- 0
if (!missing(istate)) cstat[firstid] <- istate
keep <- (firstid | (!is.na(stat2)& stat2 !=0))
newdata <- data.frame(id=id2[keep], time=time2[keep], status=cstat[keep])
newdata
}
temp1 <- rowSums(tldata[,c('diabetes', 'htn', 'dyslipid')], na.rm=TRUE)
temp2 <- with(tldata, cumevent(id, days, pmax(temp1, 4*death, na.rm=TRUE)))
state <- factor(pmin(temp2$status, 4), -1:4,
c("censor", paste0(0:3, "mc"), "death"))
tldata <- merge(tldata, data.frame(id=temp2$id, days=temp2$time, state=state),
all=TRUE)
tldata$age <- with(tldata, days/365.25 + age[match(id, id)])
check2 <- survcheck(Surv2(days, state) ~ 1, id=id, tldata)
check2$transitions
nfit2 <- coxph(list(Surv2(age, state) ~ nafld + male,
"0mc":state("1mc", "2mc", "3mc") ~ nafld+ male / common,
2:3 + 2:4 ~ nafld + male / common,
0:"death" ~ male / common),
data=tldata, id=id)
round(coef(nfit2), 3)
###################################################
### code chunk number 65: survival.Rnw:2440-2442 (eval = FALSE)
###################################################
## fit2 <- coxph(Surv(time, status) ~ trt + trt*time + celltype + karno,
## data = veteran)
###################################################
### code chunk number 66: zphcheck1
###################################################
dtime <- unique(veteran$time[veteran$status==1]) # unique times
newdata <- survSplit(Surv(time, status) ~ trt + celltype + karno,
data=veteran, cut=dtime)
nrow(veteran)
nrow(newdata)
fit0 <- coxph(Surv(time, status) ~ trt + celltype + karno, veteran)
fit1 <- coxph(Surv(tstart, time, status) ~ trt + celltype + karno,
data=newdata)
fit2 <- coxph(Surv(tstart, time, status) ~ trt + celltype + karno +
time:karno, newdata)
fit2
fit2b <- coxph(Surv(tstart, time, status) ~ trt + celltype + karno +
rank(time):karno, newdata)
###################################################
### code chunk number 67: zph2
###################################################
fit2 <- coxph(Surv(tstart, time, status) ~ trt + celltype + karno +
tt(karno), data =newdata,
tt = function(x, t,...) x*t)
###################################################
### code chunk number 68: zph2
###################################################
getOption("SweaveHooks")[["fig"]]()
zp0 <- cox.zph(fit0, transform='identity')
zp0
zp1 <- cox.zph(fit0, transform='log')
zp1
oldpar <- par(mfrow=c(2,2))
for (i in 1:3) {plot(zp1[i]); abline(0,0, lty=3)}
plot(zp0[3])
par(oldpar)
###################################################
### code chunk number 69: profile1
###################################################
getOption("SweaveHooks")[["fig"]]()
fit1 <- coxph(Surv(futime, fustat) ~ rx + age + resid.ds, ovarian)
fit1
# create the profile plot
imat <- solve(vcov(fit1)) #information matrix
acoef <- seq(0, .25, length=100)
profile <- matrix(0, 100, 2)
for (i in 1:100) {
icoef <- c(fit1$coef[1], acoef[i], fit1$coef[3])
tfit <- coxph(Surv(futime, fustat) ~ rx + age + resid.ds, ovarian,
init= icoef, iter.max=0)
profile[i,1] <- tfit$loglik[2]
delta <- c(0, acoef[i]- fit1$coef[2], 0)
profile[i,2] <- fit1$loglik[2] - delta%*% imat %*% delta/2
}
matplot(acoef, profile*2, type='l', lwd=2, lty=1, xlab="Coefficient for age",
ylab="2*loglik")
abline(h = 2*fit1$loglik[2] - qchisq(.95, 1), lty=3)
legend(.11, -58, c("Cox likelihood", "Wald approximation"), lty=1, lwd=2,
col=1:2, bty='n')
###################################################
### code chunk number 70: profile2
###################################################
myfun <- function(beta) {
icoef <- coef(fit1)
icoef[2] <- beta
tfit <- coxph(Surv(futime, fustat) ~ rx + age + resid.ds, ovarian,
init = icoef, iter.max=0)
(fit1$loglik - tfit$loglik)[2] - qchisq(.95, 1)/2
}
uniroot(myfun, c (0, .2))$root # lower
uniroot(myfun, c(.2, .5))$root # upper
###################################################
### code chunk number 71: survival.Rnw:3490-3491
###################################################
with(subset(aml, x=="Nonmaintained"), Surv(time, status))
###################################################
### code chunk number 72: coarsen
###################################################
getOption("SweaveHooks")[["fig"]]()
tdata <- subset(colon, etype==1) # progression or death
cmat <- matrix(0, 7, 6)
for( i in 1:7) {
if (i==1) scale <-1 else scale <- (i-1)*365/12
temp <- floor(tdata$time/scale)
tfit <- coxph(Surv(temp, status) ~ node4 + extent, tdata)
tfit2 <- coxph(Surv(temp, status) ~ node4 + extent, tdata,
ties='breslow')
tfit3 <- coxph(Surv(temp, status) ~ node4 + extent, tdata,
ties='exact')
cmat[i,] <- c(coef(tfit2), coef(tfit), coef(tfit3))
}
matplot(1:7, cmat[,c(1,3,5)], xaxt='n', pch='bec',
xlab="Time divisor", ylab="Coefficient for node4")
axis(1, 1:7, c(1, floor(1:6 *365/12)))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.