Bivariate/multivariate models for multinomial response variables I need to fit two categorical (potentially correlated) response variables (each has three classes) on a set of explanatory variables, while considering for the response variables' correlation. What type of models do you think I should better use (ideally implementable in R)? I have seen literature on copula-based joint models, but they seem to be mostly on the theoretic side, rather than being implementable. Also, Bayesian multivariate logistic regression models could be an option if they were developed for the multinomial case, but apparently they are not.
 A: The most obvious way to predict two correlated categorical variables is to merge them into a single variable! Instead of predicting classes $a_1$, $a_2$, $a_3$ for the first target and $b_1$, $b_2$, $b_3$ for the second target, you can predict 9 classes like $a_1b_1$, $a_1b_2$, $a_1b_3$, $a_2b_1$ etc. for the joint target. 
This structure is flexible enough to allow for any kind of correlation between the response variables. And because the resulting problem is an ordinary multiclass single-response classificaton, you can plug in any model you want, from multivariate logistic regression to ensemble of decision trees. 
A: What you are dealing with is multinomial response data, which is usually modelled using some kind of 
multinomial response model (e.g., multinomial logit regression, etc.).  You have two categorical output variables with three possible outcomes, which gives you $3 \times 3 = 9$ possible outputs in total.  As David Dale has suggested in his answer, you should combine your two outputs into a single categorical variable with these nine possible outcomes.  Modelling this categorical variable using a multinomial response model will give estimated probabilities for each of these nine outputs, which also gives you implicit estimation of the marginal probabilities and correlations for the individual output variables.
Since multinomial outputs from an exchangeable sequence of observations are, by definition, perfectly fit by the multinomial distribution, the only real question in fitting a multinomial response model is whether the particular regression and link functions you use give reasonable estimates of the probabilities.  As with any data-modelling exercise, we cannot specify the best model a priori - this must be assessed by the fit of the data to different models, as judged by diagnostic plots, etc.  However, a reasonable starting point would be a multinomial logit model (classical or Bayesian).  I would suggest you try fitting this model to your data and see how the diagnostic plots look, and whether your link seems reasonable.  Assuming this model is not falsified by the diagnostic plots, it will give you estimated probabilities of each possible outcome (out of nine outcomes), from which you can calculate marginal and conditional probabilities and correlations for the two outputs.
