Gamma Process

Walker & Mallick (1997) consider {λ_k} as independent gamma variables, i.e., λ_k ∼ Ga(α_k, β_k) independent for k = 1, 2, ... Nieto-Barajas & Walker (2002) consider a dependent process for {λ_k}. They start with λ₁ ∼ Ga(α₁, β₁) and take u_k|λ_k ∼ Po(c_kλ_k) and λ_k + 1|u_k ∼ Ga(α_k + 1 + u_k, β_k + 1 + c_k) and so on. These updates arise from the joint density

f(u, λ) = Ga(λ|α, β)Po(u|cλ) and so constitute a Gibbs type update. The difference is that they are changing the parameters (α, β, c) at each update so the chain is not stationary and marginally the {λ_k} are not gamma.

However, $$ \mathrm{E}[\lambda_{k+1}|\lambda_k] = \frac{\alpha_{k+1}+c_k\lambda_k}{\beta_{k+1}+c_k} $$

If c_k = 0, then P(u_k = 0) = 1 and hence the {λ_k} are independent gamma and we have the prior process of Walker and Mallick (1997).

An important result is that if we take α_k = α₁ and β_k = β₁ to be constant for all k, then the process {u_k} is a Poisson-gamma process with implied marginals λ_k ∼ Ga(α₁, β₁). One can note that if u₁ is Poisson distributed and λ₂|u₁ is conditionally gamma, then λ₂ is never gamma.

Beta Process

Nieto-Barajas and Walker (2002) start with π₁ ∼ Be(α₁, β₁) and take u_k|π_k ∼ Bi(c_k, π_k), π_k + 1|u_k ∼ Be(α_k + 1 + u_k, β_k + 1 + c_k − u_k) and so on. These arise from a binomial-beta conjugate set-up, from the joint density

f(u, π) = Be(π|α, β)Bi(u|c, π)

Clearly

$$ \mathrm{E}[\pi_{k+1}|\pi_k] = \frac{\alpha_{k+1}+c_k\pi_k}{\alpha_{k+1}+\beta_{k+1}+c_k} $$

As with the gamma process, if we choose c_k = 0, then P(u_k = 0) = 1 and so the {u_k} become independent beta and we obtain the prior of Hjort (1990).

A significant result is that if we take α_k = α₁ and β_k = β₁ to be constant for all k, then the process {u_k} is a binomial-beta process with marginals u_k ∼ BiBe(α₁, β₁, c_k). Moreover, the process {π_k} becomes strictly stationary and marginally π_k ∼ Be(α₁, β₁).

Gamma Process

Let T be a continuous random variable with cumulative distribution function F(t) = P(T ≤ t) on [0, ∞). Consider the time axis partition 0 = τ₀ < τ₁ < τ₂ < ⋯, and let λ_k be the hazard rate in the interval (τ_k − 1, τ_k], then the hazard function is given by $$ h (t) = \sum_{k=1}^{\infty} \lambda_kI_{(\tau_{k-1},\tau_k]}(t) $$

So, the cumulative distribution and density functions, given {λ_k}, are F(t|{λ_k}) = 1 − e^−H(t), f(t|{λ_k}) = h(t)e^−H(t), where H(t) = ∫₀^th(s)ds. We also have that f(λ_k|u_k − 1, u_k) = Ga(α_k + u_k − 1 + u_k, β_k + c_k − 1 + c_k)

Therefore, given a sample T₁, T₂, ..., T_n from f(t{λ_k}), it is straightforward to derive f(λ_k|u_k − 1, u_k, data) = Ga(α_k + u_k − 1 + u_k + n_k, β_k + c_k − 1 + c_k + m_k), where n_k = number of uncensored observations in (τ_k − 1, τ_k], m_k = ∑_ir_ki, and

