#####################################################
#### prediction on a COX or SHARED frailty model ####
#####################################################
data(readmission)
#-- here is a generated cluster (31 clusters of 13 subjects)
readmission <- transform(readmission,group=id%%31+1)
#-- we compute predictions of death
#-- we extract last row of each subject for the time of death
readmission <- aggregate(readmission,by=list(readmission$id),
FUN=function(x){x[length(x)]})[,-1]
##-- predictions on a Cox proportionnal hazard model --##
cox <- frailtyPenal(Surv(t.stop,death)~sex+dukes,
n.knots=10,kappa=10000,data=readmission)
#-- construction of the dataframe for predictions
datapred <- data.frame(sex=0,dukes=0)
datapred$sex <- as.factor(datapred$sex)
levels(datapred$sex)<- c(1,2)
datapred$dukes <- as.factor(datapred$dukes)
levels(datapred$dukes)<- c(1,2,3)
datapred[1,] <- c(1,2) # man, dukes 2
datapred[2,] <- c(2,3) # woman, dukes 3
#-- prediction of death for two patients between 100 and 100+w,
#-- with w in (50,100,...,1900)
pred.cox <- prediction(cox,datapred,100,seq(50,1900,50))
plot(pred.cox)
#-- prediction of death for two patients between t and t+400,
#-- with t in (100,150,...,1500)
pred.cox2 <- prediction(cox,datapred,seq(100,1500,50),400)
plot(pred.cox2)
##-- predictions on a shared frailty model for clustered data --##
sha <- frailtyPenal(Surv(t.stop,death)~cluster(group)+sex+dukes,
n.knots=10,kappa=10000,data=readmission)
#-- marginal prediction
pred.sha.marg <- prediction(sha,datapred,100,seq(50,1900,50))
plot(pred.sha.marg)
#-- conditional prediction, given a specific cluster (group=5)
pred.sha.cond <- prediction(sha,datapred,100,seq(50,1900,50),group=5)
plot(pred.sha.cond)
#####################################################
######## prediction on a JOINT frailty model ########
#####################################################
data(readmission)
##-- predictions of death on a joint model --##
joi <- frailtyPenal(Surv(t.start,t.stop,event)~cluster(id)
+sex+dukes+terminal(death),formula.terminalEvent=~sex
+dukes,data=readmission,n.knots=10,kappa=c(100,100),recurrentAG=TRUE)
#-- construction of the dataframe for predictions
datapredj <- data.frame(t.stop=0,event=0,id=0,sex=0,dukes=0)
datapredj$sex <- as.factor(datapredj$sex)
levels(datapredj$sex) <- c(1,2)
datapredj$dukes <- as.factor(datapredj$dukes)
levels(datapredj$dukes) <- c(1,2,3)
datapredj[1,] <- c(100,1,1,1,2)
datapredj[2,] <- c(200,1,1,1,2)
datapredj[3,] <- c(300,1,1,1,2)
datapredj[4,] <- c(400,1,1,1,2)
datapredj[5,] <- c(380,1,2,1,2)
#-- prediction of death between 100 and 100+500 given relapses
pred.joint0 <- prediction(joi,datapredj,100,500)
print(pred.joint0)
#-- prediction of death between 100 and 100+w given relapses (with confidence intervals)
pred.joint <- prediction(joi,datapredj,100,seq(50,1500,50),MC.sample=100)
plot(pred.joint,conf.bands=TRUE)
# each y-value of the plot corresponds to the prediction between [100,x]
#-- prediction of death between t and t+500 given relapses
pred.joint2 <- prediction(joi,datapredj,seq(100,1000,50),500)
plot(pred.joint2)
# each y-value of the plot corresponds to the prediction between [x,x+500], or in the next 500
Run the code above in your browser using DataLab