Imputing parameter from survival analysis (gamma distribution)

Modeling
, ,
1

Hi! I wanted to estimate time-till-event from survival analysis. My input data is from experiment where subject is housed in group. in each group, we have total hours that individual is positive. The Event of interest is the recovery of individual (positive->negative). All groups have a fixed end time (experiment is over) and that produced many right censored, that is individual did not recovered before the end of experiment. the majority of the data point is censored. I used gamma distribution where observed duration and censored duration is incorporated to find gamma’s parameter, the code is below.

model {
    vector[M] alphaInf_observed;
    vector[N] alphaInf_censored;
    alpha ~ normal(2,0.5);
    beta ~ normal(0,0.5);
    zalpha ~normal(0, 1);
    sigma ~ exponential(1);
  
    for ( i in 1:M ) {  //M is observed
      alphaInf_obs[i]  = alpha + zalpha[isol_obs[i]]*sigma;
      alphaInf_obs[i] = exp(alphaInf_obs[i] );
      target +=  gamma_lpdf(inf_obs[i]| alphaInf_obs[i],exp(beta));
    }
     for ( i in 1:N) {  //N is not censored
      alphaInf_cens[i]  = alpha + zalpha[isol_cens[i]]*sigmaA ;
      alphaInf_cens[i] = exp(alphaInf_cens[i] );
     target+= gamma_lccdf(inf_cens[i]| alphaInf_cens[i],exp(beta));
    }
}
}

It was very difficult to initialize and the calculation is also very off from raw data. That is the mean duration (alpha/beta) would be less than 5 or negative whereas the mean of raw data would be 200hours.
I tried many priors and also adjusted the structure

data {
  int<lower=0> N;
  int<lower=0> N_iso;
  vector[N] infp; // duration of positive
  vector[N] event; // censored is 0 and observed is 1
  int<lower=0> isol[N];  //index for grouping
}
parameters{
  real <lower=0> alpha;
  real beta;
  real<lower=0> sigma_alpha;
  real<lower=0> zalpha[N_iso]; 
}
model {
  vector[N] alpha_obs;
  alpha~normal(4,2);
  beta~normal(2,2);
  zalpha~normal(3,2);
  sigma_alpha~exponential(0.5);
  for (i in 1:N ){
    alpha_obs[i]  = alpha + zalpha[isol[i]]*sigma_alpha;
    alpha_obs[i] = exp(alpha_obs[i] );
    target += ( event[i] == 0)*gamma_lccdf(infp[i] | alpha_obs[i],exp(beta))+ (event[i]==1)*gamma_lpdf(infp[i]| alpha_obs[i],exp(beta) );
  }
}

I believed the problem is the censored data. I checked imputation in statistic rethinking book but I have difficulty adapt it to this Stan code. May I please have some suggestion? I also attached the simulated data for gamma distribution that will work the second code. With that simulation, I adjusted censored event to be lower but the model still produced a very low value compared to raw data.

inf_sim <- rgamma(n=100, shape=10, rate=0.05)
inf_sim<-sort(inf_sim, decreasing = F)
status_sim<-rep(c(1,0),times=c(90,10))
isolate<- rep(1:20, each=5)

sim_lst<- list(
  N=length(inf_sim),
  infp= inf_sim,
  event=status_sim,
  N_iso= 20,
  isol=isolate
)

I have been trying with these codes for a while. Would be really appreciated for your help.
Thank you in advance.