I have a dataset of around 30,000 people where each chooses one of 4 items: A, B, C or D. People are nested within the 600 areas. I want to fit a multinomial logit model, where random intercepts are calculated for each of the 600 areas for each option. i.e. we allow that people in some areas may have a greater tendency to pick A. I also want to include area level predictors and possibly individual level predictors.
What software would be suitable for fitting such a model? The mlogit package in R does not appear to allow alternative specific random intercepts and various other packages I've tried in Stata and R do not seem to have exactly this set of options. I'm open to trying out any software if it will do the job.