Introduction

bayesm’s posterior sampling function rhierMnlRwMixture permits the imposition of sign constraints on the individual-specific parameters of a hierarchical multinomial logit model. This may be desired if, for example, the researcher believes there are heterogenous effects from, say, price, but that all responses should be negative (i.e., sign-constrained). This vignette provides exposition of the model, discussion of prior specification, and an example.

Model

The model follows the hierarchical multinomial logit specification given in Example 3 of the “bayesm Overview” Vignette, but will be repeated here succinctly. Individuals are assumed to be rational economic agents that make utility-maximizing choices. Utility is modeled as the sum of deterministic and stochastic components, where the inverse-logit of the probability of chosing an alternative is linear in the parameters and the error is assumed to follow a Type I Extreme Value distribution:

$$ U_{ij} = X_{ij}\beta_i + \varepsilon_{ij} \hspace{0.8em} \text{with} \hspace{0.8em} \varepsilon_{ij}\ \sim \text{ iid Type I EV} $$

These assumptions yield choice probabilities of:

$$ \text{Pr}(y_i=j) = \frac{\exp \{x_{ij}'\beta_i\}}{\sum_{k=1}^p\exp\{x_{ik}'\beta_i\}} $$

x_i is n_i × k and i = 1, …, N. There are p alternatives, j = 1, …, p. An outside option, often denoted j = 0 can be introduced by assigning 0’s to that option’s covariate (x) values.

We impose sign constraints by defining a k-length constraint vector SignRes that takes values from the set {−1, 0, 1} to define β_ik = f(β_ik^*) where f(⋅) is as follows:

