Rsidual correlation in multivariate non gaussian model

Hello everyone,

I am currently working on a project where I need to to model two outcomes simultaneously, one of which follows a zero-inflated negative binomial distribution and the second a gamma distribution.

As I am used to work with brms, I can model these two responses as follows.

bf_y1 <- bf(y1~ x + (1|p|id), zi ~ x) + zero_inflated_negbinomial()
bf_y2 <- bf(y2~ x +  (1|p|id)) + gaussian()

mod <-  brm(
  bf_y1 + bf_y2 + set_rescor(F),
  data = echant)

The problem is that, as specified in the code above, I can’t estimate the residual correlation between the two variables (set_rescor(F)). Currently brms only allows this residual correlation to be estimated in the context of a multivariate model where the dependent variables follow a Gaussian or Student distribution.

My question is therefore whether it would be possible to estimate this type of dependence between the two outcomes while modelling them using an appropriate distribution (i.e. ZINB and Gamma). Perhaps this would be possible by coding the model directly in STAN?

Thank you in advance for your help on this subject!


I think this can be done using Gaussian copulas? It seems like @andrjohns was working on getting these implemented in brms a few years ago.