What is rescor in {brms} doing for this non-linear multivariate model

benjamin_hlina · November 23, 2025, 11:03pm

Hi everyone,

I am using {brms} to create a non-linear model to estimate trophic position based on Post 2002 and Heuvel et al. 2024 which proposes the following model:

\text{Trophic Position} = \lambda + \frac{(\delta^{15}N_c - [(\delta^{15}N_{b1} \times \alpha_r) + (\delta^{15}N_{b2} \times (1 - \alpha_r))]}{\Delta N}

This formula relies on the calculation of \alpha_r which we can do using the following:

\alpha_r = \frac{(\alpha - \alpha_{min})}{(\alpha_{max}-\alpha_{min})}

in {brms} this looks like this:

\alpha = \alpha_r \times (\alpha_{max} - \alpha_{min}) + \alpha_{min}

Then this can be used in the trophic position formula that has been rearranged to solve for \delta^{15}N_c

\delta^{15}N_c = \Delta N \times (\text{Trophic Position} - \lambda) + \delta^{15}N_{b1} \times \alpha_r + \delta^{15}N_{b2} \times (1 - \alpha_r)

I previously asked the following question as I was initially having a hard time wrapping my head around using the priors for the same variable in two models. I have since figured this out but now have a new question/problem.

A colleague pointed out that once the model runs the posterior draws for b_alpha_ar_Intercept make sense, however, draws for b_d15n_ar_Intercept are directly sampled from the priors supplied for this variable (in the current example it assumes a uninformed beta distribution; see posterior histograms and trace plots below). It is my understanding that rescor (residual correlation) is being applied in some manner to the second formula across the shared parameter hence why the posteriors for b_d15n_ar_Intercept are just random draws from the priors for this variable. Is this correct? What is rescor actually doing here. Thanks - Ben

I have have provided an example below:

# load packages 
suppressPackageStartupMessages(
  library(brms)
)
# test dataframe 
test <- structure(list(id = 1:30, common_name = structure(c(1L, 1L, 1L, 
 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), levels = "Lake Trout", class = "factor"), 
     ecoregion = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
     2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
     2L, 2L, 2L, 2L, 2L, 2L), levels = c("Anthropogenic", "Embayment"
     ), class = "factor"), d13c = c(-22.94658528, -22.51012416, 
     -22.78499464, -22.32669228, -22.45175064, -22.31867644, -22.28633518, 
     -22.45369046, -22.9138507, -22.25144698, -23.50689133, -22.32177069, 
     -22.24848683, -23.16317399, -22.80185743, -22.16461583, -21.86209349, 
     -21.82368985, -22.05221775, -20.76847322, -23.10855048, -22.16543141, 
     -22.27316888, -22.27460177, -21.77925327, -22.00708063, -22.75222408, 
     -20.80769779, -22.9086035, -21.37549007), d15n = c(15.8916573, 
     16.24570262, 16.9587455, 16.57961716, 16.5595928, 16.58978275, 
     16.60921641, 16.24422099, 16.42821653, 16.31949971, 16.07445668, 
     16.21850635, 16.40991133, 16.50936412, 17.24475461, 16.42799366, 
     16.45332704, 16.45004299, 16.32828884, 16.43756424, 17.24475899, 
     17.49526771, 17.64165936, 17.605025826512, 17.608486882222, 
     18.06772822, 17.577167139726, 16.98207199, 17.53002113, 17.463065668893
     ), d13c_b1 = c(-26.1767494183659, -26.6468331734217, -26.0146862822546, 
     -22.1256554796036, -24.2800883077986, -22.1070858870673, 
     -24.7339114652255, -26.567163857651, -24.6062345350996, -22.1070858870673, 
     -24.7339114652255, -22.1070858870673, -24.7339114652255, 
     -26.8903446551477, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
     NA, NA, NA, NA, NA, NA), d15n_b1 = c(8.44125347407797, 8.76982443827081, 
     8.05422509907782, 13.5929397010189, 6.98988123610252, 7.9508221409179, 
     7.37054834436797, 6.92934136835461, 6.9656336500113, 7.9508221409179, 
     7.37054834436797, 7.9508221409179, 7.37054834436797, 10.1514338013711, 
     NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
     NA), d13c_b2 = c(-23.41655522, -22.93946954, -22.70998529, 
     -23.41655522, -22.93946954, -22.70998529, -23.41655522, -22.93946954, 
     -22.70998529, -26.89034466, -23.45480259, -23.70341053, NA, 
     NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
     NA, NA), d15n_b2 = c(7.80516233, 7.611940133, 7.321611395, 
     7.80516233, 7.611940133, 7.321611395, 7.80516233, 7.611940133, 
     7.321611395, 10.1514338, 7.675356649, 7.639684858, NA, NA, 
     NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
     NA), c1 = c(-24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586, 
     -24.5593391261586, -24.5593391261586, -24.5593391261586), 
     n1 = c(8.27561744458162, 8.27561744458162, 8.27561744458162, 
     8.27561744458162, 8.27561744458162, 8.27561744458162, 8.27561744458162, 
     8.27561744458162, 8.27561744458162, 8.27561744458162, 8.27561744458162, 
     8.27561744458162, 8.27561744458162, 8.27561744458162, 8.27561744458162, 
     8.27561744458162, 8.27561744458162, 8.27561744458162, 8.27561744458162, 
     8.27561744458162, 8.27561744458162, 8.27561744458162, 8.27561744458162, 
     8.27561744458162, 8.27561744458162, 8.27561744458162, 8.27561744458162, 
     8.27561744458162, 8.27561744458162, 8.27561744458162), c2 = c(-23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333, -23.4372156608333, 
     -23.4372156608333, -23.4372156608333), n2 = c(7.80688474008333, 
     7.80688474008333, 7.80688474008333, 7.80688474008333, 7.80688474008333, 
     7.80688474008333, 7.80688474008333, 7.80688474008333, 7.80688474008333, 
     7.80688474008333, 7.80688474008333, 7.80688474008333, 7.80688474008333, 
     7.80688474008333, 7.80688474008333, 7.80688474008333, 7.80688474008333, 
     7.80688474008333, 7.80688474008333, 7.80688474008333, 7.80688474008333, 
     7.80688474008333, 7.80688474008333, 7.80688474008333, 7.80688474008333, 
     7.80688474008333, 7.80688474008333, 7.80688474008333, 7.80688474008333, 
     7.80688474008333), alpha = c(-0.437233865964201, -0.826193845402353, 
     -0.581238197923478, -0.989662381323953, -0.878214431196874, 
     -0.996805837679422, -1.02562731855865, -0.876485726593554, 
     -0.466405860857392, -1.05671854967365, 0.0620926941818009, 
     -0.994048342541306, -1.05935653924563, -0.244217039658719, 
     -0.566210626964422, -1.1340996513823, -1.40369773871247, 
     -1.43792181581869, -1.23426517101825, -2.3782966164598, -0.292895738293868, 
     -1.13337283296596, -1.03736069764468, -1.03608375259872, 
     -1.47752225317976, -1.27448990688269, -0.610442257024496, 
     -2.34334095319094, -0.471081995135074, -1.83734290792647), 
     min_alpha = c(-2.3782966164598, -2.3782966164598, -2.3782966164598, 
     -2.3782966164598, -2.3782966164598, -2.3782966164598, -2.3782966164598, 
     -2.3782966164598, -2.3782966164598, -2.3782966164598, -2.3782966164598, 
     -2.3782966164598, -2.3782966164598, -2.3782966164598, -2.3782966164598, 
     -2.3782966164598, -2.3782966164598, -2.3782966164598, -2.3782966164598, 
     -2.3782966164598, -2.3782966164598, -2.3782966164598, -2.3782966164598, 
     -2.3782966164598, -2.3782966164598, -2.3782966164598, -2.3782966164598, 
     -2.3782966164598, -2.3782966164598, -2.3782966164598), max_alpha = c(1.26738277086584, 
     1.26738277086584, 1.26738277086584, 1.26738277086584, 1.26738277086584, 
     1.26738277086584, 1.26738277086584, 1.26738277086584, 1.26738277086584, 
     1.26738277086584, 1.26738277086584, 1.26738277086584, 1.26738277086584, 
     1.26738277086584, 1.26738277086584, 1.26738277086584, 1.26738277086584, 
     1.26738277086584, 1.26738277086584, 1.26738277086584, 1.26738277086584, 
     1.26738277086584, 1.26738277086584, 1.26738277086584, 1.26738277086584, 
     1.26738277086584, 1.26738277086584, 1.26738277086584, 1.26738277086584, 
     1.26738277086584), lambda = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 
     2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
     2, 2)), row.names = c(NA, -30L), class = c("tbl_df", "tbl", 
 "data.frame"))


