Getting to see Stan's initial values after it rejects my own

Modeling
1

I’m trying to develop a mixed model for longitudinal disease data using a generalized logistic function to describe disease trajectories. The model is beta distributed. I’m looking to get some help identifying an issue regarding initial values. In the stan code, I specify the initial values as below:

init = function() list(theta_TMS0 = rnorm(1,.1,.01), theta_r = rnorm(1,.7,0.01), eta = rmvnorm(P,rep(0,2),sigma = .001*diag(2)), sigma_diag = rlnorm(2,log(.2)-.5), Omega = Q, tau = 1, beta = 1)

In the code above, P is the number of patients for which there is data (4754 patients). Unfortunately, I cannot post the data.

I have checked numerous times to make sure that these initial values are reasonable and do not cause any numeric issues. However, when I run the code, stan immediately says it’s rejecting the initial value, BUT, it eventually finds a set of initial values that work, and the model runs. (Before I continue with my questions, is the following a true statement? When Stan says it’s rejecting the initial value, that means it already rejected my initial values, and now it’s giving an error message for the one of the 100 attempts that stan is attempting for it’s initial value?)

I would like to know how to see what the initial values that stan found on it’s own that worked, so I can compare them against my own. Below is a typical error message that is produced:

SAMPLING FOR MODEL ‘HD CTS’ NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: model3e9459053fc2_HD_CTS_namespace::write_array: muS[7] is 4.71149, but must be less than or equal to 1 (in ‘model3e9459053fc2_HD_CTS’ at line 46)

I need to understand this issue, so I can extend the model to include covariates. If I attempt to fit the model with covariates, it never finds a set of workable initial values. If there is something blatantly obvious please mention it as I’m drawing a blank. Feel free to mention any other non-related tips as I’m still a novice to stats modeling. Thanks for the responsiveness ahead of time.


data{
  int N;
  int P;
  int IDp[N]; 
  real<lower=0> TMS[N];
  real<lower=0> TIME[N];
}

parameters{

  //TMS model
  real<lower=0,upper=1> theta_TMS0;
  real theta_r;
  
  //Shape parameters in beta distribution as fixed effects
  real<lower=0> tau;
  real<lower=0> beta;
  
  //Random effects
  vector[2] eta[P];
  
  //Standard deviations of random effects
  vector<lower=0>[2] sigma_diag;
  
  //Correlation matrix
  corr_matrix[2] Omega;
 
}

transformed parameters{
  
  vector[N] baseline_cov;
  vector[N] rate_cov;
  vector[N] TMS0;
  vector[N] r;
  vector<lower=0,upper=1>[N] muS;
  real tau_trans;
  real beta_trans;
  
  //Covariance matrix 
  cov_matrix[2] Sigma;
  Sigma = quad_form_diag(Omega,sigma_diag);
  
  tau_trans = (tau+1)*50;
  beta_trans = beta;
  
   for(i in 1:N){
    
    baseline_cov[i] = theta_TMS0;
    rate_cov[i] = theta_r;

    TMS0[i] = baseline_cov[i]*exp(eta[IDp[i]][1]);
    r[i] = rate_cov[i] + eta[IDp[i]][2];

    muS[i] = TMS0[i]/(TMS0[i]^beta_trans + (1-TMS0[i]^beta_trans)*exp(-beta_trans*r[i]*TIME[i]))^(1/beta_trans);
    
   }
  
}

model{
  
  //Priors
  theta_TMS0~cauchy(0,10);
  theta_r~cauchy(0,10);
  tau~cauchy(0,10); 
  beta~cauchy(0,10);

  //correlation priors
  sigma_diag~cauchy(0,10);
  Omega~lkj_corr(1);
  
  for(i in 1:P){
    eta[i] ~ multi_normal(rep_vector(0,2),Sigma);
  }
   
    TMS ~beta(muS*tau_trans, (1-muS)*tau_trans);
}
2

your parameters are constrained - the sampler works on the unconstrained scale.
this is described in the Stan jstatsoft paper:

see sections 4.3, 4.7

cheers,
Mitzi

3

Hi Mitzi,

That’s informative to know that stan performs transformations to unconstrained scales. It appears, although I can pin point it, that my at least one of my initial values is being transformed to + or - infinity, i.e. I must be specifying something right on the constrained boundary.

In any case, I decided to use stan’s random function for init, and just tightened up the init_r interval. Thankfully, that seems to be working which is a big relief. Less is more.

Thanks,
Jackson

1 Like