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

Description

rhierMnlDP is a MCMC algorithm for a hierarchical multinomial logit with a Dirichlet Process Prior for 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 there are 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-parameteric method in the sense that the DP prior can accomodate heterogeniety of an unknown form.

Usage

rhierMnlDP(Data, Prior, Mcmc)

Arguments

Data

list(p,lgtdata,Z) ( Z is optional)

Prior

list(deltabar,Ad,Prioralpha,lambda_hyper) (all 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

adraw

R/keep draws of hyperparm a

vdraw

R/keep draws of hyperparm v

nudraw

R/keep draws of hyperparm nu

Istardraw

R/keep draws of number of unique components

alphadraw

R/keep draws of number of DP tightness parameter

loglike

R/keep draws of log-likelihood

Details

Model: $y_i$ $\sim$ $MNL(X_i,\beta_i)$. i=1,…, length(lgtdata). $\theta_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.

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

Priors: $\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

Lists contain:

Data:

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$ by nvar design matrix for ith unit

Prior:

deltabarnz*nvar vector of prior means (def: 0)
Ad prior prec matrix for vec(D) (def: .01I)

Prioralpha:

Istarminexpected number of components at lower bound of support of alpha def(1)
Istarmaxexpected number of components at upper bound of support of alpha (def: min(50,.1*nlgt))
powerpower parameter for alpha prior (def: .8)

lambda_hyper:

alimdefines support of a distribution (def: (.01,2))
nulimdefines support of nu distribution (def: (.01,3))
vlimdefines support of v distribution (def: (.1,4))

Mcmc:

R number of mcmc draws
keep thinning parm, keep every keepth draw
nprint print the estimated time remaining for every nprint'th draw (def: 100)
maxuniq storage constraint on the number of unique components
gridsize number of discrete points for hyperparameter priors,def: 20

References

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

Examples

Run this code

##
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=20000} 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(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(.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 with a mixture of 3 normals
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(1){
## 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 Data and Mcmc lists
keep=5
Mcmc1=list(R=R,keep=keep)
Data1=list(p=p,lgtdata=simlgtdata,Z=Z)

out=rhierMnlDP(Data=Data1,Mcmc=Mcmc1)

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

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

Run the code above in your browser using DataLab