The logit choice model

set.seed(13)
# Define the softmax function
softmax <- function(x) exp(x)/sum(exp(x))
P <- 10 # Number of alternative choices
# Create a random characteristic matrix and weights matrix X and beta
X <- matrix(rnorm(P*5), P, 5)
beta <- rnorm(5)
# Draw some random demand shocks for each product
xi <- rnorm(P)
# Generate the fixed utilities
mu <- X%*%beta + xi

# illustrate the probabilities
probs <- softmax(mu)
data.frame(probs_percent = round(probs*100, 2), cumulative_probability = cumsum(probs))

##    probs_percent cumulative_probability
## 1           0.98            0.009779541
## 2           0.75            0.017297128
## 3           4.55            0.062772591
## 4           8.82            0.151005651
## 5           4.56            0.196561651
## 6           0.84            0.204980945
## 7          23.44            0.439378311
## 8           4.02            0.479612967
## 9           2.33            0.502867434
## 10         49.71            1.000000000

round(softmax(mu + 10)*100, 2)

##        [,1]
##  [1,]  0.98
##  [2,]  0.75
##  [3,]  4.55
##  [4,]  8.82
##  [5,]  4.56
##  [6,]  0.84
##  [7,] 23.44
##  [8,]  4.02
##  [9,]  2.33
## [10,] 49.71

Microdata model: conjoint studies and observational choice data

The microdata version of the logit model is what we fit to observations of individual choices. In this type of data, we observe what each person chooses as well as what they do not choose, that is, we see their full choice set as well as their choice. This data type is most commonly exploited in quantitative marketing, where such studies are known as “conjoint analysis”. In conjoint studies, participants are repeatedly shown sets of alternative choices, and asked to express their choice, either by chosing their favored chocie, or by making one choice and selecting least favored choice, or by indicating quantities of items—“allocating chips”. Such studies tend to be very high-powered, and can help to identify complex substitution patterns. They also have the advantage that the researcher can experimentally vary the choice set, and the characteristics of each choice. And so they might overcome some of the endogeneity concerns common in observational studies. Yet they do not actually observe choices; they only measure stated preferences.

Microdata models can also be applied in an observational setting, say, by using scanner or sales data. This form of data has the advantage of being a measurement of revealed preferences rather than stated preferences. Yet it can be harder to observe the complete choice set available to customers. Also, in most observational settings, the choice set and choice attributes normally aren’t varied experimentally, meaning we have to contend with the endogeneity of these attributes of the data. Later in this chapter we discuss approaches to doing so.

The microdata implementation of the Logit choice model is one that could simulate individual choices, which will occur with the probabilities described above. The only thing we need is a sampling distribution for the choices. We have this in the categorical distribution, which is a multinomial distribution—a distribution over unordered qualitative variables—of size 1. You can think of a categorical distribution as being like rolling a weighted dice.

Let’s continute to use the simulated data we specified above, though we’ll add \(\mu_{1} = 0\) to signify the fixed utility of the outside choice being zero. We’ll then simulate the choices of 100 customers who all face precisely the same choices. Once we’ve done that, we’ll derive the likelihood and write the model out in Stan. Finally we’ll see what happens when each customer faces a different set of choices.

First we simulate and plot some choices for 100 customers:

# Variable sizes
N_customers <- 100
N_choices <- (P+1) * N_customers
# Probabilities including outside good
probs <- softmax(c(0, mu))
simulated_data <- data_frame(choice = rep(0, N_choices),
                             customer_id = rep(1:N_customers, each = P+1),
                             product = rep(1:(P+1), N_customers))
# Draw a choice and indicate this with choice = 1; all unchosen 
# products remain with choice = 0
for(i in 1:N_customers) {
  choice <- sample(1:(P+1), size = 1, prob = probs)
  simulated_data[simulated_data$customer_id == i & simulated_data$product==choice, 1] <- 1
}

# Plot the choices vs the generative probabilities
simulated_data %>%
  group_by(product) %>% 
  summarise(`Proportion of choices` = sum(choice)/N_customers) %>%
  mutate(`Underlying probabilities` = probs) %>%
  ggplot(aes(x = factor(product,ordered = T), y = `Proportion of choices`)) +
  geom_bar(stat = "identity") +
  geom_point(aes(y = `Underlying probabilities`), colour = "red") +
  annotate("text", x = 4, y = 0.15, label = "True preference shares", colour = "red") +
  annotate("text", x = 9, y = 0.35, label = "Market shares", alpha = 0.7) +
  theme_hc() + labs(x = "Product", title = "With 100 customers")

