# rhierMnlRwMixture: MCMC Algorithm for Hierarchical Multinomial Logit with... In bayesm: Bayesian Inference for Marketing/Micro-Econometrics

## 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

 `1` ```rhierMnlRwMixture(Data, Prior, Mcmc) ```

## 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)

## Details

#### Model and Priors

y_i ~ MNL(X_i,β_i) with i = 1, …, length(lgtdata) and where β_i is nvar x 1

β_i = [i,] + u_i
Note: ZΔ is the matrix Z * Δ and [i,] refers to ith row of this product
Delta is an nz x nvar array

u_i ~ N(μ_{ind},Σ_{ind}) with ind ~ multinomial(pvec)

pvec ~ dirichlet(a)
delta = vec(Δ) ~ N(deltabar, A_d^{-1})
μ_j ~ N(mubar, Σ_j (x) Amu^{-1})
Σ_j ~ IW(nu, V)

Note: Z should NOT include an intercept and is centered for ease of interpretation. The mean of each of the `nlgt` β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 ith 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: nu*I if unrestricted; nu*D 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 β_ik has a sign restriction: β_ik = SignRes[k] * exp(β*_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 β has a non-zero sign restriction (i.e., it is constrained), we have β_k = SignRes[k] * exp(β*_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`.

## 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

## Note

Note: as of version 2.0-2 of `bayesm`, the fractional weight parameter has been changed to a weight between 0 and 1. w is the fractional weight on the normalized pooled likelihood. This differs from what is in Rossi et al chapter 5, i.e.

like_i^{(1-w)} x like_pooled^{((n_i/N)*w)}

Large `R` values may be required (>20,000).

## Author(s)

Peter Rossi, Anderson School, UCLA, [email protected].

## References

For further discussion, see Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1

`rmnlIndepMetrop`

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92``` ```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) } ```

### Example output

```Table of Y values pooled over all units
ypooled
1    2    3
5994 6618 2388

Starting MCMC Inference for Hierarchical Logit:
Normal Mixture with 5 components for first stage prior
3  alternatives;  3  variables in X
for  300  cross-sectional units

Prior Parms:
nu = 6
V
[,1] [,2] [,3]
[1,]    6    0    0
[2,]    0    6    0
[3,]    0    0    6
mubar
[,1] [,2] [,3]
[1,]    0    0    0
Amu
[,1]
[1,] 0.01
a
[1] 5 5 5 5 5
deltabar
[1] 0 0 0 0 0 0
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.01 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01

MCMC Parms:
s= 1.374  w=  0.1  R=  10  keep=  5  nprint=  100

initializing Metropolis candidate densities for  300  units ...
completed unit # 50
completed unit # 100
completed unit # 150
completed unit # 200
completed unit # 250
completed unit # 300
MCMC Iteration (est time to end - min)
Total Time Elapsed: 0.00
Summary of Delta draws
fewer than 100 draws submitted
Summary of Normal Mixture Distribution
fewer than 100 draws submitted
Table of Y values pooled over all units
ypooled
1    2    3
5994 6618 2388

Starting MCMC Inference for Hierarchical Logit:
Normal Mixture with 5 components for first stage prior
3  alternatives;  3  variables in X
for  300  cross-sectional units

Prior Parms:
nu = 6
V
[,1] [,2] [,3]
[1,]    6    0    0
[2,]    0    6    0
[3,]    0    0    6
mubar
[,1] [,2] [,3]
[1,]    0    0    0
Amu
[,1]
[1,] 0.01
a
[1] 5 5 5 5 5
deltabar
[1] 0 0 0 0 0 0
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.01 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01

MCMC Parms:
s= 1.374  w=  0.1  R=  10  keep=  5  nprint=  100

initializing Metropolis candidate densities for  300  units ...
completed unit # 50
completed unit # 100
completed unit # 150
completed unit # 200
completed unit # 250
completed unit # 300
MCMC Iteration (est time to end - min)
Total Time Elapsed: 0.00
Summary of Delta draws
fewer than 100 draws submitted
Summary of Normal Mixture Distribution
fewer than 100 draws submitted
```

bayesm documentation built on Dec. 21, 2018, 9:04 a.m.