---
title: "Estimating Weighted Logit Models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimating Weighted 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 = "#>"
)
```

This vignette demonstrates an example of how to use the `logitr()` function with the `weights` argument to estimate weighted logit models.

# The data

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

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

<div style="width:550px">
| Symbol | Variable                             |
|--------|--------------------------------------|
| $p$                        | The price in US dollars.       |
| $x_{j}^{\mathrm{hev}}$ | Dummy variable for HEV vehicle type |
| $x_{j}^{\mathrm{phev10}}$ | Dummy variable for PHEV10 vehicle type |
| $x_{j}^{\mathrm{phev20}}$ | Dummy variable for PHEV20 vehicle type |
| $x_{j}^{\mathrm{phev40}}$ | Dummy variable for PHEV40 vehicle type |
| $x_{j}^{\mathrm{bev75}}$ | Dummy variable for BEV75 vehicle type |
| $x_{j}^{\mathrm{bev100}}$ | Dummy variable for BEV100 vehicle type |
| $x_{j}^{\mathrm{bev150}}$ | Dummy variable for BEV150 vehicle type |
| $x_{j}^{\mathrm{phevFastcharge}}$ | Dummy variable for if the PHEV has a fast charging capability |
| $x_{j}^{\mathrm{bevFastcharge}}$ | Dummy variable for if the BEV has a fast charging capability |
| $x_{j}^{\mathrm{opCost}}$ | The vehicle operating costs (cents / mile) |
| $x_{j}^{\mathrm{accelTime}}$ | The vehicle 0-60mph acceleration time |
| $x_{j}^{\mathrm{american}}$ | Dummy variable for an American brand |
| $x_{j}^{\mathrm{japanese}}$ | Dummy variable for a Japanese brand |
| $x_{j}^{\mathrm{chinese}}$ | Dummy variable for a Chinese brand |
| $x_{j}^{\mathrm{skorean}}$ | Dummy variable for a S. Korean brand |
</div>

# The utility model

In this example, we'll estimate two versions of the following utility model in the WTP space: one without weights and one with weights. Notation is taken from Helveston et al. [-@Helveston2015]:

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

where all the $\omega$ parameters have units of dollars and $\lambda$ is the scale parameter.

# Unweighted model

Estimate the unweighted model using the `logitr()` function. In this example, I have set `robust = TRUE` since it will also be `TRUE` in the weighted model:

```{r}
library("logitr")

mnl_wtp_unweighted <- logitr(
  data    = cars_us,
  outcome = 'choice',
  obsID   = 'obsnum',
  pars    = c(
    'hev', 'phev10', 'phev20', 'phev40', 'bev75', 'bev100', 'bev150',
    'american', 'japanese', 'chinese', 'skorean', 'phevFastcharge',
    'bevFastcharge','opCost', 'accelTime'),
  scalePar   = 'price',
  robust     = TRUE,
  # Since WTP space models are non-convex, run a multistart
  numMultiStarts = 10
)
```

Print a summary of the results:

```{r}
summary(mnl_wtp_unweighted)
```

# Weighted model

To estimate the weighted model, simply add the `weights` argument to the call to  `logitr()`, referring to the column of weights that will be used to weight each choice observation. In this example, the weights used in the `weights` column range from 0.2 to 5:

```{r}
summary(cars_us$weights)

mnl_wtp_weighted <- logitr(
  data    = cars_us,
  outcome = 'choice',
  obsID   = 'obsnum',
  pars    = c(
    'hev', 'phev10', 'phev20', 'phev40', 'bev75', 'bev100', 'bev150',
    'american', 'japanese', 'chinese', 'skorean', 'phevFastcharge',
    'bevFastcharge','opCost', 'accelTime'),
  scalePar = 'price',
  weights  = 'weights', # This enables the weights
  robust   = TRUE,
  numMultiStarts = 10
)
```

Print a summary of the results:

```{r}
summary(mnl_wtp_weighted)
```

# Compare results

Here is a comparison of the coefficients between the weighted and unweighted models. All of the significant coefficients have the same sign, but the magnitudes shift some based on the differential weighting of each individual choice in the weighted model:

```{r}
data.frame(
  Unweighted = coef(mnl_wtp_unweighted),
  Weighted   = coef(mnl_wtp_weighted)
)
```

Here is a comparison of the log-likelihood for the weighted and unweighted models:

```{r}
c(
  "Unweighted" = mnl_wtp_unweighted$logLik,
  "Weighted" = mnl_wtp_weighted$logLik
)
```

# References