Reverse engineering values for mlogit logsum function applied to a choice model

Question

I have a dataset of observed choices (N = 26 alternatives) and a corresponding range of predictor variables. The data have been split into training and test subsets, Data_Train and Data_Test.

I have used Data_Train to estimate an mlogit multinomial logit model, M1, to quantify the relative preferences for each alternative at each choice occasion.

My overall objective is to now use M1 to calculate the welfare loss predicted by removing one specific alternative, A, from Data_Test and assuming that whenever alternative A was chosen, the chooser instead took the second-best alternative.

My intention was to use mlogit::logsum() to calculate the expected utility for each alternative at each choice occasion, and then sum the vectors of values across the two versions of Data_Test, one with alternative A and one without; I’ll call the latter Data_Test[-A]. My expectation is that the sum of the logsum vector will be higher for Data_Test than for Data_Test[-A], and the difference between the two will be the welfare loss.

In order to compute the logsum values for Data_Test[-A], I think I need to compute each alternative-specific probability at each choice occasion then manually compute the logsum value for the second-best alternative, so I’ve been trying to reverse engineer how logsum works.

From the mlogit::logsum R-help page, “The inclusive value, or inclusive utility, or log-sum is the log of the denominator of the probabilities of the multinomial logit model.” I interpreted this to mean that logsum takes the maximum probability (P_max) from the 26 calculated at each choice occasion, expresses P_max as a fraction, then takes the log of the denominator. But doing so does not match the value returned by logsum so I’m obviously doing something (or things) wrongly.

Here are the probabilities for each alternative from the first choice occasion (row) as generated by the following line:

Probs <- fitted(M1, newdata = Data_Test.LongFormat, type = "probabilities")

print(Probs)

9.324656e-07 2.804257e-14 1.451469e-04 7.924634e-05 1.099780e-05 1.142550e-02 1.089842e-03 2.157007e-05 3.849436e-03 9.678074e-01 4.846874e-03 3.258943e-05 7.138557e-04 9.775038e-03 1.842690e-04 5.789495e-06 5.373004e-35 1.105493e-05 3.140959e-07 7.864120e-08 3.632554e-09 6.768526e-08 4.103485e-09 4.762275e-17 8.694194e-25 6.025488e-14

P_max here is Probs[10], namely 0.9678074. So to convert to a fraction and get the denominator I use 1/P_max = 1.033263, then take the log of that, giving me 0.03272219. However, when I apply logsum to the same dataset, the result -0.1498398 is returned for this row of data. Indeed, because the logsum result here is negative, the input value must be <1 so the way I am calculating the denominator must be wrong because 1 divided by any probability will be >=1.

Can anyone tell me what I’m doing wrong here, or to suggest an alternative method for calculating welfare lost from removing one alternative from a choice set? Thanks in advance for your help. It will probably be fairly clear that I’m not an economist so apologies for mangled terminology and concepts.

Answer 1

I will take your question point by point.

I have a dataset of observed choices (N = 26 alternatives) and a corresponding range of predictor variables. The data have been split into training and test subsets, Data_Train and Data_Test.

I understand that you have observations of choice situations (we don't know how many) and 26 alternatives in total (which is a lot) presented through choice situations.

The data have been split into training and test subsets, Data_Train and Data_Test. I have used Data_Train to estimate an mlogit multinomial logit model, M1, to quantify the relative preferences for each alternative at each choice occasion.

You have separated your data into two sets (we don't really know why, in order to predict the choices in the second set from the estimates in the first?).

My overall objective is to now use M1 to calculate the welfare loss predicted by removing one specific alternative, A, from Data_Test and assuming that whenever alternative A was chosen, the chooser instead took the second-best alternative.

It seems that two objectives are described here. The first is to predict the probability of choices when an alternative is removed, based on the estimation of your M1 model. This can be done with the following code:

#This 2 lines will give you the probability of choices predicted by your model
fit <- apply(fitted(M1, outcome = FALSE), 2, mean)
fit

#Here we will remove one alternative and predict the probability of choices accross the remaining alternatives
X <- model.matrix(M1)
X <- X[-which(X[alt == "A",]),] #We remove one alternative from the model matrix
chid <- idx(M1, 1)[which(idx(M1)$alt != "A")] #We extract alternative's indexes without the alternative removed
eX <- as.numeric(exp(X %*% coef(M1))) #We compute utility for each alternatives
SeX <- as.numeric(tapply(eX, sort(chid), sum)) #We sum over the choice situations
P <- eX / SeX[sort(chid)] #We compute the probability of choices 
P <- matrix(P, ncol = 25, byrow = TRUE) 
P
MS <- rbind(c(1:25), apply(P, 2, mean)) #Here you have the mean
MS
MS <- t(MS)
sum(as.numeric(MS[,2])) #Here the sum should be 1

The second objective is to compare the welfare obtained from your initial model with the welfare predicted without alternative A. This can be done with the following code:

#For the welfare of your initial model
W_M1 <- logsum(M1)
summary(W_M1)

#For the welfare of the predicted model without one alternative
W_wo_A <- log(SeX)
summary(W_wo_A)

#If we want to compare
summary(W_wo_A - W_M1)

#If we want to compute the surplus and the loss of surplus (in $)
summary(-W_M1/coef(M1)["Price"])
summary(-W_wo_A/coef(M1)["Price"])
summary(-(W_wo_A - W_M1)/coef(M1)["Price"])

I hope it will help (for next time try to add your data so we can better adapt the code).