# model formula
tp_formula <- bf(
  alpha ~ ar * (max_alpha - min_alpha) + min_alpha,
  ar ~ 1,
  nl = TRUE
) +
  bf(
    d15n ~ dn * (tp - lambda) + n1 * ar + n2 * (1 - ar),
    ar ~ 1,
    tp ~ 1,
    dn ~ 1,
    nl = TRUE
  ) +
  brms::set_rescor()

# priors
tp_priors <- c(
  # Beta prior for 'ar'
  brms::prior(beta(1, 1), lb = 0, ub = 1, resp = "alpha", nlpar = "ar"),
  brms::prior(beta(1, 1), lb = 0, ub = 1, resp = "d15n", nlpar = "ar"),
  # 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), lb = 0, ub = 10, resp = "alpha", class = "sigma"),
  brms::prior(uniform(0, 10), lb = 0, ub = 10, resp = "d15n", class = "sigma")
)

# run model
tp_fit <- brm(
  formula = tp_formula,
  data = test,
  family = gaussian(),
  prior = tp_priors,
  iter = 4000,
  warmup = 1000,
  chains = 2,
  cores = 6,
  control = list(adapt_delta = 0.995),
)

# check summary
summary(tp_fit)
#>  Family: MV(gaussian, gaussian) 
#>   Links: mu = identity
#>          mu = identity 
#> Formula: alpha ~ ar * (max_alpha - min_alpha) + min_alpha 
#>          ar ~ 1
#>          d15n ~ dn * (tp - lambda) + n1 * ar + n2 * (1 - ar) 
#>          ar ~ 1
#>          tp ~ 1
#>          dn ~ 1
#>    Data: test (Number of observations: 30) 
#>   Draws: 2 chains, each with iter = 4000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 6000
#> 
#> Regression Coefficients:
#>                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> alpha_ar_Intercept     0.38      0.03     0.32     0.44 1.00     5592     3640
#> d15n_ar_Intercept      0.50      0.29     0.02     0.98 1.00     4929     3267
#> d15n_tp_Intercept      4.55      0.20     4.20     4.93 1.00     1989     2977
#> d15n_dn_Intercept      3.45      0.27     3.02     3.96 1.00     2040     2665
#> 
#> Further Distributional Parameters:
#>             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma_alpha     0.59      0.08     0.45     0.77 1.00     4545     4102
#> sigma_d15n      0.62      0.09     0.48     0.82 1.00     4644     3689
#> 
#> Residual Correlations: 
#>                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> rescor(alpha,d15n)    -0.13      0.18    -0.46     0.23 1.00     4149     3923
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

