BayesID: The function to implement Bayesian parametric and semi-parametric analyses for semi-competing risks data.

Description

Independent/cluster-correlated semi-competing risks data can be analyzed using hierarchical models. The priors for baseline hazard functions can be specified by either parametric (Weibull) model or non-parametric mixture of piecewise exponential models (PEM). The option to choose between parametric multivariate normal (MVN) and non-parametric Dirichlet process mixture of multivariate normals (DPM) is available for the prior of cluster-specific random effects distribution. The conditional hazard function for time to the terminal event given time to non-terminal event can be modeled based on either Markov (it does not depend on the timing of the non-terminal event) or semi-Markov assumption (it does depend on the timing).

Usage

BayesID(Y, lin.pred, data, cluster=NULL, model=c("semi-Markov", "Weibull"), 
        hyperParams, startValues, mcmc, path=NULL)

Arguments

a data.frame containing semi-competing risks outcomes from n subjects. It is of dimension $n\times 4$: the columns correspond to $y_1$, $\delta_1$, $y_2$, $\delta_2$.

lin.pred

a list containing three formula objects that correspond to $h_g$, $g$=1,2,3.

data

a data.frame in which to interpret the variables named in the formulas in lin.pred.

cluster

a vector of cluster information for n subjects. The cluster membership must be consecutive positive integers, $1:J$.

model

a character vector that specifies the type of components in a model. The first element is for the assumption on $h_3$: "semi-Markov" or "Markov". The second element is for the specification of baseline hazard functions: "Weibull" or "PEM". The third

hyperParams

a list containing lists or vectors for hyperparameter values in hierarchical models. Components include, theta (a numeric vector for hyperparameter in the prior of subject-specific frailty variance component), WB (a list cont

startValues

a list containing vectors of starting values for model parameters. It can be specified as the object returned by the function initiate.startValues.

mcmc

a list containing variables required for MCMC sampling. Components include, run (a list containing numeric values for setting for the overall run: numReps, total number of scans; thin, extent of thinning; bu

path

the name of directory where the results are saved.

Value

BayesID returns an object of class Bayes.

Details

We view the semi-competing risks data as arising from an underlying illness-death model system in which individuals may undergo one or more of three transitions: 1) from some initial condition to non-terminal event, 2) from some initial condition to terminal event, 3) from non-terminal event to terminal event. Let $t_{ji1}$, $t_{ji2}$ denote time to non-terminal and terminal event from subject $i=1,...,n_j$ in cluster $j=1,...,J$. The system of transitions is modeled via the specification of three hazard functions: $$h_1(t_{ji1} | \gamma_{ji}, x_{ji1}, V_{j1}) = \gamma_{ji} h_{01}(t_{ji1})\exp(x_{ji1}^{\top}\beta_1 +V_{j1}), t_{ji1}>0,$$ $$h_2(t_{ji2} | \gamma_{ji}, x_{ji2}, V_{j2}) = \gamma_{ji} h_{02}(t_{ji2})\exp(x_{ji2}^{\top}\beta_2 +V_{j2}), t_{ji2}>0,$$ $$h_3(t_{ji2} | t_{ji1}, \gamma_{ji}, x_{ji3}, V_{j3}) = \gamma_{ji} h_{03}(t_{ji2})\exp(x_{ji3}^{\top}\beta_3 +V_{j3}), 0

References

Lee, K. H., Haneuse, S., Schrag, D., and Dominici, F. (2015), Bayesian semiparametric analysis of semicompeting risks data: investigating hospital readmission after a pancreatic cancer diagnosis, Journal of the Royal Statistical Society: Series C, 64, 2, 253-273. Lee, K. H., Dominici, F., Schrag, D., and Haneuse, S. (2016), Hierarchical models for semicompeting risks data with application to quality of end-of-life care for pancreatic cancer, Journal of the American Statistical Association, in press.

Examples

Run this code

# loading a data set
data(scrData)
Y <- scrData[,c(1,2,3,4)]
cluster <- scrData[,5]
form1 <- as.formula( ~ x1 + x2 + x3)
form2 <- as.formula( ~ x1 + x2)
form3 <- as.formula( ~ x1 + x2)
lin.pred <- list(form1, form2, form3)

#####################
## Hyperparameters ##
#####################

## Subject-specific frailty variance component
##  - prior parameters for 1/theta
##
theta.ab <- c(0.7, 0.7)

## Weibull baseline hazard function: alphas, kappas
##
WB.ab1 <- c(0.5, 0.01) # prior parameters for alpha1
WB.ab2 <- c(0.5, 0.01) # prior parameters for alpha2
WB.ab3 <- c(0.5, 0.01) # prior parameters for alpha3
##
WB.cd1 <- c(0.5, 0.05) # prior parameters for kappa1
WB.cd2 <- c(0.5, 0.05) # prior parameters for kappa2
WB.cd3 <- c(0.5, 0.05) # prior parameters for kappa3

