**Summary:** Below I include details about an example model and how to fit it, it is just an example to illustrate a couple of points. My main question is at the end. I would like to generate simulated datasets of the same function but with various levels of the parameters ( for example, generate dose response curves with various levels of variation in the slope parameter), but it seems like it is not feasible to have a single rstan model that could handle all of them. My question is what strategies do people use when trying to do simulation studies using stan? I put a couple of ideas at the end.

**Example Model**

I have a general question that came after having a lot of problems fitting one of the parameters in a non-linear model. I will simplify the model as much as possible and only consider the troublesome parameter. I am fitting a four parameter logistic function such that:

Where A is the range, B the slope (specific for each block b), C the logEC50, and D the no-dose asymptote.

For this example I will only try to fit the B parameter. The centered model would be:

and the corresponding non-centered model would be:

This is some example data extracted from a simulated dataset:

data â data.frame(log2Dose=c(-13.83,-13.83,-11.31,-11.31,-8.8,-8.8,-6.29,-6.29,-3.77,-3.77,-1.26,-1.26,1.26,1.26,3.77,3.77,6.29,6.29,8.8,8.8,11.31,11.31,13.83,13.83),

y=c(0.24,0.25,0.21,0.29,0.23,0.32,0.26,0.3,0.37,0.35,0.56,0.55,0.86,0.89,1.01,1.07,1.1,1.19,1.12,1.2,1.15,1.23,1.14,1.21),

block=c(2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1,2,1))

Which looks like this:

And this is the model:

```
model_code <- 'functions {
real fplBPosRange(real logX, real logBase, real A, real B, real logC, real D) {
return A/(1+logBase^(-B*(logX-logC)))+D;
}
}
data {
int<lower=1> N; // size of dataset
vector[N] logX; // dose value in logarithmic scale
real y[N]; // Response values
int<lower=1> nBlocks; // number of blocks
int<lower=1,upper=nBlocks> cultureBlock[N]; //block
}
parameters {
real<lower=-1> sigmaB_raw;
real<lower=-10> meanB_raw;
vector[nBlocks] blockNoiseB_raw;
real<lower=0> sigma;
}
transformed parameters {
vector[N] B;
vector[nBlocks] blockNoiseB;
real meanB;
real sigmaB;
sigmaB = 0.2 + 0.2*sigmaB_raw;
for (i in 1:nBlocks) {
blockNoiseB[i] = blockNoiseB_raw[i]*sigmaB;
}
meanB = 1 + meanB_raw*0.1;
for (i in 1:N) {
B[i] = meanB + blockNoiseB[cultureBlock[i]];
}
}
model {
vector[N] mu_y;
sigmaB_raw ~ normal(0,1);
blockNoiseB_raw ~ normal(0,1);
meanB_raw ~ normal(0,1);
sigma ~ normal(0,1);
for (i in 1:N) {
mu_y[i] = fplBPosRange(logX[i],
2,
1,
B[i],
0.01,
0.07);
};
y ~ normal(mu_y,sigma);
}'
```

And it could run like this:

```
stanFit <- rstan::stan(
model_code = model_code,
data = list(N=nrow(data),
logX=data$log2Dose,
y=data$y,
nBlocks=2,
cultureBlock=data$block),
chains = 4, # number of Markov chains
warmup = 1000, # number of warmup iterations per chain
iter = 5000, # total number of iterations per chain
cores = 4, # number of cores (could use one per chain)
refresh = 0, # no progress shown
seed = 1,
pars=c("sigma","meanB", "sigmaB","blockNoiseB"),
save_warmup=FALSE,
control=list(adapt_delta=0.8)
)
```

Following guidelines for investigating funnels etc. from @betanalpha writing I can get this to work by either constraing the prior or increase adapt_delta. In this case I can constrain the prior on sigmaB because I know that the true value is 0.2. Illustrations are below showing log(sigmaB) as a function of meanB - with the red dots corresponding to the divergent transitions - for various parameterizations and levels of adapt_delta:

**QUESTION**

The question I have is, from my reading of @betanalpha 's document there are tradeoffs depending on how informative the likelihood is, how much data we have, and the presence of an infinity-supressing prior. We need to adjust depending on our dataset. My problem is that I want to do some simulation studies with many datasets in order to study some particular behaviors. What strategies do people use in this case? If you want to do simulation studies with various designs (e.g. create and fit simulated datasets at various levels of the sigmaB parameter above) it seems to me that a single âone-size fits allâ stan model would not work. In addition sometimes the simulated dataset might be generated from a true parameter very close to zero, hence requiring a non-centered parameterization, or the other way around (high information parameter - more informed likelihood - centered parameterization more suitable). I can think of the following strategies:

- Have several models available (centered, non-centered) and for each dataset one of those will hopefully work.
- Use my knowledge of the true values of the simulated dataset to inform and constrain the prior, this sounds a bit like I could bias the results to the true answer but I am not quite sure what else to do.

Any other suggestions?