---
title: "Estimating Mixed Logit Models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimating Mixed Logit Models}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
bibliography: "`r here::here('vignettes', 'library.bib')`"
---

```{r setup, include=FALSE, message=FALSE, warning=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.retina = 3,
  comment = "#>"
)

library(logitr)

# Read in results from already estimated models so that the
# examples aren't actually run when building this page,
# otherwise it'll take much longer to build
mxl_pref <- readRDS(here::here('inst', 'extdata', 'mxl_pref.Rds'))
mxl_wtp  <- readRDS(here::here('inst', 'extdata', 'mxl_wtp.Rds'))
mxl_pref_cor <- readRDS(here::here('inst', 'extdata', 'mxl_pref_cor.Rds'))
```

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.

# Supported distributions

The mixed logit model is a popular approach for modeling unobserved heterogeneity across individuals, which is implemented by assuming that parameters vary randomly across individuals according to a chosen distribution [@McFadden2000]. A mixed logit model is specified by setting the `randPars` argument in the `logitr()` function equal to a named vector defining parameter distributions. In the example below, we set `randPars = c(feat = 'n', brand = 'n')` so that `feat` and `brand` are normally distributed. The current package version supports the following distributions:

- Normal: `"n"`
- Log-normal: `"ln"`
- Zero-censored normal: `"cn"`

Mixed logit models will estimate a mean and standard deviation of the underlying normal distribution for each random coefficient. Note that log-normal or zero-censored normal parameters force positivity, so when using these it is often necessary to use the negative of a value (e.g. for "price", which typically has a negative coefficient). Mixed logit models in `logitr` assume uncorrelated heterogeneity covariances by default, though full covariances can be estimated using the `correlation = TRUE` argument. For WTP space models, the `scalePar` parameter can also be modeled as following a random distribution by setting the `randScale` argument equal to `"n"`, `"ln"`, or `"cn"`.

# The data

```{r, child='../man/rmdchunks/yogurtDataDescription.Rmd'}
```

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

```{r, child='../man/rmdchunks/mnlPrefExampleTable.Rmd'}
```

# Preference space model

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

```{r, child='../man/rmdchunks/mnlPrefExample.Rmd'}
```

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 a mean and standard deviation for each of these coefficients.

Note that since the `yogurt` data has a panel structure (i.e. multiple choice observations for each respondent), it is necessary to set the `panelID` argument to the `id` variable, which identifies the individual. This will use the panel version of the log-likelihood (see [Train 2009 chapter 6](https://eml.berkeley.edu/books/choice2nd/Ch06_p134-150.pdf), section 6.7 for details).

Finally, as with WTP space models, it is recommended to use a multistart search for mixed logit models as they are non-convex. This is implemented in the example below by setting `numMultiStarts = 10`:

```{r eval=FALSE}
library("logitr")

set.seed(456)

mxl_pref <- logitr(
  data     = yogurt,
  outcome  = 'choice',
  obsID    = 'obsID',
  panelID  = 'id',
  pars     = c('price', 'feat', 'brand'),
  randPars = c(feat = 'n', brand = 'n'),
  numMultiStarts = 10
)
```

Print a summary of the results:

```{r}
summary(mxl_pref)
```

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:

```{r}
wtp_mxl_pref <- wtp(mxl_pref, scalePar =  "price")
wtp_mxl_pref
```

# WTP space model

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

```{r, child='../man/rmdchunks/mnlWtpExample.Rmd'}
```

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. Note that this is a slightly different assumption than in the preference space model. In the WTP space, we are assuming that the WTP for these features is normally-distributed (as opposed to the preference space model where the utility coefficients are assumed to follow a normal distribution).

In the example below, we also use a 10-iteration multistart. We also set the starting values for the first iteration to the computed WTP from the preference space model:

```{r eval=FALSE}
set.seed(6789)

mxl_wtp <- logitr(
  data       = yogurt,
  outcome    = 'choice',
  obsID      = 'obsID',
  panelID    = 'id',
  pars       = c('feat', 'brand'),
  scalePar   = 'price',
  randPars   = c(feat = 'n', brand = 'n'),
  numMultiStarts = 10,
  startVals = wtp_mxl_pref$Estimate
)
```

Print a summary of the results:

```{r}
summary(mxl_wtp)
```

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

```{r}
wtpCompare(mxl_pref, mxl_wtp, scalePar = 'price')
```

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 [-@Train2005] and Sonnier, Ainslie, and Otter [-@Sonnier2007] for details on this topic.

# Correlated heterogeneity

By default, `logitr` assumes that mixed logit models have uncorrelated heterogeneity. However, correlated heterogeneity can be implemented by setting `correlation = TRUE` for models in either space (preference or WTP). The example below shows the results for the same mixed logit model in the preference space as above but now with correlated heterogeneity:

```{r eval=FALSE}
library("logitr")

set.seed(456)

mxl_pref_cor <- logitr(
  data     = yogurt,
  outcome  = 'choice',
  obsID    = 'obsID',
  panelID  = 'id',
  pars     = c('price', 'feat', 'brand'),
  randPars = c(feat = 'n', brand = 'n'),
  numMultiStarts = 10,
  correlation = TRUE
)
```

Print a summary of the results:

```{r}
summary(mxl_pref_cor)
```

# References