R - gmnl package

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.

simple multinomial logit model

(MNL)

Does not account for heterogeneity

scaled multinomial logit model

(SMNL)

Heterogeneity accounted for using a scale factor σ

mixed logit model

(MIXL)

β parameters vary by individual

generalized multinomial logit model

(GMNL)

Combination of scale factors and β''s varying by individual

latent-class model

(LC)

β parameters vary across discrete classes c

mixed-mixed multinomial logit model

(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.

Setting up R

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.

Multinomial logit template

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.

Mixed-logit template

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

 

 

Latent-class model template

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

 

Mixed-mixed logit template

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.

 

 

Comparison of techniques

You may have noticed that the Wizard included calculations for AIC and BIC criteria in each script.  The Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) are commonly used to compare models with different specifications.  Each rewards goodness of fit and penalizes the number of parameters.  In general, BIC punishes models for the number of free parameters more than AIC does.  In comparing models, the ones with lower BIC or AIC scores are considered better.  Of course, in evaluating any model, other considerations should be brought to bear as well.