bayesm (version 3.0-2)

rhierMnlRwMixture: 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(p,lgtdata,Z) ( Z is optional)

Prior

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

Mcmc

list(s,w,R,keep,nprint) (R required)

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 draws of betas

nmix

list of 3 components, probdraw, NULL, compdraw

loglike

log-likelihood for each kept draw (length R/keep)

Details

Model: \(y_i\) \(\sim\) \(MNL(X_i,\beta_i)\). i=1,…, length(lgtdata). \(\beta_i\) is nvar x 1.

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

\(u_i\) \(\sim\) \(N(\mu_{ind},\Sigma_{ind})\). \(ind\) \(\sim\) multinomial(pvec).

Priors: \(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)

Lists contain:

• p p is number of choice alternatives

• lgtdatalist of lists with each cross-section unit MNL data

• lgtdata[[i]]$y \(n_i\) vector of multinomial outcomes (1,…,m)

• lgtdata[[i]]$X \(n_i\)*p by nvar design matrix for ith unit

• avector of length ncomp of Dirichlet prior parms (def: rep(5,ncomp))

• deltabarnz*nvar vector of prior means (def: 0)

• Ad prior prec matrix for vec(D) (def: .01I)

• mubar nvar x 1 prior mean vector for normal comp mean (def: 0)

• Amu prior precision for normal comp mean (def: .01I)

• nu d.f. parm for IW prior on norm comp Sigma (def: nvar+3)

• V pds location parm for IW prior on norm comp Sigma (def: nuI)

• a Dirichlet prior parameter (def: 5)

• ncomp number of components used in normal mixture

• s scaling parm for RW Metropolis (def: 2.93/sqrt(nvar))

• w fractional likelihood weighting parm (def: .1)

• R number of MCMC draws

• keep MCMC thinning parm: keep every keepth draw (def: 1)

• nprint print the estimated time remaining for every nprint'th draw (def: 100)

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 5. http://www.perossi.org/home/bsm-1

See Also

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

simmnlwX= function(n,X,beta) {
  ##  simulate from MNL model conditional on X matrix
  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){
## set if(1) above to produce plots
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)

out=rhierMnlRwMixture(Data=Data1,Prior=Prior1,Mcmc=Mcmc1)

cat("Summary of Delta draws",fill=TRUE)
summary(out$Deltadraw,tvalues=as.vector(Delta))
cat("Summary of Normal Mixture Distribution",fill=TRUE)
summary(out$nmix)

if(0) {
## plotting examples
plot(out$betadraw)
plot(out$nmix)
}