A free, open-source statistics and graphics package, R has soared to prominence in the last few years, passing even commercial packages in capability. gmnl, an R library devoted to solving discrete-choice models, provides R devotees a way to estimate a variety of logit models that differ in the way they address heterogeneity.

We begin with the simple multinomial
logit model. In this and the following models we use a common notation.
The utility *U
*for
individual *n* choosing
attribute *j* among the
items in choice set *t* is given
by the expressions in the third column. *β
*is a vector of coefficients applied to each attribute, and *ε* represents the unobserved factors
that influence this specific choice. For ease of computation in
all these models, *ε* is assumed
to be distributed iid extreme value. In
the simple model, coefficients do not vary by individual, and unobserved
factors are not correlated with the coefficients.

The scaled multinomial logit model
incorporates heterogeneity by adding a scale factor *σ
* to the formulation. The presumption is that all coefficients
for a given individual are scaled proportionally up or down according
to the value of the scale factor.

By contrast, the mixed logit model
incorporates heterogeneity by presuming that selected coefficients differ
for each individual and conform to a distribution, most often the normal
distribution. This is expressed below by adding parameter *η*, which varies by individual to
coefficient *β*, the mean value
of the coefficient across all individuals.

The generalized multinomial logit model
combines the ideas underlying the scaled and mixed logit models. It
includes both a scale factor s that captures the overall intensity of
preference and a distribution parameter *η*
that describes how coefficients vary across individuals. In addition,
it includes another parameter *γ*
that estimates the proportion of individuals affected by both scale
and varying coefficients compared to those affected by varying coefficients
alone. This model is more flexible than the SMNL or MIXL models
alone.

The latent-class model approaches
heterogeneity in a completely different way. This model assumes
that each individual is a member of a class c, one of a fixed number of
classes. Each class has its own unique set of *β*
coefficients.

The mixed-mixed multinomial logit model takes this approach a step further, nesting mixed-logit models within each latent class. Members of each class have their own mixed-logit model, where coefficients vary according to a distribution for that class.

At the end of this page, we have included a comparison of these techniques applied to the shoe tutorial data set. Several studies have found that no one technique dominates the others, so any conclusion applied to this data set should not be seen as transferrable to others.

The table below summarizes this information.

(MNL) |
Does not account for heterogeneity | ||

(SMNL) |
Heterogeneity accounted for using a scale factor σ |
||

(MIXL) |
β parameters vary
by individual |
||

(GMNL) |
Combination of scale factors and β''s
varying by individual |
||

(LC) |
β parameters vary
across discrete classes c |
||

(MM) |
β parameters vary
by individual, with a unique β
distribution for each class |

gmnl supports all of these models within a common R package. As a command-driven language, R can take time to master, but the Data Wizard provides a way around this obstacle. Along with setting up R data frames, this page illustrates how to set up gmnl and use it to run models listed above.

Generating an R data frame and command scripts

We will make use of our old friend, the sample data file Shoe responses.xlsx, installed with the Data Wizard. Recall that this worksheet contains a hypothetical study of a fictional shoe market and looks like this:

We will start where Step 5 of the Data Wizard leaves off, as shown below.

If you don't feel comfortable getting to this point, review the tutorial through Step 5. At this point when you run the Data Wizard and see the following dialog box.

From here, use the drop-down box in panel A to select R (gmnl package), so that the dialog box looks like this:

Click the Create Files button in panel B. After some processing time, you will see the following pop-up dialog box:

This window contains some very useful information, so you may want to save it by clicking Copy to Word. Click OK to close the window, You should then see the Data Setup Facilitator's dialog box again.

Click Finish to close the Wizard. You can also close Excel, as you won't need it any more. All of the files you need have been created.

You can download a free copy of R from http://cran.r-project.org/bin/windows/base/. Installation is straightforward. When you run R, you should see the console.

Next, if you haven't loaded gmnl, do so by selecting Packages | Install package... from R's menu. If the packages don't appear in the list that appears,

Select gmnl and follow the accompanying instructions.

The multinomial logit template, named MNLTm.r, contains a framework for a simple logit model. The template is located in the same folder as your original data file. From R, click File | Open script ..., navigate to the folder, then click on the script. R should display it as below.

By default, the script runs a model with all attributes and no covariates included. Results get saved in a text file named MNLTm.txt. You can run this, and any other script, from the R console by selecting Edit | Run all from R's menu. Doing that generates the following results.

Model estimated on: Mon Feb 06 10:46:33 AM 2017

Call:

