rheteroMnlIndepMetrop: Independence Metropolis-Hastings Algorithm for Draws From Multinomial Distribution

Description

rheteroMnlIndepMetrop implements an Independence Metropolis algorithm with a Gibbs sampler.

Usage

rheteroMnlIndepMetrop(Data, draws, Mcmc)

Value

A list containing:

betadraw: A matrix of size $R \times nvar$ containing the drawn beta values from the Gibbs sampling procedure.

Arguments

Data: A list of s partitions where each partition includes: `p`- number of choice alternatives, `lgtdata` - An $nlgt$ size list of multinomial logistic data, and optional `Z`- matrix of unit characteristics.
draws: A list of draws returned from either rhierMnlRwMixtureParallel.
Mcmc: A list with one required parameter: `R`-number of iterations, and optional parameters: `keep` and `nprint`.

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/shards$ 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:`	A list of s partitions where each partition include $(nlgt/number of shards) \times nz$ matrix of unit characteristics
`p:`	number of choice alternatives

draws: A matrix with $R$ rows and $nlgt$ columns of beta draws.

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

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)

betadraws = parallel::mclapply(out2,FUN=drawPosteriorParallel,Z=Z, 
Prior = Prior1, Mcmc = Mcmc1, mc.cores=s,mc.set.seed = FALSE)
betadraws = combine_draws(betadraws, R)

out_indep = parallel::mclapply(Data2, FUN=rheteroMnlIndepMetrop, draws = betadraws, 
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)

betadraws = parallel::mclapply(out2,FUN=drawPosteriorParallel,Z=Z, 
Prior = Prior1, Mcmc = Mcmc1, mc.cores=s,mc.set.seed = FALSE)
betadraws = combine_draws(betadraws, R)

out_indep = parallel::mclapply(Data2, FUN=rheteroMnlIndepMetrop, draws = betadraws, 
Mcmc = Mcmc1, mc.cores = s, mc.set.seed = FALSE)
# }

Run the code above in your browser using DataLab

`lgtdata:`	A \(nlgt/shards\) 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:`	A list of s partitions where each partition include \((nlgt/number of shards) \times nz\) matrix of unit characteristics
`p:`	number of choice alternatives