demo/demo.R

library(PIcompete)

#the parameters for the general transformation function
#when trans.r1=0, the incidence model is proportional hazards model.
trans.r1<-0
trans.r2<-0

#prevalence model
#HR_status: binary (0.5=positive, -0.5: negative)
#risk_score: continuous
model0<-"C ~ HR_status + risk_score"

#incidence model for event 1
model1<-"C + L1 + R1 ~ HR_status + risk_score"

#incidence model for event 2
model2<-"L2 + R2 ~ HR_status + risk_score"

#the number of parameters for the prevalence and incidence models
n.beta<-3
n.gamma1<-2
n.gamma2<-2

#fitting the models to the example data
output <- ipw.pi.competing(Data=samp.data, p.model=model0, i.model1=model1, i.model2=model2, trans.r1=trans.r1, trans.r2=trans.r2, n.beta=n.beta, n.gamma1=n.gamma1, n.gamma2=n.gamma2)  #ipw.pi.competing

#two-phase weighted bootstrap (sequential computing) with 2 replications 
outputb <- ipw.bootstrap.ph2(n.boot=2, Data=samp.data, p.model=model0, i.model1=model1, i.model2=model2, trans.r1=trans.r1, trans.r2=trans.r2, n.beta=n.beta, n.gamma1=n.gamma1, n.gamma2=n.gamma2, iteration.limit=250, time.interval=0.1, time.list=seq(0,9,by=0.1))  #bootstrap

#parallel computing for two-phase sample weighted bootstrap with 2 replications
outputbp <- ipw.bootstrap.ph2(n.boot=2, Data=samp.data, p.model=model0, i.model1=model1, i.model2=model2, trans.r1=trans.r1, trans.r2=trans.r2, n.beta=n.beta, n.gamma1=n.gamma1, n.gamma2=n.gamma2, iteration.limit=250, time.interval=0.1, time.list=seq(0,9,by=0.1),parallel=TRUE,n.core=2)  #bootstrap parallel

#design matrix for prediction 
x <- cbind(X0=1, samp.data[,c("HR_status","risk_score")])  #add 1 as the intercept col
x <- data.matrix(x)  #design matrix
time <- c(1,3,5)  #timepoints

#summary from the empirical distributions generated by the bootstrap procedure
sum1 <- bootstrap.summary(o.input = output, b.input = outputb, p.mat = x, i.mat1 = x[,-1], i.mat2 = x[,-1], time.points = time)  #summary
sum2 <- bootstrap.summary(o.input = output, b.input = outputbp, p.mat = x, i.mat1 = x[,-1], i.mat2 = x[,-1], time.points = time)  #summary

#prediction 
pred <- pi.pred(input = output, p.mat = x, i.mat1 = x[,-1], i.mat2 = x[,-1], time.points = time)  #predict
xiaoli-mcw/PIcompete documentation built on May 20, 2020, 7:44 p.m.