 # Parameter estimate for each categorical variable

Hi all, I have a general question:
I have a non linear model y=a*x^b and a categorical variable ID. Is it possible in brms to estimate the model coefficients a and b for each categorical variable ID?

I was looking at the example of the non linear model, is this a correct form?

``````   bf(y ~ a *x^b,
a~ 1 + (1|ID),b~ 1 + (1|ID),  nl = TRUE),
data = data, family = gaussian(),
prior = c(
prior(normal(0.5, 0.2), nlpar = "a"),
prior(normal(0.5,0.2), nlpar = "b"))
)``````

I believe you should be able to extract your ID coefficients by using `ranef(fit)` assuming you called your model `fit <- brm(...)`.

Thanks A lot dylanjm.
Running ranef It is a bit confusing for me the fact that the coefficients are called intercepts. Is it assumed a given slope for my a and b coefficients? and if yes, how I can retrive it?

Thanks again and best regards

I think your current model only specifies random intercepts and not random slopes. So each ID should have its own intercept but they should all share the Population-Level coefficients. At least that’s my understanding.

Also `"a"` in this case represents your model intercept, so intercept is the estimate for `a` for each ID.

Have you considered fitting the log-transformed model (i.e., \log y \sim \log a + b * \log x)? It would be a lot easier to parameterize with group-level effects for a and b (`log_y ~ (1 + log_x | ID)`), you could back-transform the intercept into a after fitting, and it would probably be more efficient. I guess this would require all y to be positive and need a different likelihood distribution.

Thanks a lot Christopher, unfortunately my non linear model is slightly more complex than the example I made. Probably this solution is not working for my case. But thanks a lot,

best regards

Thanks again dylanjm. Are the value of ‘a’ and ‘b’ that i retrive from ranef the values I can put in my formula and I can run the model with?I.e for id=1 i have y=a1x^b1, for id=2 y=a2x^b2
, etc…?

Yes, in theory that should be the case. But like I said, if you only have a varying intercept model then you will have a1, a2, a3, … aN, but they will all share the same b value for the slope.

Thanks a lot dylanjm. In the case of my model y=a*x^b a is a non linear parameter so, I was thinking it is not an intercept. Could you please indicate that you also means for slope? in your answer? because b for my model is also a non linear parameter.

Thanks again

I suppose in the formal sense they are not necessarily intercept and slope coefficients. Here is a quick simulation of a non-linear model without using a multi-level structure.

``````library(tidyverse)
library(brms)

N <- 1200
x <- runif(N, 1, 20)
a <- 50
b <- -1.3
y_true <- a * x^b
sigma <- 3
y <- rnorm(N, y_true, sigma)

dat <- tibble(x, y)

fit <- brm(bf(y ~ a*x^b, nl = TRUE,
a ~ 1, b ~ 1),
data = dat, family = gaussian(),
prior = c(prior(normal(40, 10), nlpar = "a"),
prior(normal(-2, 1), nlpar = "b")),
cores = 24)

summary(fit)
#>  Family: gaussian
#>   Links: mu = identity; sigma = identity
#> Formula: y ~ a * x^b
#>          a ~ 1
#>          b ~ 1
#>    Data: dat (Number of observations: 1200)
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#>
#> Population-Level Effects:
#>             Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> a_Intercept    50.56      0.66    49.24    51.83       2169 1.00
#> b_Intercept    -1.33      0.02    -1.37    -1.30       2231 1.00
#>
#> Family Specific Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> sigma     3.01      0.06     2.90     3.13       3038 1.00
#>
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
#> is a crude measure of effective sample size, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

plot(marginal_effects(fit), points = TRUE)
`````` Created on 2019-07-24 by the reprex package (v0.3.0)

I think I understand better your question and the confusion of intercept in the model output. The population-level effects show an `a_Intercept` and `b_Intercept` those are just the estimates of your coefficients in your non-linear model.

