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.
2 Answers
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.
-
$\begingroup$ How does this structure allow for any kind of correlation between the response variables? Isn't IIA (Independence of Irrelevant Alternatives) the main flaw of multinomial logit models (and it seems to be likely to exist in my case, as the classes are all manually created and are clearly correlated)? $\endgroup$– FredCommented Jun 3, 2018 at 20:14
-
$\begingroup$ When you speak about correlation between categorical variables, I assume that you mean that your outputs $y_1$ and $y_2$ are dependent conditionally on the inputs $x$. By modeling their conditional distribution $P(y_1, y_2 | x)$ as a complete $3\times3$ matrix, you can capture any dependence between them. $\endgroup$ Commented Jun 4, 2018 at 6:44
-
$\begingroup$ What I mean is what if there is an actual correlation among classes of the response variable? For example, assume a response variable "weather" with the classes of "clear, rainy, humid". Clearly the two latter classes are more probable to be dependent on more similar conditions, and consequently may be correlated with a higher degree than with the first class. This seems to be the case of my 9 classes (since for example a1b1 and a1b2 are probably more correlated than a1b1 and a2b2). How do I address this? $\endgroup$– FredCommented Jun 4, 2018 at 16:57
-
1$\begingroup$ If these common conditions are reflected in $X$, then slope coefficients for "rainy" and "humid" classes will be similar. If these common conditions are not fully reflected in $X$ (that is, unobservable), then no model will be able to discriminate between this classes, and both $P(rainy|X)$ and $P(humid|X)$ will be high. Any model which is flexible enough can account for such a dependence. $\endgroup$ Commented Jun 4, 2018 at 18:13
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.
-
$\begingroup$ Thank you for your answer. How is this 9-class response variable approach better than estimating two 3-level response variable models? How can the inference between these two approaches be compared? For example, how can the differences between the estimated effects of one of the X's on Y, vs. these effects on Y1 and Y2 be interpreted? $\endgroup$– FredCommented Jun 3, 2018 at 20:11
-
$\begingroup$ If you estimate two 3-level models separately, you don't get any interaction estimates. This means that you only estimate the marginal probabilities of the two outcomes, but you can't estimate the joint probabilities (since you don't know how they're correlated). But if you combine them into a single 9-level model, you are estimating the probabilities of all nine possible outcomes. That gives you estimates of the joint probabilities of the two 3-level outcomes, which allows you to also get marginal probabilities, correlation, etc. $\endgroup$– BenCommented Jun 3, 2018 at 22:40
-
$\begingroup$ From a 9-level model, you will get an estimate that can be written as a $3 \times 3$ table of probabilities (that add up to one). You can then calculate marginal probabilities as row and column sums. Then use these to calculate variances, covariances, correlation. (If you want to know how to do this, please ask a new question: "How do I calculate moments from a $3 \times 3$ table of probabilties?") $\endgroup$– BenCommented Jun 3, 2018 at 22:43
-
$\begingroup$ Is there a statistical approach to test whether it is preferred to model the joint probability or the marginal probabilities separately? I assume if we find that the two response variables are actually not correlated, there is no reason for considering a joint probability, right (my main objective of this modeling is inference, specifically interpretation of the odds ratios, to see how the X's are affecting the Y's)? $\endgroup$– FredCommented Jun 4, 2018 at 16:50
-
1$\begingroup$ You can compare prediction accuracy for both approaches and choose the best one. For some classes of models you can also perform a likelihood ratio test, but you need to make sure they are nested. $\endgroup$ Commented Jun 4, 2018 at 18:12