Using brms to create generated quantities Stan code block

Interfaces brms
1

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.

2

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?

3

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.

...
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

a3 = (log(0.05 / 0.95) ./ -nlp_a2) / 2;
4

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.

5

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!

6

Inside bf() as any other formula just wrapped inside nlf(). See examples in ?brmsformula

1 Like
7

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.

8

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

9

Add nlf(a3 ~ a1 - a2) to your model formula.

10

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)
11

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)