Estimating Mixed Logit Models

This vignette demonstrates an example of how to use the logitr() function to estimate mixed logit (MXL) models with preference space and WTP space utility parameterizations.

The data

This example uses the yogurt data set from Jain et al. (1994). The data set contains 2,412 choice observations from a series of yogurt purchases by a panel of 100 households in Springfield, Missouri, over a roughly two-year period. The data were collected by optical scanners and contain information about the price, brand, and a “feature” variable, which identifies whether a newspaper advertisement was shown to the customer. There are four brands of yogurt: Yoplait, Dannon, Weight Watchers, and Hiland, with market shares of 34%, 40%, 23% and 3%, respectively.

In the utility models described below, the data variables are represented as follows:

Symbol Variable
\(p\) The price in US dollars.
\(x_{j}^{\mathrm{Feat}}\) Dummy variable for whether the newspaper advertisement was shown to the customer.
\(x_{j}^{\mathrm{Hiland}}\) Dummy variable for the “Highland” brand.
\(x_{j}^{\mathrm{Weight}}\) Dummy variable for the “Weight Watchers” brand.
\(x_{j}^{\mathrm{Yoplait}}\) Dummy variable for the “Yoplait” brand.

Preference space model

This example will estimate the following mixed logit model in the preference space:

\[\begin{equation} u_{j} = \alpha p_{j} + \beta_1 x_{j}^{\mathrm{Feat}} + \beta_2 x_{j}^{\mathrm{Hiland}} + \beta_3 x_{j}^{\mathrm{Weight}} + \beta_4 x_{j}^{\mathrm{Yoplait}} + \varepsilon_{j} \label{eq:mnlPrefExample} \end{equation}\]

where the parameters \(\alpha\), \(\beta_1\), \(\beta_2\), \(\beta_3\), and \(\beta_4\) have units of utility. In the example below, we will model \(\beta_1\), \(\beta_2\), \(\beta_3\), and \(\beta_4\) as normally distributed across the population. As a result, the model will estimate two parameters for each of these coefficients: a mean (par_mu) and a standard deviation (par_sigma).

Estimate the model using the logitr() function:

library("logitr")

mxl_pref <- logitr(
  data     = yogurt,
  choice   = 'choice',
  obsID    = 'obsID',
  pars     = c('price', 'feat', 'brand'),
  randPars = c(feat = 'n', brand = 'n'),
  # You should run a multistart for MXL models since they are non-convex
  numMultiStarts = 10
)
#> Running multistart...
#>   Iterations: 10
#>   Cores: 3
#> Done!

Print a summary of the results:

summary(mxl_pref)
#> =================================================
#> Call:
#> logitr(data = yogurt, outcome = "choice", obsID = "obsID", pars = c("price", 
#>     "feat", "brand"), randPars = c(feat = "n", brand = "n"), 
#>     numMultiStarts = 10)
#> 
#> Frequencies of alternatives:
#>        1        2        3        4 
#> 0.402156 0.029436 0.229270 0.339138 
#> 
#> Summary Of Multistart Runs:
#>    Log Likelihood Iterations Exit Status
#> 1       -2641.243         38           3
#> 2       -2641.501         32           3
#> 3       -2641.372         38           3
#> 4       -2641.501         29           3
#> 5       -2641.245         43           3
#> 6       -2641.502         35           3
#> 7       -2641.375         24           3
#> 8       -2641.371         28           3
#> 9       -2641.509         31           3
#> 10      -2642.057         36           3
#> 
#> Use statusCodes() to view the meaning of each status code
#> 
#> Exit Status: 3, Optimization stopped because ftol_rel or ftol_abs was reached.
#>                              
#> Model Type:       Mixed Logit
#> Model Space:       Preference
#> Model Run:            1 of 10
#> Iterations:                38
#> Elapsed Time:        0h:0m:2s
#> Algorithm:     NLOPT_LD_LBFGS
#> Weights Used?:          FALSE
#> Robust?                 FALSE
#> 
#> Model Coefficients: 
#>                     Estimate Std. Error  z-value  Pr(>|z|)    
#> price              -0.438379   0.038877 -11.2762 < 2.2e-16 ***
#> feat_mu             0.382304   0.223530   1.7103 0.0872096 .  
#> brandhiland_mu     -4.146425   0.215717 -19.2216 < 2.2e-16 ***
#> brandweight_mu     -1.306121   0.345491  -3.7805 0.0001565 ***
#> brandyoplait_mu     0.843576   0.125763   6.7076 1.978e-11 ***
#> feat_sigma          2.426846   0.506472   4.7917 1.654e-06 ***
#> brandhiland_sigma   0.022985   0.886992   0.0259 0.9793267    
#> brandweight_sigma   2.115097   0.651051   3.2487 0.0011592 ** 
#> brandyoplait_sigma  0.295287   0.594273   0.4969 0.6192679    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>                                      
#> Log-Likelihood:         -2641.2433390
#> Null Log-Likelihood:    -3343.7419990
#> AIC:                     5300.4866780
#> BIC:                     5352.5806000
#> McFadden R2:                0.2100936
#> Adj McFadden R2:            0.2074020
#> Number of Observations:  2412.0000000
#> 
#> Summary of 10k Draws for Random Coefficients: 
#>                             Min.    1st Qu.     Median       Mean    3rd Qu.
#> feat (normal)         -9.0501506 -1.2562283  0.3804493  0.3789060  2.0165682
#> brandhiland (normal)  -4.2344624 -4.1619360 -4.1464324 -4.1464463 -4.1309337
#> brandweight (normal)  -9.4453960 -2.7338212 -1.3078558 -1.3093844  0.1185107
#> brandyoplait (normal) -0.2827217  0.6439945  0.8430701  0.8428169  1.0422049
#>                            Max.
#> feat (normal)          9.146030
#> brandhiland (normal)  -4.064810
#> brandweight (normal)   6.083848
#> brandyoplait (normal)  1.846463

