Since you have two columns in your DV this will be a multivariate regression. The attachment information is I hope helpful. Best wishes David Booth

Thanks David Eugene Booth for the recommendation and sorry for my late reply. I'd like to add a random effect into the formula as it makes more sense with glmer without Bayesian. I've also tried using the brms function but it does not apply to multivariate regression. Is there any other package/function that can handle both mixed model and combined DV please?

Do you really need this to be multivariate eg are the two DVs really independent and then 2 separate regression would do. if that wouldn't work then Google search for a multivariate function that meets all of your criteria. David Booth

The MCMCglmm package should do exactly what you are looking for (https://cran.r-project.org/web/packages/MCMCglmm/index.html).

Ran Dai are you sure that multivariate models do not work with brms?? I think this formula should work in principal to get a multivariate model for the two DVs y1 and y2, with a radom intercept effect for "id":

bf(mvbind(y1, y2) ~ x1 + x2 + (1 | id))

Thank you for the corrections Rainer Duesing

It works with your formula.

Yes, David Eugene Booth the 2 variables are not independent unfortunately.

Thank you Jack Henry will try this one too! I'd also recommend this link on multivariate modelling with Bayesian:

https://cran.r-project.org/web/packages/brms/vignettes/brms_multivariate.html

Cheers - Ran

Here, you can find some code example with two response variables in MCMCglmm: https://tomhouslay.files.wordpress.com/2017/02/indivvar_mv_tutorial_mcmcglmm.pdf

Félicie Dhellemmes thank you, I will try the code. The haggis drawing looks interesting as well~!