In this case `a` is the “intercept” at the point the data begins and `b` is the curving coefficient. I believe that `brms` labels them as `x_Intercept` because you are essentially fitting two intercept only models on each separate non-linear coefficient (hence `a ~ 1` and `b ~ 1` in the formula statement). I could be wrong on this though.

In the case of applying a random-intercept model, you will get many different values for `a_Intercept` but they will all share the same `b_Intercept` value. This is my vague understanding as I don’t have much experience in multi-level modeling.

Thanks dylanjm it was very very useful for me and I hope also for the other users. quoting your text: " In the case of applying a random-intercept model, you will get many different values for `a_Intercept` but they will all share the same `b_Intercept` value. "
Is this true if my prior is a ~ 1+(1|ID), b ~ 1 whereas is not true if my prior is a ~ 1+(1|ID), b ~ 1+(1|ID) because in the latter case also b is depending on the class ID
Am I right?

I think your understanding is right.

Dear All, thanks a lot.
I have tried (with no success) to extend the toy-model of the example to the case of multi-level model. Can you please have a look to it? Thanks a lot and best regards

``````x1 <- runif(N, 1, 20)
a1 <- 50
b1 <- -1.3
sigma1 <- 3
id1=1
y_true1 <- a1 * x1^b1
y1 <- rnorm(N, y_true1, sigma1)

x2 <- runif(N, 1, 20)
a2 <- 20
b2 <- -0.3
id2=2
y_true2 <- a2 * x2^b2
sigma2=0.6
y2 <- rnorm(N, y_true2, sigma2)

x=c(x1,x2)
y=c(y1,y2)
ID=c(rep(id1,N),rep(id2,N))
dat <- tibble(x, y,ID)

fit <- brm(bf(y ~ a*x^b, nl = TRUE,
a ~ 1+(1|ID), b ~ 1+(1|ID)),
data = dat, family = gaussian(),
prior = c(prior(normal(40, 10), nlpar = "a"),
prior(normal(-2, 1), nlpar = "b")))

summary(fit)

Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ a * x^b
a ~ 1 + (1 | ID)
b ~ 1 + (1 | ID)
Data: dat (Number of observations: 2400)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000

Group-Level Effects:
~ID (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(a_Intercept)    16.94      7.15     8.01    35.22        604 1.01
sd(b_Intercept)     1.70      1.16     0.35     4.65        313 1.03

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
a_Intercept    37.22      6.91    23.93    51.19        502 1.00
b_Intercept    -1.35      0.72    -3.03    -0.12        297 1.02

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     2.15      0.03     2.09     2.21        880 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Warning message:
There were 74 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup``````

What makes your try unsuccesful? The results don’t look too stupid to me at first glance.

Hi Paul, and thanks again for your help.
I was expecting to have something around a=50 and b =-1.3 for the ID1 and a=20 and b=-0.3 for the ID2.
At this point I think is me not able yet to read the model output in the correct way, can you help me in that?

Thanks a lot and best regards

Any new suggestion for this? Thanks a lot

The output from summary just gives the population-level effects on the estimates, which are ~37 for a and ~-1.3 for b. However, what I think you’re looking for is the combined effect of the population-level effects and random effects due to grouping by ID. To see these results, do the following:

``````> coef(fit)
\$`ID`
, , a_Intercept

Estimate Est.Error     Q2.5    Q97.5
1 50.45061 0.5334878 49.40811 51.50074
2 19.88671 0.2604020 19.39130 20.38089

, , b_Intercept

Estimate   Est.Error       Q2.5     Q97.5
1 -1.2847395 0.012212069 -1.3090926 -1.260158
2 -0.2974429 0.006554002 -0.3103473 -0.284697
``````

You can see here the values of a at ~50 and b at ~-1.3 for ID1 and a at ~20 and b at ~-0.3 for ID2. If you’re interested in just the random effects without the population estimates, replace coef with ranef above.