Generalized Gamma distribution with absolute value on a parameter

Hello. I am working with the Generalized Gamma Distribution with an absolute value on a parameter.

I don’t know how well Stan can handle this type of distributions, since the derivative is difficult and may be indifferentiable at that point. Does anyone have experience with distribution that has an absolute value of a particular parameter in its pdf? I am thinking Generalized Gamma, Generalized Beta’s, etc.

This is what I did.

N = 500;
mu = 3;
sigma = 1;
Q = 0.5;
#censoring point > 1 (uniform for all subjects)
C=10;

#Simulate from GenGamma distribution with uniform censoring
rgengamma_cen <- function(N, mu, sigma, Q, C) {
  
  X=flexsurv::rgengamma(N, mu, sigma, Q);
  
  X[which(X>=C)]=C;
  
  return ( X );
  
}

t = rgengamma_cen(N, mu, sigma, Q, C);

df = data.frame( time = t,
                      status = (t<C) )


#Check to ensure simulation is ok and it is.
fs1 = flexsurvreg(Surv(t, status) ~ 1, data = df, dist = "gengamma")

#data input for Stan
data_list <- list(
  Nuc = length(which(df$status==1)), 
  yuc = unlist(df[which(df$status==1),"t"]),
  Nc = length(which(df$status==0)),
  yc = unlist(df[which(df$status==0),"t"])
)

Then I follow the parameterization used by flexsurv, which is from Cox et al (2007):

Generalized Gamma Stan code:

// Generalized Gamma as implemented in flexsurv R package
// https://search.r-project.org/CRAN/refmans/flexsurv/html/GenGamma.html
// only require 'sigma' >0
  
functions{

real generalized_gamma_lpdf_fs(real x, real mu, real sigma, real Q) {

  real Qw = Q*(log(x)-mu)/sigma;
  real Qis = inv(pow(Q,2));
  real ll;
  ll = log(abs(Q)) + Qis*(log(Qis)+Qw) - Qis*exp(Qw) - log(sigma*x) - lgamma(Qis);
  return ll;
  
}

real generalized_gamma_cdf_fs(real x, real mu, real sigma, real Q) {
  
  real w = (log(x)-mu)/sigma;
  real Qw = Q*(log(x)-mu)/sigma;
  real Qis = inv(pow(Q,2));
  real lik;
  if (Q > 0) lik=gamma_cdf(Qis*exp(Qw),Qis,1);
  else if (Q < 0) lik=1-gamma_cdf(Qis*exp(Qw),Qis,1);
  else lik=std_normal_cdf(w);  
  return lik;

 }



}

data{
  
  int Nuc;
  real yuc[Nuc];
  int Nc;
  real yc[Nc];
  
}

parameters{
  
  real mu;
  real<lower=0> sigma;
  real Q;
  
}
  
model{
  
  mu ~ normal(0,5);
  sigma ~ normal(0,1);
  Q ~ normal(0,5);
  
  for (i in 1:Nuc) {
    target += generalized_gamma_lpdf_fs(yuc[i],mu,sigma,Q);
  }
  
  for (i in 1:Nc){
    target += log1p(generalized_gamma_cdf_fs(yc[i],mu,sigma,Q));
  }
  
}

Surely, my Stan code could be vectorized and made faster. [The last run took me 50 mins.] Nevertheless, These are the warnings and issues I got:

Warning: 8 of 4000 (0.0%) transitions ended with a divergence.
This may indicate insufficient exploration of the posterior distribution.
Possible remedies include:

  • Increasing adapt_delta closer to 1 (default is 0.8)
  • Reparameterizing the model (e.g. using a non-centered parameterization)
  • Using informative or weakly informative prior distributions

1974 of 4000 (49.0%) transitions hit the maximum treedepth limit of 10 or 2^10-1 leapfrog steps.
Trajectories that are prematurely terminated due to this limit will result in slow exploration.
Increasing the max_treedepth limit can avoid this at the expense of more computation.
If increasing max_treedepth does not remove warnings, try to reparameterize the model.

The worst of all, the rhat’s are all bad, >2!

Could someone please take a look? Does my Stan code have an error? Or is there a better way to handle the absolute value? Thanks.