rhierMnlRwMixture: MCMC Algorithm for Hierarchical Multinomial Logit with Mixture-of-Normals Heterogeneity

Description

rhierMnlRwMixture is a MCMC algorithm for a hierarchical multinomial logit with a mixture of normals heterogeneity distribution. This is a hybrid Gibbs Sampler with a RW Metropolis step for the MNL coefficients for each panel unit.

Usage

rhierMnlRwMixture(Data, Prior, Mcmc)

Value

A list containing:

Deltadraw: $R/keep x nz*nvar$ matrix of draws of Delta, first row is initial value
betadraw: $nlgt x nvar x R/keep$ array of beta draws
nmix: a list containing: probdraw, zdraw, compdraw (see “nmix Details” section)
loglike: $R/keep x 1$ vector of log-likelihood for each kept draw
SignRes: $nvar x 1$ vector of sign restrictions

Arguments

Data: list(lgtdata, Z, p)
Prior: list(a, deltabar, Ad, mubar, Amu, nu, V, a, ncomp, SignRes)
Mcmc: list(R, keep, nprint, s, w)

Author

Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.

Details

Model and Priors

$y_i$ $\sim$ $MNL(X_i,\beta_i)$ with $i = 1, \ldots,$ length(lgtdata) and where $\beta_i$ is $nvar x 1$

$\beta_i$ = $Z\Delta$[i,] + $u_i$
Note: Z$\Delta$ is the matrix Z * $\Delta$ and [i,] refers to $i$th row of this product
Delta is an $nz x nvar$ array

$u_i$ $\sim$ $N(\mu_{ind},\Sigma_{ind})$ with $ind$ $\sim$ multinomial(pvec)

$pvec$ $\sim$ dirichlet(a)
$delta = vec(\Delta)$ $\sim$ $N(deltabar, A_d^{-1})$
$\mu_j$ $\sim$ $N(mubar, \Sigma_j (x) Amu^{-1})$
$\Sigma_j$ $\sim$ $IW(nu, V)$

Note: $Z$ should NOT include an intercept and is centered for ease of interpretation. The mean of each of the nlgt $\beta$s is the mean of the normal mixture. Use summary() to compute this mean from the compdraw output.

Be careful in assessing prior parameter Amu: 0.01 is too small for many applications. See chapter 5 of Rossi et al for full discussion.

Argument Details

Data = list(lgtdata, Z, p) [Z optional]

`lgtdata:`	list of `nlgt=length(lgtdata)` lists with each cross-section unit MNL data
`lgtdata[[i]]$y:`	$n_i x 1$ vector of multinomial outcomes (1, ..., p)
`lgtdata[[i]]$X:`	$n_i*p x nvar$ design matrix for $i$th unit
`Z:`	$nreg x nz$ matrix of unit chars (def: vector of ones)
`p:`	number of choice alternatives

Prior = list(a, deltabar, Ad, mubar, Amu, nu, V, a, ncomp, SignRes) [all but ncomp are optional]

`a:`	$ncomp x 1$ vector of Dirichlet prior parameters (def: `rep(5,ncomp)`)
`deltabar:`	$nz*nvar x 1$ vector of prior means (def: 0)
`Ad:`	prior precision matrix for vec(D) (def: 0.01*I)
`mubar:`	$nvar x 1$ prior mean vector for normal component mean (def: 0 if unrestricted; 2 if restricted)
`Amu:`	prior precision for normal component mean (def: 0.01 if unrestricted; 0.1 if restricted)
`nu:`	d.f. parameter for IW prior on normal component Sigma (def: nvar+3 if unrestricted; nvar+15 if restricted)
`V:`	PDS location parameter for IW prior on normal component Sigma (def: nuI if unrestricted; nuD if restricted with d_pp = 4 if unrestricted and d_pp = 0.01 if restricted)
`ncomp:`	number of components used in normal mixture
`SignRes:`	$nvar x 1$ vector of sign restrictions on the coefficient estimates (def: `rep(0,nvar)`)

Mcmc = list(R, keep, nprint, s, w) [only R required]

`R:`	number of MCMC draws
`keep:`	MCMC thinning parameter -- keep every `keep`th draw (def: 1)
`nprint:`	print the estimated time remaining for every `nprint`'th draw (def: 100, set to 0 for no print)
`s:`	scaling parameter for RW Metropolis (def: 2.93/`sqrt(nvar)`)
`w:`	fractional likelihood weighting parameter (def: 0.1)

Sign Restrictions

If $\beta_ik$ has a sign restriction: $\beta_ik = SignRes[k] * exp(\beta*_ik)$

To use sign restrictions on the coefficients, SignRes must be an $nvar x 1$ vector containing values of either 0, -1, or 1. The value 0 means there is no sign restriction, -1 ensures that the coefficient is negative, and 1 ensures that the coefficient is positive. For example, if SignRes = c(0,1,-1), the first coefficient is unconstrained, the second will be positive, and the third will be negative.

The sign restriction is implemented such that if the the $k$'th $\beta$ has a non-zero sign restriction (i.e., it is constrained), we have $\beta_k = SignRes[k] * exp(\beta*_k)$.

The sign restrictions (if used) will be reflected in the betadraw output. However, the unconstrained mixture components are available in nmix. Important: Note that draws from nmix are distributed according to the mixture of normals but not the coefficients in betadraw.

Care should be taken when selecting priors on any sign restricted coefficients. See related vignette for additional information.

`nmix` Details

nmix is a list with 3 components. Several functions in the bayesm package that involve a Dirichlet Process or mixture-of-normals return nmix. Across these functions, a common structure is used for nmix in order to utilize generic summary and plotting functions.

