# NOT RUN {
data(TRACE)
## Logit link for proportional odds model, using comp.risk to save time
#' cif <- prop.odds.subdist(Event(time,status)~vf+chf+wmi,data=TRACE,cause=9)
cif <- comp.risk(Event(time,status)~const(vf)+const(chf)+const(wmi),
data=TRACE,cause=9,model="logistic2")
sim1 <- sim.cif(cif,500,data=TRACE)
#' cc <- prop.odds.subdist(Event(time,status)~vf+chf+wmi,data=sim1,cause=1)
cc <- comp.risk(Event(time,status)~const(vf)+const(chf)+const(wmi),
data=sim1,cause=1,model="logistic2")
cbind(cif$gamma,cc$gamma)
plot(cif)
lines(cc$cum)
#################################################################
## Fine-Gray model model, using comp.risk to avoid dependcies
#################################################################
cif <- comp.risk(Event(time,status)~const(vf)+const(chf)+const(wmi),
data=TRACE,cause=9)
sim1 <- sim.cif(cif,500,data=TRACE)
#' cc <- crr
cc <- comp.risk(Event(time,status)~const(vf)+const(chf)+const(wmi),
data=sim1,cause=1)
cbind(cif$gamma,cc$gamma)
plot(cif)
lines(cc$cum)
# faster/better with mets package
# library(mets)
# scif <- cifreg(Event(time,status)~vf+chf+wmi,data=sim1,cause=1,prop=NULL)
#
# plot(scif$cum,type="l")
# lines(cif$cum,col=2)
# cbind(cif$gamma,scif$coef)
#
################################################################
# simulating several causes with specific cumulatives
################################################################
data(bmt)
cif1 <- comp.risk(Event(time,cause)~const(tcell)+const(age),
data=bmt,cause=1,model="logistic2")
cif2 <- comp.risk(Event(time,cause)~const(tcell)+const(age),
data=bmt,cause=2,model="logistic2")
## must look at same time-scale
cifs <- pre.cifs(list(cif1,cif2))
plot(cifs[[1]]$cum,type="l")
lines(cifs[[2]]$cum,col=2)
legend("topleft",c("cause1","cause2"),lty=1,col=1:2)
n <- 500
sim1 <- sim.cif(cifs[[1]],n,data=bmt)
Z <- sim1[,c("tcell","age")]
sim2 <- sim.cif(cifs[[2]],n,data=bmt,Z=Z,drawZ=FALSE)
###
rt <- rbinom(n,1,(sim1$F1tau+sim2$F1tau))
rb <- rbinom(n,1,sim1$F1tau/(sim1$F1tau+sim2$F1tau))
cause=ifelse(rb==1,1,2)
time=ifelse(cause==1,sim1$timecause,sim2$timecause)
cause <- rt*cause
time[cause==0] <- tail(cifs[[1]]$cum[,1],1)
bt <- data.frame(time=time,cause=cause,tcell=sim1$tcell,age=sim1$age)
scif1 <- comp.risk(Event(time,cause)~const(tcell)+const(age),
data=bt,cause=1,model="logistic2")
scif2 <- comp.risk(Event(time,cause)~const(tcell)+const(age),
data=bt,cause=2,model="logistic2")
plot(scif1$cum,type="l")
lines(scif2$cum,col=1,lty=2)
legend("topleft",c("cause1","cause2"),lty=1:2,col=1:1)
lines(cifs[[1]]$cum,col=2)
lines(cifs[[2]]$cum,col=2,lty=2)
# Everyhing wrapped in call assuming covariates work in the same way for two models
dd <- sim.cifs(list(cif1,cif2),2000,data=bmt)
scif1 <- comp.risk(Event(time,cause)~const(tcell)+const(age),
data=dd,cause=1,model="logistic2")
scif2 <- comp.risk(Event(time,cause)~const(tcell)+const(age),
data=dd,cause=2,model="logistic2")
plot(scif1$cum,type="l")
lines(scif2$cum,col=1,lty=2)
legend("topleft",c("cause1","cause2"),lty=1:2,col=1:1)
lines(cifs[[1]]$cum,col=2)
lines(cifs[[2]]$cum,col=2,lty=2)
# Everyhing wrapped in call assuming covariates work in the same way for two models
# but now draws cif1 to be of correct model, but model 2 is adapted
#(if needed) to make constraints satisfied F1+F2 <=1
# see doubleFG of mets package for paramtrization
# and drawns as "if not cause1" then distribute according to cause 2
# dd <- sim.cifsRestrict(list(cif1,cif2),2000,data=bmt)
# faster with mets package
# dd <- sim.cifs(list(cif1,cif2),1000,data=bmt)
# scif1 <- cifreg(Event(time,cause)~tcell+age,data=dd,cause=1)
# scif2 <- cifreg(Event(time,cause)~tcell+age,data=dd,cause=2)
#
# plot(scif1$cum,type="l")
# legend("topleft",c("cause1","cause2"),lty=1:2,col=1:1)
# lines(cifs[[1]]$cum,col=2)
# lines(cifs[[2]]$cum,col=2,lty=2)
#
# }
Run the code above in your browser using DataLab