# NOT RUN {
#############################################################
# Time to Cosmetic Deterioration of Breast Cancer Patients
#############################################################
data(deterioration)
attach(deterioration)
y <- cbind(left,right)
# Design matrix
x <- cbind(rep(1,length(trt)),trt)
xtf <- cbind(rep(1,length(trt)),trt)
colnames(x) <- c("(Intercept)","trt")
colnames(xtf) <- c("(Intercept)","trt")
# Prediction
xdenpred <- cbind(rep(1,2),c(0,1))
xtfdenpred <- cbind(rep(1,2),c(0,1))
xmedpred <- cbind(rep(1,2),c(0,1))
xtfmedpred <- cbind(rep(1,2),c(0,1))
prediction <- list(xdenpred=xdenpred,
xtfdenpred=xtfdenpred,
xmedpred=xmedpred,
xtfmedpred=xtfmedpred,
quans=c(0.03,0.50,0.97))
# Prior information
prior <- list(maxm=5,
a0=1,
b0=1,
mub=rep(0,2),
Sb=diag(1000,2),
tau1=2,002,
tau2=2.002)
# Initial state
state <- NULL
# MCMC parameters
mcmc <- list(nburn=5000,
nsave=5000,
nskip=4,
ndisplay=200)
# Fitting the model
fit1 <- LDTFPsurvival(y=y,
x=x,
xtf=xtf,
prediction=prediction,
prior=prior,
mcmc=mcmc,
state=state,
grid=seq(0.01,70,len=200),
status=TRUE,
compute.band=TRUE)
fit1
summary(fit1)
plot(fit1)
# Plotting survival functions estimates
par(cex=1.7,mar=c(4.1, 4.1, 1, 1))
x1 <- fit1$grid
y1 <- fit1$survml[1,]
x2 <- fit1$grid
y2 <- fit1$survmu[1,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(fit1$grid,fit1$survmu[1,],type="l",
xlab="months",ylab="survival",
lty=1,lwd=2,ylim=c(0,1),col="lightgray")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(fit1$grid,fit1$survmm[1,],lty=1,lwd=3)
par(cex=1.7,mar=c(4.1, 4.1, 1, 1))
x1 <- fit1$grid
y1 <- fit1$survml[2,]
x2 <- fit1$grid
y2 <- fit1$survmu[2,]
aa <- rbind(x2,y2)[, order(-x2, y2)]
x2 <- aa[1,]
y2 <- aa[2,]
plot(fit1$grid,fit1$survmu[2,],type="l",
xlab="months",ylab="survival",
lty=1,lwd=2,ylim=c(0,1),col="lightgray")
polygon(x=c(x1,x2),y=c(y1,y2),border=NA,col="lightgray")
lines(fit1$grid,fit1$survmm[2,],lty=1,lwd=3)
# }
Run the code above in your browser using DataLab