Can I set up multiple probability models that sum to 1? I am interested in consumer research.  For simplicity's sake, assume I want to know if people prefer Burger King or McDonald's.  My questionnaire will contain a question like "If you had to choose one restaurant, which one would you choose: Burger King, McDonald's, or unsure?"  Upon getting my results back, I want to fit logistic regression models to extrapolate my results from my sample to my population of interest.  
Suppose my first model is a "Burger King support" model.  I code all "Burger King" answers as 1 and all other answers as 0, then regress those data on my demographic variables to come up with a model.
I then construct a "McDonald's support" model.  I code all "McDonald" answers as 1 and all other answers as 0, then regress those data on my demographic variables to come up with a second model.
I can additionally construct an "Unsure" model where all "unsure" answers are coded as 1, other answers are coded as 0, and those data are regressed on my demographic variables.
So, my models should look like this:
$P(BK=1|geography,age) = logit(\beta_0+\beta_1*geography+\beta_2*age)$
$P(McD=1|geography,age) = logit(\beta_3+\beta_4*geography+\beta_5*age)$
$P(Unsure=1|geography,age) = logit(\beta_6+\beta_7*geography+\beta_8*age)$
(I'm having a bit of trouble importing a table in here to show a sample of what the data frame looks like; let me know if I can better clarify what I wrote above.  Ultimately, we will have an ID column, a column with the survey responses recorded or "NA" if they did not respond to the survey, a column of 1s and 0s for "BK," a column of 1s and 0s for "McD," a column of 1s and 0s for "Unsure," a column with their geographic indicator (categorical) and a column with their age.)
Upon running these models, I end up getting three scores that need not (and usually don't) sum to one.  For example, one individual's "Burger King score" (BK) might be 0.6, her "McDonald's score" (McD) might be 0.5, and her "unsure score" (U) might be 0.2.  Obviously, this sums to 1.3 and does not make probabilistic sense.
My goal is to be able to estimate the probability that one prefers McDonald's, that one prefers Burger King, and that one is unsure about what he/she prefers.
At this point, I have a few questions:


*

*Is this a valid method of looking at things?  If not, what are more appropriate methods for this kind of analysis?  (The framework I have built up is geared toward logistic regression so there is a strong preference for using logistic regression where possible.)  My goal is to not simply fit two of the three models and allocate the rest of the probability to the third option.  In the example above, obviously this would not be valid as the BK and McD models sum to 1.1.

*In looking at plots and building models in the "statistically proper" fashion, some covariates may be included in one or two models, but not all three.  Is this problematic when attempting to interpret these models?

*Is there a method of scaling such that we can get the probabilities to sum to 1?  For example, in the above example, an individual had a BK score of 0.6, a McD score of 0.5, and U score of 0.2.  Is it appropriate to divide each score by 1.3 such that they all sum to 1?  I don't want to make this "linear constraint," however, if there is a more appropriate way of scaling.
I have a reasonably solid background in statistics (Master's degree) so feel free to be heavy-handed with theory if need be, though practical solutions are always preferable!  Thanks so much for your help, and please let me know where I can clarify or elucidate my ideas!
 A: As others have mentioned in the comments, multinomial logistic regression (MLR) would be my first recommendation for what you are trying to do. It produces estimated probabilities that sum up to 1 across the categories for each person, and these probabilities can be interpreted as the probability of choosing category X given the covariates and the alternatives listed. MLR also gives you covariate information about how each covariates affects the outcome, but the interpretation of the coefficients of the covariates, as mentioned, is relative to the baseline category.
An alternative that you can consider for estimating these probabilities is using random forests. It is an ensemble of classification trees that works in the following way:


*

*Grows say, 1000 different classification trees based on 1000 bootstrap samples of your dataset. 

*For each person, use the 1000 classfication trees to classify the person into one of the categories you have. 

*For category X, the probability of being in that category is: # trees that classified the person into X / 1000


Random forests can also tell you which covariates are important, however as far as I know they do not tell you HOW the covariates affect the outcome.
Detailed information can be found here: https://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm
R package for implementing a simple random forest: https://cran.r-project.org/web/packages/randomForest/index.html