gmnl(formula = Choice ~ 0 + Traditnl + Standard + Price, data = ModelData,

model = "mnl", method = "nr")

Frequencies of categories:

1 2 3 4

0.20187 0.22250 0.31219 0.26344

The estimation took: 0h:0m:0s

Coefficients:

Estimate Std. Error z-value Pr(>|z|)

Traditnl -0.58327752 0.04049000 -14.4055 < 2.2e-16 ***

Standard -0.17935779 0.04126952 -4.3460 1.386e-05 ***

Price 0.00163300 0.00046192 3.5353 0.0004073 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Optimization of log-likelihood by Newton-Raphson maximisation

Log Likelihood: -4309.8

Number of observations: 3200

Number of iterations: 3

Exit of MLE: successive function values within tolerance limit

AIC[1] 8625.63

BIC[1] 8643.843

From the asterisks to the right of the coefficients table, you can see that all attributes are significant; however, Price has the wrong sign. The simple-logit approach is therefore not acceptable.

Scaled multinomial logit template

This
template extends the simple logit model to incorporate estimation of scale
factors. The template is located in the
same folder as your original data file. From R, click File
| Open script ..., navigate to the folder, then click on the script.
R should display it as below.

As before, the script runs a model with all attributes and no covariates included. Results get saved in a text file named SMNLTm.txt. You can run this, and any other script, from the R console by selecting Edit | Run all or by pressing Ctrl-A, then Ctrl-R.. Taking either one of these steps generates the following results:Model estimated on: Mon Feb 06 1:42:26 PM 2017

Call:

gmnl(formula = Choice ~ 0 + Traditnl + Standard + Price, data = ModelData,

model = "smnl", panel = TRUE, method = "bfgs")

Frequencies of categories:

1 2 3 4

0.20187 0.22250 0.31219 0.26344

The estimation took: 0h:1m:49s

Coefficients:

Estimate Std. Error z-value Pr(>|z|)

Traditnl -4.43060077 1.60484879 -2.7608 0.0057667 **

Standard 0.23094163 0.06991785 3.3030 0.0009564 ***

Price -0.00056868 0.00031029 -1.8328 0.0668391 .

tau 2.78524121 0.30199037 9.2229 < 2.2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Optimization of log-likelihood by BFGS maximization

Log Likelihood: -4113.2

Number of observations: 3200

Number of iterations: 65

Exit of MLE: successful convergence

Simulation based on 40 draws

AIC[1] 8234.316

BIC[1] 8258.599

The p-value for tau, the scale parameter, is virtually zero, indicating a significant degree of variance within unobserved factors. Also, the estimate for the price coefficient is negative, though not significantly different from zero. This model is clearly an improvement over the MNL version.

In addition to the above printout, the R template
has saved an Excel-importable *.csv file containing individual parameters
adjusted for scale. Here is the top portion of the file:

The Simulator Wizard can read this file and generate a working, presentable simulator based on this model.

The mixed-logit template is named MixedTm.r. Open it in R by clicking File | Open script... from R's menu. You should see your script appear, as you see here.

By default, the template includes all attributes as independent variables, ASC's and the designation of the first attribute (in this case Fashion) as a random variable having a normal distribution. You can run this script by clicking Edit | Run all from R's menu. Doing so generates the following text that appears on both the R console and in the MixedTm.txt file.

Model estimated on: Fri Feb 10 5:12:24 PM 2017

Call:

gmnl(formula = Choice ~ 0 + Traditnl + Standard + Price, data = ModelData,

model = "mixl", ranp = c(Price = "n"), R = 100, haltons = NA,

panel = TRUE, bi = "mixl.bi", method = "bfgs")

Frequencies of categories:

1 2 3 4

0.20187 0.22250 0.31219 0.26344

The estimation took: 0h:6m:14s

Coefficients:

Estimate Std. Error z-value Pr(>|z|)

Traditnl -5.8328e-01 4.0490e-02 -14.4055 < 2.2e-16 ***

Standard -1.7936e-01 4.1270e-02 -4.3460 1.386e-05 ***

Price 1.6330e-03 4.6192e-04 3.5353 0.0004073 ***

sd.Price 5.2711e-07 9.3734e-04 0.0006 0.9995513

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Optimization of log-likelihood by BFGS maximization

Log Likelihood: -4309.8

Number of observations: 3200

Number of iterations: 88

Exit of MLE: successful convergence

Simulation based on 100 draws

AIC[1] 8627.63

BIC[1] 8651.914

