# The function to implement Bayesian parametric and semi-parametric analyses for semi-competing risks data in the context of hazard regression (HReg) models. Independent/cluster-correlated semi-competing risks data can be analyzed using HReg models that have a hierarchical structure. 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). ```r BayesID_HReg(Formula, data, id=NULL, model=c("semi-Markov", "Weibull"), hyperParams, startValues, mcmcParams, na.action = "na.fail", subset=NULL, path=NULL) ``` ## Arguments - `Formula`: a `Formula` object, with the outcome on the left of a $\sim$, and covariates on the right. It is of the form, time to non-terminal event + corresponding censoring indicator | time to terminal event + corresponding censoring indicator $\sim$ covariates for $h_1$ | covariates for $h_2$ | covariates for $h_3$: i.e., $y_1$+$\delta_1$ | $y_2$+$\delta_2$ ~ $x_1$ | $x_2$ | $x_3$. - `data`: a data.frame in which to interpret the variables named in `Formula`. - `id`: 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 element needs to be set only for clustered semi-competing risks data and is for the specification of cluster-specific random effects distribution: "MVN" or "DPM". - `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 containing numeric vectors for Weibull hyperparameters: `WB.ab1`, `WB.ab2`, `WB.ab3`, `WB.cd1`, `WB.cd2`, `WB.cd3`), `PEM` (a list containing numeric vectors for PEM hyperparameters: `PEM.ab1`, `PEM.ab2`, `PEM.ab3`, `PEM.alpha1`, `PEM.alpha2`, `PEM.alpha3`). Models for clustered data require additional components, `MVN` (a list containing numeric vectors for MVN hyperparameters: `Psi_v`, `rho_v`), `DPM` (a list containing numeric vectors for DPM hyperparameters: `Psi0`, `rho0`, `aTau`, `bTau`). See Details and Examples below. - `startValues`: a list containing vectors of starting values for model parameters. It can be specified as the object returned by the function `initiate.startValues_HReg`. - `mcmcParams`: 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; `burninPerc`, the proportion of burn-in). `storage` (a list containing numeric values for storing posterior samples for subject- and cluster-specific random effects: `nGam_save`, the number of $\gamma$ to be stored; `storeV`, a vector of three logical values to determine whether all the posterior samples of $V_1$, $V_2$, $V_3$ are to be stored). `tuning` (a list containing numeric values relevant to tuning parameters for specific updates in Metropolis-Hastings-Green (MHG) algorithm: `mhProp_theta_var`, the variance of proposal density for $\theta$; `mhProp_Vg_var`, the variance of proposal density for $V_g$ in DPM models; `mhProp_alphag_var`, the variance of proposal density for $\alpha_g$ in Weibull models; `Cg`, a vector of three proportions that determine the sum of probabilities of choosing the birth and the death moves in PEM models. The sum of the three elements should not exceed 0.6; `delPertg`, the perturbation parameters in the birth update in PEM models. The values must be between 0 and 0.5; If `rj.scheme`=1, the birth update will draw the proposal time split from $1:s_{max}$. If `rj.scheme`=2, the birth update will draw the proposal time split from uniquely ordered failure times in the data. Only required for PEM models; `Kg_max`, the maximum number of splits allowed at each iteration in MHG algorithm for PEM models; `time_lambda1`, `time_lambda2`, `time_lambda3` - time points at which the log-hazard functions are calculated for `predict.Bayes_HReg`, Only required for PEM models). See Details and Examples below. - `na.action`: how NAs are treated. See `model.frame`. - `subset`: a specification of the rows to be used: defaults to all rows. See `model.frame`. - `path`: the name of directory where the results are saved. ## 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 ## Note The posterior samples of $\gamma$ and $V_g$ are saved separately in `working directory/path`. For a dataset with large $n$, `nGam_save` should be carefully specified considering the system memory and the storage capacity. ## See Also `initiate.startValues_HReg`, `print.Bayes_HReg`, `summary.Bayes_HReg`, `predict.Bayes_HReg` ## Examples ```r ## Not run: # loading a data set data(scrData) id=scrData$cluster form <- Formula(time1 + event1 | time2 + event2 ~ x1 + x2 + x3 | x1 + x2 | x1 + x2) ##################### ## 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(scrData$time1[scrData$event1 == 1]), max(scrData$time2[scrData$event1 == 0 & scrData$event2 == 1]), max(scrData$time2[scrData$event1 == 1 & scrData$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 <- initiate.startValues_HReg(form, scrData, model=myModel, nChain=2) ## fit_WB <- BayesID_HReg(form, scrData, id=NULL, model=myModel, hyperParams, startValues, mcmc.WB, path=myPath) fit_WB summ.fit_WB <- summary(fit_WB); names(summ.fit_WB) summ.fit_WB pred_WB <- predict(fit_WB, tseq=seq(from=0, to=30, by=5)) plot(pred_WB, plot.est="Haz") plot(pred_WB, plot.est="Surv") ######### ## PEM ## ######### ## myModel <- c("semi-Markov", "PEM") myPath <- "Output/02-Results-PEM/" startValues <- initiate.startValues_HReg(form, scrData, model=myModel, nChain=2) ## fit_PEM <- BayesID_HReg(form, scrData, id=NULL, model=myModel, hyperParams, startValues, mcmc.PEM, path=myPath) fit_PEM summ.fit_PEM <- summary(fit_PEM); names(summ.fit_PEM) summ.fit_PEM pred_PEM <- predict(fit_PEM) plot(pred_PEM, plot.est="Haz") plot(pred_PEM, plot.est="Surv") ################################################################# ## Analysis of Correlated Semi-Competing Risks Data ############# ################################################################# ################# ## WEIBULL-MVN ## ################# ## myModel <- c("semi-Markov", "Weibull", "MVN") myPath <- "Output/03-Results-WB_MVN/" startValues <- initiate.startValues_HReg(form, scrData, model=myModel, id, nChain=2) ## fit_WB_MVN <- BayesID_HReg(form, scrData, id, 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 pred_WB_MVN <- predict(fit_WB_MVN, tseq=seq(from=0, to=30, by=5)) plot(pred_WB_MVN, plot.est="Haz") plot(pred_WB_MVN, plot.est="Surv") ################# ## WEIBULL-DPM ## ################# ## myModel <- c("semi-Markov", "Weibull", "DPM") myPath <- "Output/04-Results-WB_DPM/" startValues <- initiate.startValues_HReg(form, scrData, model=myModel, id, nChain=2) ## fit_WB_DPM <- BayesID_HReg(form, scrData, id, 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 pred_WB_DPM <- predict(fit_WB_MVN, tseq=seq(from=0, to=30, by=5)) plot(pred_WB_DPM, plot.est="Haz") plot(pred_WB_DPM, plot.est="Surv") ############# ## PEM-MVN ## ############# ## myModel <- c("semi-Markov", "PEM", "MVN") myPath <- "Output/05-Results-PEM_MVN/" startValues <- initiate.startValues_HReg(form, scrData, model=myModel, id, nChain=2) ## fit_PEM_MVN <- BayesID_HReg(form, scrData, id, 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 pred_PEM_MVN <- predict(fit_PEM_MVN) plot(pred_PEM_MVN, plot.est="Haz") plot(pred_PEM_MVN, plot.est="Surv") ############# ## PEM-DPM ## ############# ## myModel <- c("semi-Markov", "PEM", "DPM") myPath <- "Output/06-Results-PEM_DPM/" startValues <- initiate.startValues_HReg(form, scrData, model=myModel, id, nChain=2) ## fit_PEM_DPM <- BayesID_HReg(form, scrData, id, 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 pred_PEM_DPM <- predict(fit_PEM_DPM) plot(pred_PEM_DPM, plot.est="Haz") plot(pred_PEM_DPM, plot.est="Surv") ## End(Not run) ```

BayesID_HReg function