Hi - I’m pretty new to Stan/brms so apologies if this is a relatively simple question. I’ve tried googling the answer to this to no avail.

I have the following non-linear model I’m attempting to fit:

```
brms_mod <- brm(
bf(trap_catch ~ a1 / (1 + exp(-a2 * r)),
a1 + a2 ~ 1,
nl = TRUE),
data = global_fit,
family = poisson(),
prior = c(prior(normal(0, 10), nlpar = "a1"),
prior(normal(0, 2), nlpar = "a2")),
cores = parallel::detectCores() - 2,
control = list(adapt_delta = 0.99)
)
```

Essentially, trap_catch_i \sim Poisson(\mu) and

log(\mu) = \frac{a1}{1 + e^{-a2 * r}} * \frac{1}{2}

and I’m trying to get estimates for a1 and a2. I want to use these estimates to estimate another parameter that is a function of both a1 and a2, a3 = \frac{log(a1 / (0.95 * a1) - 1}{2}. If I understand writing Stan code correctly, this is something I would put in the generated quantities block. Using `r stancode(brms_mod)`

I can get the Stan code that `r brm()`

creates. I’m not sure if it’s possible to specify a custom generated quantity within the `r brm() `

function, which is what I would like to be able to do.

I have tried to take what `r standcode() `

produces and add my own generated quantities block, as seen below.

```
data {
int<lower=1> N; // number of observations
int Y[N]; // response variable
int<lower=1> K_a1; // number of population-level effects
matrix[N, K_a1] X_a1; // population-level design matrix
int<lower=1> K_a2; // number of population-level effects
matrix[N, K_a2] X_a2; // population-level design matrix
// covariate vectors
vector[N] C_1;
int prior_only; // should the likelihood be ignored?
}
transformed data {
}
parameters {
vector[K_a1] b_a1; // population-level effects
vector[K_a2] b_a2; // population-level effects
}
transformed parameters {
}
model {
// initialize linear predictor term
vector[N] nlp_a1 = X_a1 * b_a1;
// initialize linear predictor term
vector[N] nlp_a2 = X_a2 * b_a2;
// initialize non-linear predictor term
vector[N] mu;
for (n in 1:N) {
// compute non-linear predictor values
mu[n] = nlp_a1[n] / (1 + exp( - nlp_a2[n] * C_1[n]));
}
// priors including all constants
target += normal_lpdf(b_a1 | 0, 10);
target += normal_lpdf(b_a2 | 0, 2);
// likelihood including all constants
if (!prior_only) {
target += poisson_log_lpmf(Y | mu);
}
}
generated quantities {
real a3;
a3 = (log((nlp_a1 / (0.95 * nlp_a1) - 1)) / -nlp_a2) / 2;
}
```

But this doesn’t seem to work. I get errors about the a3 variable. The error says the function defining a3 in the generated quantities block is ill-formed. I gather from googling that that has something to do with the types of quantities a3, nlp_a1 and nlp_a2 are.

It’s entirely possible I’m not approaching this correctly at all, so any advice is appreciated! Thank you for your patience and help.