multinomineq: Multinomial Models with Inequality Constraints

Daniel W. Heck

2019-01-30

The R package multinomineq facilitates the Bayesian analysis of discrete/categorical data using multinomial models with inequality constraints. Essentially, these types of models allow to analyze discrete choice frequencies while assuming that the corresponding choice probabilities \(\theta\) are constrained by linear inequality constraints (e.g., by the simple order \(\theta_{1} \leq \theta_{2}\leq \theta_{3}\leq \theta_{4}\)). In the following, this vignette shows how to define and test such models in R. The second section below provides a concise formal definition of the model class.

Example: Testing the Description-Experience Gap

The following example shows how inequality constraints can be used to test the description-experience (DE) gap in risky decision making (Hertwig et al., 2004). In the experiment, participants had to choose one of two risky gambles in each of 6 different paired comparisons (e.g., a lottery in which you win 10€ with \(p=.20\) versus a lottery in which you can win 20€ with \(p=.45\)). The description-experience gap states that people underweigh small probabilities \(p\) when making decisions based on experience. We will conduct a reanalysis of Hertwig’s data set using the inequality-constrained multinomial model in Regenwetter & Robinson (2017).

First, it is necessary to install the R package multinomineq from Github. Note that it is necessary to install the developer R tools first (see https://support.rstudio.com/hc/en-us/articles/200486498-Package-Development-Prerequisites). Afterwards, multinomineq can be installed and loaded via:

### uncomment to install packages:
# install.packages("devtools", "coda", "RcppArmadillo", "RcppProgress", "Rglpk", "quadprog")
# devtools::install_github("danheck/multinomineq")

library("multinomineq")
library("coda")
set.seed(1234)

Data structure

We reanalyze the data by Hertwig et al. (2004) who collected choice frequencies of \(n=25\) participants for 6 paired comparisons of risky gambles. The observed choice frequencies can simply be summarized by the frequencies of choosing “Option H” (a specific lottery, cf. original article) in each of the 6 paired comparisons:

Vertex (\(V\)-)representation

First, it is shown how to use the vertex (\(V\)-)representation for inequality-constrained models. Essentially, the rows of the matrix \(V\) list all patterns that are predicted by a substantive theory of interest. Each prediction either states that “Option H” is chosen (value 1) or not (value 0). For example, if the theory makes the prediction that “Option H” is chosen only in the 4th paired comparison but otherwise not, the matrix \(V\) needs to include the following pattern in one row: c(0,0,0,1,0,0). Similarly, all possible predictions need to be combined to form the matrix \(V\). Inequality-constrained models assume a mixture distribution over all these possible predicted patterns (e.g., due to heterogeneity across participants or trials). Geometrically, the corresponding parameter space thus represents the convex hull of the predicted patterns/vertices (i.e., all possible linear combinations of the rows of \(V\) when using nonnegative weights).

For the description-experience gap, Regenwetter & Robinson (2017) used a specific decision making model (Cumulative Prospect Theory) to generate all possible predictions. The matrix \(V\) lists all of these predictions when assuming that participants underweigh small probabilities \(p\) of winning money. The following matrix contains all predicted patterns for the 6 risky gambles:

Given this representation of the theoretical predictions, it is not easy to check whether a given vector of choice probabilities for the 6 gambles is inside the convex hull of the predicted patterns. As a remedy, one can use the function inside() to test this, or the function find_inside() to find a probability vector that is inside the convex hull of the predictions:

Inequality (\(Ab\)-)representation

Instead of listing the vertices (or edges) of the inequality-constrained parameter space, we may instead define a set of linear inequality constraints on the parameters explicitly. For this purpose, we specify a matrix \(A\) and a vector \(b\). The model only allows parameter vectors \(\theta\) that satisfy the inequalities defined by the matrix equation: \(A \cdot \theta \leq b\).

For the description-experience gap, one can show that instead of listing the predictions by the matrix \(V\) as shown above, one can equivalently describe the linear inequalities as follows:

For example, the third row of \(A\) and the third value of \(b\) specify that \(-1 \cdot \theta_1 + 1 \cdot \theta_2 + 0 \cdot \theta_3+ 0 \cdot \theta_4+ 0 \cdot \theta_5+ 0 \cdot \theta_6 \leq 0\) (which is equivalent to \(\theta_2 \leq \theta_1\)).

Using the matrix \(A\) and the vector \(b\), one can use standard matrix multiplication in R (%*%) to check whether a given vector of choice probabilities satisfies all of the inequality constraints:

Note that it is sometimes possible to convert the vertex (\(V\)-)representation (convex hull) to the inequality (\(Ab\)-)representation using specific software. To do this in R, it is necessary to install the package rPorta from Github:

Posterior sampling

In a Bayesian framework, inequality-constrained multinomial models can be fitted by drawing posterior samples for the parameter vector of choice probabilities \(\theta\). If the model is defined in terms of a matrix \(V\) with predicted patterns, we obtain \(M=2,000\) posterior samples as follows:

If the model is defined in terms of the inequality \(Ab\)-representation, we obtain posterior samples in a similar way by replacing the argument V=V by the two arguments A=A, b=b:

If possible, the \(Ab\)-representation should be used because the sampling is much more efficient than for the \(V\)-representation. Irrespective of how the posterior samples have been obtained, we can use standard methods in Bayesian analysis to check the convergence of the MCMC sampler and obtain summary statistics of the posterior distribution:

Similarly, we can use posterior-predictive checks to test whether an inequality-constrained model fits the data adequately (i.e., whether the constraints hold empirically). Essentially, this method computes Pearson’s \(X^2\)-statistic to measure the discrepancy both for the observed and the posterior-predicted frequencies. If the model does not fit the data, the posterior-predicted \(p\)-value ppp will become very small (and it will be close to 0.50 if the model fits). Hence, we usually apply a criterion such as ppp < .05 to reject models that do not fit the data, as in the following example:

Bayes factor

Often, one is interested in quantifying the evidence in favor of the inequality constraints. We can use Bayesian model selection to compare the inequality-constrained model \(M_0\) against the unconstrained model \(M_u\) that does not restrict the choice probabilities \(\theta\) (also called encompassing model). The Bayes factor bf_0u quantifies the evidence in favor of the inequality-constrained model \(M_0\) versus the unconstrained model \(M_u\). Similarly, bf_u0 is the inverse of bf_0u and provides the evidence for the unconstrained model. Moreover, the Bayes factor bf_00' compares the inequality-constrained model \(M_0\) against the model \(M_{0'}\) which has a parameter space defined as the complement of the parameter space of \(M_0\).

