I would like to study the following question: how a four-level nominal variable (A) depends on a five-level ordinal variable (B). I also would like to include a random term (C), which is a nominal variable. (For example, how choosing one of four products (product1, product2, product3, product4) depends on how often the person sees an advertisement (not seeing, rare, average, often, very often), and knowing that the district where the person lives affect both the product preference and seeing the advertisement.)
So the ideal model would be something like this: A ~ B + (1|C). (or with the example: product ~ frequency + (1|district).
I guess this is something that could be studied with General Linear Mixed Models (GLMM). However, I haven't found any type of GLMM that could deal with a four-level nominal dependent variable. If the random term wouldn't be important, using a simple chi-square test would be good, but the random term is really important.
I wonder would it be a solution to binarize the dependent variable, and to create four models instead of one. The new dependent variable in the first model would be choosing or not product1, in the second model it would be choosing or not product2, and so on. So the new models would look like this:
model1: Did they choose product1 (yes/no) ~ frequency +(1|district)
model2: Did they choose product2 (yes/no) ~ frequency +(1|district)
model3: Did they choose product3 (yes/no) ~ frequency +(1|district)
model4: Did they choose product4 (yes/no) ~ frequency +(1|district)
As a post hoc test, I would use Tukey test.
I'm not sure using four models instead of one is a good solution, but so far this is the best I could come up with. I would be really happy to receive any advice.
(I work with R.)
If I'm understanding you correctly, the implementation in mlogit should be what you are looking for. This package is designed to implement modeling for multinomial outcomes. Check out ?mlogit::mlogit for details and examples (the "mixed" model as it is sometimes referred to), but your specification would be:
m <- mlogit(A ~ B | district, data=mydata)