# brms: obtaining the coefficients of the conditional mean in a non-linear model

I'm using brms to fit a non-linear model to a set of data representing biexponential decay ($$y_i = a_1 \cdot e^{-k_1\cdot x_i}+a_2 \cdot e^{-k_2\cdot x_i}$$).

It seems that the parameter estimates listed in the model summary produces a curve that is different from the conditional mean of the model, and fits the data more poorly. I want to obtain the parameters of the conditional mean.

Here's an example using simulated data:

library(brms)

# True parameters
a1 <- 70
a2 <- 20
k1 <- .7
k2 <- .05

# CV of measurement errors
cv <- .15

x <- seq(0, 30, len = 6)
y <- (a1 * exp(-k1 * x) + a2 * exp(-k2 * x)) * rnorm(6, 1, cv)

dat <- data.frame(x = x, y = y, cv = cv)

priors <-
prior(lognormal(log(100), 1.2),
lb = 0, ub = 1000,
nlpar = "a1") +
prior(lognormal(log(100), 1.2),
lb = 0, ub = 1000,
nlpar = "a2") +
prior(lognormal(log(1), .75),
lb = 0, ub = 10,
nlpar = "k1") +
prior(lognormal(log(.05), .75),
lb = 1, ub = 50,
nlpar = "k2")

formula <- bf(log(y) ~ log(a1 * exp(-k1 * x) + a2 * exp(-k2 * x)),
a1 + a2 + k1 + k2 ~ 1,
nl = TRUE)

mod <- brm(formula, data = dat, prior = priors)


Here's the model summary:

> summary(mod)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: log(y) ~ log(a1 * exp(-k1 * x) + a2 * exp(-k2 * x))
a1 ~ 1
a2 ~ 1
k1 ~ 1
k2 ~ 1
Data: dat (Number of observations: 6)
Draws: 4 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 2000

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
a1_Intercept    51.46     59.64    14.76   219.88 1.01      375      644
a2_Intercept    77.66     74.03    12.00   267.62 1.01      950      993
k1_Intercept     0.09      0.04     0.04     0.19 1.01      352      779
k2_Intercept     1.20      0.21     1.00     1.79 1.00      924      566

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.68      0.51     0.17     2.10 1.01      242      254


Here's a plot of the conditional effects:

Now, I recreate the plot of the conditional effects with the response variable back-transformed into the natural scale. If I extract the estimates of the fixed effects using fixef(mod) and use those to plot a biexponential decay curve along with the plot of the conditional effects, it becomes apparent that the two curves are not identical:

f <- fixef(mod)[, 1]
xnew <- seq(0, 30, len = 100)
yfit <- exp(fitted(mod, robust = TRUE, newdata = data.frame(x = xnew, cv = cv)))
fitdata <- as.data.frame(cbind(
x = xnew,
yfit,
Estimate2 = f[1] * exp(-f[3] * xnew) +  f[2] * exp(-f[4] * xnew)))

library(ggplot2)

ggplot() +
geom_ribbon(
data = fitdata,
aes(x = x, ymin = Q2.5, ymax = Q97.5),
fill = "#00000022") +
geom_line(
data = fitdata,
aes(x = x, y = Estimate),
colour = "blue",
linewidth = 1) +
geom_line(
data = fitdata,
aes(x = x, y = Estimate2),
colour = "red",
linewidth = 1) +
geom_point(
data = dat,
aes(x = x, y = y)) +
scale_y_log10()


Plot of conditional effects (blue line) and biexponential decay curve using parameter estimates from the model (red line):

The red line based on the parameter estimates obviously fits more poorly. My questions are:

1. Why do the estimates of the population-level effects not produce the conditional mean of the model?
2. How can I obtain the coefficients of the conditional mean?

This is a well-known result of the non-linear relationships in your model. I.e. even if on some suitable scale your posterior distribution for the parameters is a set of independent normal distributions, the uncertainty around them will no longer cancel out once you apply non-linear functions to them.

A way of getting the same results as you did from conditional_effects(mod) is to take the posterior MCMC samples and to calculate the curve across values of x for the parameters in each MCMC sample (of course, that's 4000 curves across however many values of x you want to evaluate the curve at). Then, you could take the median of all curves (plus quantiles for credible intervals). E.g. with brms + tidyverse functions you could do that like this:

library(tidyverse)
mod %>%
as_draws_df() %>%
mutate(logy_pred = pmap(list(b_a1_Intercept, b_a2_Intercept, b_k1_Intercept, b_k2_Intercept),
function(a1, a2, k1, k2) tibble(x=seq(0, 30, 0.1)) %>% mutate(logy_pred = log(a1 * exp(-k1 * x) + a2 * exp(-k2 * x))))) %>%
unnest(logy_pred) %>%
group_by(x) %>%
summarize(median = median(logy_pred),
lcri = quantile(logy_pred, probs=0.025),
ucri = quantile(logy_pred, probs=0.975)) %>%
ggplot(aes(x=x, y=median, ymin=lcri, ymax=ucri)) +
geom_ribbon() +
geom_line() +
ylab("log(y)")


That exactly reproduces the plot you obtained originally.

Note that another confusing thing was that one of your plots is on the $$y$$-scale and the other on the $$log(y)$$-scale. However, you can calculate quantiles (median etc.) before and after the transformation without it making a difference in this case.

• Thank you! It is nice to know this way of reproducing the plot, but it doesn't really answer my question about how to obtain the coefficients of the conditional mean. I guess I could do a frequentist, non-linear regression of the conditional mean/median. It would be more satisfying, though, if the coefficients could be extracted directly from the brms model. Jan 5, 2023 at 8:36
• as_draws_df() did directly extract the MCMC samples for the coefficients. Look at the df & you can see them. But, because of the the joint posterior with correlation between different parameters going into a non-linear function, there's no reason that a single point estimate vector for the parameters would exist that exactly reproduces the conditional means curve when plugged into the non-linear function. Remember: the curve is giving the distribution of predictions at each x value across the curves for each MCMC sample, why would such an average be representable in the same functional form? Jan 5, 2023 at 12:02