# plot trace and posteriors histograms
plot(tp_fit)

^{Created on 2025-11-23 with reprex v2.1.1}

amynang · November 24, 2025, 5:34pm

Hi Ben,

I’m afraid my answer to that other question was misleading, because it results in estimating two parameters (in the current version of your model alpha_ar and d15n_ar), while what you want is to estimate a_r jointly from the two formulas.

I think this is currently not possible in brms but you could modify the generated Stan code to achieve it.

benjamin_hlina · November 24, 2025, 7:29pm

Thank you for the clearification - shortly after I posted this I found the following posts:

Which both suggested the approach you suggested.

If I use make_stancode() I get the following:

two_source_model_ar()
#> alpha ~ ar * (max_alpha - min_alpha) + min_alpha 
#> ar ~ 1
#> d15n ~ dn * tp - l1 + n1 * ar + n2 * (1 - ar) 
#> ar ~ 1
#> tp ~ 1
#> dn ~ 1
two_source_priors_ar()
#>                        prior class coef group  resp dpar nlpar    lb       ub
#>                   beta(a, b)     b            alpha         ar     0        1
#>                   beta(a, b)     b             d15n         ar     0        1
#>         normal(dn, dn_sigma)     b             d15n         dn  <NA>     <NA>
#>        uniform(tp_lb, tp_ub)     b             d15n         tp tp_lb    tp_ub
#>  uniform(sigma_lb, sigma_ub) sigma            alpha             <NA> sigma_ub
#>  uniform(sigma_lb, sigma_ub) sigma             d15n             <NA> sigma_ub
#>  tag source
#>        user
#>        user
#>        user
#>        user
#>        user
#>        user

