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

Description

rhierMnlRwMixtureParallel implements a MCMC algorithm for a hierarchical multinomial logit model with a mixture of normals heterogeneity distribution.

Usage

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

Value

A list of sharded partitions where each partition contains the following:

compdraw: A list (length: R/keep) where each list contains `mu` (vector, length: `ncomp`) and `rooti` (matrix, shape: ncomp $\times$ ncomp)
probdraw: A $(R/keep) \times (ncomp)$ matrix that reports the probability that each draw came from a particular component
Deltadraw: A $(R/keep) \times (nz \times nvar)$ matrix of draws of Delta, first row is initial value. This could be null if Z is NULL

Arguments

Data: A list containing: `p`- number of choice alternatives, `lgtdata` - A $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, federico.bumbaca@colorado.edu

Details

Model and Priors

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

$\beta_i$ = $Z\Delta$[i,] + $u_i$
Note: Z$\Delta$ is the matrix Z $ \times \Delta$ and [i,] refers to $i$th row of this product
Delta is an $nz \times 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:`	A $nlgt$ size list of multinominal lgtdata
`lgtdata[[i]]$y:`	$n_i \times 1$ vector of multinomial outcomes (1, ..., p) for $i$th unit
`lgtdata[[i]]$X:`	$(n_i \times p) \times nvar$ design matrix for $i$th unit
`Z:`	An $(nlgt) \times nz$ matrix of unit characteristics
`p:`	number of choice alternatives

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

`a:`	$ncomp \times 1$ vector of Dirichlet prior parameters (def: `rep(5,ncomp)`)
`deltabar:`	$(nz \times nvar) \times 1$ vector of prior means (def: 0)
`Ad:`	prior precision matrix for vec(D) (def: 0.01*I)
`mubar:`	$nvar \times 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

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)

References

Bumbaca, Federico (Rico), Sanjog Misra, and Peter E. Rossi (2020), "Scalable Target Marketing: Distributed Markov Chain Monte Carlo for Bayesian Hierarchical Models", Journal of Marketing Research, 57(6), 999-1018.

Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

Examples

Run this code


# \donttest{

R=500

######### Single Component with rhierMnlRwMixtureParallel########
##parameters
p = 3 # num of choice alterns                            
ncoef = 3  
nlgt = 1000                          
nz = 2
Z = matrix(runif(nz*nlgt),ncol=nz)
Z = t(t(Z)-apply(Z,2,mean)) # demean Z

ncomp = 1 # no of mixture components
Delta = matrix(c(1,0,1,0,1,2),ncol=2)
comps = NULL
comps[[1]] = list(mu=c(0,2,1),rooti=diag(rep(1,3)))
pvec = c(1)

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(5,nlgt) 
for (i in 1:nlgt) 
{
  if (is.null(Z))
  {
    betai=as.vector(bayesm::rmixture(1,pvec,comps)$x)
  } else {
    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)
}

## set MCMC parameters
Prior1=list(ncomp=ncomp) 
keep=1
Mcmc1=list(R=R,keep=keep) 
Data1=list(list(p=p,lgtdata=simlgtdata,Z=Z))
s = 1
Data2 = partition_data(Data1, s=s)

out2 = parallel::mclapply(Data2, FUN = rhierMnlRwMixtureParallel, Prior = Prior1,
Mcmc = Mcmc1, mc.cores = s, mc.set.seed = FALSE)

######### Multiple Components with rhierMnlRwMixtureParallel########
##parameters
R = 500
p=3 # num of choice alterns                            
ncoef=3  
nlgt=1000  # 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
Delta=matrix(c(1,0,1,0,1,2),ncol=2) 

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

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(5,nlgt) 
for (i in 1:nlgt) 
{
  if (is.null(Z))
  {
    betai=as.vector(bayesm::rmixture(1,pvec,comps)$x)
  } else {
    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)
}

## set parameters for priors and Z
Prior1=list(ncomp=ncomp) 
keep=1
Mcmc1=list(R=R,keep=keep) 
Data1=list(list(p=p,lgtdata=simlgtdata,Z=Z))
s = 1
Data2 = partition_data(Data1,s)

out2 = parallel::mclapply(Data2, FUN = rhierMnlRwMixtureParallel, Prior = Prior1,
Mcmc = Mcmc1, mc.cores = s, mc.set.seed = FALSE)
# }

Run the code above in your browser using DataLab

`lgtdata:`	A \(nlgt\) size list of multinominal lgtdata
`lgtdata[[i]]$y:`	\(n_i \times 1\) vector of multinomial outcomes (1, ..., p) for \(i\)th unit
`lgtdata[[i]]$X:`	\((n_i \times p) \times nvar\) design matrix for \(i\)th unit
`Z:`	An \((nlgt) \times nz\) matrix of unit characteristics
`p:`	number of choice alternatives

`a:`	\(ncomp \times 1\) vector of Dirichlet prior parameters (def: `rep(5,ncomp)`)
`deltabar:`	\((nz \times nvar) \times 1\) vector of prior means (def: 0)
`Ad:`	prior precision matrix for vec(D) (def: 0.01*I)
`mubar:`	\(nvar \times 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