The above summary table prints summaries of the estimated coefficients as well as a summary table of the distribution of 10,000 population draws for each normally-distributed model coefficient. The results show that the feat attribute has a significant standard deviation coefficient, suggesting that there is considerable heterogeneity across the population for this attribute. In contrast, the brand coefficients have small and insignificant standard deviation coefficients.

Compute the WTP implied from the preference space model:

wtp_mxl_pref <- wtp(mxl_pref, price =  "price")
wtp_mxl_pref
#>                     Estimate Std. Error  z-value  Pr(>|z|)    
#> lambda              0.438379   0.038944  11.2565 < 2.2e-16 ***
#> feat_mu             0.872085   0.525918   1.6582 0.0972739 .  
#> brandhiland_mu     -9.458547   0.675335 -14.0057 < 2.2e-16 ***
#> brandweight_mu     -2.979435   0.731220  -4.0746 4.609e-05 ***
#> brandyoplait_mu     1.924308   0.316067   6.0883 1.141e-09 ***
#> feat_sigma          5.535958   1.329015   4.1655 3.107e-05 ***
#> brandhiland_sigma   0.052431   2.041600   0.0257 0.9795115    
#> brandweight_sigma   4.824817   1.365110   3.5344 0.0004087 ***
#> brandyoplait_sigma  0.673590   1.422614   0.4735 0.6358654    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

WTP space model

This example will estimate the following mixed logit model in the WTP space:

\[\begin{equation} u_{j} = \lambda ( \omega_1 x_{j}^{\mathrm{Feat}} + \omega_2 x_{j}^{\mathrm{Hiland}} + \omega_3 x_{j}^{\mathrm{Weight}} + \omega_4 x_{j}^{\mathrm{Yoplait}} - p_{j}) + \varepsilon_{j} \label{eq:mnlWtpExample} \end{equation}\]

where the parameters \(\omega_1\), \(\omega_2\), \(\omega_3\), and \(\omega_4\) have units of dollars and \(\lambda\) is the scale parameter. In the example below, we will model \(\omega_1\), \(\omega_2\), \(\omega_3\), and \(\omega_4\) as normally distributed across the population. As a result, the model will estimate two parameters for each of these coefficients: a mean (par_mu) and a standard deviation (par_sigma).

Estimate the model using the logitr() function:

mxl_wtp <- logitr(
  data       = yogurt,
  choice     = 'choice',
  obsID      = 'obsID',
  pars       = c('feat', 'brand'),
  price      = 'price',
  randPars   = c(feat = 'n', brand = 'n'),
  modelSpace = 'wtp',
  # You should run a multistart for MXL models since they are non-convex
  numMultiStarts = 10,
  # Use the computed WTP from the preference space model as the starting
  # values for the first run:
  startVals = wtp_mxl_pref$Estimate
)
#> Running multistart...
#>   Iterations: 10
#>   Cores: 3
#>   NOTE: Using user-provided starting values for first iteration
#> Done!

Print a summary of the results:

