# Multivariate model with forced single (shared) intercept and variance?

Dear all,

this is my first post.

I am trying to fit a (simple?) multivariate model in `brms` with two responses coming from different distributions. Both are intercept-only models with the same grouping structure (‘random effects’, see reprex below), so the model I have estimates two intercepts and two group-level variances.

``````bf(resp1 | trials(resp1_trials) ~ 1 + (1|g)) + binomial() +
bf(resp2 ~ 1 + (1|g)) + poisson()
``````

Since I know (in the context of simulated data) that the underlying intercept and group level variance is the same for both responses, I would like to constrain the model such that it returns a single (‘shared’) intercept and a single (‘shared’) variance (instead of two each). Is it possible in `brms` to write the model formula so that it adheres to that constraint (rather than coding the model directly in `Stan`)?

I know that I could add a fixed (via a constant prior) or estimated correlation term for the group-level intercepts (e.g. here), which gets me closer to what I have in mind, but still estimates two group-level terms.

Any pointer would be highly appreciated.

``````library(brms)
set.seed(123)
g_val <- rnorm(50, mean = 0.3, sd = 1.2)
resp1 <- rbinom(n = 50, size = 100, prob = inv_logit_scaled(g_val))
resp2 <- rpois(n = 50, lambda = exp(g_val))

dat <- data.frame(resp1, resp1_trials = 100, resp2, g = seq_along(g_val))

form1 <- bf(resp1 | trials(resp1_trials) ~ 1 + (1|c|g)) + binomial()
form2 <- bf(resp2 ~ 1 + (1|c|g)) + poisson()
fit <- brm(form1 + form2 + set_rescor(FALSE), data = dat)
fit

par(mfrow = c(1, 2))
p <- ranef(fit)\$g[, , "resp1_Intercept"][, "Estimate"]
plot(g_val, p, xlab = "truth", ylab = "estimated for response 1")
p <- ranef(fit)\$g[, , "resp2_Intercept"][, "Estimate"]
plot(g_val, p, xlab = "truth", ylab = "estimated for response 2")
``````
• Operating System: macOS 12.5
• brms Version: 2.17.0

Hello, I think that nonlinear syntax will be able to achieve this. By tying together expressions that look something like this:

`bf(count_outcome ~ common_intercept, family = poisson(link = 'log'))`
`bf(binary_outcome ~ common_intercept, family = binomial(link = 'logit'))`
`bf(common_intercept ~ 1)`

The concept of a common variability parameter would depend on what response distributions you selected.

Thanks a lot for that pointer!

I have near-zero experience with nonlinear models and none with the corresponding `brms` syntax. So I only managed to partially replicate my original model using a ‘nonlinear’ model (not accounting for the shared varying intercepts), which still estimates two population-level intercepts.

``````priors <- c(prior(normal(0, 1), coef = "",          nlpar = "commonintercept", resp = "resp1"),
prior(normal(0, 1), coef = "",          nlpar = "commonintercept", resp = "resp2"),
prior(normal(0, 1), coef = "Intercept", nlpar = "commonintercept", resp = "resp1"),
prior(normal(0, 1), coef = "Intercept", nlpar = "commonintercept", resp = "resp2"))

f1 <- bf(resp1 | trials(resp1_trials) ~ commonintercept, nl = TRUE, family = binomial()) + lf(commonintercept ~ 1)
f2 <- bf(resp2                        ~ commonintercept, nl = TRUE, family = poisson())  + lf(commonintercept ~ 1)
f1 + f2 + set_rescor(FALSE)
r <- brm(f1 + f2 + set_rescor(FALSE), data = dat, prior = priors)
``````

Trying a to formulate a common intercept fails:

``````f <- bf(resp1 | trials(resp1_trials) ~ commonintercept, nl = TRUE, family = binomial()) +
bf(resp2                        ~ commonintercept, nl = TRUE, family = poisson()) +
lf(commonintercept ~ 1)
# fails with error

f <- bf(resp1 | trials(resp1_trials) ~ commonintercept, nl = TRUE, family = binomial()) +
bf(resp2                        ~ commonintercept, nl = TRUE, family = poisson()) +
lf(commonintercept ~ 1, resp = c("resp1", "resp2"))
# fails with different error
``````

And this brings me back to my starting point:

``````f <- bf(resp1 | trials(resp1_trials) ~ commonintercept, nl = TRUE, family = binomial()) +
bf(resp2                        ~ commonintercept, nl = TRUE, family = poisson()) +
lf(commonintercept ~ 1, resp = "resp1") +
lf(commonintercept ~ 1, resp = "resp2")
``````

I guess this boils down to me not understanding the nonlinear syntax. And I haven’t yet found examples for multivariate nonlinear models using `brms`.

So for now I think I will just code this model up in Stan directly.

Thanks again for your feedback!

That’s a shame, I thought this would work as a nonlinear model, but I was only able to define a nonlinear model for one response at a time. I suppose it would be simple enough to define two brms models and combine the resulting `stancode`.