Specifying multivariate non-linear model with one common parameter in brms

Interfaces brms
,
1

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)

2

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 ys are sampled from some probability distribution whose mean is given by the right-hand side?

3

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.

4

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?

5

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.

6

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