scalablebayesm (version 0.2)

rhierMnlDPParallel: MCMC Algorithm for Hierarchical Multinomial Logit with Dirichlet Process Prior Heterogeneity

Description

rhierMnlDPParallel is an MCMC algorithm for a hierarchical multinomial logit with a Dirichlet Process prior describing the distribution of heteorogeneity. A base normal model is used so that the DP can be interpreted as allowing for a mixture of normals with as many components as panel units. This is a hybrid Gibbs Sampler with a RW Metropolis step for the MNL coefficients for each panel unit. This procedure can be interpreted as a Bayesian semi-parametric method in the sense that the DP prior can accommodate heterogeneity of an unknown form.

Usage

rhierMnlDPParallel(Data, Prior, Mcmc, verbose = FALSE)

Value

A list containing the following elements:

• nmix: A list with the following components:

  • probdraw: A matrix of size (R/keep) x (ncomp, containing the probability draws at each Gibbs iteration.

  • compdraw: A list containing the drawn mixture components at each Gibbs iteration.

• Deltadraw (optional): A matrix of size (R/keep) x (nz * nvar), containing the delta draws, if Z is not NULL. If Z is NULL, this element is not included.

Arguments

Data

A list of: `p`- number of choice alternatives, `lgtdata` - An \(nlgt\) size list of multinomial logistic data, and optional `Z`- matrix of unit characteristics.

Prior

A list with one required parameter: `ncomp`, and optional parameters: `mubar`, `Amu`, `nu`, `V`, `Ad`, `deltaBar`, and `a`.

Mcmc

A list with one required parameter: 'R' - number of iterations, and optional parameters: `s`, `w`, `keep`, `nprint`, and `drawcomp`.

verbose

If TRUE, enumerates model parameters and timing information.

Author

Federico Bumbaca, Leeds School of Business, University of Colorado Boulder,

Details

Model and Priors

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

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

\(\beta_i\) \(\sim\) \(N(\mu_i, \Sigma_i)\)

\(\theta_i = (\mu_i, \Sigma_i)\) \(\sim\) \(DP(G_0(\lambda), alpha)\)

\(G_0(\lambda):\)
\(\mu_i | \Sigma_i\) \(\sim\) \(N(0, \Sigma_i (x) a^{-1})\)
\(\Sigma_i\) \(\sim\) \(IW(nu, nu*v*I)\)
\(delta = vec(\Delta)\) \(\sim\) \(N(deltabar, A_d^{-1})\)

\(\lambda(a, nu, v):\)
\(a\) \(\sim\) uniform[alim[1], alimb[2]]
\(nu\) \(\sim\) dim(data)-1 + exp(z)
\(z\) \(\sim\) uniform[dim(data)-1+nulim[1], nulim[2]]
\(v\) \(\sim\) uniform[vlim[1], vlim[2]]

\(alpha\) \(\sim\) \((1-(alpha-alphamin) / (alphamax-alphamin))^{power}\)
alpha = alphamin then expected number of components = Istarmin
alpha = alphamax then expected number of components = Istarmax

\(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.

We parameterize the prior on \(\Sigma_i\) such that \(mode(\Sigma) = nu/(nu+2) vI\). The support of nu enforces a non-degenerate IW density; \(nulim[1] > 0\).

The default choices of alim, nulim, and vlim determine the location and approximate size of candidate "atoms" or possible normal components. The defaults are sensible given a reasonable scaling of the X variables. You want to insure that alim is set for a wide enough range of values (remember a is a precision parameter) and the v is big enough to propose Sigma matrices wide enough to cover the data range.

A careful analyst should look at the posterior distribution of a, nu, v to make sure that the support is set correctly in alim, nulim, vlim. In other words, if we see the posterior bunched up at one end of these support ranges, we should widen the range and rerun.

If you want to force the procedure to use many small atoms, then set nulim to consider only large values and set vlim to consider only small scaling constants. Set alphamax to a large number. This will create a very "lumpy" density estimate somewhat like the classical Kernel density estimates. Of course, this is not advised if you have a prior belief that densities are relatively smooth.

References

Bumbaca, F. (Rico), Misra, S., & Rossi, P. E. (2020). Scalable Target Marketing: Distributed Markov Chain Monte Carlo for Bayesian Hierarchical Models. Journal of Marketing Research, 57(6), 999-1018.

Examples

# \donttest{
set.seed(66)
R = 500

p = 3                                # num of choice alterns
ncoef = 3  
nlgt = 1000                           # num of cross sectional units
nz = 2
nvar = 3
Z = matrix(runif(nz*nlgt), ncol=nz)
Z = t(t(Z)-apply(Z,2,mean))          # demean Z
ncomp = 3                            # no 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(2,3)))
comps[[2]] = list(mu=c(0,-1,-2)*2, rooti=diag(rep(2,3)))
comps[[3]] = list(mu=c(0,-1,-2)*4, rooti=diag(rep(2,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 with a mixture of 3 normals
simlgtdata = NULL
ni = rep(50,nlgt)
for (i in 1:nlgt) {
  betai = Delta%*%Z[i,] + as.vector(bayesm::rmixture(1,pvec,comps)$x)
  Xa = matrix(runif(ni[i]*p,min=-1.5,max=0), ncol=p)
  X = bayesm::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

  old.par = par(no.readonly = TRUE)
  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") }
  par(old.par)

## set Data and Mcmc lists
keep = 5
Mcmc1 = list(R=R, keep=keep)
Prior1=list(ncomp=1)
Data1 = list(p=p, lgtdata=simlgtdata, Z=Z)
Data2 = partition_data(Data = list(Data1), s = 1)

out = parallel::mclapply(Data2, FUN = rhierMnlDPParallel,
Prior = Prior1, Mcmc = Mcmc1, mc.cores = 1, mc.set.seed = FALSE)
# }