I want to specify a multivariate (three responses) non-linear model with two parameters: one is common between all responses (a) and one response specific (b_1, b_2, b_3). The model looks like this:

y_1 = a - \log_{10}(1+ \frac{x}{10^{-b_1}})

y_2 = a - \log_{10}(1+ \frac{x}{10^{-b_2}})

y_3 = a - \log_{10}(1+ \frac{x}{10^{-b_3}})

I know that a \sim N(\mu, \sigma),\; b_i \sim N(\mu_i, \sigma_i),\; i = 1,2,3 Therefore, I guess while three b_i should be estimated (per y_i), **a** can be defined as a random effect across the three responses. I wonder what is the best way of specifying this model in brms.

Attached is also an example simulated data.

exampleData.csv (461 Bytes)

Hi @Vahid_Nassiri and welcome! I think the easiest way to set up this model will be to avoid treating it as a multivariate model, and instead treat each `i`

as a different observation, with an associated grouping variable telling us what `i`

the observation belongs to. Then you can use the nonlinear syntax of `brms`

to fit this model with a random effect for the parameter `b`

. Note however that as written (with equals signs) the observation model is singular. Do you mean instead that the `y`

s are sampled from some probability distribution whose mean is given by the right-hand side?

Hi @jsocolar . Thanks for your response. Indeed, we can fit this as a univariate response, and add three indicator variables that shows the three responses. To me, I may not define parameter *b* as a random effect, since it has a different distribution per outcome (and *b_i*’s are independent from each other). I’d rather estimate one *b_i* per response (each corresponds to one of the indicator variables) and each of these *b_i*’s will also have their own prior (normal with a given mean and sd). And for *a*, I’d then define it as a random effect:

a = a_0 + u_i, u_i \sim N(0, \sigma)

But then the problem is that I’m not sure how to tell brms (or even nls) to estimate one b_i per indicator column. If it was a linear model, it could be straightforward, as we had 3 different intercepts. But I’m not sure how to formulate it in the non-linear way.

I hope I could explain my issue.

Oh, I see! I assumed that \mu, \sigma, \mu_i, and \sigma_i were all known a priori (i.e. that these were non-hierarchical priors) rather than to be estimated from data.

Let me see if I understand your model now, where i indexes which of the y’s we are looking at within a row, and j indexes which row we are looking at.

You want:

y_{ij} \sim \mathcal{N}(\mu_{ij}, \sigma)

where

\mu_{ij} = a_j - \mathrm{log}_{10}(1 + \frac{x}{10^{-b_{ij}}})

and

a_j \sim \mathcal{N}(m, s)

b_{1j} \sim \mathcal{N}(m1, s1)

b_{2j} \sim \mathcal{N}(m2, s2)

b_{3j} \sim \mathcal{N}(m3, s3)

Do I have it right now?

Yes, indeed, it’s very well right. And all those parameters are also known a priori (\mu, \sigma, etc.). So we do not need to estimate the hyper-parameters, they are known and fixed.

I actually created those dummy variables and fit the model like this:

```
brm(
bf(y|se(StErr, sigma = TRUE) ~ a- log10(1 + (x/ (10^(-(b1*y1+ b2*y2 + b3*y3))))),
a~1|Group, b1~1,b2~1,b3~1,
nl = TRUE),
data = data2fit, family = gaussian(),
prior = c(
prior(normal(7 , 0.05), nlpar = "b1"),
prior(normal(7, 0.03), nlpar = "b2"),
prior(normal(7, 0.07), nlpar = "b3"),
prior(normal(7, 10), nlpar = "a")
),
chains = 5,
iter = 20000,
control = list(adapt_delta = 0.9)
)
```

so, `y1`

, `y2`

, and `y3`

are three dummy variables. For example, `y1`

is 1 if the response is `y1`

and zero, otherwise.

seems to work! And the outcomes is pretty similar to the case one fits the values per outcome.

1 Like