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.
---
### Author:  [@Stephen\_Martin](https://discourse.mc-stan.org/u/Stephen_Martin)
#### Post date: [November 6, 2019, 12:12am UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/2 "2019-11-06T00:12:33Z")
Looks like a3 is supposed to be a vector of length N (vector[N] a3), and you want to do element-wise division (e.g., …nlp\_a1 ./ (.95 \* nlp\_a1) - 1…). No?
---
### Author:  [@nhuurre](https://discourse.mc-stan.org/u/nhuurre)
#### Post date: [November 6, 2019, 7:30am UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/3 "2019-11-06T07:30:42Z")
The `generated quantities` block cannot use any variables declared in the `model` block so you need to move `nlp_a1` and `nlp_a2` into `transformed parameters`. And as Stephen said, `a3` should be vector, not real, and you need to use element-wise division.
```stan
...
transformed parameters {
// initialize linear predictor term
vector[N] nlp_a1 = X_a1 * b_a1;
// initialize linear predictor term
vector[N] nlp_a2 = X_a2 * b_a2;
}
model {
// initialize non-linear predictor term
vector[N] mu;
...
}
generated quantities {
vector[N] a3;
a3 = (log((nlp_a1 ./ (0.95 * nlp_a1) - 1)) ./ -nlp_a2) / 2;
}
```
Btw, isn’t that last statement equivalent to just
```stan
a3 = (log(0.05 / 0.95) ./ -nlp_a2) / 2;
```
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [November 6, 2019, 7:39am UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/4 "2019-11-06T07:39:54Z")
You can add it as another non-linear parameter wrapped inside the nlf() function to indicate that the formula for it is itself non-linear. This is slightly less efficient as doing it in generated quantities but you dont have to move outside of brms in this case. Or alternitvly you compute it after model fitting in R.
---
### Author:  [@wittja01](https://discourse.mc-stan.org/u/wittja01)
#### Post date: [November 6, 2019, 5:50pm UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/5 "2019-11-06T17:50:03Z")
Thanks for your response Paul - could you clarify where in the brms() function I would put the function for the non-linear parameter? I’ve tried a few things but they don’t seem to work:
```
# attempt 1
brms_mod <- brm(
bf(trap_catch ~ a1 / (1 + exp(-a2 * r)),
a1 + a2 ~ 1,
bf(a3 ~ (log(a1 / (0.95 * a1) - 1) / -a2) * 0.5,
nl = TRUE),
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),
seed = 314
)
# attempt 2
brms_mod <- brm(
bf(trap_catch ~ a1 / (1 + exp(-a2 * r)),
a1 + a2 ~ 1,
nlf(a3 ~ (log(a1 / (0.95 * a1) - 1) / -a2) * 0.5),
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),
seed = 314
)
# attempt 3
brms_mod <- brm(
bf(trap_catch ~ a1 / (1 + exp(-a2 * r)),
a1 + a2 ~ 1,
a3 ~ (log(a1 / (0.95 * a1) - 1) / -a2) * 0.5),
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),
seed = 314
)
```
I would prefer to stay within the brms framework, rather than using the brms generated Stan code with some tweaks. Thanks for your help!
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [November 6, 2019, 5:58pm UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/6 "2019-11-06T17:58:25Z")
Inside bf() as any other formula just wrapped inside nlf(). See examples in ?brmsformula
---
### Author:  [@wittja01](https://discourse.mc-stan.org/u/wittja01)
#### Post date: [November 7, 2019, 5:50pm UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/7 "2019-11-07T17:50:32Z")
Thanks everyone for the input. Earlier commenters were helpful with editing my raw stan code. I’ve opted to use the brms function and just calculated my generated quantity afterward by pulling out the posterior estimates for a1 and a2 to calculate a3.
---
### Author:  [@Luwis\_Diya](https://discourse.mc-stan.org/u/Luwis_Diya)
#### Post date: [January 10, 2020, 8:47am UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/8 "2020-01-10T08:47:41Z")
> [@paul.buerkner](#):
>
> ?brmsformula
Hi Paul,
I am struggling with a similar example: say I want a3 = a1 - a2, say. How can I go about doing this in brms. You can advise me based on this example.
Regards,
Luwis
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [January 10, 2020, 9:27am UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/9 "2020-01-10T09:27:33Z")
Add `nlf(a3 ~ a1 - a2)` to your model formula.
---
### Author:  [@PhDemetri](https://discourse.mc-stan.org/u/PhDemetri)
#### Post date: [February 14, 2020, 2:05am UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/10 "2020-02-14T02:05:28Z")
Still a little unclear. I’ve read `?brmsformula`. In the documentation, I see something like
```
bf(y ~ eta, nl = TRUE) +
lf(eta ~ 1 + x) +
nlf(sigma ~ tau * sqrt(eta)) +
lf(tau ~ 1)
```
However, when I run
```
model = brm(bf(Concentration_scaled ~ y,
b1 + b2 + b0 ~ 1|Subject,
nl = T)+
nlf(y~2.5*exp(b0 + b1*Time + b2/Time)),
data = d,
family = Gamma(),
prior = prior)
```
The variable `y` is not accessible. This is also the case when I add the expression to my formula via
```
model = brm(bf(Concentration_scaled ~ y,
nlf(y~2.5*exp(b0 + b1*Time + b2/Time)),
b1 + b2 + b0 ~ 1|Subject,
nl = T),
data = d,
family = Gamma(),
prior = prior)
```
---
### Author:  [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)
#### Post date: [February 20, 2020, 9:50am UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/11 "2020-02-20T09:50:22Z")
The following works for me in the latest release or github version.
```
d <- data.frame(
Concentration_scaled = rgamma(100, 1, 1),
Subject = rep(1:10, each = 10),
Time = rnorm(100)
)
prior <- set_prior("normal(0,1)", nlpar = c("b0", "b1", "b2"))
model = brm(bf(Concentration_scaled ~ y,
b1 + b2 + b0 ~ 1|Subject,
nl = T)+
nlf(y~2.5*exp(b0 + b1*Time + b2/Time)),
data = d,
family = Gamma(),
prior = prior)
```
---
### Author:  [@amynang](https://discourse.mc-stan.org/u/amynang)
#### Post date: [March 23, 2025, 5:51pm UTC](https://discourse.mc-stan.org/t/using-brms-to-create-generated-quantities-stan-code-block/11766/12 "2025-03-23T17:51:40Z")
> [@paul.buerkner](#):
>
> You can add it as another non-linear parameter wrapped inside the nlf() function to indicate that the formula for it is itself non-linear. This is slightly less efficient as doing it in generated quantities but you dont have to move outside of brms in this case. Or alternitvly you compute it after model fitting in R.
I still think that this is not possible or, if it is, how exactly it should be specified remains unclear. I am returning to the OP’s example, adding some dummy data. Since this is an old thread, I also want to emphasize that a3 is a generated quantity, not present in the data.
```r
library(brms)
data = data.frame(r = runif(50, 1, 5),
trap_catch = NA)
data$trap_catch = rpois(50, 5 / (1 + exp(-2 * data$r)))
# the original model
brms_mod <- brm(
bf(trap_catch ~ a1 / (1 + exp(-a2 * r)),
a1 + a2 ~ 1,
nl = TRUE),
data = data,
family = poisson(),
prior = c(prior(normal(0, 10), nlpar = "a1"),
prior(normal(0, 2), nlpar = "a2")),
iter = 4000,
warmup = 2000,
chains = 4,
cores = 4,
backend = "cmdstanr")
# attempts to modify it to generate a3
brms_mod <- brm(
bf(trap_catch ~ a1 / (1 + exp(-a2 * r)),
a1 + a2 ~ 1,
nl = TRUE) +
nlf(a3 ~ (log(a1 / (0.95 * a1) - 1) / -a2) * 0.5),
data = data,
family = poisson(),
prior = c(prior(normal(0, 10), nlpar = "a1"),
prior(normal(0, 2), nlpar = "a2")),
iter = 4000,
warmup = 2000,
chains = 4,
cores = 4,
backend = "cmdstanr")
#Error: The parameter 'a3' is not a valid distributional or non-linear parameter.
#Did you forget to set 'nl = TRUE'?
brms_mod <- brm(
bf(trap_catch ~ a1 / (1 + exp(-a2 * r)),
a1 + a2 ~ 1,
nl = TRUE) +
nlf(a3 ~ (log(a1 / (0.95 * a1) - 1) / -a2) * 0.5, nl = TRUE),
data = data,
family = poisson(),
prior = c(prior(normal(0, 10), nlpar = "a1"),
prior(normal(0, 2), nlpar = "a2")),
iter = 4000,
warmup = 2000,
chains = 4,
cores = 4,
backend = "cmdstanr")
#Error: The parameter 'a3' is not a valid distributional or non-linear parameter.
#Did you forget to set 'nl = TRUE'?
#In addition: Warning message:
#Arguments '...' and 'flist' in nlf() will be reworked at some point.
#Please avoid using them if possible.
brms_mod <- brm(
bf(trap_catch ~ a1 / (1 + exp(-a2 * r)),
a1 + a2 ~ 1,
nlf(a3 ~ (log(a1 / (0.95 * a1) - 1) / -a2) * 0.5),
nl = TRUE),
data = data,
family = poisson(),
prior = c(prior(normal(0, 10), nlpar = "a1"),
prior(normal(0, 2), nlpar = "a2")),
iter = 4000,
warmup = 2000,
chains = 4,
cores = 4,
backend = "cmdstanr")
#Error: The parameter 'a3' is not a valid distributional or non-linear parameter.
#Did you forget to set 'nl = TRUE'?
```