1
$\begingroup$

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: 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): 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?
$\endgroup$

1 Answer 1

1
$\begingroup$

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. enter image description here

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.

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Jan 5, 2023 at 8:36
  • $\begingroup$ 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? $\endgroup$
    – Björn
    Jan 5, 2023 at 12:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.