make_stancode(two_source_model_ar(), data = test, 
              priors = two_source_priors_ar())
#> // generated with brms 2.23.0
#> functions {
#> }
#> data {
#>   int<lower=1> N;  // total number of observations
#>   int<lower=1> N_alpha;  // number of observations
#>   vector[N_alpha] Y_alpha;  // response variable
#>   int<lower=1> K_alpha_ar;  // number of population-level effects
#>   matrix[N_alpha, K_alpha_ar] X_alpha_ar;  // population-level design matrix
#>   // covariates for non-linear functions
#>   vector[N_alpha] C_alpha_1;
#>   vector[N_alpha] C_alpha_2;
#>   int<lower=1> N_d15n;  // number of observations
#>   vector[N_d15n] Y_d15n;  // response variable
#>   int<lower=1> K_d15n_ar;  // number of population-level effects
#>   matrix[N_d15n, K_d15n_ar] X_d15n_ar;  // population-level design matrix
#>   int<lower=1> K_d15n_tp;  // number of population-level effects
#>   matrix[N_d15n, K_d15n_tp] X_d15n_tp;  // population-level design matrix
#>   int<lower=1> K_d15n_dn;  // number of population-level effects
#>   matrix[N_d15n, K_d15n_dn] X_d15n_dn;  // population-level design matrix
#>   // covariates for non-linear functions
#>   vector[N_d15n] C_d15n_1;
#>   vector[N_d15n] C_d15n_2;
#>   vector[N_d15n] C_d15n_3;
#>   int<lower=1> nresp;  // number of responses
#>   int nrescor;  // number of residual correlations
#>   int prior_only;  // should the likelihood be ignored?
#> }
#> transformed data {
#>   array[N] vector[nresp] Y;  // response array
#>   for (n in 1:N) {
#>     Y[n] = transpose([Y_alpha[n], Y_d15n[n]]);
#>   }
#> }
#> parameters {
#>   vector[K_alpha_ar] b_alpha_ar;  // regression coefficients
#>   real<lower=0> sigma_alpha;  // dispersion parameter
#>   vector[K_d15n_ar] b_d15n_ar;  // regression coefficients
#>   vector[K_d15n_tp] b_d15n_tp;  // regression coefficients
#>   vector[K_d15n_dn] b_d15n_dn;  // regression coefficients
#>   real<lower=0> sigma_d15n;  // dispersion parameter
#>   cholesky_factor_corr[nresp] Lrescor;  // parameters for multivariate linear models
#> }
#> transformed parameters {
#>   // prior contributions to the log posterior
#>   real lprior = 0;
#>   lprior += student_t_lpdf(sigma_alpha | 3, 0, 2.5)
#>     - 1 * student_t_lccdf(0 | 3, 0, 2.5);
#>   lprior += student_t_lpdf(sigma_d15n | 3, 0, 2.5)
#>     - 1 * student_t_lccdf(0 | 3, 0, 2.5);
#>   lprior += lkj_corr_cholesky_lpdf(Lrescor | 1);
#> }
#> model {
#>   // likelihood including constants
#>   if (!prior_only) {
#>     // initialize linear predictor term
#>     vector[N_alpha] nlp_alpha_ar = rep_vector(0.0, N_alpha);
#>     // initialize non-linear predictor term
#>     vector[N_alpha] mu_alpha;
#>     // initialize linear predictor term
#>     vector[N_d15n] nlp_d15n_ar = rep_vector(0.0, N_d15n);
#>     // initialize linear predictor term
#>     vector[N_d15n] nlp_d15n_tp = rep_vector(0.0, N_d15n);
#>     // initialize linear predictor term
#>     vector[N_d15n] nlp_d15n_dn = rep_vector(0.0, N_d15n);
#>     // initialize non-linear predictor term
#>     vector[N_d15n] mu_d15n;
#>     // multivariate predictor array
#>     array[N] vector[nresp] Mu;
#>     vector[nresp] sigma = transpose([sigma_alpha, sigma_d15n]);
#>     // cholesky factor of residual covariance matrix
#>     matrix[nresp, nresp] LSigma = diag_pre_multiply(sigma, Lrescor);
#>     nlp_alpha_ar += X_alpha_ar * b_alpha_ar;
#>     nlp_d15n_ar += X_d15n_ar * b_d15n_ar;
#>     nlp_d15n_tp += X_d15n_tp * b_d15n_tp;
#>     nlp_d15n_dn += X_d15n_dn * b_d15n_dn;
#>     for (n in 1:N_alpha) {
#>       // compute non-linear predictor values
#>       mu_alpha[n] = (nlp_alpha_ar[n] * (C_alpha_1[n] - C_alpha_2[n]) + C_alpha_2[n]);
#>     }
#>     for (n in 1:N_d15n) {
#>       // compute non-linear predictor values
#>       mu_d15n[n] = (nlp_d15n_dn[n] * nlp_d15n_tp[n] - C_d15n_1[n] + C_d15n_2[n] * nlp_d15n_ar[n] + C_d15n_3[n] * (1 - nlp_d15n_ar[n]));
#>     }
#>     // combine univariate parameters
#>     for (n in 1:N) {
#>       Mu[n] = transpose([mu_alpha[n], mu_d15n[n]]);
#>     }
#>     target += multi_normal_cholesky_lpdf(Y | Mu, LSigma);
#>   }
#>   // priors including constants
#>   target += lprior;
#> }
#> generated quantities {
#>   // residual correlations
#>   corr_matrix[nresp] Rescor = multiply_lower_tri_self_transpose(Lrescor);
#>   vector<lower=-1,upper=1>[nrescor] rescor;
#>   // extract upper diagonal of correlation matrix
#>   for (k in 1:nresp) {
#>     for (j in 1:(k - 1)) {
#>       rescor[choose(k - 1, 2) + j] = Rescor[j, k];
#>     }
#>   }
#> }