$$ \beta_{ik} = f(\beta_{ik}^*) = \left\{ \begin{array}{lcl} \exp(\beta_{ik}^*) & \text{if} & \texttt{SignRes[k]} = 1 \\ \beta_{ik}^* & \text{if} & \texttt{SignRes[k]} = 0 \\ -\exp(\beta_{ik}^*) & \text{if} & \texttt{SignRes[k]} = -1 \\ \end{array} \right. $$

The “deep” individual-specific parameters (β_i^*) are assumed to be drawn from a mixture of M normal distributions with mean values driven by cross-sectional unit characteristics Z. That is, β_i^* = z_i′Δ + u_i where u_i has a mixture-of-normals distribution.¹

Considering β_i^* a length-k row vector, we will stack the N β_i^*’s vertically and write:

B = ZΔ + U Thus we have β_i, z_i, and u_i as the i^th rows of B, Z, and U. B is N × k, Z is N × M, Δ is M × k, and U is N × k where the distribution on U is such that:

$$ \Pr(\beta_{ik}^*) = \sum_{m=1}^M \pi_m \phi(z_i' \Delta \vert \mu_j, \Sigma_j) $$

ϕ is the normal pdf.

Priors

Natural conjugate priors are specified:

π ∼ Dirichlet(a) vec(Δ) = δ ∼ MVN(δ̄, A_δ⁻¹) μ_m ∼ MVN(μ̄, Σ_m ⊗ a_μ⁻¹) Σ_m ∼ IW(ν, V)

This specification of priors assumes that the (μ_m, Σ_m) are independent and that, conditional on the hyperparameters, the β_i’s are independent.

a implements prior beliefs on the number of normal components in the mixture with a default of 5. ν is a “tightness” parameter of the inverted-Wishart distribution and V is its location matrix. Without sign constraints, they default to ν = k + 3 and V = νI, which has the effect of centering the prior on I and making it “barely proper”. a_μ is a tightness parameter for the priors on μ, and when no sign constraints are imposed it defaults to an extremely diffuse prior of 0.01.

These defaults assume the logit coefficients (β_ik’s) are on the order of approximately 1 and, if so, are typically reasonable hyperprior values. However, when sign constraints are imposed, say, SignRes[k]=-1 such that β_ik = −exp {β_ik^*}, then these hyperprior defults pile up mass near zero — a result that follows from the nature of the exponential function and the fact that the β_ik^*’s are on the log scale. Let’s show this graphically.

# define function
drawprior <- function (mubar_betak, nvar, ncomp, a, nu, Amu, V, ndraw) {
  betakstar <- double(ndraw)
  betak     <- double(ndraw)
  otherbeta <- double(ndraw)
  mubar     <- c(rep(0, nvar-1), mubar_betak)
  
  for(i in 1:ndraw) {
    comps=list()
    for(k in 1:ncomp) {
      Sigma <- rwishart(nu,chol2inv(chol(V)))$IW
      comps[[k]] <- list(mubar + t(chol(Sigma/Amu)) %*% rnorm(nvar), 
                         backsolve(chol(Sigma), diag(1,nvar)) )
    }
    pvec         <- rdirichlet(a)
    beta         <- rmixture(1,pvec,comps)$x
    betakstar[i] <- beta[nvar]
    betak[i]     <- -exp(beta[nvar])
    otherbeta[i] <- beta[1]
  }
  
  return(list(betakstar=betakstar, betak=betak, otherbeta=otherbeta))
}
set.seed(1234)

# specify rhierMnlRwMixture defaults
mubar_betak <- 0
nvar  <- 10
ncomp <- 3
a     <- rep(5, ncomp)
nu    <- nvar + 3
Amu   <- 0.01
V     <- nu*diag(c(rep(1,nvar-1),1))
ndraw <- 10000
defaultprior <- drawprior(mubar_betak, nvar, ncomp, a, nu, Amu, V, ndraw)

# plot priors under defaults
par(mfrow=c(1,3))
trimhist <- -20
hist(defaultprior$betakstar, breaks=40, col="magenta", 
     main="Beta_k_star", xlab="", ylab="", yaxt="n")
hist(defaultprior$betak[defaultprior$betak>trimhist],
     breaks=40, col="magenta", main="Beta_k",
     xlab="", ylab="", yaxt="n", xlim=c(trimhist,0))
hist(defaultprior$otherbeta, breaks=40, col="magenta",
     main="Other Beta", xlab="", ylab="", yaxt="n")

We see that the hyperprior values for constrained logit parameters are far from uninformative. As a result, rhierMnlRwMixture implements different default priors for parameters when sign constraints are imposed. In particular, a_μ = 0.1, ν = k + 15, and V = ν * diag(d) where d_i = 4 if β_ik is unconstrained and d_i = 0.1 if β_ik is constrained. Additionally, μ̄_m = 0 if unconstrained and μ̄_m = 2 otherwise. As the following plots show, this yields substantially less informative hyperpriors on β_ik^* without significantly affecting the hyperpriors on β_ik or β_ij (j ≠ k).

# adjust priors for constraints
mubar_betak <- 2
nvar  <- 10
ncomp <- 3 
a     <- rep(5, ncomp)
nu    <- nvar + 15
Amu   <- 0.1 
V     <- nu*diag(c(rep(4,nvar-1),0.1))
ndraw <- 10000
tightprior <- drawprior(mubar_betak, nvar, ncomp, a, nu, Amu, V, ndraw)

# plot priors under adjusted values
par(mfrow=c(1,3))
trimhist <- -20
hist(tightprior$betakstar, breaks=40, col="magenta", 
     main="Beta_k_star", xlab="", ylab="", yaxt="n")
hist(tightprior$betak[tightprior$betak>trimhist],
     breaks=40, col="magenta", main="Beta_k",
     xlab="", ylab="", yaxt="n", xlim=c(trimhist,0))
hist(tightprior$otherbeta, breaks=40, col="magenta", 
     main="Other Beta", xlab="", ylab="", yaxt="n")

Example

Here we demonstrate the implementation of the hierarchical multinomial logit model with sign-constrained parameters. We return to the camera data used in Example 3 of the “bayesm Overview” Vignette. This dataset contains conjoint choice data for 332 respondents who evaluated digital cameras. The data are stored in a lists-of-lists format with one list per respondent, and each respondent’s list having two elements: a vector of choices (y) and a matrix of covariates (X). Notice the dimensions: there is one value for each choice occasion in each individual’s y vector but one row per alternative in each individual’s X matrix, making nrow(x) = 5 × length(y) because there are 5 alternatives per choice occasion.

library(bayesm)
data(camera)
length(camera)

## [1] 332

str(camera[[1]])

## List of 2
##  $ y: int [1:16] 1 2 2 4 2 2 1 1 1 2 ...
##  $ X: num [1:80, 1:10] 0 1 0 0 0 0 1 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:80] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:10] "canon" "sony" "nikon" "panasonic" ...

As shown next, the first 4 covariates are binary indicators for the brands Canon, Sony, Nikon, and Panasonic. These correspond to choice (y) values of 1, 2, 3, and 4. y can also take the value 5, indicating that the respondent chose “none”. The data include binary indicators for two levels of pixel count, zoom strength, swivel video display capability, and wifi connectivity. The last covaritate is price, recorded in hundreds of U.S. dollars so that the magnitude of the expected price coefficient is such that the default prior settings in rhierMnlRwMixture do not need to be adjusted.

str(camera[[1]]$y)

##  int [1:16] 1 2 2 4 2 2 1 1 1 2 ...

str(as.data.frame(camera[[1]]$X))

## 'data.frame':    80 obs. of  10 variables:
##  $ canon    : num  0 1 0 0 0 0 1 0 0 0 ...
##  $ sony     : num  0 0 0 1 0 0 0 1 0 0 ...
##  $ nikon    : num  1 0 0 0 0 0 0 0 1 0 ...
##  $ panasonic: num  0 0 1 0 0 1 0 0 0 0 ...
##  $ pixels   : num  0 1 0 0 0 1 1 1 1 0 ...
##  $ zoom     : num  1 1 0 1 0 0 0 0 0 0 ...
##  $ video    : num  0 0 0 0 0 0 1 1 1 0 ...
##  $ swivel   : num  1 1 0 1 0 1 0 0 1 0 ...
##  $ wifi     : num  0 1 1 0 0 0 1 1 1 0 ...
##  $ price    : num  0.79 2.29 1.29 2.79 0 2.79 0.79 1.79 1.29 0 ...

Let’s say we would like to estimate the part-worths of the various attributes of these digital cameras using a multinomial logit model. To incorporate individual-level heterogeneous effects, we elect to use a hierarchical (i.e., random coefficient) specification. Further, we believe that despite the heterogeneity, each consumer’s estimate price response (β_i, price) should be negative, which we will impose with a sign constraint. Following the above discussion, we use the default priors, which “adjust” automatically when sign constraints are imposed.

SignRes <- c(rep(0,nvar-1),-1)

data  <- list(lgtdata=camera, p=5)
prior <- list(mubar=mubar, Amu=Amu, ncomp=ncomp, a=a, nu=nu, V=V, SignRes=SignRes)
mcmc  <- list(R=1e4, nprint=0)

out <- rhierMnlRwMixture(Data=data, Prior=prior, Mcmc=mcmc)

While much can be done to analyze the output, we will focus here on the constrained parameters on price. We first plot the posterior distributions for the price parameter for individuals i = 1, 2, 3. Notice that the posterior distributions for the selected individual’s price parameters lie entirely below zero.

par(mfrow=c(1,3))
ind_hist <- function(mod, i) {
  hist(mod$betadraw[i , 10, ], breaks = seq(-14,0,0.5), 
     col = "dodgerblue3", border = "grey", yaxt = "n",
     xlim = c(-14,0), xlab = "", ylab = "", main = paste("Ind.",i))
}
ind_hist(out,1)
ind_hist(out,2)
ind_hist(out,3)

Next we plot a histogram of the posterior means for the 332 individual price paramters (β_i, price):

par(mfrow=c(1,1))
hist(apply(out$betadraw[ , 10, ], 1, mean), 
     xlim = c(-20,0), breaks = 20, 
     col = "firebrick2", border = "gray", xlab = "", ylab = "",
     main = "Posterior Means for Ind. Price Params, 
             With Sign Constraint")

As a point of comparison, we re-run the model without the sign constraint using the default priors (output omitted) and provide the same set of plots. Note now that the right tail of the posterior distribution of β₂^price extends to the right of zero.

data0  <- list(lgtdata = camera, p = 5)
prior0 <- list(ncomp = 5)
mcmc0  <- list(R = 1e4, nprint = 0)

out0 <- rhierMnlRwMixture(Data = data0, Prior = prior0, Mcmc = mcmc0)

par(mfrow=c(1,3))
ind_hist <- function(mod, i) {
  hist(mod$betadraw[i , 10, ], breaks = seq(-12,2,0.5), 
     col = "dodgerblue4", border = "grey", yaxt = "n",
     xlim = c(-12,2), xlab = "", ylab = "", main = paste("Ind.",i))
}
ind_hist(out0,1)
ind_hist(out0,2)
ind_hist(out0,3)

par(mfrow=c(1,1))
hist(apply(out0$betadraw[ , 10, ], 1, mean), 
     xlim = c(-15,5), breaks = 20, 
     col = "firebrick3", border = "gray", 
     xlab = "", ylab = "",
     main = "Posterior Means for Ind. Price Params, 
             No Sign Constraint")

_ Last updated July 2019._