Bayes factors are computed as follows:

Note that multinomineq also provides an error of approximation for the Bayes factor to account for random sampling error. The matrix that is provided as output shows the standard error se for the Bayes factor and the 90% credibility interval ci.5% and ci.95%.

Increased efficiency for Bayes factor computations

Within the function bf_binom, the encompassing Bayes factor is computed in two steps based on the unconstrained multinomial model \(M_u\):

  1. Draw prior parameter samples from \(M_u\) and count the proportion of samples \(c\) inside the inequality-constrained parameter space.
  2. Draw posterior parameter samples from \(M_u\) and count the proportion of samples \(f\) inside the inequality-constrained parameter space.

Based on these two proportions, the Bayes factor in favor of the inequality constraints is approximated as \(B_{0u}=f/c\). It might be convenient to perform these computations in separate steps in R for some models and datasets. For example, the prior constant \(c\) is sometimes known analytically. Moreover, when analyzing multiple datasets or participants, it is more efficient to estimate the prior constant \(c\) only once because \(c\) is identical for all analyses:

Since Bayes factors can be hard to compute for some models and data sets, multinomineq provides a stepwise procedure that aims at increasing the efficiency. The user can specify a minimum number of samples cmin that must satisy a subset of the inequality constraints. The algorithm keeps sampling until this threshold is reached. The argument steps determines at which rows of the matrix \(A\) the inequalities are split into subsets of increasing size. For example, steps = c(3, 5, 7) means that the first subset contains only the first 3 inequalities, the second subset the first 5 inequalities etc.:

Data format for multinomial frequencies

In the present example, we only observed binomial data, which can easily be described by the frequency of choosing one of two options. If more than two choice options are provided for a comparison, a different data format has to be used for multinomial data. In the binomial case, we ommitted the frequencies of not choosing a specific option because we knew the total number \(n\). For multinomial data, one needs to specify all choice frequencies \(k_{ij}\). Moreover, one needs to specify a vector options that defines how many choice options are available for each item type (e.g., options = c(3,3,) means that there are two scenarios, each with three possible options).

