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

Interfaces brms
,
1

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
2

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.

3

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!

4

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.