data(ebmt2)
Data <- ebmt2[c(1:5,7:8,30),] # just a small (non-representative) subset
Data$Tstart <- 0
Data$Tstop <- Data$time
Data$status[Data$status==6] <- 2 # recode 6 to 2
Data$dissub <- factor(Data$dissub) # get rid of missing level
Data.weight <- crprep(Data$Tstop, Data$status, trans=c(1,2), keep=Data$dissub)
# calculate cause-specific cumulative incidence, no truncation,
# compare with Cuminc (also from mstate)
ci <- Cuminc(Data$Tstop, Data$status)
ci
sf <- survfit(Surv(Tstart,Tstop,status==1)~1, data=Data.weight,
weight=weight.cens, subset=failcode==1)
sf <- summary(sf)
data.frame(time=sf$time, CI.1=1-sf$surv)
sf <- survfit(Surv(Tstart,Tstop,status==2)~1, data=Data.weight,
weight=weight.cens, subset=failcode==2)
sf <- summary(sf)
data.frame(time=sf$time, CI.2=1-sf$surv)
# Fine and Gray regression for cause 1
cw <- coxph(Surv(Tstart,Tstop,status==1)~dissub, data=Data.weight,
weight=weight.cens, subset=failcode==1)
cw
# This can be checked with the results of crr (cmprsk)
# crr(ftime=Data$Tstop, fstatus=Data$status, cov1=as.numeric(Data$dissub)-1)
# Create simulated data
set.seed(1234)
N <- 200
p <- 0.3
p.t <- 0.6
Z1 <- rnorm(N)
Z2 <- rnorm(N)
tmp <- runif(N)
Data <- data.frame(Tstart=runif(N,0,3)*rbinom(N,1,p.t),
Tstop=ifelse( tmp>(1-p)^(exp(0.5*Z1+0.5*Z2)),
-log(1-(1-exp(log(tmp)/exp(0.5*Z1+0.5*Z2)))/p),
-log(tmp)/exp(-0.5*Z1+0.5*Z2)),
stat=ifelse(tmp>(1-p)^(exp(0.5*Z1+0.5*Z2)),1,2),
Z1=Z1, Z2=Z2, tstat=rep(0,N))
Data$type <- Data$stat
Data$tevent <- Data$Tstop
Data$cens <- runif(N,0.5,1)
Data[Data$cens<Data$Tstop,"stat"] <- 0
Data$Tstop <- pmin(Data$Tstop,Data$cens)
Data$tstat <- ifelse(Data$Tstart > Data$Tstop, 1, 0)
Data.weight <- crprep(Data$Tstop, Data$stat, trans=c(1,2))
# calculate cause-specific cumulative incidence, no truncation,
# compare with Cuminc (also from mstate)
ci <- Cuminc(Data$Tstop, Data$stat)
plot(ci$time,ci$CI.1,type="s",lwd=3,col="black",ylim=c(0,0.5),
xlab="Time",ylab="Cumulative incidence")
lines(ci$time,ci$CI.2,type="s",lwd=3,col="red")
lines(ci$time,ci$CI.1-qnorm(0.975)*ci$seCI.1,type="s",lty=3)
lines(ci$time,ci$CI.1+qnorm(0.975)*ci$seCI.1,type="s",lty=3)
lines(survfit(Surv(Tstart,Tstop,status==1)~1,data=Data.weight,
weight=weight.cens,subset=failcode==1),
fun="event",col="lightblue",lwd=1,mark.time=FALSE)
lines(survfit(Surv(Tstart,Tstop,status==1)~1,data=Data.weight,
weight=weight.cens,subset=failcode==1),
fun="event",col="lightblue",mark.time=FALSE,conf.int="only",lty=2)
# Proportional hazards regression on subdistribution and cause-specific hazard,
# with truncation
Data.weight.trunc <- crprep("Tstop", "stat", Tstart="Tstart",
data=subset(Data,tstat==0), trans=1:2, keep=c("Z1","Z2"))
coxph(Surv(Tstart,Tstop,status==1)~Z1+Z2,data=Data.weight.trunc,
weight=weight.cens*weight.trunc,subset=failcode==1) #cause 1
# Both end points, assume effect Z2 same for both
coxph(Surv(Tstart,Tstop,status==1)~strata(failcode)*Z1+Z2,
data=Data.weight.trunc,weight=weight.cens*weight.trunc)
# Cause-specific hazard
coxph(Surv(Tstart,Tstop,stat==1)~Z1+Z2,data=subset(Data, tstat==0))
coxph(Surv(Tstart,Tstop,status==1)~Z1+Z2,data=Data.weight.trunc,
subset=failcode==1&count==1)
data(ebmt2)
ebmt2.long <- crprep("time", "status", data=ebmt2, trans=1:6,
keep=c("dissub", "match", "tcd", "year", "age"))
# ebmt2.long <- with(ebmt2, crprep(time, cod, trans=levels(cod)[-1], cens="Alive",
# keep=c("dissub", "match", "tcd", "year", "age")))
plot(survfit(Surv(time, status>0)~1, data=ebmt2, etype=cod), mark.time=FALSE,
col=2:7, fun="event",lwd=3)
lines(survfit(Surv(Tstart, Tstop, failcode==status) ~ failcode, data=ebmt2.long,
weight=weight.cens), lwd=2, lty=2, fun="event", mark.time=FALSE)Run the code above in your browser using DataLab