Multinomial Logistic Regression in R

The data on customers' choices can be analyzed to determine which features of a product (e.f., package size, brand, or flavor) are most attractive to customers and how they trade of desirable feature against price. 

Given the customer makes a choice among several options, each of which has its own set of attributes. To accommodate the unique data structure, marketers have adopted choice models, which are well suited to understanding the relationship between the attributes of products and customers' choices among sets of products. 

Choice models are used to understand how product attributes drive customers' choices. The most popular choice model in practice is the multinomial logit model. The model can be estimated using frequentist models with mlogit() or using Bayesian methods with MCMCmnl(). 

Before analyzing any choice data, it is useful to compute raw choice counts for each attribute. 

Estimating a choice model is similar to estimating simpler linear models, The key output of the estimation is a set of parameters that describe how much each attribute is associated with the observed choices. 

Conjoint analysis is a survey method where customers are asked to make choices among products with varying features and prices.

##################################################################
### Multinomial Logistic Regression
##################################################################
## http://www.ats.ucla.edu/stat/r/dae/mlogit.htm
library(foreign)
library(nnet)
ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")

with(ml, table(ses, prog))
with(ml, do.call(rbind, tapply(write, prog, function(x) c(Mean = mean(x), SD = sd(x)))))

ml$prog2 <- relevel(ml$prog, ref = "academic")
with(ml, table(ses, prog2))
with(ml, table(write, prog2))
test <- multinom(prog2 ~ ses + write, data = ml)
summary(test)
z <- summary(test)$coefficients/summary(test)$standard.errors; z

#2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1))*2; p

## extract the coefficients from the model and exponentiate
exp(coef(test))
head(pp <- fitted(test))
dses <- data.frame(ses = c("low", "middle", "high"), write = mean(ml$write))
predict(test, newdata = dses, "probs")