Description Usage Arguments Details Value Note Author(s) References See Also Examples
View source: R/rhierMnlRwMixture_rcpp.r
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.
1  rhierMnlRwMixture(Data, Prior, Mcmc)

Data 
list(lgtdata, Z, p) 
Prior 
list(a, deltabar, Ad, mubar, Amu, nu, V, a, ncomp, SignRes) 
Mcmc 
list(R, keep, nprint, s, w) 
y_i ~ MNL(X_i,β_i) with i = 1, …, length(lgtdata) and where β_i is nvar x 1
β_i = ZΔ[i,] + u_i
Note: ZΔ is the matrix Z * Δ and [i,] refers to ith row of this product
Delta is an nz x nvar array
u_i ~ N(μ_{ind},Σ_{ind}) with ind ~ multinomial(pvec)
pvec ~ dirichlet(a)
delta = vec(Δ) ~ N(deltabar, A_d^{1})
μ_j ~ N(mubar, Σ_j (x) Amu^{1})
Σ_j ~ IW(nu, V)
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
: 0.01 is too small for many applications.
See chapter 5 of Rossi et al for full discussion.
Data = list(lgtdata, Z, p)
[Z
optional]
lgtdata:  list of nlgt=length(lgtdata) lists with each crosssection unit MNL data 
lgtdata[[i]]$y:  n_i x 1 vector of multinomial outcomes (1, ..., m) 
lgtdata[[i]]$X:  n_i*p x nvar design matrix for ith unit 
Z:  nreg x nz matrix of unit chars (def: vector of ones) 
p:  number of choice alternatives 
Prior = list(a, deltabar, Ad, mubar, Amu, nu, V, a, ncomp, SignRes)
[all but ncomp
are optional]
a:  ncomp x 1 vector of Dirichlet prior parameters (def: rep(5,ncomp) ) 
deltabar:  nz*nvar x 1 vector of prior means (def: 0) 
Ad:  prior precision matrix for vec(D) (def: 0.01*I) 
mubar:  nvar x 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 
SignRes:  nvar x 1 vector of sign restrictions on the coefficient estimates (def: rep(0,nvar) )

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) 
If β_ik has a sign restriction: β_ik = SignRes[k] * exp(β*_ik)
To use sign restrictions on the coefficients, SignRes
must be an nvar x 1 vector containing values of either 0, 1, or 1. The value 0 means there is no sign restriction, 1 ensures that the coefficient is negative, and 1 ensures that the coefficient is positive. For example, if SignRes = c(0,1,1)
, the first coefficient is unconstrained, the second will be positive, and the third will be negative.
The sign restriction is implemented such that if the the k'th β has a nonzero sign restriction (i.e., it is constrained), we have β_k = SignRes[k] * exp(β*_k).
The sign restrictions (if used) will be reflected in the betadraw
output. However, the unconstrained mixture components are available in nmix
. Important: Note that draws from nmix
are distributed according to the mixture of normals but not the coefficients in betadraw
.
Care should be taken when selecting priors on any sign restricted coefficients. See related vignette for additional information.
nmix
Detailsnmix
is a list with 3 components. Several functions in the bayesm
package that involve a Dirichlet Process or mixtureofnormals return nmix
. Across these functions, a common structure is used for nmix
in order to utilize generic summary and plotting functions.
probdraw:  ncomp x R/keep matrix that reports the probability that each draw came from a particular component 
zdraw:  R/keep x nobs matrix that indicates which component each draw is assigned to (here, null) 
compdraw:  A list of R/keep lists of ncomp lists. Each of the innermost lists has 2 elemens: a vector of draws for mu and a matrix of draws for the Cholesky root of Sigma .

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 beta draws 
nmix 
a list containing: 
loglike 
R/keep x 1 vector of loglikelihood for each kept draw 
SignRes 
nvar x 1 vector of sign restrictions 
Note: as of version 2.02 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^{(1w)} x like_pooled^{((n_i/N)*w)}
Large R
values may be required (>20,000).
Peter Rossi, Anderson School, UCLA, [email protected].
For further discussion, see Chapter 5, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm1
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 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92  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 # num 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(0.4, 0.2, 0.4)
## simulate from MNL model conditional on X matrix
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(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){
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)
## fit model without sign constraints
out1 = rhierMnlRwMixture(Data=Data1, Prior=Prior1, Mcmc=Mcmc1)
cat("Summary of Delta draws", fill=TRUE)
summary(out1$Deltadraw, tvalues=as.vector(Delta))
cat("Summary of Normal Mixture Distribution", fill=TRUE)
summary(out1$nmix)
## plotting examples
if(0) {
plot(out1$betadraw)
plot(out1$nmix)
}
## fit model with constraint that beta_i2 < 0 forall i
Prior2 = list(ncomp=5, SignRes=c(0,1,0))
out2 = rhierMnlRwMixture(Data=Data1, Prior=Prior2, Mcmc=Mcmc1)
cat("Summary of Delta draws", fill=TRUE)
summary(out2$Deltadraw, tvalues=as.vector(Delta))
cat("Summary of Normal Mixture Distribution", fill=TRUE)
summary(out2$nmix)
## plotting examples
if(0) {
plot(out2$betadraw)
plot(out2$nmix)
}

Table of Y values pooled over all units
ypooled
1 2 3
5994 6618 2388
Starting MCMC Inference for Hierarchical Logit:
Normal Mixture with 5 components for first stage prior
3 alternatives; 3 variables in X
for 300 crosssectional units
Prior Parms:
nu = 6
V
[,1] [,2] [,3]
[1,] 6 0 0
[2,] 0 6 0
[3,] 0 0 6
mubar
[,1] [,2] [,3]
[1,] 0 0 0
Amu
[,1]
[1,] 0.01
a
[1] 5 5 5 5 5
deltabar
[1] 0 0 0 0 0 0
Ad
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.01 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01
MCMC Parms:
s= 1.374 w= 0.1 R= 10 keep= 5 nprint= 100
initializing Metropolis candidate densities for 300 units ...
completed unit # 50
completed unit # 100
completed unit # 150
completed unit # 200
completed unit # 250
completed unit # 300
MCMC Iteration (est time to end  min)
Total Time Elapsed: 0.00
Summary of Delta draws
fewer than 100 draws submitted
Summary of Normal Mixture Distribution
fewer than 100 draws submitted
Table of Y values pooled over all units
ypooled
1 2 3
5994 6618 2388
Starting MCMC Inference for Hierarchical Logit:
Normal Mixture with 5 components for first stage prior
3 alternatives; 3 variables in X
for 300 crosssectional units
Prior Parms:
nu = 6
V
[,1] [,2] [,3]
[1,] 6 0 0
[2,] 0 6 0
[3,] 0 0 6
mubar
[,1] [,2] [,3]
[1,] 0 0 0
Amu
[,1]
[1,] 0.01
a
[1] 5 5 5 5 5
deltabar
[1] 0 0 0 0 0 0
Ad
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.01 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 0.01 0.00 0.00 0.00 0.00
[3,] 0.00 0.00 0.01 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.01 0.00 0.00
[5,] 0.00 0.00 0.00 0.00 0.01 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.01
MCMC Parms:
s= 1.374 w= 0.1 R= 10 keep= 5 nprint= 100
initializing Metropolis candidate densities for 300 units ...
completed unit # 50
completed unit # 100
completed unit # 150
completed unit # 200
completed unit # 250
completed unit # 300
MCMC Iteration (est time to end  min)
Total Time Elapsed: 0.00
Summary of Delta draws
fewer than 100 draws submitted
Summary of Normal Mixture Distribution
fewer than 100 draws submitted
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.