## PEM baseline hazard function
##
PEM.ab1 <- c(0.7, 0.7) # prior parameters for 1/sigma_1^2
PEM.ab2 <- c(0.7, 0.7) # prior parameters for 1/sigma_2^2
PEM.ab3 <- c(0.7, 0.7) # prior parameters for 1/sigma_3^2
##
PEM.alpha1 <- 10 # prior parameters for K1
PEM.alpha2 <- 10 # prior parameters for K2
PEM.alpha3 <- 10 # prior parameters for K3

## MVN cluster-specific random effects
##
Psi_v <- diag(1, 3)
rho_v <- 100

## DPM cluster-specific random effects
##
Psi0  <- diag(1, 3)
rho0  <- 10
aTau  <- 1.5
bTau  <- 0.0125

##
hyperParams <- list(theta=theta.ab,
                WB=list(WB.ab1=WB.ab1, WB.ab2=WB.ab2, WB.ab3=WB.ab3,
                       WB.cd1=WB.cd1, WB.cd2=WB.cd2, WB.cd3=WB.cd3),
                   PEM=list(PEM.ab1=PEM.ab1, PEM.ab2=PEM.ab2, PEM.ab3=PEM.ab3,
                       PEM.alpha1=PEM.alpha1, PEM.alpha2=PEM.alpha2, PEM.alpha3=PEM.alpha3),
                   MVN=list(Psi_v=Psi_v, rho_v=rho_v),
                   DPM=list(Psi0=Psi0, rho0=rho0, aTau=aTau, bTau=bTau))
                    
###################
## MCMC SETTINGS ##
###################

## Setting for the overall run
##
numReps    <- 2000
thin       <- 10
burninPerc <- 0.5

## Settings for storage
##
nGam_save <- 0
storeV    <- rep(TRUE, 3)

## Tuning parameters for specific updates
##
##  - those common to all models
mhProp_theta_var  <- 0.05
mhProp_Vg_var     <- c(0.05, 0.05, 0.05)
##
## - those specific to the Weibull specification of the baseline hazard functions
mhProp_alphag_var <- c(0.01, 0.01, 0.01)
##
## - those specific to the PEM specification of the baseline hazard functions
Cg        <- c(0.2, 0.2, 0.2)
delPertg  <- c(0.5, 0.5, 0.5)
rj.scheme <- 1
Kg_max    <- c(50, 50, 50)
sg_max    <- c(max(Y$time1[Y$event1 == 1]),
               max(Y$time2[Y$event1 == 0 & Y$event2 == 1]),
               max(Y$time2[Y$event1 == 1 & Y$event2 == 1]))

time_lambda1 <- seq(1, sg_max[1], 1)
time_lambda2 <- seq(1, sg_max[2], 1)
time_lambda3 <- seq(1, sg_max[3], 1)               

##
mcmc.WB  <- list(run=list(numReps=numReps, thin=thin, burninPerc=burninPerc),
                storage=list(nGam_save=nGam_save, storeV=storeV),
                tuning=list(mhProp_theta_var=mhProp_theta_var,
                mhProp_Vg_var=mhProp_Vg_var, mhProp_alphag_var=mhProp_alphag_var))

##
mcmc.PEM <- list(run=list(numReps=numReps, thin=thin, burninPerc=burninPerc),
                storage=list(nGam_save=nGam_save, storeV=storeV),
                tuning=list(mhProp_theta_var=mhProp_theta_var,
                mhProp_Vg_var=mhProp_Vg_var, Cg=Cg, delPertg=delPertg,
                rj.scheme=rj.scheme, Kg_max=Kg_max,
                time_lambda1=time_lambda1, time_lambda2=time_lambda2,
                time_lambda3=time_lambda3))
    
#####################
## Starting Values ##
#####################

##
Sigma_V <- diag(0.1, 3)
Sigma_V[1,2] <- Sigma_V[2,1] <- -0.05
Sigma_V[1,3] <- Sigma_V[3,1] <- -0.06
Sigma_V[2,3] <- Sigma_V[3,2] <- 0.07

#################################################################
## Analysis of Independent Semi-competing risks data ############
#################################################################

#############
## WEIBULL ##
#############

##
myModel <- c("semi-Markov", "Weibull")
myPath  <- "Output/01-Results-WB/"

startValues      <- vector("list", 2)
startValues[[1]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel)
startValues[[2]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, theta = 0.23)
                
##
fit_WB <- BayesID(Y, lin.pred, scrData, cluster=NULL, model=myModel, 
                hyperParams, startValues, mcmc.WB, path=myPath)
				
fit_WB
summ.fit_WB <- summary(fit_WB); names(summ.fit_WB)
summ.fit_WB
plot(fit_WB, tseq=seq(from=0, to=30, by=5))
plot(fit_WB, tseq=seq(from=0, to=30, by=5), plot.est = "BH")
names(fit_WB.plot <- plot(fit_WB, tseq=seq(0, 30, 5), plot=FALSE))

