Hello!

I am new to Stan, and I am trying to fit a log binomial random intercept model. I was encouraged to use Stan by colleagues who suggested it was quite efficient, so I used the brms package in R to output the Stan code for a logistic random intercept model. I then updated the link to be log instead of logit (and added a line to directly calculate an additive interaction parameter of interest, the relative excess risk due to interaction). My code is below.

This code works well for the unadjusted model and when I adjust for one dichotomous covariate. However, when I try to adjust for a few more covariates, I get sampling errors. I believe this is due to the fact that I’m using a log link, and the risks are being pushed outside the parameter space (>1). So, I would like to add explicit constraints into the code such that the mu[n]'s are < 0 or exp(mu) < 1. However, I have not been able to figure out how to successfully incorporate that constraint into the Stan code.

Any suggestions would be much appreciated!

```
// generated with brms 1.6.1 (and then edited to change from logit to log link; and to directly calculate RERI)
functions {
}
data {
int<lower=1> N; // total number of observations
int Y[N]; // response variable
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
// data for group-level effects of ID 1
int<lower=1> J_1[N];
int<lower=1> N_1;
int<lower=1> M_1;
vector[N] Z_1_1;
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc;
matrix[N, K - 1] Xc; // centered version of X
vector[K - 1] means_X; // column means of X before centering
Kc = K - 1; // the intercept is removed from the design matrix
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b; // population-level effects
real temp_Intercept; // temporary intercept
vector<lower=0>[M_1] sd_1; // group-level standard deviations
vector[N_1] z_1[M_1]; // unscaled group-level effects
}
transformed parameters {
// group-level effects
vector[N_1] r_1_1;
r_1_1 = sd_1[1] * (z_1[1]);
}
model {
vector[N] mu;
mu = Xc * b + temp_Intercept;
for (n in 1:N) {
mu[n] = mu[n] + (r_1_1[J_1[n]]) * Z_1_1[n];
}
// prior specifications
sd_1 ~ gamma(2,0.1);
z_1[1] ~ normal(0, 1);
// likelihood contribution
if (!prior_only) {
Y ~ bernoulli(exp(mu));
}
}
generated quantities {
real b_Intercept; // population-level intercept
real RERI; // relative excess risk due to interaction
b_Intercept = temp_Intercept - dot_product(means_X, b);
RERI = exp(b[1]+b[2]+b[6]) - exp(b[1]) - exp(b[2]) + 1;
}
```

[edit: escaped code]