The means of all attributes are significant, as is the standard deviation of Fashion. As a next step, let's designate Quality as a random variable. We do that by adding Quality = 'n' to the rpar phrase in the gmnl command.

# Set up initial model.

SWData01.mxl <- gmnl(Choice ~ Fashion + Quality + Price,

ModelData,

rpar = c(Fashion = 'n', Quality = 'n'),

R = 100, halton = NA)

Running this model generates

Call:

gmnl(formula = Choice ~ Fashion + Quality + Price, data = ModelData,

rpar = c(Fashion = "n", Quality = "n"), R = 100, halton = NA)

Frequencies of alternatives:

1 2 3 4

0.20187 0.22250 0.31219 0.26344

bfgs method

16 iterations, 0h:0m:54s

g'(-H)^-1g = 5.28E-07

gradient close to zero

Coefficients :

Estimate Std. Error t-value Pr(>|t|)

2:(intercept) 0.182390 0.106912 1.7060 0.08801 .

3:(intercept) 0.048496 0.101486 0.4779 0.63275

4:(intercept) -0.147533 0.146101 -1.0098 0.31259

Fashion 2.000847 0.155008 12.9080 < 2.2e-16 ***

Quality 1.152848 0.128984 8.9379 < 2.2e-16 ***

Price -0.029541 0.003099 -9.5324 < 2.2e-16 ***

sd.Fashion 2.729048 0.476403 5.7284 1.014e-08 ***

sd.Quality 2.967424 0.545057 5.4442 5.202e-08 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -4121.8

McFadden R^2: 0.061335

Likelihood ratio test : chisq = 538.67 (p.value = < 2.22e-16)

random coefficients

Min. 1st Qu. Median Mean 3rd Qu. Max.

Fashion -Inf 0.1601321 2.000847 2.000847 3.841562 Inf

Quality -Inf -0.8486497 1.152848 1.152848 3.154345 Inf

The coefficient for the standard deviation of Quality passes its significance test, along with most other coefficients.

The exception is ASC's, all of which fail significance tests. ASC's have no meaning in the context of this model, so we can drop them from the next specification. We do this by placing 0 + in front of the independent variable list within the gmnl phrase, so our new specification looks like this.

# Set up initial model.

SWData02.mxl <- gmnl(Choice ~ 0 + Fashion + Quality + Price,

ModelData,

rpar = c(Fashion = 'n', Quality = 'n'),

R = 100, halton = NA)

Running this model generates the following output:

Call:

gmnl(formula = Choice ~ 0 + Fashion + Quality + Price, data = ModelData,

rpar = c(Fashion = "n", Quality = "n"), R = 100, halton = NA)

Frequencies of alternatives:

1 2 3 4

0.20187 0.22250 0.31219 0.26344

bfgs method

16 iterations, 0h:0m:52s

g'(-H)^-1g = 5.68E-08

gradient close to zero

Coefficients :

Estimate Std. Error t-value Pr(>|t|)

Fashion 1.9576656 0.1317055 14.8640 < 2.2e-16 ***

Quality 1.0465534 0.1331903 7.8576 3.997e-15 ***

Price -0.0259093 0.0016055 -16.1379 < 2.2e-16 ***

sd.Fashion 2.9635255 0.4664983 6.3527 2.116e-10 ***

sd.Quality 3.3654954 0.5563404 6.0493 1.454e-09 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -4123.8

random coefficients

Min. 1st Qu. Median Mean 3rd Qu. Max.

Fashion -Inf -0.04120195 1.957666 1.957666 3.956533 Inf

Quality -Inf -1.22343879 1.046553 1.046553 3.316546 Inf

All coefficients pass significance tests. Let's continue by specifying Price as a random variable with a normal distribution.

# Set up initial model.

SWData03.mxl <- gmnl(Choice ~ 0 + Fashion + Quality + Price,

ModelData,

rpar = c(Fashion = 'n', Quality = 'n', Price = 'n'),

R = 100, halton = NA)

Running this model gives us

gmnl(formula = Choice ~ 0 + Fashion + Quality + Price, data = ModelData,

rpar = c(Fashion = "n", Quality = "n", Price = "n"), R = 100,

halton = NA)

Frequencies of alternatives:

1 2 3 4

0.20187 0.22250 0.31219 0.26344

bfgs method

32 iterations, 0h:1m:53s

g'(-H)^-1g = 6.89E-07

gradient close to zero

Coefficients :

Estimate Std. Error t-value Pr(>|t|)

Fashion 1.9596142 0.1321759 14.8258 < 2.2e-16 ***