#########
## PEM ##
#########

##						
myModel <- c("semi-Markov", "PEM")
myPath  <- "Output/02-Results-PEM/"

startValues      <- vector("list", 2)
startValues[[1]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel)
startValues[[2]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, theta = 0.23)
                                                           
##
fit_PEM <- BayesID(Y, lin.pred, scrData, cluster=NULL, model=myModel, 
                 hyperParams, startValues, mcmc.PEM, path=myPath)
				
fit_PEM
summ.fit_PEM <- summary(fit_PEM); names(summ.fit_PEM)
summ.fit_PEM
plot(fit_PEM)
plot(fit_PEM, plot.est = "BH")
names(fit_PEM.plot <- plot(fit_PEM, plot=FALSE))
					
#################################################################
## Analysis of Correlated Semi-competing risks data #############
#################################################################

#################
## WEIBULL-MVN ##
#################

##
myModel <- c("semi-Markov", "Weibull", "MVN")
myPath  <- "Output/03-Results-WB_MVN/"

startValues      <- vector("list", 2)
startValues[[1]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, cluster)
startValues[[2]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, cluster,
MVN.SigmaV=Sigma_V)

##
fit_WB_MVN <- BayesID(Y, lin.pred, scrData, cluster, model=myModel,
                    hyperParams, startValues, mcmc.WB, path=myPath)
                    
fit_WB_MVN
summ.fit_WB_MVN <- summary(fit_WB_MVN); names(summ.fit_WB_MVN)
summ.fit_WB_MVN
plot(fit_WB_MVN, tseq=seq(from=0, to=30, by=5))
plot(fit_WB_MVN, tseq=seq(from=0, to=30, by=5), plot.est = "BH")
names(fit_WB_MVN.plot <- plot(fit_WB_MVN, tseq=seq(0, 30, 5), plot=FALSE))                    

#################
## WEIBULL-DPM ##
#################

##
myModel <- c("semi-Markov", "Weibull", "DPM")
myPath  <- "Output/04-Results-WB_DPM/"

startValues      <- vector("list", 2)
startValues[[1]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, cluster)
startValues[[2]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, cluster,
                    MVN.SigmaV=Sigma_V)

##
fit_WB_DPM <- BayesID(Y, lin.pred, scrData, cluster, model=myModel,
                    hyperParams, startValues, mcmc.WB, path=myPath)

fit_WB_DPM
summ.fit_WB_DPM <- summary(fit_WB_DPM); names(summ.fit_WB_DPM)
summ.fit_WB_DPM
plot(fit_WB_DPM, tseq=seq(from=0, to=30, by=5))
plot(fit_WB_DPM, tseq=seq(from=0, to=30, by=5), plot.est = "BH")
names(fit_WB_DPM.plot <- plot(fit_WB_DPM, tseq=seq(0, 30, 5), plot=FALSE))                     

#############
## PEM-MVN ##
#############

##
myModel <- c("semi-Markov", "PEM", "MVN")
myPath  <- "Output/05-Results-PEM_MVN/"

startValues      <- vector("list", 2)
startValues[[1]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, cluster)
startValues[[2]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, cluster,
                    MVN.SigmaV=Sigma_V)

##
fit_PEM_MVN <- BayesID(Y, lin.pred, scrData, cluster, model=myModel,
                    hyperParams, startValues, mcmc.PEM, path=myPath)
                    
fit_PEM_MVN
summ.fit_PEM_MVN <- summary(fit_PEM_MVN); names(summ.fit_PEM_MVN)
summ.fit_PEM_MVN
plot(fit_PEM_MVN)
plot(fit_PEM_MVN, plot.est = "BH")
names(fit_PEM_MVN.plot <- plot(fit_PEM_MVN, plot=FALSE))                    

#############
## PEM-DPM ##
#############

##
myModel <- c("semi-Markov", "PEM", "DPM")
myPath  <- "Output/06-Results-PEM_DPM/"

startValues      <- vector("list", 2)
startValues[[1]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, cluster)
startValues[[2]] <- initiate.startValues(Y, lin.pred, scrData, model=myModel, cluster,
                    MVN.SigmaV=Sigma_V)
                    
##
fit_PEM_DPM <- BayesID(Y, lin.pred, scrData, cluster, model=myModel,
                    hyperParams, startValues, mcmc.PEM, path=myPath)
                    
fit_PEM_DPM
summ.fit_PEM_DPM <- summary(fit_PEM_DPM); names(summ.fit_PEM_DPM)
summ.fit_PEM_DPM
plot(fit_PEM_DPM)
plot(fit_PEM_DPM, plot.est = "BH")
names(fit_PEM_DPM.plot <- plot(fit_PEM_DPM, plot=FALSE))

Run the code above in your browser using DataLab