MCMC Algorithm for Hierarchical Multinomial Logit with Mixture of Normals Heterogeneity

Share:

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

1
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)

Details

Model:
y_i ~ MNL(X_i,β_i). i=1,..., length(lgtdata). β_i is nvar x 1.

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

u_i ~ N(μ_{ind},Σ_{ind}). ind ~ multinomial(pvec).

Priors:
pvec ~ dirichlet (a)
delta= vec(Δ) ~ N(deltabar,A_d^{-1})
μ_j ~ N(mubar,Σ_j (x) Amu^{-1})
Σ_j ~ 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)

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)

Note

More on probdraw component of nmix list:
R/keep x ncomp matrix of draws of probs of mixture components (pvec)
More on compdraw component of return value list:

  • compdraw[[i]] the ith draw of components for mixtures

  • compdraw[[i]][[j]] ith draw of the jth normal mixture comp

  • compdraw[[i]][[j]][[1]] ith draw of jth normal mixture comp mean vector

  • compdraw[[i]][[j]][[2]] ith draw of jth normal mixture cov parm (rooti)

Note: Z should not include an intercept and is centered for ease of interpretation. The mean of each of the nlgt β 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. .01 is too small for many applications. See Rossi et al, chapter 5 for full discussion.

Note: as of version 2.0-2 of bayesm, the fractional weight parameter has been changed to a weight between 0 and 1. w is the fractional weight on the normalized pooled likelihood. This differs from what is in Rossi et al chapter 5, i.e.

like_i^{(1-w)} x like_pooled^{((n_i/N)*w)}

Large R values may be required (>20,000).

Author(s)

Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
##
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)
}

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.