Non-linear power relationship model under a lognormal distribution - realtionship to parameters

I am unfortunately very confused about the relationship between parameters and the user supplied formula in a non-linear model under a log normal distribution. I have read the Bürkner, 2017 paper “Bayesian Distributional Non-Linear Multilevel Modeling with the R Package brms” but I fear I may lack the statistical background to understand how the equations therein actually relate to a real world example. I would be extremely grateful if someone can explain in simpleton friendly terms, and/or point me towards some resources to help me get a better grounding. The Package brms is fantastic, and its applications would be extremely useful for my work, so I am trying to get my head around it. Basically, I just want to get a mathematical representation of the equation using the model output.

For example, say I have this model that follows (what I hope is the correct syntax) for a power relationship a*x^b:

b <- c(2.12, 0.35)
x <- rnorm(500)
y <- rlnorm(500, mean = b[1] * x^ b[2])
dat1 <- data.frame(x, y) %>% filter(x > 0)

prior1 <- prior(normal(0, 1), nlpar = "b1", lb = 0) + 
  prior(normal(0, 1), nlpar = "b2", lb = 0) 

formula.fit.1.gaus <- brm(bf(y ~ b1 * x ^ b2, 
                        b1 ~ 1,
                        b2 ~ 1,
                        nl = TRUE),
                     data = dat1, 
                     prior = prior1, 
                     iter = 5000, 
                     warmup = 2000,
                     thin = 5, 
                     cores = 4,
                     family = gaussian(),
                     file = "formula.fit.1.gaussian",
                     sample_prior = "yes",
                     save_all_pars = T, 
                     seed = 1234)

plot(conditional_effects(formula.fit.1.gaus), points = TRUE)

Returning output

Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: y ~ b1 * x^b2 
         b1 ~ 1
         b2 ~ 1
   Data: dat1 (Number of observations: 262) 
Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 5;
         total post-warmup samples = 2400

Population-Level Effects: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b1_Intercept     6.53      0.82     4.97     8.09 1.00     2352     2220
b2_Intercept     1.47      0.26     0.88     1.94 1.00     2101     2046

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    12.87      0.62    11.69    14.17 1.00     2377     2333

Samples were drawn 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).

My understanding is the estimates from the model output cannot simply be put back into the user supplied formula, as each non-linear parameter is a placeholder for a linear equation? If so, how then can I obtain the values that can be used to summaries the model?

However, when I try to put the estimates into the power equation (b1*x^b2) and plot them ( curve(6.53*x^1.47, from = 0, to = 3.315768) ), at a glance it appears to be fairly similar to the model described using plot(conditional_effects(formula.fit.1.gaus), points = TRUE) . Can someone explain why this is? Is it just coincidence, or can the estimates sometimes be used, or have I just completely misunderstood?

Following this, if I do the same model, but under a lognormal distribution:


prior1 <- prior(normal(0, 1), nlpar = "b1", lb = 0) + 
  prior(normal(0, 1), nlpar = "b2", lb = 0) 

formula.fit.1.log <- brm(bf(y ~ b1 * x ^ b2, 
                        b1 ~ 1,
                        b2 ~ 1,
                        nl = TRUE),
                     data = dat1, 
                     prior = prior1, 
                     iter = 5000, 
                     warmup = 2000,
                     thin = 5, 
                     cores = 4,
                     family = lognormal(),
                     #control = list(adapt_delta = 0.99,
                     #              max_treedepth = 12), 
                     file = "formula.fit.1.log",
                     sample_prior = "yes",
                     save_all_pars = T, 
                     seed = 1234)

with results

 Family: lognormal 
  Links: mu = identity; sigma = identity 
Formula: y ~ b1 * x^b2 
         b1 ~ 1
         b2 ~ 1
   Data: dat1 (Number of observations: 262) 
Samples: 4 chains, each with iter = 5000; warmup = 2000; thin = 5;
         total post-warmup samples = 2400

Population-Level Effects: 
             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b1_Intercept     2.11      0.07     1.97     2.24 1.00     2607     2353
b2_Intercept     0.44      0.05     0.34     0.54 1.00     2138     2292

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.97      0.04     0.89     1.06 1.00     2394     2290

Samples were drawn 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).

I understand I cannot simply back transform parameters using exp() - but how could I go about getting the parameters from this model to place back into the supplied formula as a mathematical approximation?

Again would appreciate any insight whatsoever!