rhierMnlRwMixture {bayesm} R Documentation

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

### 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, perossichi@gmail.com.

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

```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[] = list(mu=c(0,-1,-2),   rooti=diag(rep(1,3)))
comps[] = list(mu=c(0,-1,-2)*2, rooti=diag(rep(1,3)))
comps[] = 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)

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

## plotting examples
if(0) {
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)