Leuk: Cox regression

Several authors have discussed Bayesian inference for censored survival data where the integrated baseline hazard function is to be estimated non-parametrically Kalbfleisch (1978) ,Kalbfleisch and Prentice (1980), Clayton (1991), Clayton (1994).Clayton (1994) formulates the Cox model using the counting process notation introduced by Andersen and Gill (1982) and discusses estimation of the baseline hazard and regression parameters using MCMC methods. Although his approach may appear somewhat contrived, it forms the basis for extensions to random effect (frailty) models, time-dependent covariates, smoothed hazards, multiple events and so on. We show below how to implement this formulation of the Cox model in BUGS .

For subjects i = 1,...,n, we observe processes N _i (t) which count the number of failures which have occurred up to time t. The corresponding intensity process I _i (t) is given by

   I _i (t)dt = E(dN _i (t) | F _t- )
   
where dN _i (t) is the increment of N _i over the small time interval [t, t+dt), and F _t- represents the available data just before time t. If subject i is observed to fail during this time interval, dN _i (t) will take the value 1; otherwise dN _i (t) = 0. Hence E(dN _i (t) | F _t- ) corresponds to the probability of subject i failing in the interval [t, t+dt). As dt -> 0 (assuming time to be continuous) then this probability becomes the instantaneous hazard at time t for subject i . This is assumed to have the proportional hazards form   

   I _i (t) = Y _i (t) l ₀ (t)exp( b z _i )
   
where Y _i (t) is an observed process taking the value 1 or 0 according to whether or not subject i is observed at time t and l ₀ (t)exp( b z _i ) is the familiar Cox regression model. Thus we have observed data D = N _i (t), Y _i (t), z _i ; i = 1,..n and unknown parameters b and L ₀ (t) = Integral( l ₀ (u), u, t, 0), the latter to be estimated non-parametrically.    

The joint posterior distribution for the above model is defined by

   P( b , L ₀ () | D) ~ P(D | b , L ₀ ()) P( b ) P( L ₀ ())
   
For BUGS , we need to specify the form of the likelihood P(D | b , L ₀ ()) and prior distributions for b and L ₀ (). Under non-informative censoring, the likelihood of the data is proportional to   

    n
    P [ P I _i (t) ^{dN i (t)} ] exp(- I _i (t)dt)
    i = 1 t >= 0
   
   
This is essentially as if the counting process increments dN _i (t) in the time interval [t, t+dt) are independent Poisson random variables with means I _i (t)dt:   

   dN _i (t) ~ Poisson(I _i (t)dt)
   
We may write   

   I _i (t)dt = Y _i (t)exp( b z _i )d L ₀ (t)
   
where d L ₀ (t) = L ₀ (t)dt is the increment or jump in the integrated baseline hazard function occurring during the time interval [t, t+dt). Since the conjugate prior for the Poisson mean is the gamma distribution, it would be convenient if L ₀ () were a process in which the increments d L ₀ (t) are distributed according to gamma distributions. We assume the conjugate independent increments prior suggested by Kalbfleisch (1978), namely   

   d L ₀ (t) ~ Gamma(cd L ^* ₀ (t), c)
   
Here, d L ^* ₀ (t) can be thought of as a prior guess at the unknown hazard function, with c representing the degree of confidence in this guess. Small values of c correspond to weak prior beliefs. In the example below, we set d L ^* ₀ (t) = r dt where r is a guess at the failure rate per unit time, and dt is the size of the time interval.    
   
The above formulation is appropriate when genuine prior information exists concerning the underlying hazard function. Alternatively, if we wish to reproduce a Cox analysis but with, say, additional hierarchical structure, we may use the multinomial-Poisson trick described in the BUGS manual. This is equivalent to assuming independent increments in the cumulative `non-informative' priors. This formulation is also shown below.

The fixed effect regression coefficients b are assigned a vague prior

    b ~ Normal(0.0, 0.000001)