BayesSurv_HReg: The function to implement Bayesian parametric and...

Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/BayesSurv_HReg.R

Description

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).

Usage

1
2
BayesSurv_HReg(Formula, data, id=NULL, model="Weibull", hyperParams,
        startValues, mcmcParams, na.action = "na.fail", subset=NULL, path=NULL)

Arguments

Formula

a Formula object of the form y+δ ~ x. See Details and Examples below.

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 α 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}β + V_{j}), t_{ji}>0,

where h_0 is an unspecified baseline hazard function and β 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 σ^2. Specifically, the priors can be written as follows:

V_j \sim Normal(0, σ^2),

ζ=1/σ^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, (μ_m, σ_m^2) for m=1,...,M. Let m_j\in\{1,...,M\} denote the specific component to which the jth cluster belongs. Since the class-specific (μ_m, σ_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 jth 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(τ/M,...,τ/M) distribution; the hyperparameter, τ, 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(μ_{m_j}, σ_{m_j}^2),

(μ_m, σ_m^2) \sim G_{0},~~ for ~m=1,...,M,

m_j | p \sim Discrete(m_j| p_1,...,p_M),

p \sim Dirichlet(τ/M,...,τ/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}(μ, σ^2 | μ_0, ζ_0, a_0, b_0) = f_{Normal}(μ | 0, 1/ζ_0^2) \times f_{Gamma}(ζ=1/σ^2 | a_0, b_0).

In addition, we use Gamma(a_{τ}, b_{τ}) as the hyperprior for τ.

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<s_1<s_2<...<s_{K+1} = s_{\max}. For notational convenience, let I_k=(s_{k-1}, s_k] denote the k^{th} partition. For given a partition, s = (s_1, …, s_{K + 1}), we assume the log-baseline hazard functions is piecewise constant:

λ_{0}(t)=\log h_{0}(t) = ∑_{k=1}^{K + 1} λ_{k} I(t\in I_{k}),

where I(\cdot) is the indicator function and s_0 \equiv 0. In our proposed Bayesian framework, our prior choices are:

π(β) \propto 1,

λ | K, μ_{λ}, σ_{λ}^2 \sim MVN_{K+1}(μ_{λ}1, σ_{λ}^2Σ_{λ}),

K \sim Poisson(α),

π(s | K) \propto \frac{(2K+1)! ∏_{k=1}^{K+1}(s_k-s_{k-1})}{(s_{K+1})^{(2K+1)}},

π(μ_{λ}) \propto 1,

σ_{λ}^{-2} \sim Gamma(a, b).

Note that K and s are treated as random and the priors for K and s jointly form a time-homogeneous Poisson process prior for the partition. The number of time splits and their positions are therefore updated within our computational scheme using reversible jump MCMC.

For parametric Weibull prior specification for baseline hazard function, h_{0}(t) = α κ t^{α-1}. In our Bayesian framework, our prior choices are:

π(β) \propto 1,

π(α) \sim Gamma(a, b),

π(κ) \sim Gamma(c, d).

We provide a detailed description of the hierarchical models for cluster-correlated univariate survival data. The models for independent data can be obtained by removing cluster-specific random effects, V_j, and its corresponding prior specification from the description given above.

Value

BayesSurv_HReg returns an object of class Bayes_HReg.

Note

The posterior samples of V_g are saved separately in working directory/path.

Author(s)

Kyu Ha Lee and Sebastien Haneuse
Maintainer: Kyu Ha Lee <[email protected]>

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, 111, 515, 1075-1095.

Alvares, D., Haneuse, S., Lee, C., Lee, K. H. (2018+), SemiCompRisks: an R package for independent and cluster-correlated analyses of semi-competing risks data, submitted, arXiv:1801.03567.

See Also

initiate.startValues_HReg, print.Bayes_HReg, summary.Bayes_HReg, predict.Bayes_HReg

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
	
## 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)

SemiCompRisks documentation built on May 7, 2018, 9:04 a.m.