timereg (version 2.0.1)

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,
  n = NULL,
  entry = NULL,
  type = "cloglog",
  startcum = c(0, 0),
  ...
)

Arguments

cumhazard

matrix that specified times and values of some cumulative hazard.

rr

"relative risk" terms

n

number of simulation if rr not given

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

# 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,n=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=1000
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)
# }