Using this notation, the 6 binomial choice frequencies above can equivalently be specified by the multinomial data format as follows:

The binomial format can also be translated to the multinomial format automatically:

The multinomial data format is used as input for the analysis as follows:

Formal definition of multinomial models with inequality constraints

In the following, we define inequality-constrained multinomial models for discrete data. We use the term “item type” to refer to a category system \(i\) that is modeled by a multinomial distribution with a fixed total number of observations \(n_i\). For an item type \(i\), we label the observed frequencies as \(k_{ij}\).

Parameters and Likelihood

The parameter vector of interest refers to the choice probabilities: \[\theta = (\theta_{11},\dots,\theta_{1 (J_1-1)},\theta_{21},\dots, \theta_{I (J_I-1)})\] Note that the last probability within each item type is ommited because probabilities have to sum to one. Hence, for \(n_i=3\) choice options there are only 2 free parameters.

Using this notation, we can define a product-multinomial likelihood function: \[p(k \mid \theta) = \prod_{i=1}^I {n_i \choose k_{i1} ,\dots ,k_{i J_i} } \prod_{j=1}^{J_i - 1} \theta_{ij}^{k_{ij}} \left(1- \sum_{r=1}^{J_i-1} \theta_{ir}\right)^{k_{J_i}}\]

Linear Inequality Constraints

Inequality constraints can be defined in two different, equivalent ways. First, one may define a list of inequality constraints that have to hold for the choice probabilities, for instance: \[\begin{align*} \theta_{11} \leq \theta_{21} \, \text{ and } \, \theta_{21} \leq \theta_{31} \end{align*}\] These inequalities can be expressed in vector notation as \(A \,\theta \leq b\) by defining a matrix $ A$ and a vector $ b$ as follows: \[ A = \pmatrix{1 & -1 & 0 \\ 0 & 1 & -1} \,\, \text{ and }\,\, b = \pmatrix{0 \\ 0}\] To formalize this notation, we define the constrained parameter space \(\Omega_0\) via a matrix \(A\) and a vector \(b\) that specify a list of inequalities that have to hold: \[\Omega_0 = \left\{\theta \in \Omega \, \middle | \, A \, \theta \leq b \right\}\]

Alternatively, one may list all the vertices \(v^{(s)}\) defining the constrained parameter space in the rows of a matrix \(V\). In the example above, we have \[ v^{(1)} = \pmatrix{0\\0\\0} \,,\, v^{(2)} = \pmatrix{0\\0\\1} \,,\, v^{(3)} = \pmatrix{0\\1\\1} \,,\, v^{(1)} = \pmatrix{1\\1\\1}\] These vectors can be convenietly summarizes in a matrix $ V$ with one vertex per row: \[ V = \pmatrix{0 & 0 & 0\\ 0 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & 1 & 1}\] In general, the restricted parameter space \(\Omega_0\) is defined as the convex hull of the set of vertices \(v^{(s)}\): \[\Omega_0 = \left\{ \theta = \sum_{s=1}^S \alpha_s v^{(s)} \,\middle |\, \alpha_s \geq 0 \text{ for all } s=1,\dots, S \text{ and }\sum_{s=1}^S \alpha_s=1 \right\}\]

Prior & Posterior

As a prior distribution for the parameters \(\theta\), we assume independent Dirichlet distributions for the different item types \(i\) with shape parameters \(\alpha_{ij}\). Note that we only allow parameters that are inside the constrained parameter space \(\Omega_0\): \[p(\theta) = \frac {1}{c} \, \mathbb I_{\Omega_0}(\theta) \, \prod_{i=1}^I \prod_{j=1}^{J_i} \theta_{ij}^{\alpha_{ij}-1}\] Here, \(\mathbb I_{\Omega_0}(\theta)\) is the indicator function which is one if \(\theta\in\Omega_0\) and zero otherwise. Moreover, the normalizing constant \(c\) ensures that the prior density function integrates to one.

Finally, the analytical solution of the posterior distribution is also a product-Dirichlet distribution, truncated to the constrained parameter space: \[p(\theta \mid k) = \frac {1}{f} \, \mathbb I_{\Omega_0}(\theta) \, \prod_{i=1}^I \prod_{j=1}^{J_i} \theta_{ij}^{k_{ij} + \alpha_{ij}-1}\] Similarly as before, \(f\) is the normalizing constant of the posterior density function.

Reference

For a more detailed introduction with technical details, see: