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);
}
```