$$ r_{ki} = \left \{ \begin{array}{l l} \tau_k - \tau_{k-1} & t_i > \tau_k \\ t_i - \tau_{k-1} & $t$ \in (\tau_{k-1},\tau_k] \\ 0 & $otherwise$ \end{array} \right. $$

Additionally,

$$ P(u_k=u | \lambda_k,\lambda_{k+1},data) \propto f(u|\lambda) \propto \frac{[c_k(c_k+\beta_{k+1})\lambda_k\lambda_{k+1}]^u}{\Gamma(u+1)\Gamma(\alpha_{k+1}+u)} $$

with u = 0, 1, 2, .... Hence, with these full conditional distributions, a Gibbs sampler is straightforward to implement in order to obtain posterior summaries.

We can learn about the {c_k} by assigning an independent exponential distribution with mean ϵ for each c_k, k = 1, .., K − 1. The Gibbs sampler can be extended to include the full conditional densities for each c_k. It is not difficult to derive that a c_k from f(c_k|u, λ, data) can be taken from the density

$$ f(c_k|u_k,\lambda_k,\lambda_{k+1}) \propto (\beta_{k+1}+c_k)^{\alpha_{k+1}+u_k}c_k^{u_k}\exp\left\{ -c_k \left (\lambda_{k+1}+\lambda_k+\frac{1}{\epsilon} \right) \right\} $$ for c_k > 0.

Dependence between c_k′s can be introduced through a hierarchical model via assigning a distribution to ϵ ∼ Ga(a₀, b₀). The update would be given by: $$ f(\epsilon|\{c_k\}) = Ga\left(\epsilon |a_0 + K, b_0 + \sum_{k=1}^K c_k \right) $$ where K is the number of intervals generated by the time axis partition. This hierarchical specification of the initial distribution for c_k let us assign a better value for c_k.

Simulating from this distribution is not so straightforward. However, we construct a hybrid algorithm using a Metropolis-Hastings scheme taking advantage of the Markov chain generated by the Gibbs sampling.

Beta process

Let T be a discrete random variable taking values in the set {τ₁, τ₂, ...} with probability density function f(τ_k) = P(T = τ_k). Let π_k be the hazard rate at τ_k, then the cumulative distribution and the density functions, given {π_k}, are $F(\tau_j|\{\pi_k\}) = 1 - \prod_{k=1}^j(1-\pi_k)$ and $f(\tau_j|\{\pi_k\}) = \pi_j \prod_{k=1}^{j-1}(1-\pi_k)$

The conditional distribution of π_k is f(π_k|u_k − 1, u_k) = Be(π_k|α_k + u_k − 1 + u_k, β_k + c_k − 1 − u_k − 1 + c_k − u_k),

Thus, given a sample T₁, T₂, ..., T_n form f(⋅|{π_k}) it is straightforward to derive f(π_k|u_k − 1, u_k, data) = Be(π_k|α_k + u_k − 1 + u_k + n_k, β_k + c_k − 1 − u_k − 1 + c_k − u_k + m_k), where n_k = number of failures at τ_k, m_k = ∑_ir_ki and

$$ r_{ki} = \left \{ \begin{array}{l l} 1 & \quad t_i > \tau_k \\ 0 & \quad $otherwise$ \end{array} \right. $$

Additionally, $$ P(u_k=u | \pi_k,\pi_{k+1}, data) \propto \frac{\theta_k^u}{\Gamma(u+1)\Gamma(c_k-u+1)\Gamma(\alpha_{k+1}+u)\Gamma(\beta_{k+1}+c_k-u)} $$ with u = 0, 1, ..., c_k and $$ \theta_k = \frac{\pi_k\pi_{k+1}}{(1-\pi_k)(1-\pi_{k+1})} $$ As before, obtaining posterior summaries via Gibbs sampler is simple.

We can learn about the {c_k} via assigning each c_k an independent Poisson distribution with mean ϵ. The Gibbs sampler can be extended to include the full conditional densities for each c_k. A c_k from f(c_k|u, π, data) can be taken from the density

$$ f(c_k|u_k,\pi_k,\pi_{k+1}) \propto \frac{\Gamma(\alpha_{k+1}+\beta_{k+1}+c_k)}{\Gamma(\beta_{k+1}+c_k-u_k)\Gamma(c_k-u_k+1)} \left[\epsilon (1-\pi_{k+1})(1-\pi_k) \right]^{c_k} $$ with c_k ∈ {u_k, u_k + 1, u_k + 2, ...}.

Dependence between c_k’s can be introduced through a hierarchical model via assigning a distribution to ϵ ∼ Ga(a₀, b₀). So the update would be given by:

$$ f(\epsilon|\{c_k\}) = Ga\left(\epsilon |a_0 + K, b_0 + \sum_{k=1}^K c_k \right) $$ where K is the number of discrete values in random variable T.

Cox-gamma model

Differing from most of the previous Bayesian analysis of the proportional hazards model, Nieto-Barajas (2003) models the baseline hazard rate function with a stochastic process. Let T_i be a non negative random variable which represents the failure time of i and Z_i = (Z_i1, ..., Z_ip) the vector containing its p explanatory variables. Therefore, the hazard function for individual i is: λ_i(t) = λ₀(t)exp{Z_i′θ} where λ₀(t) is the baseline hazard rate and θ is regression’s coefficient vector. The cumulative hazard function for individual i becomes $$ H_i(t)=\sum_{k=1}^{\infty}\lambda_k W_{i,k}(t,\lambda), $$ where,

$$ W_{i,k}(t,\theta) = \left \{ \begin{array}{l l} (\tau_k - \tau_{k-1})\exp\{Z_i'\theta \} & t_i > \tau_k \\ (t_i - \tau_{k-1})\exp\{Z_i'\theta \} & t_i \in (\tau_{k-1},\tau_k] \\ 0 & otherwise \end{array} \right. $$

Given a sample of possible right-censored observations where T₁, ..., T_n are uncensored and T_nu + 1, ..., T_n are right-censored, the conditional posterior distributions for the parameters of the semi-parametric model are:

where $n_k=\sum_{i=1}^{n_u}I_{(\tau_{k-1},\tau_k]}(t_i)$ and $m_k(\theta)=\sum_{i=1}^n W_{i,k}(t_i,\theta)$.

Similarly as we pointed out in the previous two cases, we incorporate a hyper prior process for the {c_k} such that c_k ∼ Ga(1, ϵ_k). The set of full conditional posterior distributions can then be extended to include

$$ f(c_k|u_k,\lambda_k,\lambda_{k+1}) \propto (\beta_{k+1}+c_k)^{\alpha_{k+1}+u_k}c_k^{u_k}\exp\left\{ -c_k \left (\lambda_{k+1}+\lambda_k+\frac{1}{\epsilon} \right) \right\} $$

Dependence between c_k’s can be introduced through a hierarchical model via assigning a distribution to ϵ ∼ Ga(a₀, b₀). So the update would be given by: where K is the number of intervals generated by the time axis partition. This hierarchical specification of the initial distribution for c_k let us assign a better value for c_k.

Example 1. Independence case. Unitary length intervals}

We obtain with our model the Nelson-Aalen and Kaplan-Meier estimators by defining c_k as a null vector –through fixing type.c = 1–. A unitary partitioned axis is obtained by fixing type.t = 2. In Figure we show that our estimators –under independence– return the same results as the Nelson-Aalen and Kaplan-Meier estimators.

ExG1 <- GaMRes(times, delta, type.t = 2, K = 35, type.c = 1, 
               iterations = 3000)
GaPloth(ExG1, confint = FALSE)

Example 2. Introducing dependence through c. Unitary length intervals

The influence of c –or c.r– can be understood as a dependence parameter: the greater the value of each c_k, k = 0, 1, 2, ..., K − 1, the higher the dependence between intervals k and k + 1. For this example, we assign a fixed value to vector c, so we use type.c = 2. Comparing with the previous example, we see that the estimates that were zero, now have a positive value –see Figures and –. Note how this model compares to Example 5 –see Figure –. The difference between those examples is how the partition is defined.

ExG2 <- GaMRes(times, delta, type.t = 2, K = 35, type.c = 2, 
               c.r = rep(50, 34), iterations = 3000)
GaPloth(ExG2)

Additionally, we can get further detail on the Gibbs sampler with a diagnosis of the resulting Markov chain. We can run this diagnosis for each entry of λ, u, c or ϵ. In Figure we show the diagnosis for λ₆ which includes the trace, the ergodic mean, the ACF function and the histogram for the generated chain.

GaPlotDiag(ExG2, variable = "lambda", pos = 6)

Example 3. Varying c through a distribution. Unitary length intervals

As we reviewed, we can learn about the {c_k} via assigning an exponential distribution with mean ϵ. The estimates using ϵ = 1 –type.c = 3– and unitary length intervals –type.t = 2– are shown in Figure . We can compare the hazard function with the previous example, where c was fixed, and observe that because of the variability given to c, the change on the estimated values between two contiguous intervals are greater in this example. The survival function echoes the shape of the Kaplan-Meier with higher decrease rates.

Comparing this example with Example 6 –Figure –, as with the previous example, the main difference is given by the partition of the time axis. This affects the estimates as we will note later.

ExG3 <- GaMRes(times, delta, type.t = 2, K = 35, type.c = 3, epsilon = 1,
               iterations = 3000)
GaPloth(ExG3)

Example 4. Using a hierarchical model to estimate c. Unitary intervals

Previous example can be extended with a hierarchical model, assigning a distribution to ϵ  ∼ Ga(a₀, b₀) with a₀ = b₀ = 0.01. In order to set up the model we should set type.c = 4. The result displayed on Figure is a soft hazard function and a survival function that decreases faster than the Kaplan-Meier estimate.

ExG4 <- GaMRes(times, delta, type.t = 2, K = 35, type.c = 4, 
               iterations = 3000)
GaPloth(ExG4)

Example 5. Introducing dependence through c. Equally dense intervals

This example illustrates the same concept as Example 2 –how c introduces dependence–, but with a different partition of the time axis. To get this partition we set type.t = 1. Figure shows that the survival function is close to the Kaplan-Meier estimate. The fact that it gets closer to the K-M estimates does not makes it a better estimate, we only could say that this partition yields, in average, a smaller hazard rate than the Example 2 built with a unitary length partition.

ExG5 <- GaMRes(times, delta, type.t = 1, K = 8, type.c = 2, 
               c.r = rep(50, 7), iterations = 3000)
GaPloth(ExG5)

Example 6. Varying c through a distribution. Equally dense intervals

We can compare this example on Figure with Example 3 –Figure –. We use less intervals, so the survival function is smoother. Our estimate decreases faster than the Kaplan-Meier estimate.

ExG6 <- GaMRes(times, delta, type.t = 1, K = 8, type.c = 3,
               iterations=3000)
GaPloth(ExG6)

Example 7. Using a hierarchical model to estimate c. Equally dense intervals

The survival curve for this particular example results in a smoothed version of the Kaplan-Meier estimate (Figure ). As with previous examples, it can be compared with the unitary partition example (see Example 4, Figure ).

ExG7 <- GaMRes(times, delta, type.t = 1, K = 8, type.c = 4,
               iterations = 3000)
GaPloth(ExG7)

%

Example 1. Independence case

As with the Gamma Example, we obtain the Nelson-Aalen and Kaplan-Meier estimators by defining c_k as a null vector through fixing type.c = 1 –see Figure – The conclusion does not change: the independence case of our model results on the N-A and K-M estimators.

ExB1 <- BeMRes(times, delta, type.c = 1, iterations = 3000)
BePloth(ExB1, confint = FALSE)

Example 2. Introducing dependence through c

The influence of c –or c.r– can be also understood as a dependence parameter: the greater the value of each c_k, k = 0, 1, 2, ..., K − 1, the higher the dependence between intervals k and k + 1. In this example, we fix each c entry at 100. As we are defining vector c with fixed values, we should fix type.c = 2. We see on Figure that the hazard function estimate turns smoother than on the previous example. The steps on the survival function appear to be more uniform than on the Kaplan-Meier estimate.

ExB2 <- BeMRes(times, delta, type.c = 2, c.r = rep(100, 39), 
               iterations = 3000)
BePloth(ExB2)

Additionally, we can get further detail on the Gibbs sampler with a diagnosis of the resulting Markov chain. We can run this diagnosis for each entry from π, u, c or ϵ. In Figure we show the diagnosis for π₁₀ which includes plot of the trace, the ergodic mean, the ACF function and the histogram of the chain.

BePlotDiag(ExB2, variable = "Pi", pos = 6)

Example 3. Varying c through a distribution

As with the gamma model, we can learn about the {c_k} via assigning an exponential distribution with mean ϵ. The estimates using ϵ = 1 and unitary length intervals are shown on Figure . Note that the confidence intervals have widen: it is a signal that we have introduced variability to the model –because of the distribution assigned to c–.

ExB3 <- BeMRes(times, delta, type.c = 3, epsilon = 1, iterations = 3000)
BePloth(ExB3)

Example 4. Using a hierarchical model to estimate c

The previous example can be extended with a hierarchical model, assigning a distribution to ϵ  ∼ Ga(a₀, b₀), with a₀ = b₀ = 0.01. In order to set up the model we should set type.c = 4. The result displayed on Figure is a soft hazard function and a survival function that decreases faster than the Kaplan-Meier estimate.

ExB4 <- BeMRes(times, delta, type.c = 4, iterations = 3000)
BePloth(ExB4)