rhierMnlRwMixtureParallel: MCMC Algorithm for Hierarchical Multinomial Logit with...

In scalablebayesm: Distributed Markov Chain Monte Carlo for Bayesian Inference in Marketing

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)

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.

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 ith 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 \betas 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 ith unit
lgtdata[[i]]$X:	(n_i \times p) \times nvar design matrix for ith 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: 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

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

R:	number of MCMC draws
keep:	MCMC thinning parameter -- keep every keepth 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)

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

Author(s)

Federico Bumbaca, Leeds School of Business, University of Colorado Boulder, [e-mail]

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.

See Also

partition_data, drawPosteriorParallel, combine_draws, rheteroMnlIndepMetrop

Examples

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)

scalablebayesm documentation built on April 3, 2025, 7:55 p.m.