For some reasons, I am trying to fit a model with a Gaussian-shape optimum on binary data. This can be done quite easily and nicely when no random effects are involved, but I am now trying to include random effects on the parameter W_{max} that constrains the latent model between 0 and 1. This, of course, is difficult as the random effects are not individually constrained, only their sum must be constrained between 0 and 1. Worse, the aim would to include more than one effect.

But let’s focus on a single random effect for now. I figured it would be more easy on a log scale as this allow for only one boundary (namely 0) rather than two (0 and 1). The following model runs, but results in a lot of divergent transitions (roughly 25% with adapt_delta at 0.80, needless to say that increasing this value doesn’t result in total remove of such transitions):

```
functions {
}
data {
int<lower=1> Nobs; // Total number of observations
int<lower=0> W[Nobs]; // Response variable
vector[Nobs] z; // Covariate
int<lower=1> Nind; // Total number of individuals
int<lower=1,upper=Nind> J_ind[Nobs]; // Grouping indices for ID
}
parameters {
real<lower=0> omega; // Peak Width
real theta; // Values of the optimum parameter
real<upper=0> log_Wmax; // Values of the max fitness parameter
real<lower=0> ind_sigma; // Sigma of individuals effects
vector[Nind] ind_extend; // Individual effects
}
transformed parameters {
vector<upper=0>[Nind] log_w_ind; // Effect of individual on log_Wmax
vector<lower=0,upper=1>[Nobs] prob;
// Wmax with individual effects
log_w_ind = log_Wmax + ind_extend * ind_sigma;
for (n in 1:Nobs) {
prob[n] = exp(log_w_ind[J_ind[n]] -((z[n] - theta) .* (z[n] - theta)) ./ (2 * omega * omega));
}
}
model {
// Likelihood
W ~ bernoulli(prob);
// Priors
omega ~ gamma(3.36,0.78);
theta ~ normal(0,1000);
log_Wmax ~ normal(0,1);
ind_extend ~ normal(0,1);
ind_sigma ~ normal(0,1);
}
generated quantities {
}
```

Any guidance would be appreciated here. Am I being unreasonable in wanting to fit this model? Is there any trick I could use to make the boundary issue less painful?

FYI, a (unconstrained) logit-link model on the same data runs nicely.

Thanks!