Suppose we have N individuals in the dataset where the time-to-event data for individual i is given by T_{i}. Of those N individuals, N_{obs} experience the event and N_{cen} are censored at \{C_{i}\}_{i=1}^{N_{cen}}. The vector representing the true event times, say \textbf{T}_{true} , can be split into \textbf{T} _{true}= [\textbf{T}_{obs,true}, \textbf{T}_{cen,true} ] where \textbf{T}_{obs,true} denotes the vector of (true ) death times for observed individuals and \textbf{T}_{cen,true} denotes the vectorof (true ) death times for censored individuals. This implies T_{i,obs,true} = T_{i} and T_{i,cen,true} \geq C_{i} .
In this simple example, suppose that the data comes from Weibull distribution (with shape \eta and scale \lambda). Our goal is augment the missing true death times for censored individuals. In other words, at each iteration q, q \in \{ 1, ..., Q\}, we wish to draw:
\textbf{T}_{i,cen,true}^{(q+1)} \sim p(\textbf{T}_{i,cen,true}| \textbf{T}_{i,cen,true} \geq C_{i}, \textbf{T}_{i,obs,true}, \eta^{(q)}, \lambda^{(q)})
\eta^{(q+1)} , \lambda^{(q+1)}\sim p(\eta , \lambda | \textbf{T}_{i,cen,true}^{(q+1)}, \textbf{T}_{i,obs,true})
Does the following code give us what we want? Please provide your thoughts and comments. Thanks!
data {
int<lower=0> Nobs; // Number of individuals who experienced death before censoring
int<lower=0> Ncen; // Number of censored individuals
real<lower=0> ObsDeathTime[Nobs]; // Observed death times
real<lower=0> ObsCenTime[Ncen]; // Observed Censoring times
}
parameters {
real<lower = 0> eta; // shape parameter for ObsTime
real<lower = 0> lambda; // scale parameter for ObsTime
vector<lower=0>[Ncen] T_miss_raw; // vector for missing death times
}
transformed parameters{
vector<lower=0>[Ncen] T_miss;
for(j in 1:Ncen){
T_miss[j] = T_miss_raw[j] + ObsCenTime[j]; // make sure that T_miss[j] is greater than or equal to C_i
}
}
model {
// Priors
eta ~ gamma(1,1);
lambda ~ gamma(1,1);
//Likelihood
for(i in 1:Nobs){
target+= weibull_lpdf(ObsDeathTime[i]|eta,lambda); // observed death time
}
for(j in 1:Ncen){
target+= weibull_lpdf(T_miss[j]|eta,lambda); // contribution of censored death time to the likelihood once we know what their values are
}
}