# The function to implement Bayesian parametric and semi-parametric regression analyses for univariate time-to-event data in the context of hazard regression (HReg) models. Independent/cluster-correlated univariate right-censored survival data can be analyzed using hierarchical models. The prior for the baseline hazard function can be specified by either parametric (Weibull) model or non-parametric mixture of piecewise exponential models (PEM). ```r BayesSurv_HReg(Formula, data, id=NULL, model="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 event + censoring indicator $\sim$ covariates: i.e., $y$+$\delta$ ~ $x$. - `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 specification of baseline hazard functions: "Weibull" or "PEM". The second element needs to be set only for clustered data and is for the specification of cluster-specific random effects distribution: "Normal" or "DPM". - `hyperParams`: a list containing lists or vectors for hyperparameter values in hierarchical models. Components include, `WB` (a list containing a numeric vector for Weibull hyperparameters: `WB.ab`, `WB.cd`), `PEM` (a list containing numeric vectors for PEM hyperparameters: `PEM.ab`, `PEM.alpha`). Models for clustered data require additional components, `Normal` (a list containing a numeric vector for hyperparameters in a Normal prior: `Normal.ab`), `DPM` (a list containing numeric vectors for DPM hyperparameters: `DPM.ab`, `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 cluster-specific random effects: `storeV`, a logical value to determine whether all the posterior samples of $V$ are to be stored.) `tuning` (a list containing numeric values relevant to tuning parameters for specific updates in Metropolis-Hastings-Green (MHG) algorithm: `mhProp_V_var`, the variance of proposal density for $V$ in DPM models; `mhProp_alpha_var`, the variance of proposal density for $\alpha$ in Weibull models; `C`, a numeric value for the proportion that determines the sum of probabilities of choosing the birth and the death moves in PEM models. The value should not exceed 0.8; `delPert`, the perturbation parameter 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; `K_max`, the maximum number of splits allowed at each iteration in MHG algorithm for PEM models; `time_lambda` - time points at which the log-hazard function is 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 The function `BayesSurv_HReg` implements Bayesian semi-parametric (piecewise exponential mixture) and parametric (Weibull) models to univariate time-to-event data. Let $t_{ji}$ denote time to event of interest from subject $i=1,...,n_j$ in cluster $j=1,...,J$. The covariates $x_{ji}$ are incorporated via Cox proportional hazards model: $$ h(t_{ji} | x_{ji}) = h_{0}(t_{ji})\exp(x_{ji}^{\top}\beta + V_{j}), t_{ji}>0, $$ where $h_0$ is an unspecified baseline hazard function and $\beta$ is a vector of $p$ log-hazard ratio regression parameters. $V_j$'s are cluster-specific random effects. For parametric Normal prior specification for a vector of cluster-specific random effects, we assume $V$ arise as i.i.d. draws from a mean 0 Normal distribution with variance $\sigma^2$. Specifically, the priors can be written as follows: $$ V_j \sim Normal(0, \sigma^2), $$ $$ \zeta=1/\sigma^2 \sim Gamma(a_{N}, b_{N}). $$ For DPM prior specification for $V_j$, we consider non-parametric Dirichlet process mixture of Normal distributions: the $V_j$'s' are draws from a finite mixture of M Normal distributions, each with their own mean and variance, ($\mu_m$, $\sigma_m^2$) for $m=1,...,M$. Let $m_j\in\{1,...,M\}$ denote the specific component to which the $j$th cluster belongs. Since the class-specific ($\mu_m$, $\sigma_m^2$) are not known they are taken to be draws from some distribution, $G_0$, often referred to as the centering distribution. Furthermore, since the true class memberships are unknown, we denote the probability that the $j$th cluster belongs to any given class by the vector $p=(p_1,..., p_M)$ whose components add up to 1.0. In the absence of prior knowledge regarding the distribution of class memberships for the $J$ clusters across the $M$ classes, a natural prior for $p$ is the conjugate symmetric $Dirichlet(\tau/M,...,\tau/M)$ distribution; the hyperparameter, $\tau$, is often referred to as a the precision parameter. The prior can be represented as follows ($M$ goes to infinity): $$ V_j | m_j \sim Normal(\mu_{m_j}, \sigma_{m_j}^2), $$ $$ (\mu_m, \sigma_m^2) \sim G_{0},~~ for ~m=1,...,M, $$ $$ m_j | p \sim Discrete(m_j| p_1,...,p_M), $$ $$ p \sim Dirichlet(\tau/M,...,\tau/M), $$ where $G_0$ is taken to be a multivariate Normal/inverse-Gamma (NIG) distribution for which the probability density function is the following product: $$ f_{NIG}(\mu, \sigma^2 | \mu_0, \zeta_0, a_0, b_0) = f_{Normal}(\mu | 0, 1/\zeta_0^2) \times f_{Gamma}(\zeta=1/\sigma^2 | a_0, b_0). $$ In addition, we use $Gamma(a_{\tau}, b_{\tau})$ as the hyperprior for $\tau$. For non-parametric prior specification (PEM) for baseline hazard function, let $s_{\max}$ denote the largest observed event time. Then, consider the finite partition of the relevant time axis into $K + 1$ disjoint intervals: $0 ## Note The posterior samples of $V_g$ are saved separately in `working directory/path`. ## See Also `initiate.startValues_HReg`, `print.Bayes_HReg`, `summary.Bayes_HReg`, `predict.Bayes_HReg` ## Examples ```r ## Not run: # loading a data set data(survData) id=survData$cluster form <- Formula(time + event ~ cov1 + cov2) ##################### ## Hyperparameters ## ##################### ## Weibull baseline hazard function: alpha1, kappa1 ## WB.ab <- c(0.5, 0.01) # prior parameters for alpha ## WB.cd <- c(0.5, 0.05) # prior parameters for kappa ## PEM baseline hazard function: ## PEM.ab <- c(0.7, 0.7) # prior parameters for 1/sigma^2 ## PEM.alpha <- 10 # prior parameters for K ## Normal cluster-specific random effects ## Normal.ab <- c(0.5, 0.01) # prior for zeta ## DPM cluster-specific random effects ## DPM.ab <- c(0.5, 0.01) aTau <- 1.5 bTau <- 0.0125 ## hyperParams <- list(WB=list(WB.ab=WB.ab, WB.cd=WB.cd), PEM=list(PEM.ab=PEM.ab, PEM.alpha=PEM.alpha), Normal=list(Normal.ab=Normal.ab), DPM=list(DPM.ab=DPM.ab, aTau=aTau, bTau=bTau)) ################### ## MCMC SETTINGS ## ################### ## Setting for the overall run ## numReps <- 2000 thin <- 10 burninPerc <- 0.5 ## Settings for storage ## storeV <- TRUE ## Tuning parameters for specific updates ## ## - those common to all models mhProp_V_var <- 0.05 ## ## - those specific to the Weibull specification of the baseline hazard functions mhProp_alpha_var <- 0.01 ## ## - those specific to the PEM specification of the baseline hazard functions C <- 0.2 delPert <- 0.5 rj.scheme <- 1 K_max <- 50 s_max <- max(survData$time[survData$event == 1]) time_lambda <- seq(1, s_max, 1) ## mcmc.WB <- list(run=list(numReps=numReps, thin=thin, burninPerc=burninPerc), storage=list(storeV=storeV), tuning=list(mhProp_alpha_var=mhProp_alpha_var, mhProp_V_var=mhProp_V_var)) ## mcmc.PEM <- list(run=list(numReps=numReps, thin=thin, burninPerc=burninPerc), storage=list(storeV=storeV), tuning=list(mhProp_V_var=mhProp_V_var, C=C, delPert=delPert, rj.scheme=rj.scheme, K_max=K_max, time_lambda=time_lambda)) ################################################################ ## Analysis of Independent Univariate Survival Data ############ ################################################################ ############# ## WEIBULL ## ############# ## myModel <- "Weibull" myPath <- "Output/01-Results-WB/" startValues <- initiate.startValues_HReg(form, survData, model=myModel, nChain=2) ## fit_WB <- BayesSurv_HReg(form, survData, 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 <- "PEM" myPath <- "Output/02-Results-PEM/" startValues <- initiate.startValues_HReg(form, survData, model=myModel, nChain=2) ## fit_PEM <- BayesSurv_HReg(form, survData, 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 Univariate Survival Data ############ ############################################################### #################### ## WEIBULL-NORMAL ## #################### ## myModel <- c("Weibull", "Normal") myPath <- "Output/03-Results-WB_Normal/" startValues <- initiate.startValues_HReg(form, survData, model=myModel, id, nChain=2) ## fit_WB_N <- BayesSurv_HReg(form, survData, id, model=myModel, hyperParams, startValues, mcmc.WB, path=myPath) fit_WB_N summ.fit_WB_N <- summary(fit_WB_N); names(summ.fit_WB_N) summ.fit_WB_N pred_WB_N <- predict(fit_WB_N, tseq=seq(from=0, to=30, by=5)) plot(pred_WB_N, plot.est="Haz") plot(pred_WB_N, plot.est="Surv") ################# ## WEIBULL-DPM ## ################# ## myModel <- c("Weibull", "DPM") myPath <- "Output/04-Results-WB_DPM/" startValues <- initiate.startValues_HReg(form, survData, model=myModel, id, nChain=2) ## fit_WB_DPM <- BayesSurv_HReg(form, survData, 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_DPM, tseq=seq(from=0, to=30, by=5)) plot(pred_WB_DPM, plot.est="Haz") plot(pred_WB_DPM, plot.est="Surv") ################ ## PEM-NORMAL ## ################ ## myModel <- c("PEM", "Normal") myPath <- "Output/05-Results-PEM_Normal/" startValues <- initiate.startValues_HReg(form, survData, model=myModel, id, nChain=2) ## fit_PEM_N <- BayesSurv_HReg(form, survData, id, model=myModel, hyperParams, startValues, mcmc.PEM, path=myPath) fit_PEM_N summ.fit_PEM_N <- summary(fit_PEM_N); names(summ.fit_PEM_N) summ.fit_PEM_N pred_PEM_N <- predict(fit_PEM_N) plot(pred_PEM_N, plot.est="Haz") plot(pred_PEM_N, plot.est="Surv") ############# ## PEM-DPM ## ############# ## myModel <- c("PEM", "DPM") myPath <- "Output/06-Results-PEM_DPM/" startValues <- initiate.startValues_HReg(form, survData, model=myModel, id, nChain=2) ## fit_PEM_DPM <- BayesSurv_HReg(form, survData, 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) ```

BayesSurv_HReg function