summary(mxl_wtp)
#> =================================================
#> Call:
#> logitr(data = yogurt, outcome = "choice", obsID = "obsID", pars = c("feat", 
#>     "brand"), price = "price", randPars = c(feat = "n", brand = "n"), 
#>     modelSpace = "wtp", startVals = wtp_mxl_pref$Estimate, numMultiStarts = 10)
#> 
#> Frequencies of alternatives:
#>        1        2        3        4 
#> 0.402156 0.029436 0.229270 0.339138 
#> 
#> Summary Of Multistart Runs:
#>    Log Likelihood Iterations Exit Status
#> 1       -2670.953         95           3
#> 2       -2641.501        145           3
#> 3       -2654.770         94           3
#> 4       -2658.865         91           3
#> 5       -2641.243         85           3
#> 6       -2652.828         77           3
#> 7       -2643.957         87          -1
#> 8       -2641.733        112           3
#> 9       -2641.372         86           3
#> 10      -2652.458         78           3
#> 
#> Use statusCodes() to view the meaning of each status code
#> 
#> Exit Status: 3, Optimization stopped because ftol_rel or ftol_abs was reached.
#>                                  
#> Model Type:           Mixed Logit
#> Model Space:   Willingness-to-Pay
#> Model Run:                5 of 10
#> Iterations:                    85
#> Elapsed Time:            0h:0m:8s
#> Algorithm:         NLOPT_LD_LBFGS
#> Weights Used?:              FALSE
#> Robust?                     FALSE
#> 
#> Model Coefficients: 
#>                     Estimate Std. Error  z-value  Pr(>|z|)    
#> lambda              0.438390   0.038655  11.3411 < 2.2e-16 ***
#> feat_mu             0.876940   0.519094   1.6894 0.0911492 .  
#> brandhiland_mu     -9.458196   0.653132 -14.4813 < 2.2e-16 ***
#> brandweight_mu     -2.990549   0.725301  -4.1232 3.737e-05 ***
#> brandyoplait_mu     1.927919   0.310703   6.2050 5.469e-10 ***
#> feat_sigma          5.540046   1.205409   4.5960 4.307e-06 ***
#> brandhiland_sigma   0.045039   2.026861   0.0222 0.9822718    
#> brandweight_sigma   4.842951   1.351538   3.5833 0.0003393 ***
#> brandyoplait_sigma  0.656073   1.330646   0.4930 0.6219786    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>                                      
#> Log-Likelihood:         -2641.2428565
#> Null Log-Likelihood:    -3343.7419990
#> AIC:                     5300.4857129
#> BIC:                     5352.5796000
#> McFadden R2:                0.2100937
#> Adj McFadden R2:            0.2074021
#> Number of Observations:  2412.0000000
#> 
#> Summary of 10k Draws for Random Coefficients: 
#>                              Min.   1st Qu.     Median       Mean    3rd Qu.
#> feat (normal)         -20.6556311 -2.863529  0.8727069  0.8691839  4.6076677
#> brandhiland (normal)   -9.6307061 -9.488590 -9.4582103 -9.4582376 -9.4278404
#> brandweight (normal)  -21.6271020 -6.259564 -2.9945215 -2.9980215  0.2714395
#> brandyoplait (normal)  -0.5745001  1.484488  1.9267957  1.9262331  2.3692354
#>                            Max.
#> feat (normal)         20.882927
#> brandhiland (normal)  -9.298271
#> brandweight (normal)  13.930311
#> brandyoplait (normal)  4.156144

If you want to compare the WTP from the two different model spaces, use the wtpCompare() function:

wtpCompare(mxl_pref, mxl_wtp, price = 'price')
#>                              pref            wtp  difference
#> lambda                 0.43837864     0.43839031  0.00001167
#> feat_mu                0.87208545     0.87694005  0.00485460
#> brandhiland_mu        -9.45854712    -9.45819583  0.00035129
#> brandweight_mu        -2.97943538    -2.99054893 -0.01111355
#> brandyoplait_mu        1.92430825     1.92791890  0.00361065
#> feat_sigma             5.53595796     5.54004578  0.00408782
#> brandhiland_sigma      0.05243094     0.04503861 -0.00739232
#> brandweight_sigma      4.82481684     4.84295105  0.01813421
#> brandyoplait_sigma     0.67359004     0.65607279 -0.01751725
#> logLik             -2641.24333902 -2641.24285646  0.00048257

Note that the WTP will not necessarily be the same between preference space and WTP space MXL models. This is because the distributional assumptions in MXL models imply different distributions on WTP depending on the model space. See Train and Weeks (2005) and Sonnier, Ainslie, and Otter (2007) for details on this topic.

References

Jain, Dipak C, Naufel J Vilcassim, and Pradeep K Chintagunta. 1994. “A Random-Coefficients Logit Brand-Choice Model Applied to Panel Data.” Journal of Business & Economic Statistics 12 (3): 317–28.
Sonnier, Garrett, Andrew Ainslie, and Thomas Otter. 2007. “Heterogeneity distributions of willingness-to-pay in choice models.” Quant. Mark. Econ. 5 (3): 313–31. https://doi.org/10.1007/s11129-007-9024-6.
Train, Kenneth E., and Melvyn Weeks. 2005. “Discrete Choice Models in Preference and Willingness-to-Pay Space.” In Appl. Simul. Methods Environ. Resour. Econ., 1–16.