What do we notice here? First, the more popular products (those with higher values for \(\mu_{j}\)) indeed get more choices. But the variation around the probabilities is quite high. How does this change with the number of decisionmakers?

The lesson here is that the variation of the observed choice shares around the true preference shares depends strongly on the sample size (or market size, if we’re modeling markets). This shouldn’t surprise you much—more observations drive more precise estimates. And so when you are modeling markets with few sales, or few customers, you should always use a model that makes use of this information. Of the three Logit implementations described here, this includes the microdata model and sales model—not the shares model.

Constructing the likelihood of the microdata model

Recall from the likelihood post what a likelihood is in a practical sense: a way scoring whether some model parameters \(\theta\) make sense of the observed data under the assumption that the data were generated by our model. If our choice data is a sequence of choices that were made by decisionmakers, the likelihood is simply the product of the probabilities that our model would give to each of those choices being made (for a given set of parameter values). More formally, let \(y\) be the vector of all customer choices (So it might look like \(y = (\text{Good 12, Good 4, Outside good}, \, \dots)'\)). If we denote the choice probability for individual \(i\)’s choice under the model parameterized by some values \(\theta\) as \(p_{y_i}(\theta)\), then the likelihood is

\[ p(y |\, \theta) = \prod_{i=1}^{N}p_{y_i}(\theta) \] Note that each of these probabilities is between 0 and 1 by construction, so their product will be a very small number, and so will invite computational issues. This motivates taking the log likelihood

\[ \log\left(p(y |\, \theta)\right) = \sum_{i=1}^{N}\log\left(p_{y_i}(\theta)\right) \]

To construct the log likelihood, all we need is to sum up the log probabilities (under the model evaluated with a given set of parameters \(\theta\)) of each choice that was made. One straightforward way of doing that is to encode \(y\) as a binary vector which is 1 if a customer made a choice, and 0 otherwise. Note that it each decision-maker’s choice set and choice is “stacked” on top of the other. Then we can construct a vector of choice probabilities \(\pi\) for our dataset—both for the choices that were made and those that were not. The log likelihood is then

\[ \log\left(p(y |\, \theta)\right) = y' \log(\pi) \] Of course this can be evaluated only where \(\pi \neq 0\). To make this more concrete, it might help to look at the data constructed this way

Individual	Possible choice	y	pi_theta	price	Size = Large
1	Outside option	1	\(\pi_{\text{Outside option}}(\theta)\)	0	0
1	Choice 1	0	\(\pi_{\text{Choice 1}}(\theta)\)	2.50	0
1	Choice 2	0	\(\pi_{\text{Choice 2}}(\theta)\)	3.99	1
2	Outside option	0	\(\pi_{\text{Outside option}}(\theta)\)	0	0
2	Choice 1	0	\(\pi_{\text{Choice 1}}(\theta)\)	2.50	0
2	Choice 2	1	\(\pi_{\text{Choice 2}}(\theta)\)	3.99	1
3	Outside option	0	\(\pi_{\text{Outside option}}(\theta)\)	0	0
3	Choice 1	0	\(\pi_{\text{Choice 1}}(\theta)\)	2.50	0
3	Choice 2	1	\(\pi_{\text{Choice 2}}(\theta)\)	3.99	1
…	…	…	…	…	…

The log likelihood of this model is the sum of the product of columns y and the log of pi_theta. Or in linear algebra

\[ \log\left(p(y |\, \theta)\right) = y' \log(\pi) = 1 \times \log(\pi_{\text{Outside option}}(\theta)) + 1 \times \log(\pi_{\text{Choice 2}}(\theta)) + 1 \times \log(\pi_{\text{Choice 2}}(\theta)) \]

In Stan, the log likelihood is added to the log posterior density accumulator (the target) like so:

target += y ' * log(choice_probabilities)

All we need to do is show Stan how to construct the choice_probabilities vector.

Calculating the choice probabilities

We just saw how to construct the log likelihood from a vector of choices and a set of choice probabilities. And we know that the choice probabilities are just the softmax of the fixed component of utility. If all decision-makers faced precisely the same choice-set, then we could construct a common choice probability vector

\[ \pi_{i} = \pi = \text{softmax}(X\beta + \xi) \]

in Stan, like

model {
  vector[N] choice_probability = rep_vector(softmax(X*beta + xi), I);
  
  // priors on beta and xi here
  //... 
  
  // evaluate log likelihood
  target += y' * log(choice_probability);
}

Here N is the number of observations in y and I is the number of decision-makers \(X\) has as many rows as there are choices—including the outside option—as does xi. Both the first entry in xi and X’s first row are all 0s, so that X[1]*beta + xi[1] is equal to zero—the fixed utility from the outside option. Using the fake data illustrated above, \(X\) would be

Intercept column	price	Size = Large
0	0	0
1	2.50	0
1	3.99	1

Note that we include an intercept column here. In cases where the choice characteristics are multicolinear with 1 (such as brand dummies), you should not use an intercept column. Also note that the intercept entries for the outside good are equal to 0, for identifiability reasons discussed above.

The only tricky thing to do is deal with the possibility that each individual might face different choice sets. We accomplish this using indexes that indicate where in our data each choice set starts and finishes, and which choice is being evaluated in each row, then constructing the choice probabilities for each individual in a loop across choice sets.

To illustrate, say that instead our choice data is

Individual	Possible choice	y	price	Size = Large
1	Outside option	1	0	0
1	Choice 1	0	2.50	0
1	Choice 2	0	3.99	1
2	Outside option	0	0	0
2	Choice 1	0	2.50	0
2	Choice 2	1	3.99	1
3	Outside option	0	0	0
3	Choice 2	0	3.99	1
3	Choice 3	1	1.50	1
3	Choice 4	0	1.99	0
…	…	…	…	…

Now each decision-maker faces different choices, or even different numbers of choices. We can no longer describe an \(X\) matrix that applies to all individuals, so instead the \(X\) matrix contains the attributes for all individuals, stacked on top of each other, like so

Intercept column	price	Size = Large
0	0	0
1	2.50	0
1	3.99	1
0	0	0
1	2.50	0
1	3.99	1
0	0	0
1	3.99	1
1	1.50	1
1	1.99	0
…	…	…

Recall that \(\mu = X\beta + \xi\). It doesn’t matter how many rows \(X\) has, each row of \(X\beta\) will be the same. But we’d ideally want to construct mu in Stan like

vector[N] mu = X * beta + xi;

yet xi has only as many elements as there are choices, while X*beta has as many rows as there are possible choices made. The trick is to use an indexing vector choice_index which tells us which possible choice is being considered in each row. Then we can define

vector[N] mu = X * beta + xi[choice_index];

This gives us the fixed component of utility for each decision-maker, but we still need to convert it into probabilities. Note that you must not just take the softmax of mu; we need to take the softmax of mu within each choice set. If each individual has only one choice set, then we can loop across these with indexing vectors that tell us on which row each choice set starts and ends. Let’s call these indexes start and end. Our model block is therefore

model {
  vector[N] mu = X * beta + xi[choice_index];
  vector[N] choice_probability;
  
  // fill in the choice probabilities for each choice set
  for(i in 1:I) {
    choice_probability[start[i]:end[i]] = softmax(mu[start[i]:end[i]]);
  }
  
  // priors go here
  // ...
  
  // Evaluate the log likelihood
  target += y' * log(choice_probabilities);
}

So how best to build up the indexes start, end and choice_index? There are many ways, but we find that this is easy:

library(tidyverse)
# Note that we name "Outside option" "1. Outside option"
my_data <- data_frame(Individual = c(rep(1:3, each = 3), 3), 
           `Possible choice` = c(rep(c("1. Outside option", "Choice 1", "Choice 2"), 2), 
                                 c("1. Outside option", "Choice 2", "Choice 3", "Choice 4")),
           y = c(1, 0, 0, 0, 0, 1, 0, 0, 1, 0),
           price = c(rep(c(0, "2.50", "3.99"), 2), 0, "3.99", "1.50", "1.99"),
           `Size = Large` = c(0, 0, 1, 0, 0, 1, 0, 1, 1, 0))

# Ordered will order the factor alphabetically, and we want the outside option to have index = 1
choice_index <- as.numeric(ordered(my_data$`Possible choice`))

# Now construct start and end
tmp <- my_data %>% 
  mutate(row = 1:n()) %>% 
  group_by(Individual) %>%
  summarise(start = first(row),
            end = last(row))
  
start <- tmp$start
end <- tmp$end