Quality 1.0445642 0.1342316 7.7818 7.105e-15 ***

Price -0.0259201 0.0016094 -16.1057 < 2.2e-16 ***

sd.Fashion 2.9643739 0.4676684 6.3386 2.318e-10 ***

sd.Quality 3.3811353 0.5598058 6.0398 1.543e-09 ***

sd.Price 0.0016672 0.0184995 0.0901 0.9282

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -4123.8

random coefficients

Min. 1st Qu. Median Mean 3rd Qu. Max.

Fashion -Inf -0.03982563 1.95961418 1.95961418 3.95905398 Inf

Quality -Inf -1.23597691 1.04456423 1.04456423 3.32510536 Inf

Price -Inf -0.02704459 -0.02592011 -0.02592011 -0.02479563 Inf

The coefficient for the standard deviation of Price roundly fails its significance test, so we reject the hypothesis that Price is a random variable. Our best model is therefore the prior one. You can submit the output file from the prior model to StatWizards' Simulator Builder Wizard to construct a full-featured simulator in Excel.

Generalized multinomial logit model

For the generalized multinomial logit model,
the Data Wizard generates a template named GMNLTm.r.
Open it in R by clicking File |
Open script..., then selecting the name from R's menu. The
template for the shoe model looks like this:

Pressing Ctrl-A then Ctrl-R runs this script, yielding the following results.

Model estimated on: Fri Feb 10 5:29:23 PM 2017

Call:

gmnl(formula = Choice ~ 0 + Traditnl + Standard + Price, data = ModelData,

model = "gmnl", ranp = c(Price = "n"), R = 100, method = "bfgs")

Frequencies of categories:

1 2 3 4

0.20187 0.22250 0.31219 0.26344

The estimation took: 0h:11m:27s

Coefficients:

Estimate Std. Error z-value Pr(>|z|)

Traditnl -3.4587511 1.1100399 -3.1159 0.001834 **

Standard 10.1388772 2.5519411 3.9730 7.097e-05 ***

Price -0.0431272 0.0297845 -1.4480 0.147625

sd.Price 0.0053192 NA NA NA

tau 21.1141916 NA NA NA

gamma -2.1878966 NA NA NA

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Optimization of log-likelihood by BFGS maximization

Log Likelihood: -4431.1

Number of observations: 3200

Number of iterations: 104

Exit of MLE: successful convergence

Simulation based on 100 draws

AIC[1] 8874.118

BIC[1] 8910.544

For the latent-class model, the Data Wizard generates a template named LCTm.r. Open it in R by clicking File | Open script..., then selecting the name from R's menu. The template for the shoe model looks like this:

Pressing Ctrl-A then Ctrl-R runs this script, yielding the following results.

Model estimated on: Fri Feb 10 5:30:18 PM 2017

Call:

gmnl(formula = Choice ~ 0 + Traditnl + Standard + Price | 0 |

0 | 0 | 1, data = ModelData, model = "lc", Q = 2, panel = TRUE,

method = "bfgs")

Frequencies of categories:

1 2 3 4

0.20187 0.22250 0.31219 0.26344

The estimation took: 0h:0m:2s

Coefficients:

Estimate Std. Error z-value Pr(>|z|)

class.1.Traditnl -2.04860251 0.10383216 -19.7299 < 2.2e-16 ***

class.1.Standard 0.63025947 0.06193638 10.1759 < 2.2e-16 ***

class.1.Price -0.00156352 0.00061907 -2.5256 0.01155 *

class.2.Traditnl 0.80788380 0.07358715 10.9786 < 2.2e-16 ***

class.2.Standard -1.41166627 0.09712171 -14.5350 < 2.2e-16 ***

class.2.Price 0.00365128 0.00092884 3.9310 8.458e-05 ***

(class)2 -0.46246341 0.04116283 -11.2350 < 2.2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Optimization of log-likelihood by BFGS maximization

Log Likelihood: -3835.2

Number of observations: 3200

Number of iterations: 88

Exit of MLE: successful convergence

AIC[1] 7684.31

BIC[1] 7726.806

For the latent-class model, the Data
Wizard generates a template named MMTm.r.
Open it in R by clicking File | Open script..., then
selecting the name from R's menu. The template for the shoe model
looks like this:

As of this writing, a bug in gmnl prevents the script from running. As soon as a fix becomes available, we will present the results.

You may have
noticed that the Wizard included calculations for AIC and
BIC criteria in each script.
The ** Akaike
Information Criterion** (AIC) and