Brms - non-linear model estimate parameter used within another model

Interfaces brms
,
1

Hi everyone,

I am trying to create a non-linear model to estimate trophic position based on Post 2002 and Heuvel et al. 2024. I can write this model in JAGS and use {rjags} to estimate ac and tp but would like to move over to stan and use {brms}. However, I consistently run into the following error message Error: The parameter 'd15n' is not a valid distributional or non-linear parameter. Did you forget to set 'nl = TRUE'? . From what I can understand {brms} has issues with circular phased estimates considering I’m first estimating ac which is then to be used in the formula for tp. Is that correct? If so how can I get around this? I’ve thought of just inserting the formula for ac into the tp but this doesn’t seem to doing what I expect. Would it be best to write this in {rstan} instead?

If this question has been asked and answered before, sorry for taking up your time. Thanks

suppressPackageStartupMessages(
  library(brms)
)

test <- structure(list(common_name = c("Lake Trout", "Lake Trout", "Lake Trout", 
"Lake Trout", "Lake Trout", "Lake Trout", "Lake Trout", "Lake Trout", 
"Lake Trout", "Lake Trout", "Lake Trout", "Lake Trout", "Lake Trout", 
"Lake Trout", "Lake Trout"), d13c = c(-22.24432219, -22.19363157, 
-23.12710892, -21.94797704, -22.09868526, -23.15772812, -22.46564577, 
-22.06822162, -22.58020397, -22.78587869, -23.2490471, -22.09275059, 
-22.27354166, -21.77119478, -22.02081495), d15n = c(18.0137765, 
17.9057604, 17.04653059, 18.08178664, 17.44761289, 17.20065531, 
17.98916764, 17.602950331349, 16.84128291, 17.14235167, 16.65215117, 
17.7093308, 16.77823155, 18.284430048708, 17.12016153), n1 = c(8.823124606, 
8.823124606, 8.823124606, 8.823124606, 8.823124606, 8.823124606, 
8.823124606, 8.823124606, 8.823124606, 8.823124606, 8.823124606, 
8.823124606, 8.823124606, 8.823124606, 8.823124606), n2 = c(12.629882536375, 
12.629882536375, 12.629882536375, 12.629882536375, 12.629882536375, 
12.629882536375, 12.629882536375, 12.629882536375, 12.629882536375, 
12.629882536375, 12.629882536375, 12.629882536375, 12.629882536375, 
12.629882536375, 12.629882536375), lambda = c(2, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2), alpha = c(1.23643632363876, 1.24970560263278, 
1.00534932149529, 1.31401056664222, 1.27455968957911, 0.99733413629192, 
1.17850047020125, 1.28253415385541, 1.14851258059064, 1.0946731278142, 
0.973429575190285, 1.27611320760606, 1.22878754563037, 1.36028684278827, 
1.29494379355616), amin = c(0.885791396936097, 0.885791396936097, 
0.885791396936097, 0.885791396936097, 0.885791396936097, 0.885791396936097, 
0.885791396936097, 0.885791396936097, 0.885791396936097, 0.885791396936097, 
0.885791396936097, 0.885791396936097, 0.885791396936097, 0.885791396936097, 
0.885791396936097), amax = c(1.53574192008444, 1.53574192008444, 
1.53574192008444, 1.53574192008444, 1.53574192008444, 1.53574192008444, 
1.53574192008444, 1.53574192008444, 1.53574192008444, 1.53574192008444, 
1.53574192008444, 1.53574192008444, 1.53574192008444, 1.53574192008444, 
1.53574192008444)), row.names = c(NA, -15L), spec = structure(list(
    cols = list(common_name = structure(list(), class = c("collector_character", 
    "collector")), d13c = structure(list(), class = c("collector_double", 
    "collector")), d15n = structure(list(), class = c("collector_double", 
    "collector")), n1 = structure(list(), class = c("collector_double", 
    "collector")), n2 = structure(list(), class = c("collector_double", 
    "collector")), lambda = structure(list(), class = c("collector_double", 
    "collector")), alpha = structure(list(), class = c("collector_double", 
    "collector")), amin = structure(list(), class = c("collector_double", 
    "collector")), amax = structure(list(), class = c("collector_double", 
    "collector"))), default = structure(list(), class = c("collector_guess", 
    "collector")), delim = ","), class = "col_spec"), class = c("spec_tbl_df", 
"tbl_df", "tbl", "data.frame"))