`probdraw:`	$ncomp x R/keep$ matrix that reports the probability that each draw came from a particular component
`zdraw:`	$R/keep x nobs$ matrix that indicates which component each draw is assigned to (here, null)
`compdraw:`	A list of $R/keep$ lists of $ncomp$ lists. Each of the inner-most lists has 2 elemens: a vector of draws for `mu` and a matrix of draws for the Cholesky root of `Sigma`.

References

For further discussion, see Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

Examples

Run this code

if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=10000} else {R=10}
set.seed(66)

p = 3                                # num of choice alterns
ncoef = 3  
nlgt = 300                           # num of cross sectional units
nz = 2
Z = matrix(runif(nz*nlgt),ncol=nz)
Z = t(t(Z) - apply(Z,2,mean))        # demean Z
ncomp = 3                            # num of mixture components
Delta = matrix(c(1,0,1,0,1,2),ncol=2)

comps=NULL
comps[[1]] = list(mu=c(0,-1,-2),   rooti=diag(rep(1,3)))
comps[[2]] = list(mu=c(0,-1,-2)*2, rooti=diag(rep(1,3)))
comps[[3]] = list(mu=c(0,-1,-2)*4, rooti=diag(rep(1,3)))
pvec = c(0.4, 0.2, 0.4)

##  simulate from MNL model conditional on X matrix
simmnlwX= function(n,X,beta) {
  k = length(beta)
  Xbeta = X%*%beta
  j = nrow(Xbeta) / n
  Xbeta = matrix(Xbeta, byrow=TRUE, ncol=j)
  Prob = exp(Xbeta)
  iota = c(rep(1,j))
  denom = Prob%*%iota
  Prob = Prob/as.vector(denom)
  y = vector("double",n)
  ind = 1:j
  for (i in 1:n) { 
    yvec = rmultinom(1, 1, Prob[i,])
    y[i] = ind%*%yvec
  }
  return(list(y=y, X=X, beta=beta, prob=Prob))
}

## simulate data
simlgtdata = NULL
ni = rep(50, 300)
for (i in 1:nlgt) {
  betai = Delta%*%Z[i,] + as.vector(rmixture(1,pvec,comps)$x)
   Xa = matrix(runif(ni[i]*p,min=-1.5,max=0), ncol=p)
   X = createX(p, na=1, nd=NULL, Xa=Xa, Xd=NULL, base=1)
   outa = simmnlwX(ni[i], X, betai)
   simlgtdata[[i]] = list(y=outa$y, X=X, beta=betai)
}

## plot betas
if(0){
  bmat = matrix(0, nlgt, ncoef)
  for(i in 1:nlgt) {bmat[i,] = simlgtdata[[i]]$beta}
  par(mfrow = c(ncoef,1))
  for(i in 1:ncoef) { hist(bmat[,i], breaks=30, col="magenta") }
}

## set parms for priors and Z
Prior1 = list(ncomp=5)
keep = 5
Mcmc1 = list(R=R, keep=keep)
Data1 = list(p=p, lgtdata=simlgtdata, Z=Z)

## fit model without sign constraints
out1 = rhierMnlRwMixture(Data=Data1, Prior=Prior1, Mcmc=Mcmc1)

cat("Summary of Delta draws", fill=TRUE)
summary(out1$Deltadraw, tvalues=as.vector(Delta))

cat("Summary of Normal Mixture Distribution", fill=TRUE)
summary(out1$nmix)

## plotting examples
if(0) {
  plot(out1$betadraw)
  plot(out1$nmix)
}

## fit model with constraint that beta_i2 < 0 forall i
Prior2 = list(ncomp=5, SignRes=c(0,-1,0))
out2 = rhierMnlRwMixture(Data=Data1, Prior=Prior2, Mcmc=Mcmc1)

cat("Summary of Delta draws", fill=TRUE)
summary(out2$Deltadraw, tvalues=as.vector(Delta))

cat("Summary of Normal Mixture Distribution", fill=TRUE)
summary(out2$nmix)

## plotting examples
if(0) {
  plot(out2$betadraw)
  plot(out2$nmix)
}

Run the code above in your browser using DataLab

`a:`	\(ncomp x 1\) vector of Dirichlet prior parameters (def: `rep(5,ncomp)`)
`deltabar:`	\(nz*nvar x 1\) vector of prior means (def: 0)
`Ad:`	prior precision matrix for vec(D) (def: 0.01*I)
`mubar:`	\(nvar x 1\) prior mean vector for normal component mean (def: 0 if unrestricted; 2 if restricted)
`Amu:`	prior precision for normal component mean (def: 0.01 if unrestricted; 0.1 if restricted)
`nu:`	d.f. parameter for IW prior on normal component Sigma (def: nvar+3 if unrestricted; nvar+15 if restricted)
`V:`	PDS location parameter for IW prior on normal component Sigma (def: nuI if unrestricted; nuD if restricted with d_pp = 4 if unrestricted and d_pp = 0.01 if restricted)
`ncomp:`	number of components used in normal mixture
`SignRes:`	\(nvar x 1\) vector of sign restrictions on the coefficient estimates (def: `rep(0,nvar)`)

`probdraw:`	\(ncomp x R/keep\) matrix that reports the probability that each draw came from a particular component
`zdraw:`	\(R/keep x nobs\) matrix that indicates which component each draw is assigned to (here, null)
`compdraw:`	A list of \(R/keep\) lists of \(ncomp\) lists. Each of the inner-most lists has 2 elemens: a vector of draws for `mu` and a matrix of draws for the Cholesky root of `Sigma`.