timereg (version 1.9.3)

simsubdist: Simulation from subdistribution function assuming piecwise linearity

Description

Simulation from subdistribution function assuming piecwise linearity for Fine-Gray or logistic link.

Usage

simsubdist(cumhazard, rr, entry = NULL, type = "cloglog",
  startcum = c(0, 0), ...)

Arguments

cumhazard

matrix that specified times and values of some cumulative hazard.

rr

"relative risk" terms

entry

not implemented yet

type

either cloglog or logistic

startcum

c(0,0) to make cumulativ start in 0 with 0 cumhazard.

...

further arguments

Examples

Run this code
# NOT RUN {
data(sTRACE)

cif <- comp.risk(Event(time,status)~const(vf),data=sTRACE,cause=9,model="logistic2")
cumhaz <- cif$cum

## 1000 logistic without covariates with baseline from model fit  
sim1 <- simsubdist(cumhaz,1000,type="logistic")
###
cifs <- comp.risk(Event(time,status)~+1,data=sim1,cause=1,model="logistic2")
###
plot(cifs)
lines(cifs$cum,col=2)

## 1000 logistic with covariates with baseline from model fit  
x <- rbinom(1000,1,0.5)
rr <- exp(x*0.3)
sim1 <- simsubdist(cumhaz,rr,type="logistic")
sim1$x <- x

cifs <- comp.risk(Event(time,status)~+const(x),data=sim1,cause=1,model="logistic2")
###
cifs$gamma
plot(cifs)
lines(cumhaz,col=2)

##################################################################
### simulation of cumulative incidence with specified functions
##################################################################
F1logit<-function(t,lam0=0.2,beta=0.3,x=0) 
{ pt <- t*lam0; rr <- exp(x*beta); return(pt*rr/(1+pt*rr)); } 

F1p<-function(t,lam0=0.4,beta=0.3,x=0) # proportional version 
{ return( 1 - exp(-(t*lam0)*exp(x*beta))) }

n=10000
tt=seq(0,3,by=.01)
tt=seq(0,3,by=.01)
t1 <- invsubdist(cbind(tt,F1p(tt)),runif(n))
t2 <- invsubdist(cbind(tt,F1p(tt,lam0=0.1)),runif(n))
rt <- rbinom(n,1,(F1p(3)+F1p(3,lam0=0.1)))
rb <- rbinom(n,1,F1p(3)/(F1p(3)+F1p(3,lam0=0.1)))
cause=ifelse(rb==1,1,2)
time=ifelse(cause==1,t1$time,t2$time)
cause <- rt*cause
time[cause==0] <- 3

datC=data.frame(time=time,cause=cause)
p1=comp.risk(Event(time,cause)~+1,data=datC,cause=1)
p2=comp.risk(Event(time,cause)~+1,data=datC,cause=2)
pp1=predict(p1,X=1,se=0)
pp2=predict(p2,X=1,se=0)
par(mfrow=c(1,2))
plot(pp1)
lines(tt,F1p(tt),col=2)
plot(pp2)
lines(tt,F1p(tt,lam0=0.1),col=2)

#to avoid dependencies when checking 
#library(prodlim)
#pp=prodlim(Hist(time,cause)~+1)
#par(mfrow=c(1,2))
#plot(pp,cause="1")
#lines(tt,F1p(tt),col=2)
#plot(pp,cause="2")
#lines(tt,F1p(tt,lam0=0.1),col=2)
# }

Run the code above in your browser using DataCamp Workspace