tp_formula <- bf(
  alpha ~ ac * (amax - amin) + amin,
  d15n ~ dn * (tp - (lambda + n1 * ac + n2 * (1 - ac))),
  ac ~ 1,
  n1 ~ 1,
  n2 ~ 1,
  tp ~ 1,
  dn ~ 1,
  nl = TRUE
)

tp_priors <- c(
  # Beta prior for 'ac'
  brms::prior(beta(1, 1), class = "b", lb = 0, ub = 1, nlpar = "ac"),
  # Baseline (n1)
  brms::prior(normal(8, 1), class = "b", coef = "Intercept", nlpar = "n1"),
  # Baseline (n2)
  brms::prior(normal(10, 1), class = "b", coef = "Intercept", nlpar = "n2"),
  # Trophic enrichment factor (ΔN)
  brms::prior(normal(3.4, 0.5), class = "b", lb = 3, ub = 4, nlpar = "dn"),
  # Trophic Position (tp)
  brms::prior(uniform(2, 10), class = "b", nlpar = "tp", lb = 2, ub = 10),
  # Standard deviation prior
  brms::prior(uniform(0, 10), class = "sigma")
)


tp_fit <- brm(
  formula = tp_formula,
  data = test,
  family = gaussian(),
  prior = tp_priors,
  iter = 4000,
  warmup = 1000,
  chains = 2,
  cores = 4,
  control = list(adapt_delta = 0.95)
)
#> Error: The parameter 'd15n' is not a valid distributional or non-linear parameter. Did you forget to set 'nl = TRUE'?

Created on 2025-02-12 with reprex v2.1.1

2

Hey @benjamin_hlina , welcome!

This worked: I separated your formula into two distinct formulas and specified the response in the priors. (I also upped the iterations to 40000.)

tp_formula <- bf(
  alpha ~ ac * (amax - amin) + amin,
  ac ~ 1,
  nl = TRUE
) + 
  bf(
    d15n ~ dn * (tp - (lambda + n1 * ac + n2 * (1 - ac))),
    ac ~ 1,
    n1 ~ 1,
    n2 ~ 1,
    tp ~ 1,
    dn ~ 1,
    nl = TRUE
  )

tp_priors <- c(
  # Beta prior for 'ac'
  brms::prior(beta(1, 1), lb = 0, ub = 1, resp = "alpha", nlpar = "ac"),
  brms::prior(beta(1, 1), lb = 0, ub = 1, resp = "d15n", nlpar = "ac"),
  # Baseline (n1)
  brms::prior(normal(8, 1), resp = "d15n", coef = "Intercept", nlpar = "n1"),
  # Baseline (n2)
  brms::prior(normal(10, 1), resp = "d15n", coef = "Intercept", nlpar = "n2"),
  # Trophic enrichment factor (ΔN)
  brms::prior(normal(3.4, 0.5), resp = "d15n", lb = 3, ub = 4, nlpar = "dn"),
  # Trophic Position (tp)
  brms::prior(uniform(2, 10), resp = "d15n", nlpar = "tp", lb = 2, ub = 10),
  # Standard deviation prior
  brms::prior(uniform(0, 10), resp = "alpha", class = "sigma"),
  brms::prior(uniform(0, 10), resp = "d15n", class = "sigma")
)

(I had to comment out problems = < pointer:0x120723a80 > , for your test data to work)

1 Like
3

Thank you so so much! This makes a lot of sense, I just kept getting stuck in having it make it all in one.