^{Created on 2025-11-24 with reprex v2.1.1}

Would I want to modify for instance K_d15n_ar in data to K_alpha_ar and then the same in priorsfound in themodel` section of the Stan code?

amynang · November 26, 2025, 9:26am

I think the main changes would be in the parameters block:

//vector[K_alpha_ar] b_alpha_ar;
//vector[K_d15n_ar] b_d15n_ar;
vector[K_alpha_ar] b_ar

and the model block:

//nlp_alpha_ar += X_alpha_ar * b_alpha_ar;
//nlp_d15n_ar += X_d15n_ar * b_d15n_ar;
nlp_alpha_ar += X_alpha_ar * b_ar;
nlp_d15n_ar += X_d15n_ar * b_ar;

But it would be great if someone more experienced with hacking brms-generated Stan code could double-check this.

benjamin_hlina · November 26, 2025, 1:13pm

Thanks for your help with this - this definitely points me in the right direction.

Topic		Replies	Views
Autoregression on non-linear model parameters in brms brms	16	1635	September 20, 2018
Brms - non-linear model estimate parameter used within another model brms fitting-issues , brms	2	86	February 12, 2025
Problems fitting multivariate model with option set_rescor brms brms	4	107	June 17, 2025
Fit multivariate non-linear model with shared parameters Modeling rstan , r , brms	2	202	August 22, 2024
Brms automatically scaling response variables in a multivariate model? brms fitting-issues , brms	5	870	October 6, 2022

What is rescor in {brms} doing for this non-linear multivariate model

Related topics