Using delta method to get confidence intervals for multinomial logit?

Question

I am working with some choice modeling data and am interested in trying to potentially use the delta method with the multinomial logit model that I'm analyzing the data with. Here's an example:

First, I uploaded some data in R:

# First upload and format the data
cbc.df <- read.csv("http://goo.gl/5xQObB",
                   colClasses = c(seat = "factor", price = "factor",
                                  choice="integer"))
cbc.df$eng <- factor(cbc.df$eng, levels=c("gas", "hyb", "elec"))
cbc.df$carpool <- factor(cbc.df$carpool, levels=c("yes", "no"))
summary(cbc.df)

Format data for mlogit package depending on version:

if (packageVersion("mlogit") < "1.1") {
  # for mlogit up through 1.0.3
  cbc.mlogit <- mlogit.data(data=cbc.df, choice="choice", shape="long", 
                            varying=3:6, alt.levels=paste("pos", 1:3), 
                            id.var="resp.id")
} else {
  # for mlogit starting with version 1.1
  library(dfidx)      # install if needed
  
  # add a column with unique question numbers, as needed in mlogit 1.1+
  cbc.df$chid <- rep(1:(nrow(cbc.df)/3), each=3)
  
  # shape the data for mlogit
  cbc.mlogit <- dfidx(cbc.df, choice="choice", 
                      idx=list(c("chid", "resp.id"), "alt" ))
}

I think fit a basic multinomial logit model:

# fit the models
library(mlogit)
m1 <- mlogit(choice ~ 0 + seat + cargo + eng + price, data = cbc.mlogit)

Next, I want to predict the probability of making a choice given different combinations of the predictors. I found this function to do that:

# Predicting shares
predict.mnl <- function(model, data) {
  # Function for predicting shares from a multinomial logit model 
  # model: mlogit object returned by mlogit()
  # data: a data frame containing the set of designs for which you want to 
  #       predict shares.  Same format at the data used to estimate model. 
  data.model <- model.matrix(update(model$formula, 0 ~ .), data = data)[ , -1]
  utility <- data.model%*%model$coef
  share <- exp(utility)/sum(exp(utility))
  cbind(share, data)
}

So then I generated some combinations of predictors to generate predictions for:

> # and set up attributes needed later
> attrib <- list(seat  = c("6", "7", "8"),
+                cargo = c("2ft", "3ft"),
+                eng   = c("gas", "hyb", "elec"), 
+                price = c("30", "35", "40"))
> 
> 
> 
> (new.data <- expand.grid(attrib)[c(8, 1, 3, 41, 49, 26), ]) # find attrib at top
   seat cargo  eng price
8     7   2ft  hyb    30
1     6   2ft  gas    30
3     8   2ft  gas    30
41    7   3ft  gas    40
49    6   2ft elec    40
26    7   2ft  hyb    35

And finally, get the predictions for these combinations:

> predict.mnl(m1, new.data)
        share seat cargo  eng price
8  0.11273356    7   2ft  hyb    30
1  0.43336911    6   2ft  gas    30
3  0.31917819    8   2ft  gas    30
41 0.07281396    7   3ft  gas    40
49 0.01669280    6   2ft elec    40
26 0.04521237    7   2ft  hyb    35

This share variable is great, but what I'm really wanting is something like a 95% CI for it and I've heard that the delta method could potentially help. I just think that knowing that the first combination has an 11% share is meaningless without knowing the variability in the estimate.

Can anyone help?

Answer 1

Sort of not answering your question, but the brms R package covers the multinomial model. E.g. like this brm(bf(y | trials(size) ~ 1, muy2 ~ x), data = mydata, family = multinomial(), prior = myprior) (see e.g. here). The great thing - assuming you are happy with a Bayesian approach - is that you can then get predictions based on each MCMC sample.

One of the reasons why this gets a lot easier in a Bayesian context is that if you have MCMC samples for some parameters, then to get MCMC samples for some transformation (including adding, subtracting, multiplying etc. the parameters) just requires you to calculate the transformation for each MCMC sample. Then, voilà, you've got pseudo-random samples for the transformation. If you then want to predict for new data, you just draw pseudo-random samples from the outcome distribution for each of the samples from the combination of parameters.

In short, what I'm suggesting is that reflecting uncertainties gets a lot easier, if you take a Bayesian approach. You don't need to derive any results about delta-methods etc., but of course this comes with the downside of people potentially criticizing your prior choices and MCMC sampling taking a bit longer than fitting a model using maximum likelihood.